From 6339ca3d1896ef2538247398e34801d4fcbe4af9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20Meyer-M=C3=B6lleringhof?=
 <mixmax199989@gmail.com>
Date: Wed, 27 Nov 2024 11:55:38 +0900
Subject: [PATCH 1/4] First version of nm mcsolve notebook

- replicated plots from paper
- added introduction and explanations in own words
- added general structure for notbeook such as about, etc.
---
 .../time-evolution/025_v5_paper-nm_mcsolve.md | 358 ++++++++++++++++++
 1 file changed, 358 insertions(+)
 create mode 100644 tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md

diff --git a/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md b/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md
new file mode 100644
index 00000000..7ee46257
--- /dev/null
+++ b/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md
@@ -0,0 +1,358 @@
+---
+jupyter:
+  jupytext:
+    text_representation:
+      extension: .md
+      format_name: markdown
+      format_version: '1.3'
+      jupytext_version: 1.13.8
+  kernelspec:
+    display_name: Python 3
+    language: python
+    name: python3
+---
+
+# QuTiPv5 Paper Example: Monte Carlo Solver for non-Markovian Baths
+
+Authors: Maximilian Meyer-Mölleringhof (m.meyermoelleringhof@gmail.com), Neill Lambert (nwlambert@gmail.com)
+
+## Introduction
+
+When a quantum system experiences non-Markovian effects, it can no longer be described using the standard Lindblad formalism.
+To compute the time evolution of the density matrix, we can however use time-convolutionless (TCL) projection operators that lead to a differential equation in time-local form:
+
+$\dot{\rho} (t) = - \dfrac{i}{\hbar} [H(t), \rho(t)] + \sum_n \gamma_n(t) \mathcal{D}_n(t) [\rho(t)]$
+
+with
+
+$\mathcal{D}_n[\rho(t)] = A_n \rho(t) A^\dagger_n - \dfrac{1}{2} [A^\dagger_n A_n \rho(t) + \rho(t) A_n^\dagger A_n]$.
+
+These equations include the system Hamiltonian $H(t)$ and jump operators $A_n$.
+Contrary to a Lindblad equation, the coupling rates $\gamma_n(t)$ may be negative here.
+
+In QuTiPv5, the non-Markovian Monte Carlo solver is introduced.
+It enables the mapping of the general master equation given above, to a Lindblad equation on the same Hilbert space.
+This is achieved by the introduction of so called "influence martingale" which act as trajectory weighting [\[1, 2, 3\]](#References):
+
+$\mu (t) = \exp \left[ \alpha \int_0^t s(\tau) d \tau \right] \Pi_k \dfrac{\gamma_{n_k} (t_k)}{\Gamma_{n_k} (t_k)}$.
+
+Here, the product runs over all jump operators on the trajectory with jump channels $n_k$ and jump times $t_k < t$.
+
+The jump operators are required to fulfill the completeness relation $\sum_n A_n^\dagger A_n = \alpha \mathbb{1}$ for $\alpha > 0$.
+This condition is automatically taken care off by QuTiP's `nm_mcsolve()` function.
+
+To finally arrive at the Lindblad form, the shift function
+
+$s(t) = 2 | \min \{ 0, \gamma_1(t), \gamma_2(t), ... \} |$,
+
+such that the shifted rates $\Gamma_n (t) = \gamma_n(t) + s(t)$ non-negative.
+Using the regular MCWF method, we obtain the completely positive Lindblad equation
+
+$\dot{\rho}'(t) = - \dfrac{i}{\hbar} [ H(t), \rho'(t) ] + \sum_n \Gamma(t) \mathcal{D}_n[\rho'(t)]$,
+
+so that $\rho'(t) = \mathcal{E} \{\ket{\psi(t)} \bra{\psi(t)}\}$ for the generated trajecotries $\ket{\psi (t)}$.
+
+
+
+
+```python
+import matplotlib.pyplot as plt
+import numpy as np
+import qutip as qt
+from qutip import (about, basis, brmesolve, expect, lindblad_dissipator,
+                   liouvillian, mesolve, nm_mcsolve, qeye, sigmam, sigmap,
+                   sigmax, sigmay, sigmaz)
+from scipy.interpolate import CubicSpline
+from scipy.optimize import root_scalar
+```
+
+## Damped Jaynes-Cumming Model
+
+To illustrate the application of `nm_mcsolve()`, we want to look at the damped Jaynes-Cumming model.
+It describes a two-level atom coupled to a damped cavity mode.
+Such a model can be accurately model as a two-level system coupled to an environment with the power spectrum [\[4\]](#References):
+
+$S(\omega) = \dfrac{\lambda \Gamma^2}{(\omega_0 - \Delta - \omega)^2 + \Gamma^2}$,
+
+where $\lambda$ is the atom-cavity coupling strength, $\omega_0$ is the transition frequency of the atom, $\Delta$ the detuning of the cavity and $\Gamma$ the spectral width.
+At zero temperature and after performing the rotating wave approximation, the dynamics of the two-level atom can be described by the master equation
+
+$\dot{\rho}(t) = \dfrac{A(t)}{2i\hbar} [ \sigma_+ \sigma_-, \rho(t) ] + \gamma (t) \mathcal{D}_- [\rho (t)]$,
+
+with the state $\rho(t)$ in the interaction picture, the the ladder operators $\sigma_\pm$, the dissipator $\mathcal{D}_-$ for the Lindblad operator $\sigma_-$, and lastly $\gamma (t)$ and $A(t)$ being the real and imaginary parts of
+
+$\gamma(t) + i A(t) = \dfrac{2 \lambda \Gamma \sinh (\delta t / 2)}{\delta \cosh (\delta t / 2) + (\Gamma - i \Delta) \sinh (\delta t/2)}$,
+
+with $\delta = [(\Gamma - i \Delta)^2 - 2 \lambda \Gamma]^{1/2}$.
+
+```python
+H = sigmap() * sigmam() / 2
+initial_state = (basis(2, 0) + basis(2, 1)).unit()
+tlist = np.linspace(0, 5, 500)
+```
+
+```python
+# Constants
+gamma0 = 1
+lamb = 0.3 * gamma0
+Delta = 8 * lamb
+
+# Derived Quantities
+delta = np.sqrt((lamb - 1j * Delta) ** 2 - 2 * gamma0 * lamb)
+deltaR = np.real(delta)
+deltaI = np.imag(delta)
+deltaSq = deltaR**2 + deltaI**2
+```
+
+```python
+# calculate gamma and A
+def prefac(t):
+    return (
+        2
+        * gamma0
+        * lamb
+        / (
+            (lamb**2 + Delta**2 - deltaSq) * np.cos(deltaI * t)
+            - (lamb**2 + Delta**2 + deltaSq) * np.cosh(deltaR * t)
+            - 2 * (Delta * deltaR + lamb * deltaI) * np.sin(deltaI * t)
+            + 2 * (Delta * deltaI - lamb * deltaR) * np.sinh(deltaR * t)
+        )
+    )
+
+
+def cgamma(t):
+    return prefac(t) * (
+        lamb * np.cos(deltaI * t)
+        - lamb * np.cosh(deltaR * t)
+        - deltaI * np.sin(deltaI * t)
+        - deltaR * np.sinh(deltaR * t)
+    )
+
+
+def cA(t):
+    return prefac(t) * (
+        Delta * np.cos(deltaI * t)
+        - Delta * np.cosh(deltaR * t)
+        - deltaR * np.sin(deltaI * t)
+        + deltaI * np.sinh(deltaR * t)
+    )
+```
+
+```python
+_gamma = np.zeros_like(tlist)
+_A = np.zeros_like(tlist)
+for i in range(len(tlist)):
+    _gamma[i] = cgamma(tlist[i])
+    _A[i] = cA(tlist[i])
+
+gamma = CubicSpline(tlist, np.complex128(_gamma))
+A = CubicSpline(tlist, np.complex128(_A))
+```
+
+```python
+unitary_gen = liouvillian(H)
+dissipator = lindblad_dissipator(sigmam())
+me_sol = mesolve([[unitary_gen, A], [dissipator, gamma]], initial_state, tlist)
+```
+
+```python
+mc_sol = nm_mcsolve(
+    [[H, A]],
+    initial_state,
+    tlist,
+    ops_and_rates=[(sigmam(), gamma)],
+    ntraj=1_000,
+    options={"map": "parallel"},
+    seeds=0,
+)
+```
+
+## Comparison to other methods
+
+We want to compare this computation with the HEOM and the Bloch-Refield solver.
+For these methods we directly apply a spin-boson model and the free reservoir auto-correlation function
+
+$C(t) = \dfrac{\lambda \Gamma}{2} e^{-i (\omega - \Delta) t - \lambda |t|}$
+
+as well as the same power specture as we used above.
+We use the Hamiltonian $H = \omega_0 \sigma_+ \sigma_-$ in the Schrödinger picture and the coupling operator $Q = \sigma_+ + \sigma_-$.
+Here, we chose $\omega_0 \gg \Delta$ to ensure validity of the rotating wave approximation.
+
+```python
+omega_c = 100
+omega_0 = omega_c + Delta
+
+H = omega_0 * qt.sigmap() * qt.sigmam()
+Q = qt.sigmap() + qt.sigmam()
+
+
+def power_spectrum(w):
+    return gamma0 * lamb**2 / ((omega_c - w) ** 2 + lamb**2)
+```
+
+Firstly, the HEOM solver can directly be applied after separating the auto-correlation function into its real and imaginary parts:
+
+```python
+ck_real = [gamma0 * lamb / 4] * 2
+vk_real = [lamb - 1j * omega_c, lamb + 1j * omega_c]
+ck_imag = np.array([1j, -1j]) * gamma0 * lamb / 4
+vk_imag = vk_real
+```
+
+```python
+heom_bath = qt.heom.BosonicBath(Q, ck_real, vk_real, ck_imag, vk_imag)
+heom_sol = qt.heom.heomsolve(H, heom_bath, 10, qt.ket2dm(initial_state), tlist)
+```
+
+And secondly, for the Bloch-Redfield solver we can directly use the power spectrum as input:
+
+```python
+br_sol = brmesolve(H, initial_state, tlist, a_ops=[(sigmax(), power_spectrum)])
+```
+
+Finally, in order to compare these results with the `nm_mesolve()` method, we transform them to the interaction picture:
+
+```python
+Us = [(-1j * H * t).expm() for t in tlist]
+heom_states = [U * s * U.dag() for (U, s) in zip(Us, heom_sol.states)]
+br_states = [U * s * U.dag() for (U, s) in zip(Us, br_sol.states)]
+```
+
+### Plotting the Time Evolution
+
+We choose to plot three comparisons to highlight the different solvers' computational results.
+In all plots, the gray areas show when $\gamma(t)$ is negative.
+
+```python
+root1 = root_scalar(lambda t: cgamma(t), method="bisect", bracket=(1, 2)).root
+root2 = root_scalar(lambda t: cgamma(t), method="bisect", bracket=(2, 3)).root
+root3 = root_scalar(lambda t: cgamma(t), method="bisect", bracket=(3, 4)).root
+root4 = root_scalar(lambda t: cgamma(t), method="bisect", bracket=(4, 5)).root
+```
+
+```python
+projector = (sigmaz() + qeye(2)) / 2
+
+rho11_me = expect(projector, me_sol.states)
+rho11_mc = expect(projector, mc_sol.states[::10])
+rho11_br = expect(projector, br_sol.states)
+rho11_heom = expect(projector, heom_states)
+
+plt.plot(tlist, rho11_me, "-", color="orange", label="mesolve")
+plt.plot(
+    tlist[::10],
+    rho11_mc,
+    "x",
+    color="blue",
+    label="nm_mcsolve",
+)
+plt.plot(tlist, rho11_br, "-.", color="gray", label="brmesolve")
+plt.plot(tlist, rho11_heom, "--", color="green", label="heomsolve")
+
+plt.xlabel(r"$t\, /\, \lambda^{-1}$")
+plt.xlim((-0.2, 5.2))
+plt.xticks([0, 2.5, 5], labels=["0", "2.5", "5"])
+plt.title(r"$\rho_{11}$")
+plt.ylim((0.4376, 0.5024))
+plt.yticks([0.44, 0.46, 0.48, 0.5], labels=["0.44", "0.46", "0.48", "0.50"])
+
+plt.axvspan(root1, root2, color="gray", alpha=0.08, zorder=0)
+plt.axvspan(root3, root4, color="gray", alpha=0.08, zorder=0)
+
+plt.legend()
+plt.show()
+```
+
+```python
+me_x = expect(sigmax(), me_sol.states)
+mc_x = expect(sigmax(), mc_sol.states[::10])
+heom_x = expect(sigmax(), heom_states)
+br_x = expect(sigmax(), br_states)
+
+me_y = expect(sigmay(), me_sol.states)
+mc_y = expect(sigmay(), mc_sol.states[::10])
+heom_y = expect(sigmay(), heom_states)
+br_y = expect(sigmay(), br_states)
+```
+
+```python
+# We smooth the HEOM result because it oscillates quickly and gets hard to see
+heom_plot = heom_x * heom_x + heom_y * heom_y
+heom_plot = np.convolve(heom_plot, np.array([1 / 11] * 11), mode="valid")
+heom_tlist = tlist[5:-5]
+```
+
+```python
+rho01_me = me_x * me_x + me_y * me_y
+rho01_mc = mc_x * mc_x + mc_y * mc_y
+rho01_br = br_x * br_x + br_y * br_y
+
+plt.plot(tlist, rho01_me, "-", color="orange", label=r"mesolve")
+plt.plot(tlist[::10], rho01_mc, "x", color="blue", label=r"nm_mcsolve")
+plt.plot(heom_tlist, heom_plot, "--", color="green", label=r"heomsolve")
+plt.plot(tlist, rho01_br, "-.", color="gray", label=r"brmesolve")
+
+plt.xlabel(r"$t\, /\, \lambda^{-1}$")
+plt.xlim((-0.2, 5.2))
+plt.xticks([0, 2.5, 5], labels=["0", "2.5", "5"])
+plt.title(r"$| \rho_{01} |^2$")
+plt.ylim((0.8752, 1.0048))
+plt.yticks(
+    [0.88, 0.9, 0.92, 0.94, 0.96, 0.98, 1],
+    labels=["0.88", "0.90", "0.92", "0.94", "0.96", "0.98", "1"],
+)
+
+plt.axvspan(root1, root2, color="gray", alpha=0.08, zorder=0)
+plt.axvspan(root3, root4, color="gray", alpha=0.08, zorder=0)
+
+plt.legend()
+plt.show()
+```
+
+```python
+plt.plot(tlist, np.zeros_like(tlist), "-", color="orange", label=r"Zero")
+plt.plot(
+    tlist[::10],
+    1000 * (mc_sol.trace[::10] - 1),
+    "x",
+    color="blue",
+    label=r"nm_mcsolve",
+)
+
+plt.xlabel(r"$t\, /\, \lambda^{-1}$")
+plt.xlim((-0.2, 5.2))
+plt.xticks([0, 2.5, 5], labels=["0", "2.5", "5"])
+plt.title(r"$(\mu - 1)\, /\, 10^{-3}$")
+plt.ylim((-5.8, 15.8))
+plt.yticks([-5, 0, 5, 10, 15])
+
+plt.axvspan(root1, root2, color="gray", alpha=0.08, zorder=0)
+plt.axvspan(root3, root4, color="gray", alpha=0.08, zorder=0)
+
+plt.legend()
+plt.show()
+```
+
+## References
+
+\[1\] [Donvil and Muratore-Ginanneschi. Nat Commun (2022).](https://www.nature.com/articles/s41467-022-31533-8)
+
+\[2\] [Donvil and Muratore-Ginanneschi. New J. Phys. (2023).](https://dx.doi.org/10.1088/1367-2630/acd4dc)
+
+\[3\] [Donvil and Muratore-Ginanneschi. *Open Systems & Information Dynamics*.](https://www.worldscientific.com/worldscinet/osid)
+
+\[4\] [Breuer and Petruccione *The Theory of Open Quantum Systems*.](https://doi.org/10.1093/acprof:oso/9780199213900.001.0001)
+
+
+
+## About
+
+```python
+about()
+```
+
+## Testing
+
+```python
+# TODO
+```

From cb35d3cb3ff305626c84f226e5213b18a350978a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20Meyer-M=C3=B6lleringhof?=
 <mixmax199989@gmail.com>
Date: Thu, 28 Nov 2024 15:10:31 +0900
Subject: [PATCH 2/4] Added tests, refined instructions

- added tests relying on proximity of the solutions
- rewrote instructions so they are less like the paper and more clear
- fixed typos and improved code structure
---
 .../time-evolution/025_v5_paper-nm_mcsolve.md | 69 +++++++++++++------
 1 file changed, 49 insertions(+), 20 deletions(-)

diff --git a/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md b/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md
index 7ee46257..598f7e99 100644
--- a/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md
+++ b/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md
@@ -45,15 +45,13 @@ To finally arrive at the Lindblad form, the shift function
 
 $s(t) = 2 | \min \{ 0, \gamma_1(t), \gamma_2(t), ... \} |$,
 
-such that the shifted rates $\Gamma_n (t) = \gamma_n(t) + s(t)$ non-negative.
+such that the shifted rates $\Gamma_n (t) = \gamma_n(t) + s(t)$ are non-negative.
 Using the regular MCWF method, we obtain the completely positive Lindblad equation
 
 $\dot{\rho}'(t) = - \dfrac{i}{\hbar} [ H(t), \rho'(t) ] + \sum_n \Gamma(t) \mathcal{D}_n[\rho'(t)]$,
 
-so that $\rho'(t) = \mathcal{E} \{\ket{\psi(t)} \bra{\psi(t)}\}$ for the generated trajecotries $\ket{\psi (t)}$.
-
-
-
+so that $\rho'(t) = \mathbb{E} \{\ket{\psi(t)} \bra{\psi(t)}\}$ for the generated trajecotries $\ket{\psi (t)}$, where $\mathbb{E}$ is the ensemble average over the trajectories.
+This can furthermore be used to finally reconstruct the original states via $\rho(t) = \mathbb{E}\{\mu(t) \ket{\psi(t)} \bra{\psi(t)}\}$.
 
 ```python
 import matplotlib.pyplot as plt
@@ -152,7 +150,6 @@ A = CubicSpline(tlist, np.complex128(_A))
 ```python
 unitary_gen = liouvillian(H)
 dissipator = lindblad_dissipator(sigmam())
-me_sol = mesolve([[unitary_gen, A], [dissipator, gamma]], initial_state, tlist)
 ```
 
 ```python
@@ -169,12 +166,18 @@ mc_sol = nm_mcsolve(
 
 ## Comparison to other methods
 
-We want to compare this computation with the HEOM and the Bloch-Refield solver.
-For these methods we directly apply a spin-boson model and the free reservoir auto-correlation function
+We want to compare this computation with the standard `mesolve()`, the HEOM and the Bloch-Refield solver.
+For the `mesolve()` we can directly use the operators we have created for `nm_mcsolve()`:
+
+```python
+me_sol = mesolve([[unitary_gen, A], [dissipator, gamma]], initial_state, tlist)
+```
+
+For the other methods we directly apply a spin-boson model and the free reservoir auto-correlation function
 
 $C(t) = \dfrac{\lambda \Gamma}{2} e^{-i (\omega - \Delta) t - \lambda |t|}$
 
-as well as the same power specture as we used above.
+as well as the same power specturm as we used above.
 We use the Hamiltonian $H = \omega_0 \sigma_+ \sigma_-$ in the Schrödinger picture and the coupling operator $Q = \sigma_+ + \sigma_-$.
 Here, we chose $\omega_0 \gg \Delta$ to ensure validity of the rotating wave approximation.
 
@@ -204,7 +207,7 @@ heom_bath = qt.heom.BosonicBath(Q, ck_real, vk_real, ck_imag, vk_imag)
 heom_sol = qt.heom.heomsolve(H, heom_bath, 10, qt.ket2dm(initial_state), tlist)
 ```
 
-And secondly, for the Bloch-Redfield solver we can directly use the power spectrum as input:
+Secondly, for the Bloch-Redfield solver we can directly use the power spectrum as input:
 
 ```python
 br_sol = brmesolve(H, initial_state, tlist, a_ops=[(sigmax(), power_spectrum)])
@@ -234,14 +237,14 @@ root4 = root_scalar(lambda t: cgamma(t), method="bisect", bracket=(4, 5)).root
 projector = (sigmaz() + qeye(2)) / 2
 
 rho11_me = expect(projector, me_sol.states)
-rho11_mc = expect(projector, mc_sol.states[::10])
+rho11_mc = expect(projector, mc_sol.states)
 rho11_br = expect(projector, br_sol.states)
 rho11_heom = expect(projector, heom_states)
 
 plt.plot(tlist, rho11_me, "-", color="orange", label="mesolve")
 plt.plot(
     tlist[::10],
-    rho11_mc,
+    rho11_mc[::10],
     "x",
     color="blue",
     label="nm_mcsolve",
@@ -265,20 +268,20 @@ plt.show()
 
 ```python
 me_x = expect(sigmax(), me_sol.states)
-mc_x = expect(sigmax(), mc_sol.states[::10])
+mc_x = expect(sigmax(), mc_sol.states)
 heom_x = expect(sigmax(), heom_states)
 br_x = expect(sigmax(), br_states)
 
 me_y = expect(sigmay(), me_sol.states)
-mc_y = expect(sigmay(), mc_sol.states[::10])
+mc_y = expect(sigmay(), mc_sol.states)
 heom_y = expect(sigmay(), heom_states)
 br_y = expect(sigmay(), br_states)
 ```
 
 ```python
 # We smooth the HEOM result because it oscillates quickly and gets hard to see
-heom_plot = heom_x * heom_x + heom_y * heom_y
-heom_plot = np.convolve(heom_plot, np.array([1 / 11] * 11), mode="valid")
+rho01_heom = heom_x * heom_x + heom_y * heom_y
+rho01_heom = np.convolve(rho01_heom, np.array([1 / 11] * 11), mode="valid")
 heom_tlist = tlist[5:-5]
 ```
 
@@ -288,8 +291,8 @@ rho01_mc = mc_x * mc_x + mc_y * mc_y
 rho01_br = br_x * br_x + br_y * br_y
 
 plt.plot(tlist, rho01_me, "-", color="orange", label=r"mesolve")
-plt.plot(tlist[::10], rho01_mc, "x", color="blue", label=r"nm_mcsolve")
-plt.plot(heom_tlist, heom_plot, "--", color="green", label=r"heomsolve")
+plt.plot(tlist[::10], rho01_mc[::10], "x", color="blue", label=r"nm_mcsolve")
+plt.plot(heom_tlist, rho01_heom, "--", color="green", label=r"heomsolve")
 plt.plot(tlist, rho01_br, "-.", color="gray", label=r"brmesolve")
 
 plt.xlabel(r"$t\, /\, \lambda^{-1}$")
@@ -309,11 +312,15 @@ plt.legend()
 plt.show()
 ```
 
+```python
+mart_dev = mc_sol.trace - 1
+```
+
 ```python
 plt.plot(tlist, np.zeros_like(tlist), "-", color="orange", label=r"Zero")
 plt.plot(
     tlist[::10],
-    1000 * (mc_sol.trace[::10] - 1),
+    1000 * mart_dev[::10],
     "x",
     color="blue",
     label=r"nm_mcsolve",
@@ -333,6 +340,14 @@ plt.legend()
 plt.show()
 ```
 
+In these plots we notice two things.
+First, the result from the Bloch-Redfield equation deviates greatly from the others, showing us that non-Markovian effects are strongly influencing the dynamics.
+Second, in the grey areas - when $\gamma(t)$ is negative - the atom state restores coherence.
+In the last plot especially we see that during these times, the average influence martingale fluctuates, but is constant otherwise.
+Its deviation from unity tells us how well the simulation has converged.
+
+
+
 ## References
 
 \[1\] [Donvil and Muratore-Ginanneschi. Nat Commun (2022).](https://www.nature.com/articles/s41467-022-31533-8)
@@ -354,5 +369,19 @@ about()
 ## Testing
 
 ```python
-# TODO
+assert np.allclose(
+    rho11_me, rho11_heom, atol=1e-3
+), "rho11 of mesolve and heomsolve do not agree."
+assert np.allclose(
+    rho11_me, rho11_mc, atol=1e-2
+), "rho11 of nm_mcsolve deviates from mesolve too much."
+assert np.allclose(
+    rho01_me[5:-5], rho01_heom, atol=1e-3
+), "|rho01|^2 of mesolve and heomsolve do not agree."
+assert np.allclose(
+    rho01_me, rho01_mc, atol=1e-1
+), "|rho01|^2 of nm_mcsolve deviates from mesolve too much."
+assert (
+    np.max(mart_dev) < 1e-1
+), "MC Simulation has not converged well enough. Average infl. mart. > 1e-1"
 ```

From 78279c16cfe3305c2c9792d81289cc4399b04d12 Mon Sep 17 00:00:00 2001
From: Paul Menczel <paul@menczel.net>
Date: Mon, 2 Dec 2024 16:35:50 +0900
Subject: [PATCH 3/4] Minor text improvements, use environment API

---
 .../time-evolution/025_v5_paper-nm_mcsolve.md | 42 ++++++++-----------
 1 file changed, 18 insertions(+), 24 deletions(-)

diff --git a/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md b/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md
index 598f7e99..5c016f94 100644
--- a/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md
+++ b/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md
@@ -5,16 +5,16 @@ jupyter:
       extension: .md
       format_name: markdown
       format_version: '1.3'
-      jupytext_version: 1.13.8
+      jupytext_version: 1.16.4
   kernelspec:
-    display_name: Python 3
+    display_name: Python 3 (ipykernel)
     language: python
     name: python3
 ---
 
 # QuTiPv5 Paper Example: Monte Carlo Solver for non-Markovian Baths
 
-Authors: Maximilian Meyer-Mölleringhof (m.meyermoelleringhof@gmail.com), Neill Lambert (nwlambert@gmail.com)
+Authors: Maximilian Meyer-Mölleringhof (m.meyermoelleringhof@gmail.com), Paul Menczel (paul@menczel.net), Neill Lambert (nwlambert@gmail.com)
 
 ## Introduction
 
@@ -30,29 +30,29 @@ $\mathcal{D}_n[\rho(t)] = A_n \rho(t) A^\dagger_n - \dfrac{1}{2} [A^\dagger_n A_
 These equations include the system Hamiltonian $H(t)$ and jump operators $A_n$.
 Contrary to a Lindblad equation, the coupling rates $\gamma_n(t)$ may be negative here.
 
-In QuTiPv5, the non-Markovian Monte Carlo solver is introduced.
+In QuTiP v5, the non-Markovian Monte Carlo solver is introduced.
 It enables the mapping of the general master equation given above, to a Lindblad equation on the same Hilbert space.
-This is achieved by the introduction of so called "influence martingale" which act as trajectory weighting [\[1, 2, 3\]](#References):
+This is achieved by the introduction of the so called "influence martingale" which acts as trajectory weighting [\[1, 2, 3\]](#References):
 
 $\mu (t) = \exp \left[ \alpha \int_0^t s(\tau) d \tau \right] \Pi_k \dfrac{\gamma_{n_k} (t_k)}{\Gamma_{n_k} (t_k)}$.
 
 Here, the product runs over all jump operators on the trajectory with jump channels $n_k$ and jump times $t_k < t$.
-
-The jump operators are required to fulfill the completeness relation $\sum_n A_n^\dagger A_n = \alpha \mathbb{1}$ for $\alpha > 0$.
-This condition is automatically taken care off by QuTiP's `nm_mcsolve()` function.
-
 To finally arrive at the Lindblad form, the shift function
 
-$s(t) = 2 | \min \{ 0, \gamma_1(t), \gamma_2(t), ... \} |$,
+$s(t) = 2 \left| \min \{ 0, \gamma_1(t), \gamma_2(t), ... \} \right|$
 
-such that the shifted rates $\Gamma_n (t) = \gamma_n(t) + s(t)$ are non-negative.
-Using the regular MCWF method, we obtain the completely positive Lindblad equation
+is applied, such that the shifted rates $\Gamma_n (t) = \gamma_n(t) + s(t)$ are non-negative.
+We obtain the completely positive Lindblad equation
 
 $\dot{\rho}'(t) = - \dfrac{i}{\hbar} [ H(t), \rho'(t) ] + \sum_n \Gamma(t) \mathcal{D}_n[\rho'(t)]$,
 
-so that $\rho'(t) = \mathbb{E} \{\ket{\psi(t)} \bra{\psi(t)}\}$ for the generated trajecotries $\ket{\psi (t)}$, where $\mathbb{E}$ is the ensemble average over the trajectories.
+and $\rho'(t) = \mathbb{E} \{\ket{\psi(t)} \bra{\psi(t)}\}$ using the regular MCWF method.
+Here, $\ket{\psi (t)}$ are the generated trajectories and $\mathbb{E}$ is the ensemble average over the trajectories.
 This can furthermore be used to finally reconstruct the original states via $\rho(t) = \mathbb{E}\{\mu(t) \ket{\psi(t)} \bra{\psi(t)}\}$.
 
+Note that, for this technique, the jump operators are required to fulfill the completeness relation $\sum_n A_n^\dagger A_n = \alpha \mathbb{1}$ for $\alpha > 0$.
+This condition is automatically taken care of by QuTiP's `nm_mcsolve()` function.
+
 ```python
 import matplotlib.pyplot as plt
 import numpy as np
@@ -106,10 +106,7 @@ deltaSq = deltaR**2 + deltaI**2
 # calculate gamma and A
 def prefac(t):
     return (
-        2
-        * gamma0
-        * lamb
-        / (
+        2 * gamma0 * lamb / (
             (lamb**2 + Delta**2 - deltaSq) * np.cos(deltaI * t)
             - (lamb**2 + Delta**2 + deltaSq) * np.cosh(deltaR * t)
             - 2 * (Delta * deltaR + lamb * deltaI) * np.sin(deltaI * t)
@@ -117,7 +114,6 @@ def prefac(t):
         )
     )
 
-
 def cgamma(t):
     return prefac(t) * (
         lamb * np.cos(deltaI * t)
@@ -126,7 +122,6 @@ def cgamma(t):
         - deltaR * np.sinh(deltaR * t)
     )
 
-
 def cA(t):
     return prefac(t) * (
         Delta * np.cos(deltaI * t)
@@ -175,9 +170,9 @@ me_sol = mesolve([[unitary_gen, A], [dissipator, gamma]], initial_state, tlist)
 
 For the other methods we directly apply a spin-boson model and the free reservoir auto-correlation function
 
-$C(t) = \dfrac{\lambda \Gamma}{2} e^{-i (\omega - \Delta) t - \lambda |t|}$
+$C(t) = \dfrac{\lambda \Gamma}{2} e^{-i (\omega - \Delta) t - \lambda |t|}$,
 
-as well as the same power specturm as we used above.
+which corresponds to the power spectrum defined above.
 We use the Hamiltonian $H = \omega_0 \sigma_+ \sigma_-$ in the Schrödinger picture and the coupling operator $Q = \sigma_+ + \sigma_-$.
 Here, we chose $\omega_0 \gg \Delta$ to ensure validity of the rotating wave approximation.
 
@@ -188,7 +183,6 @@ omega_0 = omega_c + Delta
 H = omega_0 * qt.sigmap() * qt.sigmam()
 Q = qt.sigmap() + qt.sigmam()
 
-
 def power_spectrum(w):
     return gamma0 * lamb**2 / ((omega_c - w) ** 2 + lamb**2)
 ```
@@ -203,8 +197,8 @@ vk_imag = vk_real
 ```
 
 ```python
-heom_bath = qt.heom.BosonicBath(Q, ck_real, vk_real, ck_imag, vk_imag)
-heom_sol = qt.heom.heomsolve(H, heom_bath, 10, qt.ket2dm(initial_state), tlist)
+heom_env = qt.ExponentialBosonicEnvironment(ck_real, vk_real, ck_imag, vk_imag)
+heom_sol = qt.heom.heomsolve(H, (heom_env, Q), 10, qt.ket2dm(initial_state), tlist)
 ```
 
 Secondly, for the Bloch-Redfield solver we can directly use the power spectrum as input:

From 5f7d28ee6d4a9e6499809f9eff431c124d159a77 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20Meyer-M=C3=B6lleringhof?=
 <mixmax199989@gmail.com>
Date: Mon, 9 Dec 2024 14:05:11 +0900
Subject: [PATCH 4/4] Added v5 paper reference, fixed formatting

---
 .../time-evolution/025_v5_paper-nm_mcsolve.md | 30 ++++++++++++-------
 1 file changed, 19 insertions(+), 11 deletions(-)

diff --git a/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md b/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md
index 5c016f94..93d23f69 100644
--- a/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md
+++ b/tutorials-v5/time-evolution/025_v5_paper-nm_mcsolve.md
@@ -5,9 +5,9 @@ jupyter:
       extension: .md
       format_name: markdown
       format_version: '1.3'
-      jupytext_version: 1.16.4
+      jupytext_version: 1.13.8
   kernelspec:
-    display_name: Python 3 (ipykernel)
+    display_name: qutip-tutorials-v5
     language: python
     name: python3
 ---
@@ -56,10 +56,10 @@ This condition is automatically taken care of by QuTiP's `nm_mcsolve()` function
 ```python
 import matplotlib.pyplot as plt
 import numpy as np
-import qutip as qt
-from qutip import (about, basis, brmesolve, expect, lindblad_dissipator,
-                   liouvillian, mesolve, nm_mcsolve, qeye, sigmam, sigmap,
-                   sigmax, sigmay, sigmaz)
+from qutip import (ExponentialBosonicEnvironment, about, basis, brmesolve,
+                   expect, heom, ket2dm, lindblad_dissipator, liouvillian,
+                   mesolve, nm_mcsolve, qeye, sigmam, sigmap, sigmax, sigmay,
+                   sigmaz)
 from scipy.interpolate import CubicSpline
 from scipy.optimize import root_scalar
 ```
@@ -106,7 +106,10 @@ deltaSq = deltaR**2 + deltaI**2
 # calculate gamma and A
 def prefac(t):
     return (
-        2 * gamma0 * lamb / (
+        2
+        * gamma0
+        * lamb
+        / (
             (lamb**2 + Delta**2 - deltaSq) * np.cos(deltaI * t)
             - (lamb**2 + Delta**2 + deltaSq) * np.cosh(deltaR * t)
             - 2 * (Delta * deltaR + lamb * deltaI) * np.sin(deltaI * t)
@@ -114,6 +117,7 @@ def prefac(t):
         )
     )
 
+
 def cgamma(t):
     return prefac(t) * (
         lamb * np.cos(deltaI * t)
@@ -122,6 +126,7 @@ def cgamma(t):
         - deltaR * np.sinh(deltaR * t)
     )
 
+
 def cA(t):
     return prefac(t) * (
         Delta * np.cos(deltaI * t)
@@ -180,8 +185,9 @@ Here, we chose $\omega_0 \gg \Delta$ to ensure validity of the rotating wave app
 omega_c = 100
 omega_0 = omega_c + Delta
 
-H = omega_0 * qt.sigmap() * qt.sigmam()
-Q = qt.sigmap() + qt.sigmam()
+H = omega_0 * sigmap() * sigmam()
+Q = sigmap() + sigmam()
+
 
 def power_spectrum(w):
     return gamma0 * lamb**2 / ((omega_c - w) ** 2 + lamb**2)
@@ -197,8 +203,8 @@ vk_imag = vk_real
 ```
 
 ```python
-heom_env = qt.ExponentialBosonicEnvironment(ck_real, vk_real, ck_imag, vk_imag)
-heom_sol = qt.heom.heomsolve(H, (heom_env, Q), 10, qt.ket2dm(initial_state), tlist)
+heom_env = ExponentialBosonicEnvironment(ck_real, vk_real, ck_imag, vk_imag)
+heom_sol = heom.heomsolve(H, (heom_env, Q), 10, ket2dm(initial_state), tlist)
 ```
 
 Secondly, for the Bloch-Redfield solver we can directly use the power spectrum as input:
@@ -352,6 +358,8 @@ Its deviation from unity tells us how well the simulation has converged.
 
 \[4\] [Breuer and Petruccione *The Theory of Open Quantum Systems*.](https://doi.org/10.1093/acprof:oso/9780199213900.001.0001)
 
+\[5\] [QuTiP 5: The Quantum Toolbox in Python](https://arxiv.org/abs/2412.04705)
+
 
 
 ## About