diff --git a/cmake/custom/fetch_mrcpp.cmake b/cmake/custom/fetch_mrcpp.cmake
index 517e320..b9a2b0a 100644
--- a/cmake/custom/fetch_mrcpp.cmake
+++ b/cmake/custom/fetch_mrcpp.cmake
@@ -18,7 +18,7 @@ else()
     GIT_REPOSITORY
       https://github.com/MRChemSoft/mrcpp.git
     GIT_TAG
-      8107aabe28d6e75f04d66c95a94c157731484eae
+      e2a8eea2b4c3470aaef16e4156feaa1b0570b5d8
     )
 
   set(CMAKE_CXX_COMPILER ${CMAKE_CXX_COMPILER})
diff --git a/src/vampyr/operators/convolutions.h b/src/vampyr/operators/convolutions.h
index 72d8bb6..b127b81 100644
--- a/src/vampyr/operators/convolutions.h
+++ b/src/vampyr/operators/convolutions.h
@@ -7,6 +7,7 @@
 #include <MRCPP/operators/IdentityConvolution.h>
 #include <MRCPP/operators/PoissonOperator.h>
 #include <MRCPP/operators/TimeEvolutionOperator.h>
+#include <MRCPP/operators/HeatOperator.h>
 #include <MRCPP/treebuilders/apply.h>
 
 namespace vampyr {
@@ -15,6 +16,7 @@ void cartesian_convolution(pybind11::module &);
 void helmholtz_operator(pybind11::module &);
 void poisson_operator(pybind11::module &);
 void time_evolution_operator(pybind11::module &m);
+void heat_operator(pybind11::module &m);
 
 template <int D> void convolutions(pybind11::module &m) {
     namespace py = pybind11;
@@ -53,6 +55,7 @@ template <int D> void convolutions(pybind11::module &m) {
     if constexpr (D == 3) helmholtz_operator(m);
     if constexpr (D == 3) poisson_operator(m);
     if constexpr (D == 1) time_evolution_operator(m);
+    if constexpr (D == 1) heat_operator(m);
 }
 
 void cartesian_convolution(pybind11::module &m) {
@@ -135,6 +138,12 @@ void time_evolution_operator(pybind11::module &m)
              "finest_scale"_a,
              "imaginary"_a,
              "max_Jpower"_a = 20)
+        .def(py::init<const MultiResolutionAnalysis<1> &, double, double, bool, int>(),
+             "mra"_a,
+             "prec"_a,
+             "time"_a,
+             "imaginary"_a,
+             "max_Jpower"_a = 40)
         .def(
             "__call__",
             [](TimeEvolutionOperator<1> &T, FunctionTree<1> *inp) {
@@ -145,4 +154,26 @@ void time_evolution_operator(pybind11::module &m)
             "inp"_a);
 }
 
+
+void heat_operator(pybind11::module &m)
+{
+    namespace py = pybind11;
+    using namespace mrcpp;
+    using namespace pybind11::literals;
+
+    py::class_<HeatOperator<1>, ConvolutionOperator<1>>(m, "HeatOperator")
+        .def(py::init<const MultiResolutionAnalysis<1> &, double, double>(),
+             "mra"_a,
+             "time"_a,
+             "prec"_a)
+        .def(
+            "__call__",
+            [](HeatOperator<1> &T, FunctionTree<1> *inp) {
+                auto out = std::make_unique<FunctionTree<1>>(inp->getMRA());
+                apply<1>(T.getBuildPrec(), *out, T, *inp);
+                return out;
+            },
+            "inp"_a);
+}
+
 } // namespace vampyr