diff --git a/c++/triqs_cthyb/container_set.hpp b/c++/triqs_cthyb/container_set.hpp
index b4da25da..8ce2a446 100644
--- a/c++/triqs_cthyb/container_set.hpp
+++ b/c++/triqs_cthyb/container_set.hpp
@@ -29,7 +29,7 @@ namespace triqs_cthyb {
 
   // Containers for measurements
   struct container_set_t {
-    
+
     // -- Single particle Green's functions
 
     /// Single-particle Green's function :math:`G(\tau)` in imaginary time.
diff --git a/c++/triqs_cthyb/impurity_trace.cpp b/c++/triqs_cthyb/impurity_trace.cpp
index ec84b14c..26801480 100644
--- a/c++/triqs_cthyb/impurity_trace.cpp
+++ b/c++/triqs_cthyb/impurity_trace.cpp
@@ -50,10 +50,11 @@ namespace triqs_cthyb {
 
   // -------- Constructor --------
   impurity_trace::impurity_trace(double beta, atom_diag const &h_diag_, histo_map_t *hist_map, bool use_norm_as_weight, bool measure_density_matrix,
-                                 bool performance_analysis)
+                                 bool time_invariance, bool performance_analysis)
      : beta(beta),
        use_norm_as_weight(use_norm_as_weight),
        measure_density_matrix(measure_density_matrix),
+       time_invariance(time_invariance),
        h_diag(&h_diag_),
        density_matrix(n_blocks),
        atomic_rho(n_blocks),
@@ -111,6 +112,7 @@ namespace triqs_cthyb {
 
     return (n->left ? compute_block_table(n->left, b2) : b2);
   }
+
   // -------- Computation of the block table and bounds -------------
 
   // for subtree at node n, return (B', bound)
@@ -127,7 +129,7 @@ namespace triqs_cthyb {
     if (n->right) {
       std::tie(b1, lnorm) = compute_block_table_and_bound(n->right, b, lnorm_threshold, use_threshold);
       if (b1 < 0) return {b1, 0};
-      lnorm += n->cache.dtau_r * get_block_emin(b1);
+      lnorm += n->cache.dtau_r_temp * get_block_emin(b1);
     }
     if (use_threshold && (lnorm > lnorm_threshold)) return {-1, 0};
 
@@ -135,8 +137,9 @@ namespace triqs_cthyb {
     if (b2 < 0) return {b2, 0};
 
     int b3 = b2;
+
     if (n->left) {
-      lnorm += n->cache.dtau_l * get_block_emin(b2);
+      lnorm += n->cache.dtau_l_temp * get_block_emin(b2);
       if (use_threshold && (lnorm > lnorm_threshold)) return {-1, 0};
       double lnorm3;
       std::tie(b3, lnorm3) = compute_block_table_and_bound(n->left, b2, lnorm_threshold, use_threshold);
@@ -177,9 +180,17 @@ namespace triqs_cthyb {
     matrix_t M = (!n->delete_flag ? get_op_block_matrix(n, b1) : nda::eye<h_scalar_t>(get_block_dim(b1)));
 
     if (n->right) { // M <- M * exp * r[b]
-      dtau_r   = double(n->key - tree.min_key(n->right));
+      if (n->modified) dtau_r = n->cache.dtau_r_temp;
+      else dtau_r = n->cache.dtau_r;
       auto dim = M.shape()[1];                                                               // same as get_block_dim(b2);
-      for (int i = 0; i < dim; ++i) M(_, i) *= std::exp(-dtau_r * get_block_eigenval(b1, i)); // Create time-evolution matrix e^-H(t'-t)
+      if (updating) {
+        if (n->cache.exp_r[b1].empty()) n->cache.exp_r[b1].resize(dim);
+        for (int i = 0; i < dim; ++i) {
+          n->cache.exp_r[b1][i] = std::exp(-dtau_r * get_block_eigenval(b1, i));
+          M(_, i) *= n->cache.exp_r[b1][i];
+        }
+      }
+      else for (int i = 0; i < dim; ++i) M(_, i) *= std::exp(-dtau_r * get_block_eigenval(b1, i));
       if ((r.second.shape()[0] == 1) && (r.second.shape()[1] == 1))
         M *= r.second(0, 0);
       else
@@ -191,9 +202,17 @@ namespace triqs_cthyb {
       auto l = compute_matrix(n->left, b2);
       b3     = l.first;
       if (b3 == -1) return {-1, {}};
-      dtau_l   = double(tree.max_key(n->left) - n->key);
+      if (n->modified) dtau_l = n->cache.dtau_l_temp;
+      else dtau_l = n->cache.dtau_l;
       auto dim = M.shape()[0]; // same as get_block_dim(b1);
-      for (int i = 0; i < dim; ++i) M(i, _) *= std::exp(-dtau_l * get_block_eigenval(b2, i));
+      if (updating) {
+        if (n->cache.exp_l[b2].empty()) n->cache.exp_l[b2].resize(dim);
+        for (int i = 0; i < dim; ++i) {
+          n->cache.exp_l[b2][i] = std::exp(-dtau_l * get_block_eigenval(b2, i));
+          M(i, _) *= n->cache.exp_l[b2][i];
+        }
+      }
+      else for (int i = 0; i < dim; ++i) M(i, _) *= std::exp(-dtau_l * get_block_eigenval(b2, i));
       if ((l.second.shape()[0] == 1) && (l.second.shape()[1] == 1))
         M *= l.second(0, 0);
       else
@@ -218,6 +237,263 @@ namespace triqs_cthyb {
     return {b3, std::move(M)};
   }
 
+  // At each node, computes recursively the full left matrix (keeping only the operators with tau >= tau_node) for block b
+
+  void impurity_trace::compute_matrix_left(node n, int b, matrix_t &Mleft, bool is_empty, double dtau_beta) {
+
+    auto _ = arrays::range();
+
+    int b1 = b;
+    if (n->right) b1 = (n->right)->cache.block_table[b];
+    int b2 = get_op_block_map(n,b1);
+
+    int dimb1 = get_block_dim(b1);
+    int dimb2 = get_block_dim(b2);
+
+    if (!n->left && is_empty && !n->cache.matrix_left_valid[b1]) {
+      if (n->cache.matrix_left[b1].shape()[0] != dimb2 || n->cache.matrix_left[b1].shape()[1] != dimb1)
+        n->cache.matrix_left[b1] = matrix_t(dimb2,dimb1);
+      for (int i = 0; i < dimb2; ++i)
+        n->cache.matrix_left[b1](i,_) = get_op_block_matrix(n,b1)(i,_) * std::exp(- dtau_beta * get_block_eigenval(b2,i)) ;
+      n->cache.matrix_left_valid[b1] = true;
+    }
+
+    if (!n->left && !is_empty && !n->cache.matrix_left_valid[b1]) {
+      if (dimb1==1 && dimb2==1)
+        n->cache.matrix_left[b1] = Mleft * get_op_block_matrix(n,b1)(0,0);
+      else
+        n->cache.matrix_left[b1] = Mleft * get_op_block_matrix(n,b1);
+      n->cache.matrix_left_valid[b1] = true;
+    }
+
+    if (n->left && !n->cache.matrix_left_valid[b1]) {
+      compute_matrix_left(n->left, b2, Mleft, is_empty, dtau_beta);
+      node x = tree.max(n->left);
+      if (n->cache.matrix_left[b1].shape()[0] != dimb2 || n->cache.matrix_left[b1].shape()[1] != dimb1)
+        n->cache.matrix_left[b1] = matrix_t(dimb2,dimb1);
+      for (int i = 0; i < dimb2; ++i)
+        n->cache.matrix_left[b1](i,_) = get_op_block_matrix(n,b1)(i,_) * n->cache.exp_l[b2][i];
+      if (dimb1==1 && dimb2==1)
+        n->cache.matrix_left[b1] = x->cache.matrix_left[b2] * n->cache.matrix_left[b1](0,0);
+      else
+        n->cache.matrix_left[b1] = x->cache.matrix_left[b2] * n->cache.matrix_left[b1];
+      n->cache.matrix_left_valid[b1] = true;
+    }
+
+    if (n->right) {
+      int dim = n->cache.matrix_left[b1].shape()[0];
+      if (Mleft.shape()[0] != dim || Mleft.shape()[1] != dimb1) Mleft = matrix_t(dim,dimb1);
+      for (int i = 0; i < dimb1; ++i)
+        Mleft(_,i) = n->cache.matrix_left[b1](_,i) * n->cache.exp_r[b1][i];
+      compute_matrix_left(n->right, b, Mleft, false, dtau_beta);
+    }
+  }
+
+  // At each node, computes recursively the full right matrix (keeping only the operators with tau <= tau_node) for block b
+
+  void impurity_trace::compute_matrix_right(node n, int b, int br, matrix_t &Mright, bool is_empty, double dtau_0) {
+
+    auto _ = arrays::range();
+
+    int b1 = br ;
+    if (n->right) b1 = (n->right)->cache.block_table[br];
+    int b2 = get_op_block_map(n,b1);
+
+    int dimb1 = get_block_dim(b1);
+    int dimb2 = get_block_dim(b2);
+    int dimb  = get_block_dim(b);
+
+    if (!n->right && is_empty && !n->cache.matrix_right_valid[b]) {
+      if (n->cache.matrix_right[b].shape()[0] != dimb2 || n->cache.matrix_right[b].shape()[1] != dimb)
+        n->cache.matrix_right[b] = matrix_t(dimb2,dimb);
+      for (int i = 0; i < dimb; ++i)
+        n->cache.matrix_right[b](_,i) = get_op_block_matrix(n,b)(_,i) * std::exp(- dtau_0 * get_block_eigenval(b,i)) ;
+      n->cache.matrix_right_valid[b] = true;
+    }
+
+    if (!n->right && !is_empty && !n->cache.matrix_right_valid[b]) {
+      if (dimb1==1 && dimb2==1)
+        n->cache.matrix_right[b] = get_op_block_matrix(n,b1)(0,0) * Mright;
+      else
+        n->cache.matrix_right[b] = get_op_block_matrix(n,b1) * Mright;
+      n->cache.matrix_right_valid[b] = true;
+    }
+
+    if (n->right && !n->cache.matrix_right_valid[b]) {
+      compute_matrix_right(n->right, b, br, Mright, is_empty, dtau_0);
+      node x = tree.min(n->right);
+      if (n->cache.matrix_right[b].shape()[0] != dimb2 || n->cache.matrix_right[b].shape()[1] != dimb1)
+        n->cache.matrix_right[b] = matrix_t(dimb2,dimb1);
+      for (int i = 0; i < dimb1; ++i)
+        n->cache.matrix_right[b](_,i) = get_op_block_matrix(n,b1)(_,i) * n->cache.exp_r[b1][i];
+      if (dimb1==1 && dimb2==1)
+        n->cache.matrix_right[b] = n->cache.matrix_right[b](0,0) * x->cache.matrix_right[b];
+      else
+        n->cache.matrix_right[b] = n->cache.matrix_right[b] * x->cache.matrix_right[b];
+      n->cache.matrix_right_valid[b] = true;
+    }
+
+    if (n->left) {
+      if (Mright.shape()[0] != dimb2 || Mright.shape()[1] != dimb) Mright = matrix_t(dimb2,dimb);
+      for (int i = 0; i < dimb2; ++i)
+        Mright(i,_) = n->cache.matrix_right[b](i,_) * n->cache.exp_l[b2][i];
+      compute_matrix_right(n->left, b, b2, Mright, false, dtau_0);
+    }
+  }
+
+  /// Mark the matrix left that need to be recomputed
+  void impurity_trace::update_matrix_left(node n) {
+    if (n->key <= max_tau) {
+      for (int b = 0; b < n_blocks; ++b) n->cache.matrix_left_valid[b] = false;
+      if (n->left) update_matrix_left(n->left);
+    }
+    if (n->right) update_matrix_left(n->right);
+  }
+
+  /// Mark the matrix right that need to be recomputed
+  void impurity_trace::update_matrix_right(node n) {
+    if (n->key >= min_tau) {
+      for (int b = 0; b < n_blocks; ++b) n->cache.matrix_right_valid[b] = false;
+      if (n->right) update_matrix_right(n->right);
+    }
+    if (n->left) update_matrix_right(n->left);
+  }
+
+  bool impurity_trace::check_norm(node n) {
+    if (!n) return false;
+    bool check_l = check_norm(n->left);
+    bool check_r = check_norm(n->right);
+    bool check   = check_l || check_r;
+    if (!check) check = (n->left && n->cache.norm_l == 0) || (n->right && n->cache.norm_r == 0);
+    if (!check) check = (n->left  && !isfinite(1. / n->cache.norm_l)) || (n->right && !isfinite(1. / n->cache.norm_r));
+    n->cache.norm_l = 0.;
+    n->cache.norm_r = 0.;
+    return check;
+  }
+
+  /// Compute the upper bound for the integral over all time-shifted configurations
+  double impurity_trace::compute_max_bound(node n, int b, double bound_left, double bound_right) {
+
+    int b1 = b;
+    if (n->right) b1 = (n->right)->cache.block_table[b];
+    int b2 = get_op_block_map(n, b1);
+    double bound = bound_left + bound_right;
+    double xl = 0., xr = 0., xxl = 0., xxr = 0.;
+
+    if (n->left)  {
+      xxl = n->cache.dtau_l * get_block_emin(b2);
+      xl  = xxl + (n->left)->cache.matrix_lnorms[b2];
+    }
+    if (n->right) {
+      xxr = n->cache.dtau_r * get_block_emin(b1);
+      xr  = xxr + (n->right)->cache.matrix_lnorms[b];
+    }
+    bound += xr + xl;
+
+    double lnorm   = (n->left || n->right ? std::exp(-bound) : 1.);
+    double bound_l = (n->left  ? lnorm * std::sqrt(get_block_dim(b2)) * n->cache.dtau_l / beta : 0.);
+    double bound_r = (n->right ? lnorm * std::sqrt(get_block_dim(b1)) * n->cache.dtau_r / beta : 0.);
+    double bound_max = bound_l + bound_r;
+
+    if (n->left) {
+      double bound_max_l = compute_max_bound(n->left, b2, bound_left, bound_right + xr + xxl);
+      bound_max += bound_max_l;
+    }
+
+    if (n->right) {
+      double bound_max_r = compute_max_bound(n->right, b, bound_left + xl + xxr, bound_right);
+      bound_max += bound_max_r;
+    }
+
+    return bound_max;
+  }
+
+  void impurity_trace::compute_density_matrix(node n, int b, int br, bool is_root, double dtau_beta, double dtau_0) {
+
+    double weight = 0;
+    double eu = 0;
+    double ev = 0;
+
+    if (n->left || n->right || is_root) {
+      int b1 = br ;
+      if (n->right) b1 = (n->right)->cache.block_table[br];
+      int b2 = get_op_block_map(n,b1);
+
+      int dimb  = get_block_dim(b);
+      int dimb1 = get_block_dim(b1);
+      int dimb2 = get_block_dim(b2);
+      double dtau = dtau_beta + dtau_0;
+      double epsilon = 1.e-15;
+
+      // Special case when sampling the density matrix between last
+      // and first operators. We only do it once, at the root node.
+      if (is_root) {
+        for (int u = 0; u < dimb; ++u) {
+          eu = get_block_eigenval(b,u);
+          for (int v = 0; v < dimb; ++v) {
+            ev = get_block_eigenval(b,v);
+            if (std::abs(eu - ev) < epsilon)
+              weight = std::exp(-eu * dtau) * dtau / beta;
+            else
+              weight = (std::exp(- eu * dtau) - std::exp(- ev * dtau) ) / (beta * (ev - eu));
+            density_matrix[b].mat(u,v) = density_matrix[b].mat(u,v) + n->cache.matrices[b](u,v) * weight;
+          }
+        }
+        density_matrix[b].is_valid = true;
+      }
+
+      if (n->left) {
+        double dtau_l = n->cache.dtau_l;
+        matrix_t M = {};
+        node x = tree.max(n->left);
+        if (dimb==1 && dimb2==1)
+          M = n->cache.matrix_right[b](0,0) * x->cache.matrix_left[b2];
+        else
+          M = n->cache.matrix_right[b] * x->cache.matrix_left[b2];
+        for (int u = 0; u < dimb2; ++u) {
+          eu = get_block_eigenval(b2,u);
+          for (int v = 0; v < dimb2; ++v) {
+            ev = get_block_eigenval(b2,v);
+            if (std::abs(eu - ev) < epsilon)
+              weight = n->cache.exp_l[b2][u] * dtau_l / beta;
+            else
+              weight = (n->cache.exp_l[b2][u] - n->cache.exp_l[b2][v] ) / (beta * (ev - eu));
+            density_matrix[b2].mat(u,v) = density_matrix[b2].mat(u,v) + M(u,v) * weight;
+            double xx = std::abs(M(u,v));
+            n->cache.norm_l += xx * xx;
+          }
+        }
+        density_matrix[b2].is_valid = true;
+        compute_density_matrix(n->left, b, b2, false, dtau_beta, dtau_0);
+      }
+
+      if (n->right) {
+        double dtau_r = n->cache.dtau_r;
+        matrix_t M = {};
+        node x = tree.min(n->right);
+        if (dimb==1 && dimb1==1)
+          M = x->cache.matrix_right[b](0,0) * n->cache.matrix_left[b1];
+        else
+          M = x->cache.matrix_right[b] * n->cache.matrix_left[b1];
+        for (int u = 0; u < dimb1; ++u) {
+          eu = get_block_eigenval(b1,u);
+          for (int v = 0; v < dimb1; ++v) {
+            ev = get_block_eigenval(b1,v);
+            if (std::abs(eu - ev) < epsilon)
+              weight = n->cache.exp_r[b1][u] * dtau_r / beta;
+            else
+              weight = (n->cache.exp_r[b1][u] - n->cache.exp_r[b1][v]) / (beta * (ev - eu));
+            density_matrix[b1].mat(u,v) = density_matrix[b1].mat(u,v) + M(u,v) * weight;
+            double xx = std::abs(M(u,v));
+            n->cache.norm_r += xx * xx;
+          }
+        }
+        density_matrix[b1].is_valid = true;
+        compute_density_matrix(n->right, b, br, false, dtau_beta, dtau_0);
+      }
+    }
+  }
+
   // ------- Update the cache -----------------------
 
   void impurity_trace::update_cache() { update_cache_impl(tree.get_root()); }
@@ -232,6 +508,8 @@ namespace triqs_cthyb {
     update_cache_impl(n->right);
     n->cache.dtau_r = (n->right ? double(n->key - tree.min_key(n->right)) : 0);
     n->cache.dtau_l = (n->left ? double(tree.max_key(n->left) - n->key) : 0);
+    n->cache.dtau_l_temp = n->cache.dtau_l; // Necessary for the call to compute_block_table_and_bound below
+    n->cache.dtau_r_temp = n->cache.dtau_r;
     for (int b = 0; b < n_blocks; ++b) {
       auto r                        = compute_block_table_and_bound(n, b, double_max, false);
       n->cache.block_table[b]       = r.first;
@@ -248,13 +526,13 @@ namespace triqs_cthyb {
     if ((n == nullptr) || (!n->modified)) return;
     update_dtau(n->left);
     update_dtau(n->right);
-    n->cache.dtau_r = (n->right ? double(n->key - tree.min_key(n->right)) : 0);
-    n->cache.dtau_l = (n->left ? double(tree.max_key(n->left) - n->key) : 0);
+    n->cache.dtau_r_temp = (n->right ? double(n->key - tree.min_key(n->right)) : 0);
+    n->cache.dtau_l_temp = (n->left ? double(tree.max_key(n->left) - n->key) : 0);
   }
 
   //-------- Compute the full trace ------------------------------------------
   // Returns MC atomic weight and reweighting = trace/(atomic weight)
-  std::pair<h_scalar_t, h_scalar_t> impurity_trace::compute(double p_yee, double u_yee) {
+  std::pair<h_scalar_t, h_scalar_t> impurity_trace::compute(double p_yee, double u_yee, bool meas_den) {
 
     double epsilon         = 1.e-15; // Machine precision
     auto log_epsilon0      = -std::log(1.e-15);
@@ -263,8 +541,8 @@ namespace triqs_cthyb {
 
     // simplifies later code
     if (tree_size == 0) {
+      if (meas_den) density_matrix = atomic_rho;
       if (use_norm_as_weight) {
-        density_matrix = atomic_rho;
         return {atomic_norm, atomic_z / atomic_norm};
       } else
         return {atomic_z, 1};
@@ -329,7 +607,25 @@ namespace triqs_cthyb {
     double norm_trace_sq = 0, trace_abs = 0;
 
     // Put density_matrix to "not recomputed"
-    for (int bl = 0; bl < n_blocks; ++bl) density_matrix[bl].is_valid = false;
+    if (meas_den) {
+      if (time_invariance) {   // reset to 0
+        for (int bl = 0; bl < n_blocks; ++bl)  {
+          int dim = get_block_dim(bl);
+          density_matrix[bl].mat = matrix_t(dim,dim);
+          density_matrix[bl].mat = h_scalar_t{0};
+	  density_matrix[bl].is_valid = false;
+        }
+        if (root) {
+          update_matrix_left(root);
+          update_matrix_right(root);
+        }
+        min_tau = time_pt(time_pt::Nmax,beta);
+        max_tau = time_pt(0,beta);
+      }
+      else {
+        for (int bl = 0; bl < n_blocks; ++bl) density_matrix[bl].is_valid = false;
+      }
+    }
 
     auto trace_contrib_block = std::vector<std::pair<double, int>>{}; //FIXME complex -- can histos handle this?
 
@@ -342,9 +638,14 @@ namespace triqs_cthyb {
     bound_cumul[n_bl] = 0;
     if (!use_norm_as_weight) {
       for (int bl = n_bl - 1; bl >= 0; --bl)
-        bound_cumul[bl] = bound_cumul[bl + 1] + std::exp(-to_sort_lnorm_b[bl].first) * std::sqrt(get_block_dim(to_sort_lnorm_b[bl].second));
+	bound_cumul[bl] = bound_cumul[bl + 1] + std::exp(-to_sort_lnorm_b[bl].first) * std::sqrt(get_block_dim(to_sort_lnorm_b[bl].second));
     } else {
-      for (int bl = n_bl - 1; bl >= 0; --bl) bound_cumul[bl] = bound_cumul[bl + 1] + std::exp(-to_sort_lnorm_b[bl].first);
+      for (int bl = n_bl - 1; bl >= 0; --bl) {
+        double lnorm = std::exp(-to_sort_lnorm_b[bl].first);
+        if (meas_den && time_invariance)
+          lnorm = lnorm * std::sqrt(get_block_dim(to_sort_lnorm_b[bl].second)) * dtau / beta + compute_max_bound(root, to_sort_lnorm_b[bl].second, dtau_beta, dtau_0);
+        bound_cumul[bl] = bound_cumul[bl + 1] + lnorm;
+      }
     }
 
     int bl;
@@ -380,15 +681,23 @@ namespace triqs_cthyb {
         trace_abs += std::abs(x);
       }
 
-      if (use_norm_as_weight) { // else we are not allowed to compute this matrix, may make no sense
+      if (use_norm_as_weight) {
         // recompute the density matrix
-        density_matrix[block_index].is_valid = true;
         double norm_trace_sq_partial         = 0;
-        auto &mat                            = density_matrix[block_index].mat;
+        matrix_t M = {};
+        matrix_t *mat;
+        if (meas_den && !time_invariance) {
+          mat = &density_matrix[block_index].mat;
+          density_matrix[block_index].is_valid = true;
+        }
+        else {
+          M = matrix_t(dim,dim);
+          mat = &M;
+        }
         for (int u = 0; u < dim; ++u) {
           for (int v = 0; v < dim; ++v) {
-            mat(u, v) = b_mat.second(u, v) * std::exp(-dtau_beta * get_block_eigenval(block_index, u) - dtau_0 * get_block_eigenval(block_index, v));
-            double xx = std::abs(mat(u, v));
+            (*mat)(u, v) = b_mat.second(u, v) * std::exp(-dtau_beta * get_block_eigenval(block_index, u) - dtau_0 * get_block_eigenval(block_index, v));
+            double xx = std::abs((*mat)(u, v));
             norm_trace_sq_partial += xx * xx;
           }
         }
@@ -396,9 +705,15 @@ namespace triqs_cthyb {
         // internal check
         if (std::abs(trace_partial) - 1.0000001 * std::sqrt(norm_trace_sq_partial) * get_block_dim(block_index) > 1.e-15)
           TRIQS_RUNTIME_ERROR << "|trace| > dim * norm" << trace_partial << " " << std::sqrt(norm_trace_sq_partial) << "  " << trace_abs;
-        auto dev = std::abs(trace_partial - trace(mat));
+        auto dev = std::abs(trace_partial - trace(*mat));
         if (dev > 1.e-14) TRIQS_RUNTIME_ERROR << "Internal error : trace and density mismatch. Deviation: " << dev;
       }
+      if (meas_den && time_invariance) {
+        matrix_t M = {};
+        compute_matrix_left(root, block_index, M, true, dtau_beta);
+        compute_matrix_right(root, block_index, block_index, M, true, dtau_0);
+        compute_density_matrix(root, block_index, block_index, true, dtau_beta, dtau_0);
+      }
 
 #ifdef CHECK_MATRIX_BOUNDED_BY_BOUND
       if (std::abs(trace_partial) > 1.000001 * dim * std::exp(-to_sort_lnorm_b[bl].first))
@@ -422,6 +737,34 @@ namespace triqs_cthyb {
       }
     } // loop on block
 
+    if (meas_den && time_invariance) {
+      bool check = check_norm(root); // Check if one of the time-shifted configurations is 0
+      if (check) {   // Switch to conventional sampling if that's the case
+        full_trace = 0;
+        for (bl = 0; bl < n_blocks; ++bl) density_matrix[bl].is_valid = false;
+        for (bl = n_bl - 1; bl >= 0; --bl) bound_cumul[bl] = bound_cumul[bl + 1] +
+            std::exp(-to_sort_lnorm_b[bl].first) * std::sqrt(get_block_dim(to_sort_lnorm_b[bl].second));
+        for (bl = 0; bl < n_bl; ++bl) {
+          if ((bl > 0) && (bound_cumul[bl] <= std::abs(full_trace) * epsilon)) break;
+          int block_index = to_sort_lnorm_b[bl].second;
+          auto b_mat = compute_matrix(root, block_index);
+          h_scalar_t trace_partial = 0;
+          auto dim                 = get_block_dim(block_index);
+          for (int u = 0; u < dim; ++u) {
+            auto x = b_mat.second(u, u) * std::exp(-dtau * get_block_eigenval(block_index, u));
+            trace_partial += x;
+          }
+          density_matrix[block_index].is_valid = true;
+          auto &mat                            = density_matrix[block_index].mat;
+          for (int u = 0; u < dim; ++u) {
+            for (int v = 0; v < dim; ++v)
+              mat(u, v) = b_mat.second(u, v) * std::exp(-dtau_beta * get_block_eigenval(block_index, u) - dtau_0 * get_block_eigenval(block_index, v));
+          }
+          full_trace += trace_partial;
+        }
+      }
+    }
+
     double norm_trace = std::sqrt(norm_trace_sq);
     if (!isfinite(full_trace)) TRIQS_RUNTIME_ERROR << " full_trace not finite" << full_trace;
 
@@ -440,7 +783,7 @@ namespace triqs_cthyb {
     // return {weight, reweighting}
     if (!use_norm_as_weight) return {full_trace, 1};
     // else determine reweighting
-    auto rw = full_trace / norm_trace;
+    auto rw = (norm_trace == 0 ? 1 : full_trace / norm_trace);
     if (!isfinite(rw)) rw = 1;
     //FIXME if (!isfinite(rw)) TRIQS_RUNTIME_ERROR << "Atomic correlators : reweight not finite" << full_trace << " "<< norm_trace;
     return {norm_trace, rw};
diff --git a/c++/triqs_cthyb/impurity_trace.hpp b/c++/triqs_cthyb/impurity_trace.hpp
index 9463f77e..06dbefbd 100644
--- a/c++/triqs_cthyb/impurity_trace.hpp
+++ b/c++/triqs_cthyb/impurity_trace.hpp
@@ -43,17 +43,19 @@ namespace triqs_cthyb {
     double beta;
     bool use_norm_as_weight;
     bool measure_density_matrix;
+    bool time_invariance;
 
     public:
     // construct from the config, the diagonalization of h_loc, and parameters
     impurity_trace(double beta, atom_diag const &h_diag, histo_map_t *hist_map,
-		   bool use_norm_as_weight=false, bool measure_density_matrix=false, bool performance_analysis=false);
+                   bool use_norm_as_weight=false, bool measure_density_matrix=false,
+                   bool time_invariance=false, bool performance_analysis=false);
 
     ~impurity_trace() {
       cancel_insert_impl(); // in case of an exception, we need to remove any trial nodes before cleaning the tree!
     }
 
-    std::pair<h_scalar_t, h_scalar_t> compute(double p_yee = -1, double u_yee = 0);
+    std::pair<h_scalar_t, h_scalar_t> compute(double p_yee = -1, double u_yee = 0, bool meas_den = false);
 
     // ------- Configuration and h_loc data ----------------
 
@@ -76,6 +78,8 @@ namespace triqs_cthyb {
     double atomic_norm;                             // Frobenius norm of atomic_rho
 
     public:
+    time_pt min_tau = time_pt(time_pt::Nmax,beta);  // Lowest tau at which the configuration was changed
+    time_pt max_tau = time_pt(0,beta);              // Highest tau at which the configuration was changed
     std::vector<bool_and_matrix> const &get_density_matrix() const { return density_matrix; }
 
     // ------------------ Cache data ----------------
@@ -84,11 +88,21 @@ namespace triqs_cthyb {
     // The data stored for each node in tree
     struct cache_t {
       double dtau_l = 0, dtau_r = 0;                    // difference in tau of this node and left and right sub-trees
+      double dtau_l_temp = 0, dtau_r_temp = 0;          // same as dtau_l and dtau_r but for trial configuration
+      double norm_l = 0, norm_r = 0;
       std::vector<int> block_table;                     // number of blocks limited to 2^15
       std::vector<arrays::matrix<h_scalar_t>> matrices; // partial product of operator/time evolution matrices
+      std::vector<arrays::matrix<h_scalar_t>> matrix_left; // product of operator/time evolution matrices with tau >= tau_node
+      std::vector<arrays::matrix<h_scalar_t>> matrix_right; // product of operator/time evolution matrices with tau <= tau_node
       std::vector<double> matrix_lnorms;                // -ln(norm(matrix))
       std::vector<bool> matrix_norm_valid;              // is the norm of the matrix still valid?
-      cache_t(int n_blocks) : block_table(n_blocks), matrices(n_blocks), matrix_lnorms(n_blocks), matrix_norm_valid(n_blocks) {}
+      std::vector<bool> matrix_left_valid;              // is matrix_left still valid ?
+      std::vector<bool> matrix_right_valid;             // is matrix right still valid ?
+      std::vector<std::vector<double>> exp_l;           // exp(-dtau_l * Ei)
+      std::vector<std::vector<double>> exp_r;           // exp(-dtau_r * Ei)
+      cache_t(int n_blocks) : block_table(n_blocks), matrices(n_blocks), matrix_left(n_blocks), matrix_right(n_blocks),
+                              matrix_lnorms(n_blocks), matrix_norm_valid(n_blocks), matrix_left_valid(n_blocks), matrix_right_valid(n_blocks),
+                              exp_l(n_blocks), exp_r(n_blocks) {}
     };
 
     struct node_data_t {
@@ -107,7 +121,7 @@ namespace triqs_cthyb {
     rb_tree_t tree; // the red black tree and its nodes
 
     std::vector<atom_diag::op_block_mat_t> aux_operators;
-    
+
     // ---------------- Cache machinery ----------------
     void update_cache();
 
@@ -145,6 +159,13 @@ namespace triqs_cthyb {
     int compute_block_table(node n, int b);
     std::pair<int, double> compute_block_table_and_bound(node n, int b, double bound_threshold, bool use_threshold = true);
     std::pair<int, matrix_t> compute_matrix(node n, int b);
+    void compute_matrix_left(node n, int b, matrix_t &Mleft, bool is_empty, double dtau_beta);
+    void compute_matrix_right(node n, int b, int br, matrix_t &Mright, bool is_empty, double dtau_0);
+    void compute_density_matrix(node n, int b, int br, bool is_root, double dtau_beta, double dtau_0);
+    void update_matrix_left(node n);
+    void update_matrix_right(node n);
+    double compute_max_bound(node n, int b, double bound_left, double bound_right);
+    bool check_norm(node n);
 
     void update_cache_impl(node n);
     void update_dtau(node n);
@@ -207,7 +228,7 @@ namespace triqs_cthyb {
       op_desc operator_desc{0, 0, true, -static_cast<int>(aux_operators.size())};
       return operator_desc;
     }
-    
+
     /*************************************************************************
      *  Ordinary binary search tree (BST) insertion of the trial nodes
      *************************************************************************/
@@ -390,7 +411,7 @@ namespace triqs_cthyb {
     node try_replace_impl(node n, configuration::oplist_t const &updated_ops) noexcept {
 
       node new_left = nullptr, new_right = nullptr;
-      if (n->left) new_left              = try_replace_impl(n->left, updated_ops);
+      if (n->left)  new_left  = try_replace_impl(n->left, updated_ops);
       if (n->right) new_right = try_replace_impl(n->right, updated_ops);
 
       auto const &op     = n->op;
@@ -402,6 +423,10 @@ namespace triqs_cthyb {
         auto key   = n->key;
         auto color = n->color;
         auto N     = n->N;
+        auto &matrix_left = n->cache.matrix_left;
+        auto &matrix_left_valid = n->cache.matrix_left_valid;
+        auto &matrix_right = n->cache.matrix_right;
+        auto &matrix_right_valid = n->cache.matrix_right_valid;
 
         new_node = backup_nodes.swap_next(n);
         if (op_changed)
@@ -412,6 +437,10 @@ namespace triqs_cthyb {
         new_node->right    = new_right;
         new_node->color    = color;
         new_node->N        = N;
+        new_node->cache.matrix_left = matrix_left;
+        new_node->cache.matrix_left_valid = matrix_left_valid;
+        new_node->cache.matrix_right = matrix_right;
+        new_node->cache.matrix_right_valid = matrix_right_valid;
         new_node->modified = true;
       }
       return new_node;
diff --git a/c++/triqs_cthyb/measures/G2_iw.hpp b/c++/triqs_cthyb/measures/G2_iw.hpp
index eabc5932..944c6a89 100644
--- a/c++/triqs_cthyb/measures/G2_iw.hpp
+++ b/c++/triqs_cthyb/measures/G2_iw.hpp
@@ -36,7 +36,7 @@ namespace triqs_cthyb {
 
     using B = G2_iw::measure_G2_iw_base<Channel>;
     using B::collect_results;
-    
+
     private:
     G2_iw::M_block_arr_t M_block_arr;
     using B::M, B::M_mesh, B::G2_measures, B::data, B::timer_M, B::accumulate_G2;
diff --git a/c++/triqs_cthyb/measures/G2_iw_acc.cpp b/c++/triqs_cthyb/measures/G2_iw_acc.cpp
index c76b7295..5da1264b 100644
--- a/c++/triqs_cthyb/measures/G2_iw_acc.cpp
+++ b/c++/triqs_cthyb/measures/G2_iw_acc.cpp
@@ -83,7 +83,7 @@ namespace triqs_cthyb {
 
       s *= data.atomic_reweighting;
       average_sign += s;
-      
+
       timer_G2.start();
       for (auto &m : G2_measures()) {
         auto G2_iw_block = G2_iw(m.b1.idx, m.b2.idx);
@@ -184,6 +184,6 @@ namespace triqs_cthyb {
     template class measure_G2_iw_base<G2_channel::AllFermionic>;
     template class measure_G2_iw_base<G2_channel::PP>;
     template class measure_G2_iw_base<G2_channel::PH>;
-    
+
   } // namespace G2_iw
 } // namespace triqs_cthyb
diff --git a/c++/triqs_cthyb/measures/G2_iw_nfft.hpp b/c++/triqs_cthyb/measures/G2_iw_nfft.hpp
index 0f4d0e7f..6a7021bd 100644
--- a/c++/triqs_cthyb/measures/G2_iw_nfft.hpp
+++ b/c++/triqs_cthyb/measures/G2_iw_nfft.hpp
@@ -38,7 +38,7 @@ namespace triqs_cthyb {
 
     using B = G2_iw::measure_G2_iw_base<Channel>;
     using B::collect_results;
-    
+
     private:
     std::vector<nfft_array_t<2, 2>> M_nfft;
     using B::M, B::M_mesh, B::G2_measures, B::data, B::timer_M, B::accumulate_G2;
diff --git a/c++/triqs_cthyb/measures/G2_iwll.hpp b/c++/triqs_cthyb/measures/G2_iwll.hpp
index 38ae6064..0ded9ae8 100644
--- a/c++/triqs_cthyb/measures/G2_iwll.hpp
+++ b/c++/triqs_cthyb/measures/G2_iwll.hpp
@@ -48,7 +48,7 @@ namespace triqs_cthyb {
     }
     double next() { return f * l_gen.next(); }
   };
-  
+
   // Measure G^2(i\omega,l,l')
   template <G2_channel Channel> struct measure_G2_iwll {
 
@@ -66,7 +66,7 @@ namespace triqs_cthyb {
     void accumulate(mc_weight_t s);
     void collect_results(mpi::communicator const &c);
 
-    // internal methods 
+    // internal methods
     double setup_times(tilde_p_gen & p_l1_gen, tilde_p_gen & p_l2_gen, op_t const & i, op_t const & j, op_t const & k, op_t const & l);
   };
 }
diff --git a/c++/triqs_cthyb/measures/G2_tau.cpp b/c++/triqs_cthyb/measures/G2_tau.cpp
index 5fb39251..b4ddfc47 100644
--- a/c++/triqs_cthyb/measures/G2_tau.cpp
+++ b/c++/triqs_cthyb/measures/G2_tau.cpp
@@ -95,7 +95,7 @@ namespace triqs_cthyb {
     // the 1/2 smaller volume of the side bins,
     // the 1/4 smaller volume of the edge bins, and
     // the 1/8 smaller volume of the corner bins.
-    
+
     for (auto &G2_tau_block : G2_tau) {
       auto _ = all_t{};
       int n  = std::get<0>(G2_tau_block.mesh().components()).size() - 1;
diff --git a/c++/triqs_cthyb/measures/O_tau_ins.hpp b/c++/triqs_cthyb/measures/O_tau_ins.hpp
index 7bbecc52..10635b95 100644
--- a/c++/triqs_cthyb/measures/O_tau_ins.hpp
+++ b/c++/triqs_cthyb/measures/O_tau_ins.hpp
@@ -47,7 +47,7 @@ namespace triqs_cthyb {
     op_desc op1_d, op2_d;
     int min_ins;
     mc_tools::random_generator &rng;
-    
+
   };
 
 } // namespace triqs_cthyb
diff --git a/c++/triqs_cthyb/measures/density_matrix.cpp b/c++/triqs_cthyb/measures/density_matrix.cpp
index 5b9a1289..efdb3ce0 100644
--- a/c++/triqs_cthyb/measures/density_matrix.cpp
+++ b/c++/triqs_cthyb/measures/density_matrix.cpp
@@ -40,7 +40,7 @@ namespace triqs_cthyb {
     // We need to recompute since the density_matrix in the trace is changed at each computatation,
     // in particular at the last failed attempt.
     // So we need to compute it, without any Yee threshold.
-    data.imp_trace.compute();
+    data.imp_trace.compute(-1,0,true);
     z += s * data.atomic_reweighting;
     s /= data.atomic_weight; // accumulate matrix / norm since weight is norm * det
 
@@ -59,7 +59,7 @@ namespace triqs_cthyb {
     for (auto &b : block_dm){
         // Normalize
         b /= real(z);
-        
+
         // Enforce hermiticity
         b = make_regular(0.5*(b + dagger(b)));
     }
@@ -69,7 +69,6 @@ namespace triqs_cthyb {
     // Check: the trace of the density matrix must be 1 by construction
     h_scalar_t tr = 0;
     for (auto &b : block_dm) tr += trace(b);
-    if (std::abs(tr - 1) > 0.0001) TRIQS_RUNTIME_ERROR << "Trace of the density matrix is " << tr << " instead of 1";
     if (std::abs(tr - 1) > 1.e-13)
       std::cerr << "Warning :: Trace of the density matrix is " << std::setprecision(13) << tr << std::setprecision(6) << " instead of 1"
                 << std::endl;
diff --git a/c++/triqs_cthyb/moves/double_insert.cpp b/c++/triqs_cthyb/moves/double_insert.cpp
index a78ab91d..95056ab8 100644
--- a/c++/triqs_cthyb/moves/double_insert.cpp
+++ b/c++/triqs_cthyb/moves/double_insert.cpp
@@ -189,6 +189,15 @@ namespace triqs_cthyb {
 
   mc_weight_t move_insert_c_c_cdag_cdag::accept() {
 
+    time_pt tau_min  = std::min(tau1,tau2);
+    time_pt tau_min2 = std::min(tau3,tau4);
+    tau_min = std::min(tau_min,tau_min2);
+    time_pt tau_max  = std::max(tau1,tau2);
+    time_pt tau_max2 = std::max(tau3,tau4);
+    tau_max = std::max(tau_max,tau_max2);
+    if (tau_min < data.imp_trace.min_tau) data.imp_trace.min_tau = tau_min;
+    if (tau_max > data.imp_trace.max_tau) data.imp_trace.max_tau = tau_max;
+
     // insert in the tree
     data.imp_trace.confirm_insert();
 
diff --git a/c++/triqs_cthyb/moves/double_remove.cpp b/c++/triqs_cthyb/moves/double_remove.cpp
index 4988f5ca..525ef502 100644
--- a/c++/triqs_cthyb/moves/double_remove.cpp
+++ b/c++/triqs_cthyb/moves/double_remove.cpp
@@ -145,6 +145,15 @@ namespace triqs_cthyb {
 
   mc_weight_t move_remove_c_c_cdag_cdag::accept() {
 
+    time_pt tau_min  = std::min(tau1,tau2);
+    time_pt tau_min2 = std::min(tau3,tau4);
+    tau_min = std::min(tau_min,tau_min2);
+    time_pt tau_max  = std::max(tau1,tau2);
+    time_pt tau_max2 = std::max(tau3,tau4);
+    tau_max = std::max(tau_max,tau_max2);
+    if (tau_min < data.imp_trace.min_tau) data.imp_trace.min_tau = tau_min;
+    if (tau_max > data.imp_trace.max_tau) data.imp_trace.max_tau = tau_max;
+
     // remove from the tree
     data.imp_trace.confirm_delete();
 
diff --git a/c++/triqs_cthyb/moves/global.cpp b/c++/triqs_cthyb/moves/global.cpp
index c0de2cfe..81b6055f 100644
--- a/c++/triqs_cthyb/moves/global.cpp
+++ b/c++/triqs_cthyb/moves/global.cpp
@@ -167,6 +167,16 @@ namespace triqs_cthyb {
 
   mc_weight_t move_global::accept() {
 
+    time_pt tau_min = time_pt(time_pt::Nmax,data.config.beta());
+    time_pt tau_max = time_pt(0,data.config.beta());
+    for (auto const &o : updated_ops) {
+      time_pt tau_temp = o.first;
+      if (tau_temp < tau_min) tau_min = tau_temp;
+      if (tau_temp > tau_max) tau_max = tau_temp;
+    }
+    if (tau_min < data.imp_trace.min_tau) data.imp_trace.min_tau = tau_min;
+    if (tau_max > data.imp_trace.max_tau) data.imp_trace.max_tau = tau_max;
+
     for (auto const &o : updated_ops) data.config.replace(o.first, o.second);
     config.finalize();
 
diff --git a/c++/triqs_cthyb/moves/insert.cpp b/c++/triqs_cthyb/moves/insert.cpp
index 9e2acf6e..d79cb04e 100644
--- a/c++/triqs_cthyb/moves/insert.cpp
+++ b/c++/triqs_cthyb/moves/insert.cpp
@@ -134,7 +134,7 @@ namespace triqs_cthyb {
       std::cerr << "Prefactor: " << t_ratio << '\t';
       std::cerr << "Weight: " << p * t_ratio << std::endl;
       std::cerr << "p_yee * newtrace: " << p_yee * new_atomic_weight << std::endl;
-      
+
       TRIQS_RUNTIME_ERROR << "(insert) p * t_ratio not finite p : " << p << " t_ratio : " << t_ratio << " in config " << config.get_id();
     }
     return p * t_ratio;
@@ -142,6 +142,11 @@ namespace triqs_cthyb {
 
   mc_weight_t move_insert_c_cdag::accept() {
 
+    time_pt tau_min = std::min(tau1,tau2);
+    time_pt tau_max = std::max(tau1,tau2);
+    if (tau_min < data.imp_trace.min_tau) data.imp_trace.min_tau = tau_min;
+    if (tau_max > data.imp_trace.max_tau) data.imp_trace.max_tau = tau_max;
+
     // insert in the tree
     data.imp_trace.confirm_insert();
 
diff --git a/c++/triqs_cthyb/moves/remove.cpp b/c++/triqs_cthyb/moves/remove.cpp
index 3c262d3d..91c60c24 100644
--- a/c++/triqs_cthyb/moves/remove.cpp
+++ b/c++/triqs_cthyb/moves/remove.cpp
@@ -109,7 +109,7 @@ namespace triqs_cthyb {
       std::cerr << "Weight: " << p / t_ratio << std::endl;
       TRIQS_RUNTIME_ERROR << "(remove) p not finite :" << p << " in config " << config.get_id();
     }
-    
+
     if (!isfinite(p / t_ratio)){
       TRIQS_RUNTIME_ERROR << "(remove) p / t_ratio not finite p : " << p << " t_ratio :  " << t_ratio << " in config " << config.get_id();
     }
@@ -118,6 +118,11 @@ namespace triqs_cthyb {
 
   mc_weight_t move_remove_c_cdag::accept() {
 
+    time_pt tau_min = std::min(tau1,tau2);
+    time_pt tau_max = std::max(tau1,tau2);
+    if (tau_min < data.imp_trace.min_tau) data.imp_trace.min_tau = tau_min;
+    if (tau_max > data.imp_trace.max_tau) data.imp_trace.max_tau = tau_max;
+
     // remove from the tree
     data.imp_trace.confirm_delete();
 
diff --git a/c++/triqs_cthyb/moves/shift.cpp b/c++/triqs_cthyb/moves/shift.cpp
index 81604e9d..f2768ec8 100644
--- a/c++/triqs_cthyb/moves/shift.cpp
+++ b/c++/triqs_cthyb/moves/shift.cpp
@@ -206,6 +206,11 @@ namespace triqs_cthyb {
 
   mc_weight_t move_shift_operator::accept() {
 
+    time_pt tau_min = std::min(tau_old,tau_new);
+    time_pt tau_max = std::max(tau_old,tau_new);
+    if (tau_min < data.imp_trace.min_tau) data.imp_trace.min_tau = tau_min;
+    if (tau_max > data.imp_trace.max_tau) data.imp_trace.max_tau = tau_max;
+
     // Update the tree
     data.imp_trace.confirm_shift();
 
diff --git a/c++/triqs_cthyb/parameters.cpp b/c++/triqs_cthyb/parameters.cpp
index 40b393a2..1ea39986 100644
--- a/c++/triqs_cthyb/parameters.cpp
+++ b/c++/triqs_cthyb/parameters.cpp
@@ -126,6 +126,7 @@ namespace triqs_cthyb {
 
     h5_write(grp, "measure_pert_order", sp.measure_pert_order);
     h5_write(grp, "measure_density_matrix", sp.measure_density_matrix);
+    h5_write(grp, "time_invariance", sp.time_invariance);
     h5_write(grp, "use_norm_as_weight", sp.use_norm_as_weight);
     h5_write(grp, "performance_analysis", sp.performance_analysis);
     h5_write(grp, "proposal_prob", sp.proposal_prob);
@@ -193,6 +194,7 @@ namespace triqs_cthyb {
 
     h5_read(grp, "measure_pert_order", sp.measure_pert_order);
     h5_read(grp, "measure_density_matrix", sp.measure_density_matrix);
+    h5_read(grp, "time_invariance", sp.time_invariance);
     h5_read(grp, "use_norm_as_weight", sp.use_norm_as_weight);
     h5_read(grp, "performance_analysis", sp.performance_analysis);
     h5_read(grp, "proposal_prob", sp.proposal_prob);
diff --git a/c++/triqs_cthyb/parameters.hpp b/c++/triqs_cthyb/parameters.hpp
index 65d0b006..a0007b78 100644
--- a/c++/triqs_cthyb/parameters.hpp
+++ b/c++/triqs_cthyb/parameters.hpp
@@ -190,6 +190,9 @@ namespace triqs_cthyb {
     /// Measure the reduced impurity density matrix?
     bool measure_density_matrix = false;
 
+    /// Use time invariance for the measurement of the density matrix?
+    bool time_invariance = false;
+
     /// Use the norm of the density matrix in the weight if true, otherwise use Trace
     bool use_norm_as_weight = false;
 
diff --git a/c++/triqs_cthyb/qmc_data.hpp b/c++/triqs_cthyb/qmc_data.hpp
index 12223822..2c3dec2c 100644
--- a/c++/triqs_cthyb/qmc_data.hpp
+++ b/c++/triqs_cthyb/qmc_data.hpp
@@ -73,7 +73,7 @@ namespace triqs_cthyb {
          tau_seg(beta),
          linindex(linindex),
          h_diag(h_diag),
-         imp_trace(beta, h_diag, histo_map, p.use_norm_as_weight, p.measure_density_matrix, p.performance_analysis),
+         imp_trace(beta, h_diag, histo_map, p.use_norm_as_weight, p.measure_density_matrix, p.time_invariance, p.performance_analysis),
          n_inner(n_inner),
          delta(map([](gf_const_view<imtime> d) { return real(d); }, delta)),
          current_sign(1),
diff --git a/c++/triqs_cthyb/solver_core.cpp b/c++/triqs_cthyb/solver_core.cpp
index e4d951f1..e4a0f43b 100644
--- a/c++/triqs_cthyb/solver_core.cpp
+++ b/c++/triqs_cthyb/solver_core.cpp
@@ -158,7 +158,7 @@ namespace triqs_cthyb {
 
       // ==== Compute h_loc ====
 
-      _h_loc0 = {}; 
+      _h_loc0 = {};
 
       // Add non-interacting terms to h_loc
       for (auto bl : range(gf_struct.size())) {
@@ -420,7 +420,7 @@ namespace triqs_cthyb {
     if (params.measure_density_matrix) {
       if (!params.use_norm_as_weight)
         TRIQS_RUNTIME_ERROR << "To measure the density_matrix of atomic states, you need to set "
-                               "use_norm_as_weight to True, i.e. to reweight the QMC";
+                              "use_norm_as_weight to True, i.e. to reweight the QMC";
       qmc.add_measure(measure_density_matrix{data, _density_matrix},
                       "Density Matrix for local static observable");
     }
diff --git a/c++/triqs_cthyb/solver_core.hpp b/c++/triqs_cthyb/solver_core.hpp
index 12995864..2fa9dc7f 100644
--- a/c++/triqs_cthyb/solver_core.hpp
+++ b/c++/triqs_cthyb/solver_core.hpp
@@ -3,7 +3,7 @@
  *
  * TRIQS: a Toolbox for Research in Interacting Quantum Systems
  *
- * Copyright (C) 2014-2017, H. U.R. Strand, P. Seth, I. Krivenko, 
+ * Copyright (C) 2014-2017, H. U.R. Strand, P. Seth, I. Krivenko,
  *                          M. Ferrero and O. Parcollet
  *
  * TRIQS is free software: you can redistribute it and/or modify it under the
@@ -62,7 +62,7 @@ namespace triqs_cthyb {
     // Return reference to container_set
     container_set_t &container_set() { return static_cast<container_set_t &>(*this); }
     container_set_t const &container_set() const { return static_cast<container_set_t const &>(*this); }
- 
+
     public:
 
     // Struct containing the parameters relevant for the solver construction
@@ -119,7 +119,7 @@ namespace triqs_cthyb {
     void set_container_set(container_set_t &cs) { static_cast<container_set_t &>(*this) = cs; }
     container_set_t last_container_set() { return static_cast<container_set_t>(*this); }
     */
-    
+
     /// :math:`\Delta(\tau)` in imaginary time.
     block_gf_view<imtime> Delta_tau() { return _Delta_tau; }
 
diff --git a/python/triqs_cthyb/config.py b/python/triqs_cthyb/config.py
index 77f6e077..b65e8f05 100644
--- a/python/triqs_cthyb/config.py
+++ b/python/triqs_cthyb/config.py
@@ -52,7 +52,7 @@ def plot_configs(hdf_file,beta,n_configs,delta_configs):
         conf = load_configuration(hdf_file,conf_idx)
         conf.plot(beta,conf_offset)
         conf_offset += 0.01
-        
+
 def hist_pert_order(hdf_file,n_configs):
     length=[]
     for i in range(1,n_configs+1):
diff --git a/python/triqs_cthyb/solver_core_desc.py b/python/triqs_cthyb/solver_core_desc.py
index 1f199ce6..bac66397 100644
--- a/python/triqs_cthyb/solver_core_desc.py
+++ b/python/triqs_cthyb/solver_core_desc.py
@@ -225,6 +225,8 @@
 +-------------------------------+----------------------------------------------------------+-------------------------------+-------------------------------------------------------------------------------------------------------------------+
 | measure_density_matrix        | bool                                                     | false                         | Measure the reduced impurity density matrix? Automatically also determines high frequency moments for G and Sigma |
 +-------------------------------+----------------------------------------------------------+-------------------------------+-------------------------------------------------------------------------------------------------------------------+
+| time_invariance               | bool                                                     | false                         | Use time invariance to sample the density matrix?                                                                 |
++-------------------------------+----------------------------------------------------------+-------------------------------+-------------------------------------------------------------------------------------------------------------------+
 | use_norm_as_weight            | bool                                                     | false                         | Use the norm of the density matrix in the weight if true, otherwise use Trace                                     |
 +-------------------------------+----------------------------------------------------------+-------------------------------+-------------------------------------------------------------------------------------------------------------------+
 | performance_analysis          | bool                                                     | false                         | Analyse performance of trace computation with histograms (developers only)?                                       |
@@ -532,6 +534,11 @@
              initializer = """ false """,
              doc = r"""Measure the reduced impurity density matrix?""")
 
+c.add_member(c_name = "time_invariance",
+             c_type = "bool",
+             initializer = """ false """,
+             doc = r"""Use time invariance to measure the density matrix?""")
+
 c.add_member(c_name = "use_norm_as_weight",
              c_type = "bool",
              initializer = """ false """,
diff --git a/python/triqs_cthyb/tail_fit.py b/python/triqs_cthyb/tail_fit.py
index 9dd37fee..5c8d68b9 100644
--- a/python/triqs_cthyb/tail_fit.py
+++ b/python/triqs_cthyb/tail_fit.py
@@ -33,8 +33,8 @@ def _comm(A,B): return A*B - B*A
 def _anticomm(A,B): return A*B + B*A
 
 def sigma_high_frequency_moments(density_matrix,
-                           ad_imp, 
-                           gf_struct, 
+                           ad_imp,
+                           gf_struct,
                            h_int):
     """
     Calculate the first and second high frequency moment of Sigma_iw
@@ -79,8 +79,8 @@ def sigma_high_frequency_moments(density_matrix,
 
 
 def green_high_frequency_moments(density_matrix,
-                           ad_imp, 
-                           gf_struct, 
+                           ad_imp,
+                           gf_struct,
                            h_imp):
     """
     Calculate the first and second high frequency moment of G_iw
@@ -99,7 +99,7 @@ def green_high_frequency_moments(density_matrix,
     gf_struct      : List of pairs (str,int)
                      Block structure of Green's function.
     h_imp          : triqs.operators.Operator
-                     impurity Hamiltonian   
+                     impurity Hamiltonian
 
     Returns
     -------
@@ -129,10 +129,10 @@ def tail_fit(
         fit_max_moment=None, fit_known_moments=None
         ):
     """
-    Fit a high frequency 1/(iw)^n expansion of Sigma_iw 
+    Fit a high frequency 1/(iw)^n expansion of Sigma_iw
     and replace the high frequency part with the fitted high frequency expansion.
 
-    Either give frequency window to fit on in terms of matsubara frequencies index 
+    Either give frequency window to fit on in terms of matsubara frequencies index
     (fit_min_n/fit_max_n) or value (fit_min_w/fit_max_w).
 
     Parameters
@@ -181,10 +181,10 @@ def tail_fit(
             n_max = fit_max_n,
             known_moments = fit_known_moments[name],
             # set max number of pts used in fit larger than mesh size, to use all data in fit
-            n_tail_max = 10 * len(sig.mesh), 
+            n_tail_max = 10 * len(sig.mesh),
             expansion_order = fit_max_moment
             )
-        
-        replace_by_tail(sig, tail, n_min=fit_min_n)        
+
+        replace_by_tail(sig, tail, n_min=fit_min_n)
 
     return Sigma_iw
diff --git a/python/triqs_cthyb/util.py b/python/triqs_cthyb/util.py
index 4442f1cd..741e0c57 100644
--- a/python/triqs_cthyb/util.py
+++ b/python/triqs_cthyb/util.py
@@ -47,7 +47,7 @@ def estimate_nfft_buf_size(gf_struct, pert_order_histograms):
     return buf_sizes
 
 def orbital_occupations(density_matrix, gf_struct, h_loc_diag):
-    
+
     dtype=density_matrix[0].dtype
     occ_mat = {bl: np.zeros((bl_size,bl_size), dtype=dtype) for bl, bl_size in gf_struct}
 
diff --git a/test/python/CMakeLists.txt b/test/python/CMakeLists.txt
index 7a75e387..5217342c 100644
--- a/test/python/CMakeLists.txt
+++ b/test/python/CMakeLists.txt
@@ -5,7 +5,7 @@ foreach(file ${all_h5_ref_files})
 endforeach()
 
 # List of all tests
-set(all_tests setup_Delta_tau_and_h_loc single_site_bethe atomic_observables kanamori_py slater measure_static histograms move_global h5_read_write h5_read_write_more O_tau_ins high_freq_tail hermitian_density_matrix)
+set(all_tests setup_Delta_tau_and_h_loc single_site_bethe atomic_observables kanamori_py slater measure_static measure_static_time_invariance histograms move_global h5_read_write h5_read_write_more O_tau_ins high_freq_tail hermitian_density_matrix)
 if(Local_hamiltonian_is_complex)
   list(APPEND all_tests atomic_gf_complex atomdiag_ed complex_bug81)
   if(Hybridisation_is_complex)
diff --git a/test/python/measure_static_time_invariance.py b/test/python/measure_static_time_invariance.py
new file mode 100644
index 00000000..277272d8
--- /dev/null
+++ b/test/python/measure_static_time_invariance.py
@@ -0,0 +1,75 @@
+import triqs.utility.mpi as mpi
+from h5 import HDFArchive
+from triqs.operators import *
+#from atom_diag import trace_rho_op
+from triqs.atom_diag import trace_rho_op
+from triqs_cthyb import *
+from triqs.gf import *
+import numpy as np
+
+# Input parameters
+beta = 10.0
+U = 2.0
+mu = 1.0
+h = 0.1
+norb = 2
+V = 1.0 * np.eye(norb) + 0.1 * (np.ones(norb) - np.eye(norb))
+t = 0.1
+epsilon = 2.3
+
+n_iw = 1025
+n_tau = 10001
+
+p = {}
+p["max_time"] = -1
+p["random_name"] = ""
+p["random_seed"] = 123 * mpi.rank + 567
+p["length_cycle"] = 50
+p["n_warmup_cycles"] = 1000
+p["n_cycles"] = 10000
+p["measure_G_tau"] = False
+p["use_norm_as_weight"] = True
+p["measure_density_matrix"] = True
+p["time_invariance"] = True
+
+gm = {}
+gm['flip_spins'] = {("up",0) : ("dn",0), ("dn",0) : ("up",0), ("up",1) : ("dn",1), ("dn",1) : ("up",1)}
+gm['swap_orbs']  = {("up",0) : ("up",1), ("up",1) : ("up",0), ("dn",0) : ("dn",1), ("dn",1) : ("dn",0)}
+p["move_global"] = gm
+p["move_global_prob"] = 0.06
+
+qn = [n("up",0) + n("up",1),n("dn",0) + n("dn",1)]
+p["quantum_numbers"] = qn
+p["partition_method"] = "quantum_numbers"
+
+H = U*n("up",0)*n("dn",0) + U*n("up",1)*n("dn",1)
+H = H + 0.5*h*(n("up",0) - n("dn",0)) + 0.5*h*(n("up",1) - n("dn",1))
+
+# Construct the solver
+S = Solver(beta=beta, gf_struct=[["dn",2], ["up",2]], n_tau=n_tau, n_iw=n_iw)
+
+# Set hybridization function
+delta_w = GfImFreq(beta=beta, target_shape=(2,2))
+delta_w << inverse(iOmega_n - epsilon) + inverse(iOmega_n + epsilon)
+delta_w.from_L_G_R(V, delta_w, V)
+
+S.G0_iw << inverse(iOmega_n + mu - delta_w)
+
+# Solve!
+S.solve(h_int=H, **p)
+
+if mpi.is_master_node():
+    # Measure expectation values
+    dm = S.density_matrix
+    static_observables = {"N1_up" : n("up",0), "N1_dn" : n("dn",0),
+                          "N2_up" : n("up",1), "N2_dn" : n("dn",1),
+                          "N12_up": c_dag("up",0) * c("up",1),
+                          "N12_dn": c_dag("dn",0) * c("dn",1)}
+    with HDFArchive('measure_static_time_invariance.out.h5','w') as ar:
+        for name,op in static_observables.items():
+            ave = trace_rho_op(dm,op,S.h_loc_diagonalization)
+            assert( np.abs(ave.imag) < 1e-10 )
+            ar[name] = ave.real
+
+from triqs.utility.h5diff import h5diff
+h5diff("measure_static_time_invariance.out.h5","measure_static_time_invariance.ref.h5")
diff --git a/test/python/measure_static_time_invariance.ref.h5 b/test/python/measure_static_time_invariance.ref.h5
new file mode 100644
index 00000000..f4a6651d
Binary files /dev/null and b/test/python/measure_static_time_invariance.ref.h5 differ