diff --git a/python/energy_flow_proper/11_new_ho_comparison/finite_t.py b/python/energy_flow_proper/11_new_ho_comparison/finite_t.py
index ce176be..5ab8823 100644
--- a/python/energy_flow_proper/11_new_ho_comparison/finite_t.py
+++ b/python/energy_flow_proper/11_new_ho_comparison/finite_t.py
@@ -89,7 +89,7 @@ def ho_zero(max_HO_level=15, wc=2, s=1, bcf_terms=7, k_max=5, bcf_scale=1, T=1):
 
     return params
 
-model_keys = [dict(bcf_scale=1, wc=ω_c) for ω_c in [1, 2]] + [dict(bcf_scale=s, wc=2) for s in [.5, 2]]
+model_keys = [dict(bcf_scale=1, wc=ω_c) for ω_c in [1, 3]] + [dict(bcf_scale=s, wc=2) for s in [2, 3]]
 
 model_keys = model_keys
 multi_params = [ho_zero(**keys) for keys in model_keys]
@@ -103,7 +103,7 @@ supervisors = []
 for params in multi_params:
       supervisor = HOPSSupervisor(
             params,
-            100_000,
+            150_000,
             data_path="ho_data",
             data_name="finite_t",
         )
@@ -133,10 +133,17 @@ for params, flow, keys in zip(multi_params, flow_hops, model_keys):
         -1 * flow.for_bath(0),
         ax=ax,
         linewidth=.8,
-        label=rf"$α(0)={params.SysP.g[0].sum().real:.2f}$, $ω_c={keys['wc']}$",
+        label=rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$",
     )
 ax.legend()
 
+for params, flow, ex_flow, keys, sup in zip(
+    multi_params, flow_hops, exact_flows, model_keys, supervisors[:1]
+):
+    system_energy = util.operator_expectation_from_data(
+        sup.get_data(True), params.SysP.H_sys, save="system_energy"
+    )
+
 from hopsflow import gaussflow as gf
 
 exact_flows = []
@@ -157,19 +164,16 @@ for params, keys in zip(multi_params, model_keys):
 
     α_0 = gf.BCF(
         params.IntP.t_max,
-        # factors=params.SysP.g[0],
-        # exponents=params.SysP.w[0],
-        hops.util.bcf.OhmicBCF_zeroTemp(
-            s=1,
-            eta=params.SysP.bcf_scale[0],
-            w_c=keys["wc"],
-            normed=True,
-        ),
-        # hops_bcf,
-        # num_terms=12,
+        factors=params.SysP.g[0]/params.SysP.bcf_scale[0],
+        exponents=params.SysP.w[0],
+        # hops.util.bcf.OhmicBCF_zeroTemp(
+        #     s=1,
+        #     eta=params.SysP.bcf_scale[0],
+        #     w_c=keys["wc"],
+        #     normed=True,
+        # ),
+        # num_terms=7,
         # resolution=0.01,
-        num_terms=7,
-        resolution=0.01,
     )
 
     α_0_dot = gf.BCF(
@@ -186,53 +190,119 @@ for params, keys in zip(multi_params, model_keys):
     flow_τ = flow(params.IntP.t)
     exact_flows.append(flow_τ)
 
-_tt = np.linspace(0, 4, 1000)
+_tt = np.linspace(0, 10, 1000)
 fig, ax = pu.plot_complex(_tt, (α(_tt)- α.approx(_tt)), label="finite", absolute=True)
-pu.plot_complex(_tt, α_0(_tt)- α_0.approx(_tt), label="finite", ax=ax, absolute=True)
+pu.plot_complex(
+    _tt,
+    hops.util.bcf.OhmicBCF_zeroTemp(
+        s=1,
+        eta=1,
+        normed=True,
+        w_c=keys["wc"],
+    )(_tt)
+    - α_0.approx(_tt),
+    label="zero",
+    ax=ax,
+    absolute=True,
+    linestyle="dashed"
+)
 plt.yscale("log")
-  # ut.plot_complex(result.τ, α(result.τ)- α.approx(result.τ), label="finite", ax=ax)
-  # ut.plot_complex(result.τ, α_0_dot(result.τ)- α_0_dot.approx(result.τ), label="finite", ax=ax)
+# ut.plot_complex(result.τ, α(result.τ)- α.approx(result.τ), label="finite", ax=ax)
+# ut.plot_complex(result.τ, α_0_dot(result.τ)- α_0_dot.approx(result.τ), label="finite", ax=ax)
 
-fig, ax = plt.subplots()
+tt = np.logspace(-1, 2, 1000)
+plt.plot(tt, np.abs(α(tt)), label="$T=1$")
+plt.plot(tt, np.abs(hops.util.bcf.OhmicBCF_zeroTemp(
+            s=1,
+            eta=1 / keys["wc"] ** (1 + 1),
+            w_c=keys["wc"],
+        )(tt)), label="$T=0$")
+plt.yscale("log")
+plt.xscale("log")
+plt.suptitle(rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$")
+plt.legend()
+plt.ylabel(r"$|\alpha(\tau)|$")
+plt.xlabel(r"$\tau$")
+plt.tight_layout()
+fs.export_fig("bcf_decay", y_scaling=.4)
+
+fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1, sharey=True)
 for params, flow, ex_flow, keys in zip(multi_params, flow_hops, exact_flows, model_keys):
-    consistency = (flow[-1].for_bath(0)).consistency(-ex_flow)
+    consistency = (-1 * flow.for_bath(0)).consistency(ex_flow)
     pu.plot_with_σ(
         params.IntP.t,
         -1 * flow,
         bath=0,
         ax=ax,
         linewidth=.8,
-        label=rf"$α(0)={params.SysP.g[0].sum().real:.2f}$ $ω_c={keys['wc']}$ | (${consistency}\%)$",
+        label=rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$, (${consistency:.1f}\%)$",
     )
-    ax.plot(params.IntP.t, ex_flow, linestyle="dotted", color="black", linewidth=.8)
+    ax.plot(params.IntP.t, ex_flow, color="black", linestyle="dashed", linewidth=.8)
 
-    plt.xlabel("$τ$")
-    plt.ylabel("$-J$")
-
-# ax.plot(multi_params[-1].IntP.t, flow_τ, label="Analytic", linestyle="dotted", color="black")
-ax.legend()
-fs.export_fig("comparison_finite")
-
-for params, flow, ex_flow, keys in zip(multi_params, flow_hops, exact_flows, model_keys):
-    fig, ax = plt.subplots()
-    consistency = (flow).consistency(-ex_flow)
-    pu.plot_convergence(
+    pu.plot_with_σ(
         params.IntP.t,
-        (flow.for_bath(0)[9::10]
-        +ex_flow) * (1/abs(ex_flow).max()),
-        reference=np.zeros_like(ex_flow),
-        ax=ax
-    )
-
-    ax.axhline(0, color="black")
-    ax.set_xlabel("$τ$")
-    ax.set_ylabel(r"$(J-J_{\mathrm{ref}})/J_\mathrm{max}$")
-    ax.legend()
-
-for params, flow, ex_flow, keys in zip(multi_params, flow_hops, exact_flows, model_keys):
-    fig, axs = pu.plot_consistency_development(
         -1 * flow,
-        ex_flow,
+        bath=0,
+        ax=ax2,
+        linewidth=.8,
+    )
+    ax2.plot(params.IntP.t, ex_flow, color="black", linestyle="dashed", linewidth=.8)
+
+
+ax.set_xlabel(r"$\tau$")
+ax.set_ylabel("$-J$")
+ax2.set_xlabel(r"$\tau$")
+ax2.set_ylabel("$-J$")
+fig.tight_layout()
+ax.legend()
+ax2.set_xlim((-.2, 10))
+fs.export_fig("flow_comp_nonzero", y_scaling=.4)
+
+for params, flow, ex_flow, keys, i in zip(
+    multi_params, flow_hops, exact_flows, model_keys, range(len(model_keys))
+):
+    subfig = fig.add_subfigure(grid[int(i / 2), i % 2])
+    f, ax = pu.plot_consistency_development(-1 * flow, ex_flow)
+
+    f.suptitle(
+        rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$"
     )
 
-    fig.suptitle(rf"$α(0)={params.SysP.g[0].sum().real:.2f}$ $ω_c={keys['wc']}$")
+    fs.export_fig(f"consistency_development_{i}", x_scaling=.48, y_scaling=.6, tikz=False)
+
+import hopsflow.util as util
+import hops.util.utilities as hops_utils
+
+fig, (ax, ax2) = plt.subplots(nrows=1, ncols=2)
+
+for params, flow, ex_flow, keys, sup in zip(
+    multi_params, flow_hops, exact_flows, model_keys, supervisors[:1]
+):
+    entropy, σ_entropy = hops_utils.entropy(sup.get_data(True))
+
+    system_energy = util.operator_expectation_from_data(
+        sup.get_data(True), params.SysP.H_sys, save="system_energy"
+    )
+
+    pu.plot_with_σ(
+        params.IntP.t,
+        system_energy,
+        ax=ax2,
+        linewidth=0.5,
+        label=rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$",
+    )
+
+    pu.plot_with_σ(
+        params.IntP.t,
+        util.EnsembleValue((entropy, σ_entropy)),
+        ax=ax,
+        linewidth=0.5,
+        label=rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$",
+    )
+
+ax.legend()
+ax.set_ylabel(r"$S(\rho)$")
+ax.set_xlabel(r"$\tau$")
+ax2.set_ylabel(r"$\langle H_\mathrm{S}\rangle$")
+ax2.set_xlabel(r"$\tau$")
+fs.export_fig("entropy_nonzero", fig, y_scaling=0.4)
diff --git a/python/energy_flow_proper/11_new_ho_comparison/ho_analytic_finite_temperature.org b/python/energy_flow_proper/11_new_ho_comparison/ho_analytic_finite_temperature.org
index 2877ca8..1baee22 100644
--- a/python/energy_flow_proper/11_new_ho_comparison/ho_analytic_finite_temperature.org
+++ b/python/energy_flow_proper/11_new_ho_comparison/ho_analytic_finite_temperature.org
@@ -94,7 +94,7 @@
 
       return params
 
-  model_keys = [dict(bcf_scale=1, wc=ω_c) for ω_c in [1, 2]] + [dict(bcf_scale=s, wc=2) for s in [.5, 2]]
+  model_keys = [dict(bcf_scale=1, wc=ω_c) for ω_c in [1, 3]] + [dict(bcf_scale=s, wc=2) for s in [2, 3]]
 
   model_keys = model_keys
   multi_params = [ho_zero(**keys) for keys in model_keys]
@@ -129,7 +129,7 @@
 #+end_src
 
 #+RESULTS:
-: RayContext(dashboard_url='', python_version='3.9.13', ray_version='1.13.0', ray_commit='e4ce38d001dbbe09cd21c497fedd03d692b2be3e', address_info={'node_ip_address': '141.30.17.225', 'raylet_ip_address': '141.30.17.225', 'redis_address': None, 'object_store_address': '/tmp/ray/session_2022-08-03_11-40-28_689081_82382/sockets/plasma_store', 'raylet_socket_name': '/tmp/ray/session_2022-08-03_11-40-28_689081_82382/sockets/raylet', 'webui_url': '', 'session_dir': '/tmp/ray/session_2022-08-03_11-40-28_689081_82382', 'metrics_export_port': 55919, 'gcs_address': '141.30.17.225:60072', 'address': '141.30.17.225:60072', 'node_id': '5167888344737f0b671c14167855fd8cac5e414985e3455d3d80962e'})
+: RayContext(dashboard_url='', python_version='3.9.13', ray_version='1.13.0', ray_commit='e4ce38d001dbbe09cd21c497fedd03d692b2be3e', address_info={'node_ip_address': '141.30.17.225', 'raylet_ip_address': '141.30.17.225', 'redis_address': None, 'object_store_address': '/tmp/ray/session_2022-08-05_11-07-29_852147_42277/sockets/plasma_store', 'raylet_socket_name': '/tmp/ray/session_2022-08-05_11-07-29_852147_42277/sockets/raylet', 'webui_url': '', 'session_dir': '/tmp/ray/session_2022-08-05_11-07-29_852147_42277', 'metrics_export_port': 64875, 'gcs_address': '141.30.17.225:60838', 'address': '141.30.17.225:60838', 'node_id': 'b15c3fa20dbefeb5f2e89cf630db49c5ed4f4a089c290f6414bd6a52'})
 
 
 #+begin_src jupyter-python
@@ -137,7 +137,7 @@
   for params in multi_params:
         supervisor = HOPSSupervisor(
               params,
-              100_000,
+              150_000,
               data_path="ho_data",
               data_name="finite_t",
           )
@@ -178,16 +178,24 @@
           -1 * flow.for_bath(0),
           ax=ax,
           linewidth=.8,
-          label=rf"$α(0)={params.SysP.g[0].sum().real:.2f}$, $ω_c={keys['wc']}$",
+          label=rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$",
       )
   ax.legend()
 #+end_src
 
 #+RESULTS:
 :RESULTS:
-: <matplotlib.legend.Legend at 0x7ff2dc5dd040>
-[[file:./.ob-jupyter/2c176c77b2e065a86242b68f6d12af3d080852b9.svg]]
+: <matplotlib.legend.Legend at 0x7f9d90647a60>
+[[file:./.ob-jupyter/95ec71cf26259aa1cd2837925e745d6d742fe7f6.svg]]
 :END:
+#+begin_src jupyter-python
+  for params, flow, ex_flow, keys, sup in zip(
+      multi_params, flow_hops, exact_flows, model_keys, supervisors[:1]
+  ):
+      system_energy = util.operator_expectation_from_data(
+          sup.get_data(True), params.SysP.H_sys, save="system_energy"
+      )
+#+end_src
 
 * Analytic
 #+begin_src jupyter-python :results none
@@ -213,19 +221,16 @@
 
       α_0 = gf.BCF(
           params.IntP.t_max,
-          # factors=params.SysP.g[0],
-          # exponents=params.SysP.w[0],
-          hops.util.bcf.OhmicBCF_zeroTemp(
-              s=1,
-              eta=params.SysP.bcf_scale[0],
-              w_c=keys["wc"],
-              normed=True,
-          ),
-          # hops_bcf,
-          # num_terms=12,
+          factors=params.SysP.g[0]/params.SysP.bcf_scale[0],
+          exponents=params.SysP.w[0],
+          # hops.util.bcf.OhmicBCF_zeroTemp(
+          #     s=1,
+          #     eta=params.SysP.bcf_scale[0],
+          #     w_c=keys["wc"],
+          #     normed=True,
+          # ),
+          # num_terms=7,
           # resolution=0.01,
-          num_terms=7,
-          resolution=0.01,
       )
 
       α_0_dot = gf.BCF(
@@ -244,83 +249,217 @@
 #+end_src
 
 #+begin_src jupyter-python
-  _tt = np.linspace(0, 4, 1000)
+  _tt = np.linspace(0, 10, 1000)
   fig, ax = pu.plot_complex(_tt, (α(_tt)- α.approx(_tt)), label="finite", absolute=True)
-  pu.plot_complex(_tt, α_0(_tt)- α_0.approx(_tt), label="finite", ax=ax, absolute=True)
+  pu.plot_complex(
+      _tt,
+      hops.util.bcf.OhmicBCF_zeroTemp(
+          s=1,
+          eta=1,
+          normed=True,
+          w_c=keys["wc"],
+      )(_tt)
+      - α_0.approx(_tt),
+      label="zero",
+      ax=ax,
+      absolute=True,
+      linestyle="dashed"
+  )
   plt.yscale("log")
-    # ut.plot_complex(result.τ, α(result.τ)- α.approx(result.τ), label="finite", ax=ax)
-    # ut.plot_complex(result.τ, α_0_dot(result.τ)- α_0_dot.approx(result.τ), label="finite", ax=ax)
+  # ut.plot_complex(result.τ, α(result.τ)- α.approx(result.τ), label="finite", ax=ax)
+  # ut.plot_complex(result.τ, α_0_dot(result.τ)- α_0_dot.approx(result.τ), label="finite", ax=ax)
 #+end_src
 
 #+RESULTS:
-[[file:./.ob-jupyter/2ca86b3de7d37e9b4fb413cf5c0d684e979fae18.svg]]
+[[file:./.ob-jupyter/27bfece24c66de5059f65e5c8a56a32987aada6b.svg]]
 
 #+begin_src jupyter-python
-  fig, ax = plt.subplots()
+  tt = np.logspace(-1, 2, 1000)
+  plt.plot(tt, np.abs(α(tt)), label="$T=1$")
+  plt.plot(tt, np.abs(hops.util.bcf.OhmicBCF_zeroTemp(
+              s=1,
+              eta=1 / keys["wc"] ** (1 + 1),
+              w_c=keys["wc"],
+          )(tt)), label="$T=0$")
+  plt.yscale("log")
+  plt.xscale("log")
+  plt.suptitle(rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$")
+  plt.legend()
+  plt.ylabel(r"$|\alpha(\tau)|$")
+  plt.xlabel(r"$\tau$")
+  plt.tight_layout()
+  fs.export_fig("bcf_decay", y_scaling=.4)
+#+end_src
+
+#+RESULTS:
+[[file:./.ob-jupyter/712a5f7581d046696cdffc6b0b8bc00be5959a20.svg]]
+
+
+#+begin_src jupyter-python
+  fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1, sharey=True)
   for params, flow, ex_flow, keys in zip(multi_params, flow_hops, exact_flows, model_keys):
-      consistency = (flow[-1].for_bath(0)).consistency(-ex_flow)
+      consistency = (-1 * flow.for_bath(0)).consistency(ex_flow)
       pu.plot_with_σ(
           params.IntP.t,
           -1 * flow,
           bath=0,
           ax=ax,
           linewidth=.8,
-          label=rf"$α(0)={params.SysP.g[0].sum().real:.2f}$ $ω_c={keys['wc']}$ | (${consistency}\%)$",
+          label=rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$, (${consistency:.1f}\%)$",
       )
-      ax.plot(params.IntP.t, ex_flow, linestyle="dotted", color="black", linewidth=.8)
+      ax.plot(params.IntP.t, ex_flow, color="black", linestyle="dashed", linewidth=.8)
 
-      plt.xlabel("$τ$")
-      plt.ylabel("$-J$")
-
-  # ax.plot(multi_params[-1].IntP.t, flow_τ, label="Analytic", linestyle="dotted", color="black")
-  ax.legend()
-  fs.export_fig("comparison_finite")
-#+end_src
-
-#+RESULTS:
-[[file:./.ob-jupyter/fff75a11b82859bb9354720b1fff2914e78c5809.svg]]
-
-#+begin_src jupyter-python
-  for params, flow, ex_flow, keys in zip(multi_params, flow_hops, exact_flows, model_keys):
-      fig, ax = plt.subplots()
-      consistency = (flow).consistency(-ex_flow)
-      pu.plot_convergence(
+      pu.plot_with_σ(
           params.IntP.t,
-          (flow.for_bath(0)[9::10]
-          +ex_flow) * (1/abs(ex_flow).max()),
-          reference=np.zeros_like(ex_flow),
-          ax=ax
+          -1 * flow,
+          bath=0,
+          ax=ax2,
+          linewidth=.8,
       )
+      ax2.plot(params.IntP.t, ex_flow, color="black", linestyle="dashed", linewidth=.8)
 
-      ax.axhline(0, color="black")
-      ax.set_xlabel("$τ$")
-      ax.set_ylabel(r"$(J-J_{\mathrm{ref}})/J_\mathrm{max}$")
-      ax.legend()
+
+  ax.set_xlabel(r"$\tau$")
+  ax.set_ylabel("$-J$")
+  ax2.set_xlabel(r"$\tau$")
+  ax2.set_ylabel("$-J$")
+  fig.tight_layout()
+  ax.legend()
+  ax2.set_xlim((-.2, 10))
+  fs.export_fig("flow_comp_nonzero", y_scaling=.4)
 #+end_src
 
 #+RESULTS:
-:RESULTS:
-[[file:./.ob-jupyter/8f7c7578463e60a2cc5cb55609c5e2350b31da18.svg]]
-[[file:./.ob-jupyter/4ac7a887a6d9b366edcb0f341c4badece1496baf.svg]]
-[[file:./.ob-jupyter/607fb3109397e4e84d6b771261fb7663416e4592.svg]]
-[[file:./.ob-jupyter/f79ccad4a11fade3d8c4ba214a0d0a89de26a2b8.svg]]
-:END:
+[[file:./.ob-jupyter/7448cfeaed2a113113be03e718a5cfb83f6fe64e.svg]]
+
 
 
 #+begin_src jupyter-python
-  for params, flow, ex_flow, keys in zip(multi_params, flow_hops, exact_flows, model_keys):
-      fig, axs = pu.plot_consistency_development(
-          -1 * flow,
-          ex_flow,
+  for params, flow, ex_flow, keys, i in zip(
+      multi_params, flow_hops, exact_flows, model_keys, range(len(model_keys))
+  ):
+      subfig = fig.add_subfigure(grid[int(i / 2), i % 2])
+      f, ax = pu.plot_consistency_development(-1 * flow, ex_flow)
+
+      f.suptitle(
+          rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$"
       )
 
-      fig.suptitle(rf"$α(0)={params.SysP.g[0].sum().real:.2f}$ $ω_c={keys['wc']}$")
+      fs.export_fig(f"consistency_development_{i}", x_scaling=.48, y_scaling=.6, tikz=False)
 #+end_src
 
 #+RESULTS:
 :RESULTS:
-[[file:./.ob-jupyter/4ed1f49960958709001508ae704fd85be1ee6510.svg]]
-[[file:./.ob-jupyter/ee332a2555dd751744bd435694ebf2d34b1e52ec.svg]]
-[[file:./.ob-jupyter/5860d96da073f680bb8196fae95791d766b82ba6.svg]]
-[[file:./.ob-jupyter/5b7b03d93715fd4534608e346cfca442ccb5b922.svg]]
+[[file:./.ob-jupyter/caff0b1a90552564d17910f32315989eb71934cd.svg]]
+[[file:./.ob-jupyter/8b8a13eadaca8a02418bcbbbcd9466465bd67b52.svg]]
+[[file:./.ob-jupyter/9f6fa17894b524c01e6506a459592d69b6026e24.svg]]
+[[file:./.ob-jupyter/b16c07d9f9fad57cd15543c684174c4c6142e3d1.svg]]
+:END:
+
+
+* Steady State Purity?
+#+begin_src jupyter-python
+  import hopsflow.util as util
+  import hops.util.utilities as hops_utils
+
+  fig, (ax, ax2) = plt.subplots(nrows=1, ncols=2)
+
+  for params, flow, ex_flow, keys, sup in zip(
+      multi_params, flow_hops, exact_flows, model_keys, supervisors[:1]
+  ):
+      entropy, σ_entropy = hops_utils.entropy(sup.get_data(True))
+
+      system_energy = util.operator_expectation_from_data(
+          sup.get_data(True), params.SysP.H_sys, save="system_energy"
+      )
+
+      pu.plot_with_σ(
+          params.IntP.t,
+          system_energy,
+          ax=ax2,
+          linewidth=0.5,
+          label=rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$",
+      )
+
+      pu.plot_with_σ(
+          params.IntP.t,
+          util.EnsembleValue((entropy, σ_entropy)),
+          ax=ax,
+          linewidth=0.5,
+          label=rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$",
+      )
+
+  ax.legend()
+  ax.set_ylabel(r"$S(\rho)$")
+  ax.set_xlabel(r"$\tau$")
+  ax2.set_ylabel(r"$\langle H_\mathrm{S}\rangle$")
+  ax2.set_xlabel(r"$\tau$")
+  fs.export_fig("entropy_nonzero", fig, y_scaling=0.4)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+: /nix/store/fgn7scpk3q7zcbva77zdph197h8ma5g9-python3-3.9.13-env/lib/python3.9/site-packages/scipy/linalg/_matfuncs_inv_ssq.py:827: LogmExactlySingularWarning: The logm input matrix is exactly singular.
+:   warnings.warn(exact_singularity_msg, LogmExactlySingularWarning)
+: /nix/store/fgn7scpk3q7zcbva77zdph197h8ma5g9-python3-3.9.13-env/lib/python3.9/site-packages/scipy/linalg/_matfuncs_inv_ssq.py:836: LogmNearlySingularWarning: The logm input matrix may be nearly singular.
+:   warnings.warn(near_singularity_msg, LogmNearlySingularWarning)
+: /home/hiro/src/hops/hops/util/utilities.py:43: RuntimeWarning: invalid value encountered in sqrt
+:   np.sqrt(np.trace(_square_mat(Δρ @ (logm2(ρ) + 1))).real)
+: Loading:   0% 40/16667 [00:15<1:53:52,  2.43it/s]
+# [goto error]
+#+begin_example
+  ---------------------------------------------------------------------------
+  KeyboardInterrupt                         Traceback (most recent call last)
+  Input In [65], in <cell line: 6>()
+        6 for params, flow, ex_flow, keys, sup in zip(
+        7     multi_params, flow_hops, exact_flows, model_keys, supervisors[:1]
+        8 ):
+        9     entropy, σ_entropy = hops_utils.entropy(sup.get_data(True))
+  ---> 11     system_energy = util.operator_expectation_from_data(
+       12         sup.get_data(True), params.SysP.H_sys, save="system_energy"
+       13     )
+       15     pu.plot_with_σ(
+       16         params.IntP.t,
+       17         system_energy,
+     (...)
+       20         label=rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$",
+       21     )
+       23     pu.plot_with_σ(
+       24         params.IntP.t,
+       25         util.EnsembleValue((entropy, σ_entropy)),
+     (...)
+       28         label=rf"$\alpha(0)={params.SysP.g[0].sum().real:.2f}$, $\omega_c={keys['wc']}$",
+       29     )
+
+  File ~/src/hopsflow/hopsflow/util.py:1031, in operator_expectation_from_data(data, op, **kwargs)
+     1028 if "save" in kwargs:
+     1029     kwargs["save"] += "_" + data.get_hi_key_hash()
+  -> 1031 return operator_expectation_ensemble(
+     1032     ψs=d.valid_sample_iterator(d.stoc_traj),
+     1033     op=op,
+     1034     t=d.get_time(),
+     1035     normalize=d.get_hi_key().HiP.nonlinear,
+     1036     **(dict(N=d.samples) | kwargs),
+     1037 )
+
+  File ~/src/hopsflow/hopsflow/util.py:542, in operator_expectation_ensemble(ψs, op, t, normalize, real, **kwargs)
+      538         sandwhiches = sandwhiches.real
+      540     return sandwhiches
+  --> 542 return ensemble_mean(ψs, op_exp_task, **kwargs)
+
+  File ~/src/hopsflow/hopsflow/util.py:827, in ensemble_mean(arg_iter, function, N, every, save, overwrite_cache, chunk_size, in_flight, gc_sleep)
+      817                 results.append(
+      818                     (
+      819                         aggregate.n,
+     (...)
+      822                     )
+      823                 )
+      825         highest_index += 1
+  --> 827     gc.collect()
+      828     time.sleep(gc_sleep)
+      830 if next_val:
+
+  KeyboardInterrupt:
+#+end_example
+[[file:./.ob-jupyter/c87c221267c7b2ed19d0590660a0140251c0437c.svg]]
 :END: