diff --git a/tutorials-v5/heom/heom-1a-spin-bath-model-basic.md b/tutorials-v5/heom/heom-1a-spin-bath-model-basic.md
index 4aa6a779..86f7dea7 100644
--- a/tutorials-v5/heom/heom-1a-spin-bath-model-basic.md
+++ b/tutorials-v5/heom/heom-1a-spin-bath-model-basic.md
@@ -5,7 +5,7 @@ jupyter:
       extension: .md
       format_name: markdown
       format_version: '1.3'
-      jupytext_version: 1.14.5
+      jupytext_version: 1.16.1
   kernelspec:
     display_name: Python 3 (ipykernel)
     language: python
diff --git a/tutorials-v5/heom/heom-1b-spin-bath-model-very-strong-coupling.md b/tutorials-v5/heom/heom-1b-spin-bath-model-very-strong-coupling.md
index e4851516..b15f24a2 100644
--- a/tutorials-v5/heom/heom-1b-spin-bath-model-very-strong-coupling.md
+++ b/tutorials-v5/heom/heom-1b-spin-bath-model-very-strong-coupling.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.16.1
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -81,7 +81,7 @@ Note that in the above, and the following, we set $\hbar = k_\mathrm{B} = 1$.
 
 ## Setup
 
-```{code-cell} ipython3
+```{code-cell}
 import contextlib
 import time
 
@@ -113,13 +113,13 @@ from qutip.solver.heom import (
 
 Let's define some helper functions for calculating correlation function expansions, plotting results and timing how long operations take:
 
-```{code-cell} ipython3
+```{code-cell}
 def cot(x):
     """ Vectorized cotangent of x. """
     return 1. / np.tan(x)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 @contextlib.contextmanager
 def timer(label):
     """ Simple utility for timing functions:
@@ -133,7 +133,7 @@ def timer(label):
     print(f"{label}: {end - start}")
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Solver options:
 
 options = {
@@ -150,19 +150,19 @@ options = {
 
 And let us set up the system Hamiltonian, bath and system measurement operators:
 
-```{code-cell} ipython3
+```{code-cell}
 # Defining the system Hamiltonian
 eps = .0     # Energy of the 2-level system.
 Del = .2     # Tunnelling term
 Hsys = 0.5 * eps * sigmaz() + 0.5 * Del * sigmax()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Initial state of the system.
 rho0 = basis(2, 0) * basis(2, 0).dag()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # System-bath coupling (Drude-Lorentz spectral density)
 Q = sigmaz()  # coupling operator
 
@@ -185,7 +185,7 @@ NC = 13
 tlist = np.linspace(0, np.pi / Del, 600)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Define some operators with which we will measure the system
 # 1,1 element of density matrix - corresonding to groundstate
 P11p = basis(2, 0) * basis(2, 0).dag()
@@ -198,7 +198,7 @@ P12p = basis(2, 0) * basis(2, 1).dag()
 
 Let us briefly inspect the spectral density.
 
-```{code-cell} ipython3
+```{code-cell}
 w = np.linspace(0, 5, 1000)
 J = w * 2 * lam * gamma / ((gamma**2 + w**2))
 
@@ -211,7 +211,7 @@ axes.set_ylabel(r'J', fontsize=28);
 
 ## Simulation 1: Matsubara decomposition, not using Ishizaki-Tanimura terminator
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     bath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
     HEOMMats = HEOMSolver(Hsys, bath, NC, options=options)
@@ -222,7 +222,7 @@ with timer("ODE solver time"):
 
 ## Simulation 2: Matsubara decomposition (including terminator)
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     bath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
     _, terminator = bath.terminator()
@@ -233,7 +233,7 @@ with timer("ODE solver time"):
     resultMatsT = HEOMMatsT.run(rho0, tlist)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Plot the results
 fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 8))
 
@@ -256,7 +256,7 @@ axes.legend(loc=0, fontsize=12);
 
 ## Simulation 3: Pade decomposition
 
-```{code-cell} ipython3
+```{code-cell}
 # First, compare Matsubara and Pade decompositions
 matsBath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
 padeBath = DrudeLorentzPadeBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
@@ -324,7 +324,7 @@ ax2.set_xlabel(r't', fontsize=28)
 ax2.legend(loc=0, fontsize=12);
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     bath = DrudeLorentzPadeBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
     HEOMPade = HEOMSolver(Hsys, bath, NC, options=options)
@@ -333,7 +333,7 @@ with timer("ODE solver time"):
     resultPade = HEOMPade.run(rho0, tlist)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Plot the results
 fig, axes = plt.subplots(figsize=(8, 8))
 
@@ -358,7 +358,7 @@ axes.legend(loc=0, fontsize=12);
 
 ## Simulation 4: Fitting approach
 
-```{code-cell} ipython3
+```{code-cell}
 def wrapper_fit_func(x, N, args):
     """ Fit function wrapper that unpacks its arguments. """
     x = np.array(x)
@@ -410,7 +410,7 @@ def fitter(ans, tlist, k):
     return (a, b)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Fitting the real part of the correlation function:
 
 # Correlation function values to fit:
@@ -425,7 +425,7 @@ with timer("Correlation (real) fitting time"):
         poptR.append(fitter(corrRana, tlist_fit, i + 1))
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 plt.plot(tlist_fit, corrRana, label="Analytic")
 
 for i in range(kR):
@@ -437,7 +437,7 @@ plt.legend()
 plt.show()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Set the exponential coefficients from the fit parameters
 
 ckAR1 = poptR[-1][0]
@@ -452,7 +452,7 @@ ckAI = [lam * gamma * (-1.0) + 0j]
 vkAI = [gamma + 0j]
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     bath = BosonicBath(Q, ckAR, vkAR, ckAI, vkAI)
     # We reduce NC slightly here for speed of execution because we retain
@@ -466,7 +466,7 @@ with timer("ODE solver time"):
 
 ## Simulation 5: Bloch-Redfield
 
-```{code-cell} ipython3
+```{code-cell}
 DL = (
     "2 * pi * 2.0 * {lam} / (pi * {gamma} * {beta}) if (w==0) "
     "else 2 * pi * (2.0 * {lam} * {gamma} * w / (pi * (w**2 + {gamma}**2))) "
@@ -484,7 +484,7 @@ with timer("ODE solver time"):
 
 Finally, let's plot all of our different results to see how they shape up against each other.
 
-```{code-cell} ipython3
+```{code-cell}
 # Calculate expectation values in the bases:
 P11_mats = np.real(expect(resultMats.states, P11p))
 P11_matsT = np.real(expect(resultMatsT.states, P11p))
@@ -493,7 +493,7 @@ P11_fit = np.real(expect(resultFit.states, P11p))
 P11_br = np.real(expect(resultBR.states, P11p))
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 rcParams = {
     "axes.titlesize": 25,
     "axes.labelsize": 30,
@@ -510,7 +510,7 @@ rcParams = {
 }
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 fig, axes = plt.subplots(1, 1, sharex=True, figsize=(12, 7))
 
 with plt.rc_context(rcParams):
@@ -550,7 +550,7 @@ with plt.rc_context(rcParams):
 
 ## About
 
-```{code-cell} ipython3
+```{code-cell}
 qutip.about()
 ```
 
@@ -558,7 +558,7 @@ qutip.about()
 
 This section can include some tests to verify that the expected outputs are generated within the notebook. We put this section at the end of the notebook, so it's not interfering with the user experience. Please, define the tests using assert, so that the cell execution fails if a wrong output is generated.
 
-```{code-cell} ipython3
+```{code-cell}
 assert np.allclose(P11_matsT, P11_pade, rtol=1e-3)
 assert np.allclose(P11_matsT, P11_fit, rtol=1e-3)
 ```
diff --git a/tutorials-v5/heom/heom-1c-spin-bath-model-underdamped-sd.md b/tutorials-v5/heom/heom-1c-spin-bath-model-underdamped-sd.md
index a20dd73b..faf448f3 100644
--- a/tutorials-v5/heom/heom-1c-spin-bath-model-underdamped-sd.md
+++ b/tutorials-v5/heom/heom-1c-spin-bath-model-underdamped-sd.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.16.1
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -76,7 +76,7 @@ Note that in the above, and the following, we set $\hbar = k_\mathrm{B} = 1$.
 
 ## Setup
 
-```{code-cell} ipython3
+```{code-cell}
 import contextlib
 import time
 
@@ -107,19 +107,19 @@ from qutip.solver.heom import (
 
 Let's define some helper functions for calculating correlation function expansions, plotting results and timing how long operations take:
 
-```{code-cell} ipython3
+```{code-cell}
 def cot(x):
     """ Vectorized cotangent of x. """
     return 1. / np.tan(x)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def coth(x):
     """ Vectorized hyperbolic cotangent of x. """
     return 1. / np.tanh(x)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def underdamped_matsubara_params(lam, gamma, T, nk):
     """ Calculation of the real and imaginary expansions of the
         underdamped correlation functions.
@@ -161,7 +161,7 @@ def underdamped_matsubara_params(lam, gamma, T, nk):
     return ckAR, vkAR, ckAI, vkAI
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def plot_result_expectations(plots, axes=None):
     """ Plot the expectation values of operators as functions of time.
 
@@ -189,7 +189,7 @@ def plot_result_expectations(plots, axes=None):
     return fig
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 @contextlib.contextmanager
 def timer(label):
     """ Simple utility for timing functions:
@@ -203,7 +203,7 @@ def timer(label):
     print(f"{label}: {end - start}")
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Solver options:
 
 options = {
@@ -220,19 +220,19 @@ options = {
 
 And let us set up the system Hamiltonian, bath and system measurement operators:
 
-```{code-cell} ipython3
+```{code-cell}
 # Defining the system Hamiltonian
 eps = .5     # Energy of the 2-level system.
 Del = 1.0    # Tunnelling term
 Hsys = 0.5 * eps * sigmaz() + 0.5 * Del * sigmax()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Initial state of the system.
 rho0 = basis(2, 0) * basis(2, 0).dag()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # System-bath coupling (underdamed spectral density)
 Q = sigmaz()  # coupling operator
 
@@ -256,7 +256,7 @@ NC = 10
 tlist = np.linspace(0, 50, 1000)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Define some operators with which we will measure the system
 # 1,1 element of density matrix - corresonding to groundstate
 P11p = basis(2, 0) * basis(2, 0).dag()
@@ -267,7 +267,7 @@ P12p = basis(2, 0) * basis(2, 1).dag()
 
 ### First let us look at what the underdamped spectral density looks like:
 
-```{code-cell} ipython3
+```{code-cell}
 def plot_spectral_density():
     """ Plot the underdamped spectral density """
     w = np.linspace(0, 5, 1000)
@@ -288,7 +288,7 @@ The correlation functions are now very oscillatory, because of the Lorentzian pe
 
 ### So next, let us plot the correlation functions themselves:
 
-```{code-cell} ipython3
+```{code-cell}
 def Mk(t, k, gamma, w0, beta):
     """ Calculate the Matsubara terms for a given t and k. """
     Om = np.sqrt(w0**2 - (gamma / 2)**2)
@@ -340,7 +340,7 @@ plot_correlation_function()
 
 It is useful to look at what the Matsubara contributions do to this spectral density. We see that they modify the real part around $t=0$:
 
-```{code-cell} ipython3
+```{code-cell}
 def plot_matsubara_correlation_function_contributions():
     """ Plot the underdamped correlation function. """
     t = np.linspace(0, 20, 1000)
@@ -374,7 +374,7 @@ Next we calculate the exponents using the Matsubara decompositions. Here we spli
 
 The HEOM code will optimize these, and reduce the number of exponents when real and imaginary parts have the same exponent. This is clearly the case for the first term in the vkAI and vkAR lists.
 
-```{code-cell} ipython3
+```{code-cell}
 ckAR, vkAR, ckAI, vkAI = underdamped_matsubara_params(
     lam=lam, gamma=gamma, T=T, nk=Nk,
 )
@@ -386,7 +386,7 @@ The solver constructs the "right hand side" (RHS) determinining how the system a
 
 Below we create the bath and solver and then solve for the dynamics by calling `.run(rho0, tlist)`.
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     bath = BosonicBath(Q, ckAR, vkAR, ckAI, vkAI)
     HEOMMats = HEOMSolver(Hsys, bath, NC, options=options)
@@ -395,7 +395,7 @@ with timer("ODE solver time"):
     resultMats = HEOMMats.run(rho0, tlist)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 plot_result_expectations([
     (resultMats, P11p, 'b', "P11 Mats"),
     (resultMats, P12p, 'r', "P12 Mats"),
@@ -408,7 +408,7 @@ to perform this expansion will allow you to construct your own baths for other s
 
 Below we show how to use this built-in functionality:
 
-```{code-cell} ipython3
+```{code-cell}
 # Compare to built-in under-damped bath:
 
 with timer("RHS construction time"):
@@ -419,7 +419,7 @@ with timer("ODE solver time"):
     result_udbath = HEOM_udbath.run(rho0, tlist)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 plot_result_expectations([
     (result_udbath, P11p, 'b', "P11 (UnderDampedBath)"),
     (result_udbath, P12p, 'r', "P12 (UnderDampedBath)"),
@@ -432,7 +432,7 @@ plot_result_expectations([
 
 ### We can compare these results to those of the Bloch-Redfield solver in QuTiP:
 
-```{code-cell} ipython3
+```{code-cell}
 UD = (
     f"2 * {lam}**2 * {gamma} / ( {w0}**4 * {beta}) if (w==0)"
     " else "
@@ -447,7 +447,7 @@ with timer("ODE solver time"):
     )
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 plot_result_expectations([
     (resultMats, P11p, 'b', "P11 Mats"),
     (resultMats, P12p, 'r', "P12 Mats"),
@@ -462,7 +462,7 @@ plot_result_expectations([
 
 The thermal state of a reaction coordinate (treating the environment as a single damped mode) should, at high temperatures and small gamma, tell us the steady-state:
 
-```{code-cell} ipython3
+```{code-cell}
 dot_energy, dot_state = Hsys.eigenstates()
 deltaE = dot_energy[1] - dot_energy[0]
 
@@ -495,7 +495,7 @@ P11RC = tensor(qeye(NRC), basis(2, 0) * basis(2, 0).dag())
 P11RC = expect(rhoss, P11RC)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 rcParams = {
     "axes.titlesize": 25,
     "axes.labelsize": 30,
@@ -512,7 +512,7 @@ rcParams = {
 }
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 fig, axes = plt.subplots(1, 1, sharex=True, figsize=(12, 7))
 
 with plt.rc_context(rcParams):
@@ -541,7 +541,7 @@ with plt.rc_context(rcParams):
 
 ## About
 
-```{code-cell} ipython3
+```{code-cell}
 qutip.about()
 ```
 
@@ -549,7 +549,7 @@ qutip.about()
 
 This section can include some tests to verify that the expected outputs are generated within the notebook. We put this section at the end of the notebook, so it's not interfering with the user experience. Please, define the tests using assert, so that the cell execution fails if a wrong output is generated.
 
-```{code-cell} ipython3
+```{code-cell}
 assert np.allclose(
     expect(P11p, resultMats.states[-100:]), P11RC, rtol=1e-2,
 )
diff --git a/tutorials-v5/heom/heom-1d-spin-bath-model-ohmic-fitting.md b/tutorials-v5/heom/heom-1d-spin-bath-model-ohmic-fitting.md
index 591e2d51..3742dcd5 100644
--- a/tutorials-v5/heom/heom-1d-spin-bath-model-ohmic-fitting.md
+++ b/tutorials-v5/heom/heom-1d-spin-bath-model-ohmic-fitting.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.16.1
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -39,7 +39,7 @@ In each case we will use the fit parameters to determine the correlation functio
 
 ## Setup
 
-```{code-cell} ipython3
+```{code-cell}
 import contextlib
 import dataclasses
 import time
@@ -78,13 +78,13 @@ mp.pretty = True
 
 Let's define some helper functions for plotting results and timing how long operations take:
 
-```{code-cell} ipython3
+```{code-cell}
 def coth(x):
     """ Vectorized hyperbolic cotangent of x. """
     return 1. / np.tanh(x)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def plot_result_expectations(plots, axes=None):
     """ Plot the expectation values of operators as functions of time.
 
@@ -118,7 +118,7 @@ def plot_result_expectations(plots, axes=None):
     return fig
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 @contextlib.contextmanager
 def timer(label):
     """ Simple utility for timing functions:
@@ -132,7 +132,7 @@ def timer(label):
     print(f"{label}: {end - start}")
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Solver options:
 
 options = {
@@ -153,21 +153,21 @@ And let us set up the system Hamiltonian, bath and system measurement operators:
 
 ### System Hamiltonian
 
-```{code-cell} ipython3
+```{code-cell}
 # Defining the system Hamiltonian
 eps = 0.0    # Energy of the 2-level system.
 Del = 0.2    # Tunnelling term
 Hsys = 0.5 * eps * sigmaz() + 0.5 * Del * sigmax()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Initial state of the system.
 rho0 = basis(2, 0) * basis(2, 0).dag()
 ```
 
 ### System measurement operators
 
-```{code-cell} ipython3
+```{code-cell}
 # Define some operators with which we will measure the system
 # 1,1 element of density matrix - corresonding to groundstate
 P11p = basis(2, 0) * basis(2, 0).dag()
@@ -203,7 +203,7 @@ The corresponding spectral density for the Ohmic case is:
 J(\omega) = \omega \alpha e^{- \frac{\omega}{\omega_c}}
 \end{equation}
 
-```{code-cell} ipython3
+```{code-cell}
 def ohmic_correlation(t, alpha, wc, beta, s=1):
     """ The Ohmic bath correlation function as a function of t
         (and the bath parameters).
@@ -221,7 +221,7 @@ def ohmic_correlation(t, alpha, wc, beta, s=1):
     ], dtype=np.complex128)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def ohmic_spectral_density(w, alpha, wc):
     """ The Ohmic bath spectral density as a function of w
         (and the bath parameters).
@@ -229,7 +229,7 @@ def ohmic_spectral_density(w, alpha, wc):
     return w * alpha * np.e**(-w / wc)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def ohmic_power_spectrum(w, alpha, wc, beta):
     """ The Ohmic bath power spectrum as a function of w
         (and the bath parameters).
@@ -246,7 +246,7 @@ def ohmic_power_spectrum(w, alpha, wc, beta):
 
 Finally, let's set the bath parameters we will work with and write down some measurement operators:
 
-```{code-cell} ipython3
+```{code-cell}
 # Bath parameters:
 
 @dataclasses.dataclass
@@ -270,7 +270,7 @@ obp = OhmicBathParameters()
 
 And set the cut-off for the HEOM hierarchy:
 
-```{code-cell} ipython3
+```{code-cell}
 # HEOM parameters:
 
 # The max_depth defaults to 5 so that the notebook executes more
@@ -290,7 +290,7 @@ J_{\mathrm approx}(\omega; a, b, c) = \sum_{i=0}^{k-1} \frac{2 a_i b_i w}{((w +
 
 where $a, b$ and $c$ are the fit parameters and each is a vector of length $k$.
 
-```{code-cell} ipython3
+```{code-cell}
 # Helper functions for packing the paramters a, b and c into a single numpy
 # array as required by SciPy's curve_fit:
 
@@ -308,7 +308,7 @@ def unpack(params):
     return a, b, c
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # The approximate spectral density and a helper for fitting the approximate
 # spectral density to values calculated from the analytical formula:
 
@@ -349,7 +349,7 @@ def fit_spectral_density(J, w, alpha, wc, N):
 
 With the spectral density approximation $J_{\mathrm approx}(w; a, b, c)$ implemented above, we can now perform the fit and examine the results.
 
-```{code-cell} ipython3
+```{code-cell}
 w = np.linspace(0, 25, 20000)
 J = ohmic_spectral_density(w, alpha=obp.alpha, wc=obp.wc)
 
@@ -361,7 +361,7 @@ params_k = [
 
 Let's plot the fit for each $k$ and examine how it improves with an increasing number of terms:
 
-```{code-cell} ipython3
+```{code-cell}
 for k, params in enumerate(params_k):
     lam, gamma, w0 = params
     y = spectral_density_approx(w, lam, gamma, w0)
@@ -372,7 +372,7 @@ for k, params in enumerate(params_k):
 
 The fit with four terms looks good. Let's take a closer look at it by plotting the contribution of each term of the fit:
 
-```{code-cell} ipython3
+```{code-cell}
 # The parameters for the fit with four terms:
 
 lam, gamma, w0 = params_k[-1]
@@ -380,7 +380,7 @@ lam, gamma, w0 = params_k[-1]
 print(f"Parameters [k={len(params_k) - 1}]: lam={lam}; gamma={gamma}; w0={w0}")
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Plot the components of the fit separately:
 
 def spectral_density_ith_component(w, i, lam, gamma, w0):
@@ -414,7 +414,7 @@ plot_spectral_density_fit_components(J, w, lam, gamma, w0);
 
 And let's also compare the power spectrum of the fit and the analytical spectral density:
 
-```{code-cell} ipython3
+```{code-cell}
 def plot_power_spectrum(alpha, wc, beta, lam, gamma, w0, save=True):
     """ Plot the power spectrum of a fit against the actual power spectrum. """
     w = np.linspace(-10, 10, 50000)
@@ -442,7 +442,7 @@ plot_power_spectrum(obp.alpha, obp.wc, obp.beta, lam, gamma, w0, save=False)
 
 Now that we have a good fit to the spectral density, we can calculate the Matsubara expansion terms for the `BosonicBath` from them. At the same time we will calculate the Matsubara terminator for this expansion.
 
-```{code-cell} ipython3
+```{code-cell}
 def matsubara_coefficients_from_spectral_fit(lam, gamma, w0, beta, Q, Nk):
     """ Calculate the Matsubara co-efficients for a fit to the spectral
         density.
@@ -511,7 +511,7 @@ def matsubara_coefficients_from_spectral_fit(lam, gamma, w0, beta, Q, Nk):
     return ckAR, vkAR, ckAI, vkAI, terminator
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def generate_spectrum_results(obp, params, Nk, max_depth):
     """ Run the HEOM with the given bath parameters and
         and return the results of the evolution.
@@ -539,7 +539,7 @@ def generate_spectrum_results(obp, params, Nk, max_depth):
 
 Below we generate results for different convergence parameters (number of terms in the fit, number of matsubara terms, and depth of the hierarchy).  For the parameter choices here, we need a relatively large depth of around '11', which can be a little slow.
 
-```{code-cell} ipython3
+```{code-cell}
 # Generate results for different number of lorentzians in fit:
 
 results_spectral_fit_pk = [
@@ -556,7 +556,7 @@ plot_result_expectations([
 ]);
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # generate results for different number of Matsubara terms per Lorentzian
 # for max number of Lorentzians:
 
@@ -575,7 +575,7 @@ plot_result_expectations([
 ]);
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Generate results for different depths:
 
 Nc_list = range(2, max_depth)
@@ -595,7 +595,7 @@ plot_result_expectations([
 
 We now combine the fitting and correlation function data into one large plot.
 
-```{code-cell} ipython3
+```{code-cell}
 def correlation_approx_matsubara(t, ck, vk):
     """ Calculate the approximate real or imaginary part of the
         correlation function from the matsubara expansion co-efficients.
@@ -605,7 +605,7 @@ def correlation_approx_matsubara(t, ck, vk):
     return np.sum(ck[:, None] * np.exp(-vk[:, None] * t), axis=0)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def plot_cr_fit_vs_actual(t, ckAR, vkAR, C, axes):
     """ Plot the C_R(t) fit. """
     yR = correlation_approx_matsubara(t, ckAR, vkAR)
@@ -732,7 +732,7 @@ def plot_matsubara_spectrum_fit_vs_actual(
     return fig
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 t = np.linspace(0, 15, 100)
 C = ohmic_correlation(t, alpha=obp.alpha, wc=obp.wc, beta=obp.beta)
 
@@ -763,7 +763,7 @@ $$C_R^F(t) = \sum_{i=1}^{k_R} c_R^ie^{-\gamma_R^i t}\cos(\omega_R^i t)$$
 
 $$C_I^F(t) = \sum_{i=1}^{k_I} c_I^ie^{-\gamma_I^i t}\sin(\omega_I^i t)$$
 
-```{code-cell} ipython3
+```{code-cell}
 # The approximate correlation functions and a helper for fitting
 # the approximate correlation function to values calculated from
 # the analytical formula:
@@ -836,7 +836,7 @@ def fit_correlation_imag(C, t, wc, N):
     return unpack(params)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 t = np.linspace(0, 15, 15000)
 C = ohmic_correlation(t, alpha=obp.alpha, wc=obp.wc, beta=obp.beta)
 
@@ -851,7 +851,7 @@ params_k_imag = [
 ]
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 for k, params in enumerate(params_k_real):
     lam, gamma, w0 = params
     y = correlation_approx_real(t, lam, gamma, w0)
@@ -862,7 +862,7 @@ for k, params in enumerate(params_k_real):
     plt.show()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 for k, params in enumerate(params_k_imag):
     lam, gamma, w0 = params
     y = correlation_approx_imag(t, lam, gamma, w0)
@@ -875,7 +875,7 @@ for k, params in enumerate(params_k_imag):
 
 Now we construct the `BosonicBath` co-efficients and frequencies from the fit to the correlation function:
 
-```{code-cell} ipython3
+```{code-cell}
 def matsubara_coefficients_from_corr_fit_real(lam, gamma, w0):
     """ Return the matsubara coefficients for the imaginary part
         of the correlation function.
@@ -904,12 +904,12 @@ def matsubara_coefficients_from_corr_fit_imag(lam, gamma, w0):
     return ckAI, vkAI
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 ckAR, vkAR = matsubara_coefficients_from_corr_fit_real(*params_k_real[-1])
 ckAI, vkAI = matsubara_coefficients_from_corr_fit_imag(*params_k_imag[-1])
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def corr_spectrum_approx(w, ckAR, vkAR, ckAI, vkAI):
     """ Calculates the approximate power spectrum from ck and vk. """
     S = np.zeros(len(w), dtype=np.complex128)
@@ -926,7 +926,7 @@ def corr_spectrum_approx(w, ckAR, vkAR, ckAI, vkAI):
     return S
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def plot_jw_correlation_fit_vs_actual(matsubara_fit, obp, axes):
     """ Plot J(w) from the correlation fit. """
     [ckAR, vkAR, ckAI, vkAI] = matsubara_fit
@@ -1012,7 +1012,7 @@ def plot_matsubara_correlation_fit_vs_actual(t, C, matsubara_fit, obp):
     )
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 t = np.linspace(0, 15, 100)
 C = ohmic_correlation(t, alpha=obp.alpha, wc=obp.wc, beta=obp.beta)
 
@@ -1023,7 +1023,7 @@ plot_matsubara_correlation_fit_vs_actual(
 )
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def generate_corr_results(params_real, params_imag, max_depth):
     ckAR, vkAR = matsubara_coefficients_from_corr_fit_real(
         *params_real
@@ -1056,7 +1056,7 @@ results_corr_fit_pk = [
 ]
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 plot_result_expectations([
     (
         result, P11p, 'rand',
@@ -1066,7 +1066,7 @@ plot_result_expectations([
 ]);
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 fig, axes = plt.subplots(1, 1, sharex=True, figsize=(12, 7))
 
 plot_result_expectations([
@@ -1091,7 +1091,7 @@ axes.legend(loc=0, fontsize=20);
 
 ## About
 
-```{code-cell} ipython3
+```{code-cell}
 qutip.about()
 ```
 
@@ -1099,7 +1099,7 @@ qutip.about()
 
 This section can include some tests to verify that the expected outputs are generated within the notebook. We put this section at the end of the notebook, so it's not interfering with the user experience. Please, define the tests using assert, so that the cell execution fails if a wrong output is generated.
 
-```{code-cell} ipython3
+```{code-cell}
 assert np.allclose(
     expect(P11p, results_spectral_fit_pk[2].states),
     expect(P11p, results_spectral_fit_pk[3].states),
diff --git a/tutorials-v5/heom/heom-1e-spin-bath-model-pure-dephasing.md b/tutorials-v5/heom/heom-1e-spin-bath-model-pure-dephasing.md
index 83103511..d8dbfdc6 100644
--- a/tutorials-v5/heom/heom-1e-spin-bath-model-pure-dephasing.md
+++ b/tutorials-v5/heom/heom-1e-spin-bath-model-pure-dephasing.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.16.1
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -66,7 +66,7 @@ Note that in the above, and the following, we set $\hbar = k_\mathrm{B} = 1$.
 
 ## Setup
 
-```{code-cell} ipython3
+```{code-cell}
 import contextlib
 import time
 
@@ -97,7 +97,7 @@ from qutip.solver.heom import (
 
 Let's define some helper functions for calculating correlation function expansions, plotting results and timing how long operations take:
 
-```{code-cell} ipython3
+```{code-cell}
 def cot(x):
     """ Vectorized cotangent of x. """
     return 1. / np.tan(x)
@@ -108,7 +108,7 @@ def coth(x):
     return 1. / np.tanh(x)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def plot_result_expectations(plots, axes=None):
     """ Plot the expectation values of operators as functions of time.
 
@@ -140,7 +140,7 @@ def plot_result_expectations(plots, axes=None):
     return fig
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 @contextlib.contextmanager
 def timer(label):
     """ Simple utility for timing functions:
@@ -154,7 +154,7 @@ def timer(label):
     print(f"{label}: {end - start}")
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Solver options:
 
 options = {
@@ -175,14 +175,14 @@ And let us set up the system Hamiltonian, bath and system measurement operators:
 
 Here we set $H_{sys}=0$, which means the interaction Hamiltonian and the system Hamiltonian commute, and we can compare the numerical results to a known analytical one.  We could in principle keep $\epsilon \neq 0$, but it just introduces fast system oscillations, so it is more convenient to set it to zero.
 
-```{code-cell} ipython3
+```{code-cell}
 # Defining the system Hamiltonian
 eps = 0.0  # Energy of the 2-level system.
 Del = 0.0  # Tunnelling term
 Hsys = 0.5 * eps * sigmaz() + 0.5 * Del * sigmax()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # System-bath coupling (Drude-Lorentz spectral density)
 Q = sigmaz()  # coupling operator
 
@@ -203,7 +203,7 @@ Nk = 3
 tlist = np.linspace(0, 50, 1000)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Define some operators with which we will measure the system
 # 1,1 element of density matrix - corresonding to groundstate
 P11p = basis(2, 0) * basis(2, 0).dag()
@@ -214,7 +214,7 @@ P12p = basis(2, 0) * basis(2, 1).dag()
 
 To get a non-trivial result we prepare the initial state in a superposition, and see how the bath destroys the coherence.
 
-```{code-cell} ipython3
+```{code-cell}
 # Initial state of the system.
 psi = (basis(2, 0) + basis(2, 1)).unit()
 rho0 = psi * psi.dag()
@@ -222,7 +222,7 @@ rho0 = psi * psi.dag()
 
 ## Simulation 1: Matsubara decomposition, not using Ishizaki-Tanimura terminator
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     bath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
     HEOMMats = HEOMSolver(Hsys, bath, NC, options=options)
@@ -231,7 +231,7 @@ with timer("ODE solver time"):
     resultMats = HEOMMats.run(rho0, tlist)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Plot the results so far
 plot_result_expectations([
     (resultMats, P11p, 'b', "P11 Matsubara"),
@@ -241,7 +241,7 @@ plot_result_expectations([
 
 ## Simulation 2: Matsubara decomposition (including terminator)
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     bath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
     _, terminator = bath.terminator()
@@ -252,7 +252,7 @@ with timer("ODE solver time"):
     resultMatsT = HEOMMatsT.run(rho0, tlist)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Plot the results
 plot_result_expectations([
     (resultMats, P11p, 'b', "P11 Matsubara"),
@@ -266,7 +266,7 @@ plot_result_expectations([
 
 As in example 1a, we can compare to Pade and Fitting approaches.
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     bath = DrudeLorentzPadeBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
     HEOMPade = HEOMSolver(Hsys, bath, NC, options=options)
@@ -275,7 +275,7 @@ with timer("ODE solver time"):
     resultPade = HEOMPade.run(rho0, tlist)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Plot the results
 plot_result_expectations([
     (resultMatsT, P11p, 'b', "P11 Matsubara (+term)"),
@@ -287,7 +287,7 @@ plot_result_expectations([
 
 ## Simulation 4: Fitting approach
 
-```{code-cell} ipython3
+```{code-cell}
 def c(t, Nk):
     """ Calculates real and imag. parts of the correlation function
         using Nk Matsubara terms.
@@ -315,7 +315,7 @@ corr_ana = c(tlist_fit, lmaxmats)
 corrRana, corrIana = np.real(corr_ana), np.imag(corr_ana)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def wrapper_fit_func(x, N, *args):
     """ Wrapper for fitting function. """
     a, b = args[0][:N], args[0][N:2*N]
@@ -383,7 +383,7 @@ for i in range(1):
     plt.show()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Set the exponential coefficients from the fit parameters
 
 ckAR = popt1[-1][0]
@@ -399,7 +399,7 @@ vkAI = -1 * popt2[-1][1]
 # vkAI = [complex(gamma)]
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     bath = BosonicBath(Q, ckAR, vkAR, ckAI, vkAI)
     HEOMFit = HEOMSolver(Hsys, bath, NC, options=options)
@@ -410,7 +410,7 @@ with timer("ODE solver time"):
 
 ## Analytic calculations
 
-```{code-cell} ipython3
+```{code-cell}
 def pure_dephasing_evolution_analytical(tlist, wq, ck, vk):
     """
     Computes the propagating function appearing in the pure dephasing model.
@@ -487,7 +487,7 @@ def correlation_integral(t, ck, vk):
 
 For the pure dephasing analytics, we just sum up as many matsubara terms as we can:
 
-```{code-cell} ipython3
+```{code-cell}
 lmaxmats2 = 15000
 
 vk = [complex(-gamma)]
@@ -509,7 +509,7 @@ P12_ana = 0.5 * pure_dephasing_evolution_analytical(
 
 Alternatively, we can just do the integral of the propagator directly, without using the correlation functions at all
 
-```{code-cell} ipython3
+```{code-cell}
 def JDL(omega, lamc, omega_c):
     return 2. * lamc * omega * omega_c / (omega_c**2 + omega**2)
 
@@ -532,7 +532,7 @@ P12_ana2 = [
 
 ## Compare results
 
-```{code-cell} ipython3
+```{code-cell}
 plot_result_expectations([
     (resultMats, P12p, 'r', "P12 Mats"),
     (resultMatsT, P12p, 'r--', "P12 Mats + Term"),
@@ -545,7 +545,7 @@ plot_result_expectations([
 
 We can't see much difference in the plot above, so let's do a log plot instead:
 
-```{code-cell} ipython3
+```{code-cell}
 fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 8))
 
 plot_result_expectations([
@@ -563,7 +563,7 @@ axes.legend(loc=0, fontsize=12);
 
 ## About
 
-```{code-cell} ipython3
+```{code-cell}
 qutip.about()
 ```
 
@@ -571,7 +571,7 @@ qutip.about()
 
 This section can include some tests to verify that the expected outputs are generated within the notebook. We put this section at the end of the notebook, so it's not interfering with the user experience. Please, define the tests using assert, so that the cell execution fails if a wrong output is generated.
 
-```{code-cell} ipython3
+```{code-cell}
 assert np.allclose(
     expect(P12p, resultMats.states[:15]), np.real(P12_ana)[:15],
     rtol=1e-2,
diff --git a/tutorials-v5/heom/heom-2-fmo-example.md b/tutorials-v5/heom/heom-2-fmo-example.md
index 7a301108..4bebe233 100644
--- a/tutorials-v5/heom/heom-2-fmo-example.md
+++ b/tutorials-v5/heom/heom-2-fmo-example.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.16.1
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -31,7 +31,7 @@ quantum environment reduces the effect of pure dephasing.
 
 ## Setup
 
-```{code-cell} ipython3
+```{code-cell}
 import contextlib
 import time
 
@@ -59,19 +59,19 @@ from qutip.solver.heom import (
 
 Let's define some helper functions for calculating correlation functions, spectral densities, thermal energy level occupations, and for plotting results and timing how long operations take:
 
-```{code-cell} ipython3
+```{code-cell}
 def cot(x):
     """ Vectorized cotangent of x. """
     return 1 / np.tan(x)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def J0(energy):
     """ Under-damped brownian oscillator spectral density. """
     return 2 * lam * gamma * energy / (energy**2 + gamma**2)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def J0_dephasing():
     """ Under-damped brownian oscillator dephasing probability.
 
@@ -80,13 +80,13 @@ def J0_dephasing():
     return 2 * lam * gamma / gamma**2
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def n_th(energy, T):
     """ The average occupation of a given energy level at temperature T. """
     return 1 / (np.exp(energy / T) - 1)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def dl_corr_approx(t, nk):
     """ Drude-Lorenz correlation function approximation.
 
@@ -99,7 +99,7 @@ def dl_corr_approx(t, nk):
     return c
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 @contextlib.contextmanager
 def timer(label):
     """ Simple utility for timing functions:
@@ -113,7 +113,7 @@ def timer(label):
     print(f"{label}: {end - start}")
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Solver options:
 
 options = {
@@ -131,7 +131,7 @@ options = {
 
 And let us set up the system Hamiltonian and bath parameters:
 
-```{code-cell} ipython3
+```{code-cell}
 # System Hamiltonian:
 #
 # We use the Hamiltonian employed in
@@ -149,7 +149,7 @@ Hsys = 3e10 * 2 * np.pi * Qobj([
 ])
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Bath parameters
 
 lam = 35 * 3e10 * 2 * np.pi
@@ -162,7 +162,7 @@ beta = 1 / T
 
 Let's quickly plot the spectral density and environment correlation functions so that we can see what they look like.
 
-```{code-cell} ipython3
+```{code-cell}
 wlist = np.linspace(0, 200 * 3e10 * 2 * np.pi, 100)
 tlist = np.linspace(0, 1e-12, 1000)
 
@@ -198,7 +198,7 @@ axes[1].legend();
 
 Now let us solve for the evolution of this system using the HEOM.
 
-```{code-cell} ipython3
+```{code-cell}
 # We start the excitation at site 1:
 rho0 = basis(7, 0) * basis(7, 0).dag()
 
@@ -226,7 +226,7 @@ for m in range(7):
     Ltot += terminator
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     HEOMMats = HEOMSolver(Hsys, baths, NC, options=options)
 
@@ -234,7 +234,7 @@ with timer("ODE solver time"):
     outputFMO_HEOM = HEOMMats.run(rho0, tlist)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 fig, axes = plt.subplots(1, 1, figsize=(12, 8))
 
 colors = ['r', 'g', 'b', 'y', 'c', 'm', 'k']
@@ -270,7 +270,7 @@ Now let us solve the same problem using the Bloch-Redfield solver. We will see t
 
 In the next section, we will examine the role of pure dephasing in the evolution to understand why this happens.
 
-```{code-cell} ipython3
+```{code-cell}
 DL = (
     f"2 * pi * 2.0 * {lam} / (pi * {gamma} * {beta}) if (w == 0) else "
     f"2 * pi * (2.0*{lam}*{gamma} *w /(pi*(w**2+{gamma}**2))) * "
@@ -287,7 +287,7 @@ with timer("BR ODE solver time"):
 
 And now let's plot the Bloch-Redfield solver results:
 
-```{code-cell} ipython3
+```{code-cell}
 fig, axes = plt.subplots(1, 1, figsize=(12, 8))
 
 for m, Q in enumerate(Q_list):
@@ -315,7 +315,7 @@ It is useful to construct the various parts of the Bloch-Redfield master equatio
 
 First we will write a function to return the list of collapse operators for a given system, either with or without the dephasing operators:
 
-```{code-cell} ipython3
+```{code-cell}
 def get_collapse(H, T, dephasing=1):
     """ Calculate collapse operators for a given system H and
         temperature T.
@@ -380,7 +380,7 @@ Now we are able to switch the pure dephasing tersms on and off.
 
 Let us starting by including the dephasing operators. We expect to see the same behaviour that we saw when using the Bloch-Redfield solver.
 
-```{code-cell} ipython3
+```{code-cell}
 # dephasing terms on, we recover the full BR solution:
 
 with timer("Building the collapse operators"):
@@ -390,7 +390,7 @@ with timer("ME ODE solver"):
     outputFMO_ME = mesolve(Hsys, rho0, tlist, collapse_list)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 fig, axes = plt.subplots(1, 1, figsize=(12, 8))
 
 for m, Q in enumerate(Q_list):
@@ -409,7 +409,7 @@ We see similar results to before.
 
 Now let us examine what happens when we remove the dephasing collapse operators:
 
-```{code-cell} ipython3
+```{code-cell}
 # dephasing terms off
 
 with timer("Building the collapse operators"):
@@ -419,7 +419,7 @@ with timer("ME ODE solver"):
     outputFMO_ME_nodephase = mesolve(Hsys, rho0, tlist, collapse_list)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 fig, axes = plt.subplots(1, 1, figsize=(12, 8))
 for m, Q in enumerate(Q_list):
     axes.plot(
@@ -443,7 +443,7 @@ And now we see that without the dephasing, the oscillations reappear. The full d
 
 ## About
 
-```{code-cell} ipython3
+```{code-cell}
 qutip.about()
 ```
 
@@ -451,7 +451,7 @@ qutip.about()
 
 This section can include some tests to verify that the expected outputs are generated within the notebook. We put this section at the end of the notebook, so it's not interfering with the user experience. Please, define the tests using assert, so that the cell execution fails if a wrong output is generated.
 
-```{code-cell} ipython3
+```{code-cell}
 assert np.allclose(
     expect(outputFMO_BR.states, Q_list[0]),
     expect(outputFMO_ME.states, Q_list[0]),
diff --git a/tutorials-v5/heom/heom-3-quantum-heat-transport.md b/tutorials-v5/heom/heom-3-quantum-heat-transport.md
index 8e956d6c..73918687 100644
--- a/tutorials-v5/heom/heom-3-quantum-heat-transport.md
+++ b/tutorials-v5/heom/heom-3-quantum-heat-transport.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.16.1
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -50,7 +50,7 @@ References:
 
 ## Setup
 
-```{code-cell} ipython3
+```{code-cell}
 import dataclasses
 
 import numpy as np
@@ -71,7 +71,7 @@ from IPython.display import display
 
 ## Helpers
 
-```{code-cell} ipython3
+```{code-cell}
 # Solver options:
 
 options = {
@@ -87,7 +87,7 @@ options = {
 
 ## System and bath definition
 
-```{code-cell} ipython3
+```{code-cell}
 @dataclasses.dataclass
 class SystemParams:
     """ System parameters and Hamiltonian. """
@@ -116,7 +116,7 @@ class SystemParams:
         return dataclasses.replace(self, **kw)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 @dataclasses.dataclass
 class BathParams:
     """ Bath parameters. """
@@ -169,7 +169,7 @@ In the expression for the bath heat currents, we left out terms involving $[Q_1,
 
 In QuTiP, these currents can be conveniently calculated as follows:
 
-```{code-cell} ipython3
+```{code-cell}
 def bath_heat_current(bath_tag, ado_state, hamiltonian, coupling_op, delta=0):
     """
     Bath heat current from the system into the heat bath with the given tag.
@@ -264,7 +264,7 @@ Note that at long times, we expect $j_{\text{B}}^1 = -j_{\text{B}}^2$ and $j_{\t
 
 For our simulations, we will represent the bath spectral densities using the first term of their Padé decompositions, and we will use $7$ levels of the HEOM hierarchy.
 
-```{code-cell} ipython3
+```{code-cell}
 Nk = 1
 NC = 7
 ```
@@ -274,7 +274,7 @@ NC = 7
 We fix $J_{12} = 0.1 \epsilon$ (as in Fig. 3(a-ii) of Ref. \[2\]) and choose the fixed coupling strength $\lambda_1 = \lambda_2 = J_{12}\, /\, (2\epsilon)$ (corresponding to $\bar\zeta = 1$ in Ref. \[2\]).
 Using these values, we will study the time evolution of the system state and the heat currents.
 
-```{code-cell} ipython3
+```{code-cell}
 # fix qubit-qubit and qubit-bath coupling strengths
 sys = SystemParams(J12=0.1)
 bath_p1 = BathParams(qubit=0, sign="+", lam=sys.J12 / 2)
@@ -287,7 +287,7 @@ rho0 = qt.tensor(qt.identity(2), qt.identity(2)) / 4
 tlist = np.linspace(0, 50, 250)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 H = sys.H()
 
 bath1 = bath_p1.bath(Nk, tag='bath 1')
@@ -316,7 +316,7 @@ result = solver.run(rho0, tlist, e_ops=[
 
 We first plot $\langle \sigma_z^1 \rangle$ to see the time evolution of the system state:
 
-```{code-cell} ipython3
+```{code-cell}
 fig, axes = plt.subplots(figsize=(8, 8))
 axes.plot(tlist, result.expect[0], 'r', linewidth=2)
 axes.set_xlabel('t', fontsize=28)
@@ -325,7 +325,7 @@ axes.set_ylabel(r"$\langle \sigma_z^1 \rangle$", fontsize=28);
 
 We find a rather quick thermalization of the system state. For the heat currents, however, it takes a somewhat longer time until they converge to their long-time values:
 
-```{code-cell} ipython3
+```{code-cell}
 fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(16, 8))
 
 ax1.plot(
@@ -377,7 +377,7 @@ ax2.legend(loc=0, fontsize=12);
 
 Here, we try to reproduce the HEOM curves in Fig. 3(a) of Ref. \[1\] by varying the coupling strength and finding the steady state for each coupling strength.
 
-```{code-cell} ipython3
+```{code-cell}
 def heat_currents(sys, bath_p1, bath_p2, Nk, NC, options):
     """ Calculate the steady sate heat currents for the given system and
         bath.
@@ -408,7 +408,7 @@ def heat_currents(sys, bath_p1, bath_p2, Nk, NC, options):
     )
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Define number of points to use for the plot
 plot_points = 10  # use 100 for a smoother curve
 
@@ -461,7 +461,7 @@ j3s = [
 
 ## Create Plot
 
-```{code-cell} ipython3
+```{code-cell}
 fig, axes = plt.subplots(figsize=(12, 7))
 
 axes.plot(
@@ -492,7 +492,7 @@ axes.legend(loc=0);
 
 ## About
 
-```{code-cell} ipython3
+```{code-cell}
 qt.about()
 ```
 
@@ -500,6 +500,6 @@ qt.about()
 
 This section can include some tests to verify that the expected outputs are generated within the notebook. We put this section at the end of the notebook, so it's not interfering with the user experience. Please, define the tests using assert, so that the cell execution fails if a wrong output is generated.
 
-```{code-cell} ipython3
+```{code-cell}
 assert 1 == 1
 ```
diff --git a/tutorials-v5/heom/heom-4-dynamical-decoupling.md b/tutorials-v5/heom/heom-4-dynamical-decoupling.md
index 952a44e9..c061ee74 100644
--- a/tutorials-v5/heom/heom-4-dynamical-decoupling.md
+++ b/tutorials-v5/heom/heom-4-dynamical-decoupling.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.16.1
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -27,7 +27,7 @@ We first show the standard example of equally spaced pulses, and then consider t
 
 ## Setup
 
-```{code-cell} ipython3
+```{code-cell}
 import numpy as np
 import matplotlib.pyplot as plt
 
@@ -55,14 +55,14 @@ from IPython.display import display
 
 Let's define some helper functions for calculating the spectral density:
 
-```{code-cell} ipython3
+```{code-cell}
 def dl_spectrum(w, lam, gamma):
     """ Return the Drude-Lorentz spectral density. """
     J = w * 2 * lam * gamma / (gamma**2 + w**2)
     return J
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Solver options:
 
 # The max_step must be set to a short time than the
@@ -84,7 +84,7 @@ options = {
 
 Now we define the system and bath properties and the HEOM parameters. The system is a single stationary qubit with $H = 0$ and the bath is a bosonic bath with a Drude-Lorentz spectrum.
 
-```{code-cell} ipython3
+```{code-cell}
 # Define the system Hamlitonian.
 #
 # The system isn't evolving by itself, so the Hamiltonian is 0 (with the
@@ -93,7 +93,7 @@ Now we define the system and bath properties and the HEOM parameters. The system
 H_sys = 0 * sigmaz()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Define some operators with which we will measure the system
 # 1,1 element of density matrix - corresonding to groundstate
 P11p = basis(2, 0) * basis(2, 0).dag()
@@ -102,7 +102,7 @@ P22p = basis(2, 1) * basis(2, 1).dag()
 P12p = basis(2, 0) * basis(2, 1).dag()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Properties for the Drude-Lorentz bath
 
 lam = 0.0005
@@ -118,7 +118,7 @@ Nk = 3
 bath = DrudeLorentzPadeBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # HEOM parameters
 
 # number of layers to keep in the hierarchy:
@@ -129,7 +129,7 @@ To perform the dynamic decoupling from the environment, we will drive the system
 
 Below we define a function that returns the pulse (which is itself a function):
 
-```{code-cell} ipython3
+```{code-cell}
 def drive(amplitude, delay, integral):
     """ Coefficient of the drive as a function of time.
 
@@ -166,7 +166,7 @@ H_drive = sigmax()
 
 Let's start by plotting the spectral density of our Drude-Lorentz bath:
 
-```{code-cell} ipython3
+```{code-cell}
 wlist = np.linspace(0, 0.5, 1000)
 J = dl_spectrum(wlist, lam, gamma)
 
@@ -184,7 +184,7 @@ First we will drive the system with fast, large amplitude pulses. Then we will d
 
 Let's start by simulating the fast pulses:
 
-```{code-cell} ipython3
+```{code-cell}
 # Fast driving (quick, large amplitude pulses)
 
 tlist = np.linspace(0, 400, 1000)
@@ -207,7 +207,7 @@ outputDD = hsolver.run(rho0, tlist)
 
 And now the longer slower pulses:
 
-```{code-cell} ipython3
+```{code-cell}
 # Slow driving (longer, small amplitude pulses)
 
 # without pulses
@@ -224,7 +224,7 @@ outputDDslow = hsolver.run(rho0, tlist)
 
 Now let's plot all of the results and the shapes of the pulses:
 
-```{code-cell} ipython3
+```{code-cell}
 def plot_dd_results(outputnoDD, outputDD, outputDDslow):
     fig, axes = plt.subplots(2, 1, sharex=False, figsize=(12, 12))
 
@@ -290,7 +290,7 @@ def plot_dd_results(outputnoDD, outputDD, outputDDslow):
     fig.tight_layout()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 plot_dd_results(outputnoDD, outputDD, outputDDslow)
 ```
 
@@ -310,7 +310,7 @@ $$
 
 This is just a convenient way to describe the varying delay. We could have chosen another monotonically increasing function to represent the cummulative delay (although it might not be as effective).
 
-```{code-cell} ipython3
+```{code-cell}
 def cummulative_delay_fractions(N):
     """ Return an array of N + 1 cummulative delay
         fractions.
@@ -364,7 +364,7 @@ Let's plot the cummulative delays and see what they look like. Note that the cum
 
 On the same axes we plot the individual $j^{th}$ delays as a fraction of the average delay.
 
-```{code-cell} ipython3
+```{code-cell}
 def plot_cummulative_delay_fractions(N):
     cummulative = cummulative_delay_fractions(N)
     individual = (cummulative[1:] - cummulative[:-1]) * N
@@ -380,7 +380,7 @@ plot_cummulative_delay_fractions(100)
 
 And now let us plot the first ten even and optimally spaced pulses together to compare them:
 
-```{code-cell} ipython3
+```{code-cell}
 def plot_even_and_optimally_spaced_pulses():
     amplitude = 10.0
     integral = np.pi / 2
@@ -408,7 +408,7 @@ Now let's simulate the effectiveness of the two sets of delays by comparing how
 
 We'll perform the simulation over a range of lambdas and gammas to show how the non-evenly spaced delays become optimal as the width of the bath spectral function increases.
 
-```{code-cell} ipython3
+```{code-cell}
 # Bath parameters to simulate over:
 
 # We use only two lambdas and two gammas so that the notebook executes
@@ -488,7 +488,7 @@ P12_results = [
 
 Now that we have the expectation values of $\rho_{01}$ let's plot them as a function of gamma for each lambda. Note how in each case the non-evenly spaced pulses become optimal once gamma is sufficiently small:
 
-```{code-cell} ipython3
+```{code-cell}
 fig, axes = plt.subplots(1, 1, sharex=False, figsize=(10, 7))
 colors = ["green", "red", "blue"]
 
@@ -518,7 +518,7 @@ And now you know about dynamically decoupling a qubit from its environment!
 
 ## About
 
-```{code-cell} ipython3
+```{code-cell}
 qutip.about()
 ```
 
@@ -526,6 +526,6 @@ qutip.about()
 
 This section can include some tests to verify that the expected outputs are generated within the notebook. We put this section at the end of the notebook, so it's not interfering with the user experience. Please, define the tests using assert, so that the cell execution fails if a wrong output is generated.
 
-```{code-cell} ipython3
+```{code-cell}
 assert 1 == 1
 ```
diff --git a/tutorials-v5/heom/heom-5a-fermions-single-impurity-model.md b/tutorials-v5/heom/heom-5a-fermions-single-impurity-model.md
index e07bf9ed..ed91e6e4 100644
--- a/tutorials-v5/heom/heom-5a-fermions-single-impurity-model.md
+++ b/tutorials-v5/heom/heom-5a-fermions-single-impurity-model.md
@@ -5,9 +5,9 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.16.1
 kernelspec:
-  display_name: Python 3 (ipykernel)
+  display_name: qutip-dev
   language: python
   name: python3
 ---
@@ -145,52 +145,11 @@ H = e1 * d1.dag() * d1
 ```
 
 ```{code-cell} ipython3
-# Define parameters for left and right fermionic baths.
-# Each bath is a lead (i.e. a wire held at a potential)
-# with temperature T and chemical potential mu.
+from qutip.core.environment import LorentzianEnvironment
 
-@dataclasses.dataclass
-class LorentzianBathParameters:
-    lead: str
-    Q: object  # coupling operator
-    gamma: float = 0.01  # coupling strength
-    W: float = 1.0  # cut-off
-    T: float = 0.025851991  # temperature
-    theta: float = 2.0  # bias
 
-    def __post_init__(self):
-        assert self.lead in ("L", "R")
-        self.beta = 1 / self.T
-        if self.lead == "L":
-            self.mu = self.theta / 2.0
-        else:
-            self.mu = - self.theta / 2.0
-
-    def J(self, w):
-        """ Spectral density. """
-        return self.gamma * self.W**2 / ((w - self.mu)**2 + self.W**2)
-
-    def fF(self, w, sign=1.0):
-        """ Fermi distribution for this bath. """
-        x = sign * self.beta * (w - self.mu)
-        return fF(x)
-
-    def lamshift(self, w):
-        """ Return the lamshift. """
-        return 0.5 * (w - self.mu) * self.J(w) / self.W
-
-    def replace(self, **kw):
-        return dataclasses.replace(self, **kw)
-
-
-def fF(x):
-    """ Return the Fermi distribution. """
-    # in units where kB = 1.0
-    return 1 / (np.exp(x) + 1)
-
-
-bath_L = LorentzianBathParameters(Q=d1, lead="L")
-bath_R = LorentzianBathParameters(Q=d1, lead="R")
+envL=LorentzianEnvironment(T= 0.025851991,W=1,mu=1,gamma=0.01,Nk=10)
+envR=LorentzianEnvironment(T= 0.025851991,W=1,mu=-1,gamma=0.01,Nk=10)
 ```
 
 ## Spectral density
@@ -198,12 +157,12 @@ bath_R = LorentzianBathParameters(Q=d1, lead="R")
 Let's plot the spectral density.
 
 ```{code-cell} ipython3
-w_list = np.linspace(-2, 2, 100)
+w_list = np.linspace(-80, 80, 2000)
 
 fig, ax = plt.subplots(figsize=(12, 7))
 
-spec_L = bath_L.J(w_list)
-spec_R = bath_R.J(w_list)
+spec_L = envL.spectral_density(w_list)
+spec_R = envR.spectral_density(w_list)
 
 ax.plot(
     w_list, spec_L,
@@ -221,159 +180,233 @@ ax.set_ylabel(r"$J(\omega)$")
 ax.legend();
 ```
 
-## Emission and absorption by the leads
-
-Next let's plot the emission and absorption by the leads.
+```{code-cell} ipython3
+from qutip.core.environment import FermionicEnvironment
+```
 
 ```{code-cell} ipython3
-w_list = np.linspace(-2, 2, 100)
+w_list = np.linspace(-30, 30, 2000)
+
+fenvL=FermionicEnvironment.from_spectral_density(envL.spectral_density(w_list),w_list,T=0.025851991,mu=1)
+fenvR=FermionicEnvironment.from_spectral_density(envR.spectral_density(w_list),w_list,T=0.025851991,mu=-1)
+# fenvL=FermionicEnvironment.from_power_spectra(envL.power_spectrum_minus(w_list),w_list,T=0.025851991,mu=1,sigma=-1)
+# fenvR=FermionicEnvironment.from_power_spectra(envR.power_spectrum_minus(w_list),w_list,T=0.025851991,mu=-1,sigma=-1)
+tk=np.linspace(0,60,2000)
+# fenvL=FermionicEnvironment.from_correlation_functions(envL.correlation_function_plus(tk),tk,T=0.025851991,mu=1,sigma=1)
+# fenvR=FermionicEnvironment.from_correlation_functions(envR.correlation_function_plus(tk),tk,T=0.025851991,mu=-1,sigma=1)
+fpenvL=fenvL._approx_by_prony(method="espira-I",tlist=tk,Np=10,Nm=10,tag="L")
+fpenvR=fenvR._approx_by_prony(method="espira-I",tlist=tk,Np=10,Nm=10,tag="R")
+```
 
+```{code-cell} ipython3
 fig, ax = plt.subplots(figsize=(12, 7))
-
-# Left lead emission and absorption
-
-gam_L_in = bath_L.J(w_list) * bath_L.fF(w_list, sign=1.0)
-gam_L_out = bath_L.J(w_list) * bath_L.fF(w_list, sign=-1.0)
-
 ax.plot(
-    w_list, gam_L_in,
-    "b--", linewidth=3,
-    label=r"S_L(w) input (absorption)",
+    w_list,envR.spectral_density(w_list) ,
+    "r", linewidth=3,
+    label=r"J_R(w)",
 )
 ax.plot(
-    w_list, gam_L_out,
-    "r--", linewidth=3,
-    label=r"S_L(w) output (emission)",
+    w_list, fenvR.spectral_density(w_list),
+    "g--", linewidth=3,
+    label=r"J_R_env(w)",
 )
-
-# Right lead emission and absorption
-
-gam_R_in = bath_R.J(w_list) * bath_R.fF(w_list, sign=1.0)
-gam_R_out = bath_R.J(w_list) * bath_R.fF(w_list, sign=-1.0)
-
 ax.plot(
-    w_list, gam_R_in,
-    "b", linewidth=3,
-    label=r"S_R(w) input (absorption)",
+    w_list, fpenvR.spectral_density(w_list),
+    "ko", linewidth=3,
+    label=r"J_R_env(w)",
 )
+ax.set_xlabel("w")
+ax.set_ylabel(r"$J(\omega)$")
+ax.legend();
+```
+
+```{code-cell} ipython3
+fig, ax = plt.subplots(figsize=(12, 7))
 ax.plot(
-    w_list, gam_R_out,
+    w_list,envL.correlation_function_plus(w_list).imag ,
     "r", linewidth=3,
-    label=r"S_R(w) output (emission)",
+    label=r"J_R(w)",
+)
+ax.plot(
+    w_list, fenvL.correlation_function_plus(w_list).imag,
+    "g--", linewidth=3,
+    label=r"J_R_env(w)",
+)
+ax.plot(
+    w_list, fpenvL.correlation_function_plus(w_list).imag,
+    "ko", linewidth=3,
+    label=r"J_R_env(w)",
 )
-
 ax.set_xlabel("w")
-ax.set_ylabel(r"$S(\omega)$")
+ax.set_ylabel(r"$J(\omega)$")
 ax.legend();
 ```
 
-## Comparing the Matsubara and Pade approximations
-
-Let's start by solving for the evolution using a Pade expansion of the correlation function of the Lorentzian spectral density:
+```{code-cell} ipython3
+from scipy.signal.windows import *
+```
 
 ```{code-cell} ipython3
-# HEOM dynamics using the Pade approximation:
+# def fg(w,x0=1):
+#     mask=w-x0<=0
+#     result=np.exp(-(w-x0)**2 / 0.03)
+#     result[mask]=1
+#     return result
+    
+# plt.plot(w_list,fg(w_list))
+```
 
-# Times to solve for and initial system state:
-tlist = np.linspace(0, 100, 1000)
-rho0 = basis(2, 0) * basis(2, 0).dag()
+```{code-cell} ipython3
+# from qutip.utilities import fermi_dirac
+# from scipy.fft import fftfreq
+
+# result=fenvL.power_spectrum_plus(w_list)
+# mask=w_list-fenvL.mu <0
+
+# result[~mask] = result[~mask]*fg(w_list[~mask])#gaussian(len(w_list),std=10*fenvL.mu)[~mask]
+# f = fermi_dirac(w_list, 1/envL.T, envL.mu)
+# ff=1/f
+# ff -= 1
+# result2=ff* result
+# problematic_indices = ~np.isfinite(result2)
+# # Only if there are problematic values, recalculate those specific points
+# if np.any(problematic_indices):
+#     # For problematic points, use an alternative calculation with a small epsilon
+#     ff_safe = 1/(fermi_dirac(w_list[problematic_indices],
+#                                 1/fenvL.T, fenvL.mu) + 1e-10)
+#     ff_safe -= 1
+# result2[problematic_indices]=ff_safe*result[problematic_indices]
+```
 
-Nk = 10  # Number of exponents to retain in the expansion of each bath
+```{code-cell} ipython3
+# w_list = np.linspace(-20, 20, 1600)
+
+# fig, ax = plt.subplots(figsize=(12, 7))
+
+# # Left lead emission and absorption
+
+# gam_L_in = envL.power_spectrum_plus(w_list)
+
+# # ax.plot(
+# #     w_list, gam_L_in,
+# #     "b", linewidth=3,
+# #     label=r"S_L(w) input (absorption)",
+# # )
+# ax.plot(
+#     w_list, gam_L_out,
+#     "r", linewidth=3,
+#     label=r"S_L(w) output (emission)",
+# )
+
+# # Right lead emission and absorption
+
+# # ax.plot(
+# #     w_list, result,
+# #     "r--", linewidth=3,
+# #     label=r"S_R(w) input (absorption)",
+# # )
+# ax.plot(
+#     w_list, result2,
+#     "b--", linewidth=3,
+#     label=r"S_R(w) output (emission)",
+# )
+# #ax.axvline(x=1)
+# ax.set_ylim(0,0.015)
+# ax.set_xlim(0,2)
+
+# ax.set_xlabel("w")
+# ax.set_ylabel(r"$S(\omega)$")
+# ax.legend();
+```
 
-bathL = LorentzianPadeBath(
-    bath_L.Q, bath_L.gamma, bath_L.W, bath_L.mu, bath_L.T,
-    Nk, tag="L",
-)
-bathR = LorentzianPadeBath(
-    bath_R.Q, bath_R.gamma, bath_R.W, bath_R.mu, bath_R.T,
-    Nk, tag="R",
-)
+```{code-cell} ipython3
+# plt.plot(w_list,envL.power_spectrum_plus(w_list)-result)
+# plt.ylim(-1e-8,1e-8)
+```
 
-with timer("RHS construction time"):
-    solver_pade = HEOMSolver(H, [bathL, bathR], max_depth=2, options=options)
+```{code-cell} ipython3
+tk = np.linspace(-60, 60, 1000)
 
-with timer("ODE solver time"):
-    result_pade = solver_pade.run(rho0, tlist)
 
-with timer("Steady state solver time"):
-    rho_ss_pade, ado_ss_pade = solver_pade.steady_state()
-```
+ig, ax = plt.subplots(figsize=(12, 7))
 
-Now let us plot the result which shows the decay of the initially excited impurity. This is not very illuminating, but we will compare it with the Matsubara expansion and analytic solution sortly:
+# Left lead emission and absorption
 
-```{code-cell} ipython3
-# Plot the Pade results
-fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 8))
+gam_L_in = envL.correlation_function_plus(tk) .imag
+gam_L_out = envL.correlation_function_minus(tk)
 
-axes.plot(
-    tlist, expect(result_pade.states, rho0),
-    'r--', linewidth=2,
-    label="P11 (Pade)",
-)
-axes.axhline(
-    expect(rho_ss_pade, rho0),
-    color='r', linestyle="dotted", linewidth=1,
-    label="P11 (Pade steady state)",
+ax.plot(
+    tk, gam_L_in,
+    "b", linewidth=3,
+    label=r"S_L(w) input (absorption)",
 )
+# ax.plot(
+#     w_list, gam_L_out,
+#     "r", linewidth=3,
+#     label=r"S_L(w) output (emission)",
+# )
 
-axes.set_xlabel('t', fontsize=28)
-axes.legend(fontsize=12);
-```
+# Right lead emission and absorption
 
-Now let us do the same for the Matsubara expansion:
+gam_R_in = fpenvL.correlation_function_plus(tk).imag
+# gam_R_out = fenvL.correlation_function_minus(w_list)
 
-```{code-cell} ipython3
-# HEOM dynamics using the Matsubara approximation:
-
-bathL = LorentzianBath(
-    bath_L.Q, bath_L.gamma, bath_L.W, bath_L.mu, bath_L.T,
-    Nk, tag="L",
-)
-bathR = LorentzianBath(
-    bath_R.Q, bath_R.gamma, bath_R.W, bath_R.mu, bath_R.T,
-    Nk, tag="R",
+ax.plot(
+    tk, gam_R_in,
+    "r--", linewidth=3,
+    label=r"S_R(w) input (absorption)",
 )
+# ax.plot(
+#     w_list, gam_R_out,
+#     "b--", linewidth=3,
+#     label=r"S_R(w) output (emission)",
+# )
 
-with timer("RHS construction time"):
-    solver_mats = HEOMSolver(H, [bathL, bathR], max_depth=2, options=options)
+ax.set_xlabel("w")
+ax.set_ylabel(r"$S(\omega)$")
+ax.legend();
+```
 
-with timer("ODE solver time"):
-    result_mats = solver_mats.run(rho0, tlist)
+```{code-cell} ipython3
+w_list = np.linspace(0, 60, 1000)
 
-with timer("Steady state solver time"):
-    rho_ss_mats, ado_ss_mats = solver_mats.steady_state()
-```
+fig, ax = plt.subplots(figsize=(12, 7))
 
-We see a marked difference in the Matsubara vs Pade results:
+# Left lead emission and absorption
 
-```{code-cell} ipython3
-# Plot the Pade results
-fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 8))
+gam_L_in = envL.correlation_function_plus(w_list).imag 
+gam_L_out = envL.correlation_function_minus(w_list).imag
 
-axes.plot(
-    tlist, expect(result_pade.states, rho0),
-    'r--', linewidth=2,
-    label="P11 (Pade)",
-)
-axes.axhline(
-    expect(rho_ss_pade, rho0),
-    color='r', linestyle="dotted", linewidth=1,
-    label="P11 (Pade steady state)",
+ax.plot(
+    w_list, gam_L_in,
+    "b", linewidth=3,
+    label=r"S_L(w) input (absorption)",
 )
+# ax.plot(
+#     w_list, gam_L_out,
+#     "r", linewidth=3,
+#     label=r"S_L(w) output (emission)",
+# )
 
-axes.plot(
-    tlist, expect(result_mats.states, rho0),
-    'b--', linewidth=2,
-    label="P11 (Mats)",
-)
-axes.axhline(
-    expect(rho_ss_mats, rho0),
-    color='b', linestyle="dotted", linewidth=1,
-    label="P11 (Mats steady state)",
+# Right lead emission and absorption
+
+gam_R_in = fpenvL.correlation_function_plus(w_list).imag
+#gam_R_out = fenvL.correlation_function_minus(w_list).imag
+
+ax.plot(
+    w_list, gam_R_in,
+    "r--", linewidth=3,
+    label=r"S_R(w) input (absorption)",
 )
+# ax.plot(
+#     w_list, gam_R_out,
+#     "b--", linewidth=3,
+#     label=r"S_R(w) output (emission)",
+# )
 
-axes.set_xlabel('t', fontsize=28)
-axes.legend(fontsize=12);
+ax.set_xlabel("w")
+ax.set_ylabel(r"$S(\omega)$")
+ax.legend();
 ```
 
 But which is more correct? The Matsubara or the Pade result?
@@ -383,15 +416,19 @@ One advantage of this simple model is that the steady state current to the baths
 See the [QuTiP-BoFiN paper](https://arxiv.org/abs/2010.10806) for a detailed description and references for the analytic result. Below we just perform the required integration numerically.
 
 ```{code-cell} ipython3
+from qutip.utilities import fermi_dirac
+def lamshift(self, w):
+    """ Return the lamshift. """
+    return 0.5 * (w - self.mu) * self.spectral_density(w) / self.W
 def analytical_steady_state_current(bath_L, bath_R, e1):
     """ Calculate the analytical steady state current. """
 
     def integrand(w):
         return (2 / np.pi) * (
-            bath_L.J(w) * bath_R.J(w) * (bath_L.fF(w) - bath_R.fF(w)) /
+            bath_L.spectral_density(w) * bath_R.spectral_density(w) * (fermi_dirac(w,1/envL.T,envL.mu) - fermi_dirac(w,1/envR.T,envR.mu)) /
             (
-                (bath_L.J(w) + bath_R.J(w))**2 +
-                4*(w - e1 - bath_L.lamshift(w) - bath_R.lamshift(w))**2
+                (bath_L.spectral_density(w) + bath_R.spectral_density(w))**2 +
+                4*(w - e1 - lamshift(envL,w) - lamshift(envL,w))**2
             )
         )
 
@@ -403,15 +440,15 @@ def analytical_steady_state_current(bath_L, bath_R, e1):
 
     # in principle the bounds for the integral should be rechecked if
     # bath or system parameters are changed substantially:
-    bounds = [-10, 10]
+    bounds = [-100, 100]
 
-    real_integral, _ = quad(real_part, *bounds)
-    imag_integral, _ = quad(imag_part, *bounds)
+    real_integral, _ = quad(real_part, *bounds,epsrel=1e-10,epsabs=1e-10)
+    imag_integral, _ = quad(imag_part, *bounds,epsrel=1e-10,epsabs=1e-10)
 
     return real_integral + 1.0j * imag_integral
 
 
-curr_ss_analytic = analytical_steady_state_current(bath_L, bath_R, e1)
+curr_ss_analytic = analytical_steady_state_current(envL, envR, e1)
 
 print(f"Analytical steady state current: {curr_ss_analytic}")
 ```
@@ -444,6 +481,24 @@ def state_current(ado_state, bath_tag):
 
 Now we can calculate the steady state currents from the Pade and Matsubara HEOM results:
 
+```{code-cell} ipython3
+# Times to solve for and initial system state:
+tlist = np.linspace(0, 100, 1000)
+rho0 = basis(2, 0) * basis(2, 0).dag()
+Nk=10
+penvL=envL.approx_by_pade(Nk=Nk,tag="L")
+penvR=envR.approx_by_pade(Nk=Nk,tag="R")
+
+with timer("RHS construction time"):
+    solver_pade = HEOMSolver(H, [(penvL,d1), (penvR,d1)], max_depth=2, options=options)
+
+with timer("ODE solver time"):
+    result_pade = solver_pade.run(rho0, tlist)
+
+with timer("Steady state solver time"):
+    rho_ss_pade, ado_ss_pade = solver_pade.steady_state()
+```
+
 ```{code-cell} ipython3
 curr_ss_pade_L = state_current(ado_ss_pade, "L")
 curr_ss_pade_R = state_current(ado_ss_pade, "R")
@@ -452,6 +507,24 @@ print(f"Pade steady state current (L): {curr_ss_pade_L}")
 print(f"Pade steady state current (R): {curr_ss_pade_R}")
 ```
 
+```{code-cell} ipython3
+# Times to solve for and initial system state:
+tlist = np.linspace(0, 100, 1000)
+rho0 = basis(2, 0) * basis(2, 0).dag()
+Nk=10
+menvL=envL.approx_by_matsubara(Nk=10,tag="L")
+menvR=envR.approx_by_matsubara(Nk=10,tag="R")
+
+with timer("RHS construction time"):
+    solver_mats = HEOMSolver(H, [(menvL,d1), (menvR,d1)], max_depth=2, options=options)
+
+with timer("ODE solver time"):
+    result_mats = solver_mats.run(rho0, tlist)
+
+with timer("Steady state solver time"):
+    rho_ss_mats, ado_ss_mats= solver_mats.steady_state()
+```
+
 ```{code-cell} ipython3
 curr_ss_mats_L = state_current(ado_ss_mats, "L")
 curr_ss_mats_R = state_current(ado_ss_mats, "R")
@@ -460,6 +533,81 @@ print(f"Matsubara steady state current (L): {curr_ss_mats_L}")
 print(f"Matsubara steady state current (R): {curr_ss_mats_R}")
 ```
 
+```{code-cell} ipython3
+# Times to solve for and initial system state:
+tlist = np.linspace(0, 100, 1000)
+rho0 = basis(2, 0) * basis(2, 0).dag()
+#tk=np.linspace(0,300,1_000)
+tk=np.linspace(0,300,2_000)
+#tk=np.linspace(0,5000,6000)
+# wk=np.concatenate((-np.logspace(1,-8,1000),np.logspace(-8,1,1000)))
+k=11
+fpenvL=envL._approx_by_prony(method="espira-I",tlist=tk,Np=k,Nm=k,tag="L")
+fpenvR=envR._approx_by_prony(method="espira-I",tlist=tk,Np=k,Nm=k,tag="R")
+# fpenvL=envL._approx_by_aaa(method="aaa",wlist=wk,Np_max=k,Nm_max=k,tag="L")
+# fpenvR=envR._approx_by_aaa(method="aaa",wlist=wk,Np_max=k,Nm_max=k,tag="R")
+with timer("RHS construction time"):
+    solver_fit = HEOMSolver(H, [(fpenvL,d1), (fpenvR,d1)], max_depth=2, options=options)
+
+with timer("ODE solver time"):
+    result_fit = solver_fit.run(rho0, tlist)
+
+with timer("Steady state solver time"):
+    rho_ss_fit, ado_ss_fit = solver_fit.steady_state()
+```
+
+```{code-cell} ipython3
+rho_ss_fit
+```
+
+```{code-cell} ipython3
+rho_ss_pade
+```
+
+```{code-cell} ipython3
+fig, ax = plt.subplots(figsize=(12, 7))
+ax.plot(
+    w_list,envL.correlation_function_plus(w_list).imag - fpenvL.correlation_function_plus(w_list).imag,
+    "r", linewidth=3,
+    label=r"J_R(w)",
+)
+ax.plot(
+    w_list,envL.correlation_function_plus(w_list).imag - penvL.correlation_function_plus(w_list).imag,
+    "b", linewidth=3,
+    label=r"J_R(w)",
+)
+
+ax.set_xlabel("w")
+ax.set_ylabel(r"$J(\omega)$")
+ax.legend();
+```
+
+```{code-cell} ipython3
+fig, ax = plt.subplots(figsize=(12, 7))
+ax.plot(
+    w_list,envL.correlation_function_plus(w_list) - fpenvL.correlation_function_plus(w_list),
+    "r", linewidth=3,
+    label=r"J_R(w)",
+)
+ax.plot(
+    w_list,envL.correlation_function_plus(w_list)  - penvL.correlation_function_plus(w_list),
+    "b", linewidth=3,
+    label=r"J_R(w)",
+)
+
+ax.set_xlabel("w")
+ax.set_ylabel(r"$J(\omega)$")
+ax.legend();
+```
+
+```{code-cell} ipython3
+curr_ss_p_L = state_current(ado_ss_fit, "L")
+curr_ss_p_R = state_current(ado_ss_fit, "R")
+
+print(f"Prony steady state current (L): {curr_ss_p_L}")
+print(f"Prony steady state current (R): {curr_ss_p_R}")
+```
+
 Note that the currents from each bath balance as is required by the steady state, but the value of the current is different for the Pade and Matsubara results.
 
 Now let's compare all three:
@@ -467,22 +615,56 @@ Now let's compare all three:
 ```{code-cell} ipython3
 print(f"Pade current (R): {curr_ss_pade_R}")
 print(f"Matsubara current (R): {curr_ss_mats_R}")
+print(f"Fit (R): {curr_ss_p_R}")
+
 print(f"Analytical curernt: {curr_ss_analytic}")
 ```
 
-In this case we observe that the Pade approximation has converged more closely to the analytical current than the Matsubara.
+```{code-cell} ipython3
+# Plot the Pade results
+fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 8))
 
-The Matsubara result could be improved by increasing the number of terms retained in the Matsubara expansion (i.e. increasing `Nk`).
+axes.plot(
+    tlist, expect(result_pade.states, rho0),
+    'r', linewidth=2,
+    label="P11 (Pade)",
+)
+axes.plot(
+    tlist, expect(result_mats.states, rho0),
+    'g--', linewidth=2,
+    label="P11 (Mats)",
+)
 
-+++
+axes.plot(
+    tlist, expect(result_fit.states, rho0),
+    'b--', linewidth=2,
+    label="P11 (Mats)",
+)
+axes.axhline(
+    expect(rho_ss_mats, rho0),
+    color='g', linestyle="dashdot", linewidth=1,
+    label="P11 (Mats steady state)",
+)
+axes.axhline(
+    expect(rho_ss_pade, rho0),
+    color='r',linewidth=2,
+    label="P11 (Pade steady state)",
+)
+axes.axhline(
+    expect(rho_ss_fit, rho0),
+    color='b', linestyle="-.", linewidth=2,
+    label="P11 (Fit steady state)",
+)
 
-## Current as a function of bias voltage
 
-+++
 
-Now lets plot the current as a function of bias voltage (the bias voltage is the parameter `theta` for the two baths).
+axes.set_xlabel('t', fontsize=28)
+axes.legend(fontsize=12);
+```
 
-We will calculate the steady state current for each `theta` both analytically and using the HEOM with the Pade correlation expansion approximation.
+```{code-cell} ipython3
+assert 1==0
+```
 
 ```{code-cell} ipython3
 # Theta (bias voltages)
@@ -496,33 +678,61 @@ display(progress)
 
 # Calculate the current for the list of thetas
 
+def lamshift(w,theta,self):
+    """ Return the lamshift. """
+    return 0.5 * (w - theta) * self.spectral_density(w) / self.W
+def analytical_steady_state_current(bath_L, bath_R, e1,theta=2):
+    """ Calculate the analytical steady state current. """
 
-def current_analytic_for_theta(e1, bath_L, bath_R, theta):
+    def integrand(w):
+        return (2 / np.pi) * (
+            bath_L.spectral_density(w) * bath_R.spectral_density(w) * (fermi_dirac(w,1/envL.T,theta/2) - fermi_dirac(w,1/envR.T,-theta/2)) /
+            (
+                (bath_L.spectral_density(w) + bath_R.spectral_density(w))**2 +
+                4*(w - e1 - lamshift(w,theta/2,envR) - lamshift(w,-theta/2,envL))**2
+            )
+        )
+
+    def real_part(x):
+        return np.real(integrand(x))
+
+    def imag_part(x):
+        return np.imag(integrand(x))
+
+    # in principle the bounds for the integral should be rechecked if
+    # bath or system parameters are changed substantially:
+    bounds = [-10, 10]
+
+    real_integral, _ = quad(real_part, *bounds)
+    imag_integral, _ = quad(imag_part, *bounds)
+
+    return real_integral + 1.0j * imag_integral
+
+def current_analytic_for_theta(e1, theta):
     """ Return the analytic current for a given theta. """
+    envL=LorentzianEnvironment(T= 0.025851991,W=1,mu=theta/2,gamma=0.01,Nk=20)
+    envR=LorentzianEnvironment(T= 0.025851991,W=1,mu=-theta/2,gamma=0.01,Nk=20)
     current = analytical_steady_state_current(
-        bath_L.replace(theta=theta),
-        bath_R.replace(theta=theta),
+        envL,
+        envR,
         e1,
+        theta
     )
     progress.value += 1
     return np.real(current)
 
 
-def current_pade_for_theta(H, bath_L, bath_R, theta, Nk):
+def current_pade_for_theta(H, theta, Nk,method):
     """ Return the steady state current using the Pade approximation. """
-    bath_L = bath_L.replace(theta=theta)
-    bath_R = bath_R.replace(theta=theta)
-
-    bathL = LorentzianPadeBath(
-        bath_L.Q, bath_L.gamma, bath_L.W, bath_L.mu, bath_L.T,
-        Nk, tag="L",
-    )
-    bathR = LorentzianPadeBath(
-        bath_R.Q, bath_R.gamma, bath_R.W, bath_R.mu, bath_R.T,
-        Nk, tag="R",
-    )
-
-    solver_pade = HEOMSolver(H, [bathL, bathR], max_depth=2, options=options)
+    envL=LorentzianEnvironment(T= 0.025851991,W=1,mu=theta/2,gamma=0.01,Nk=20)
+    envR=LorentzianEnvironment(T= 0.025851991,W=1,mu=-theta/2,gamma=0.01,Nk=20)
+    if method=="pade":
+        bathL=envL.approx_by_pade(Nk=Nk,tag="L")
+        bathR=envR.approx_by_pade(Nk=Nk,tag="R")
+    else:
+        bathL=envL._approx_by_prony(method="espira-I",tlist=tk,Np=Nk,Nm=Nk,tag="L")
+        bathR=envR._approx_by_prony(method="espira-I",tlist=tk,Np=Nk,Nm=Nk,tag="R")
+    solver_pade = HEOMSolver(H, [(bathL,d1), (bathR,d1)], max_depth=2, options=options)
     rho_ss_pade, ado_ss_pade = solver_pade.steady_state()
     current = state_current(ado_ss_pade, bath_tag="R")
 
@@ -531,23 +741,58 @@ def current_pade_for_theta(H, bath_L, bath_R, theta, Nk):
 
 
 curr_ss_analytic_thetas = [
-    current_analytic_for_theta(e1, bath_L, bath_R, theta)
+    current_analytic_for_theta(e1, theta)
     for theta in thetas
 ]
 
 # The number of expansion terms has been dropped to Nk=6 to speed
 # up notebook execution. Increase to Nk=10 for more accurate results.
+# curr_ss_pade_theta = [
+#     current_pade_for_theta(H, theta, Nk=6,method="pade")
+#     for theta in thetas
+# ]
+```
+
+```{code-cell} ipython3
 curr_ss_pade_theta = [
-    current_pade_for_theta(H, bath_L, bath_R, theta, Nk=6)
+    current_pade_for_theta(H, theta, Nk=6,method="pade")
     for theta in thetas
 ]
 ```
 
+```{code-cell} ipython3
+from tqdm import tqdm
+```
+
+```{code-cell} ipython3
+curr_ss_fit_theta = [
+    current_pade_for_theta(H, theta, Nk=11,method="pde")
+    for theta in tqdm(thetas)
+]
+```
+
+In this case we observe that the Pade approximation has converged more closely to the analytical current than the Matsubara.
+
+The Matsubara result could be improved by increasing the number of terms retained in the Matsubara expansion (i.e. increasing `Nk`).
+
++++
+
+## Current as a function of bias voltage
+
++++
+
+Now lets plot the current as a function of bias voltage (the bias voltage is the parameter `theta` for the two baths).
+
+We will calculate the steady state current for each `theta` both analytically and using the HEOM with the Pade correlation expansion approximation.
+
++++
+
 Below we plot the results and see that even with `Nk=6`, the HEOM Pade approximation gives good results for the steady state current. Increasing `Nk` to `10` gives very accurate results.
 
 ```{code-cell} ipython3
 fig, ax = plt.subplots(figsize=(12, 7))
 
+
 ax.plot(
     thetas, 2.434e-4 * 1e6 * np.array(curr_ss_analytic_thetas),
     color="black", linewidth=3,
@@ -556,10 +801,205 @@ ax.plot(
 ax.plot(
     thetas, 2.434e-4 * 1e6 * np.array(curr_ss_pade_theta),
     'r--', linewidth=3,
-    label=r"HEOM Pade $N_k=10$, $n_{\mathrm{max}}=2$",
+    label=rf"HEOM Pade $N_k={6}$,"+r"$n_{\mathrm{max}}=2$",
+)
+ax.plot(
+    thetas, 2.434e-4 * 1e6 * np.array(curr_ss_fit_theta),
+    'b--', linewidth=3,
+    label=rf"HEOM Fit $N_k={6}$,"+r"$n_{\mathrm{max}}=2$",
 )
 
+ax.locator_params(axis='y', nbins=4)
+ax.locator_params(axis='x', nbins=4)
+
+ax.set_xticks([-2.5, 0, 2.5])
+ax.set_xticklabels([-2.5, 0, 2.5])
+ax.set_xlabel(r"Bias voltage $\Delta \mu$ ($V$)", fontsize=28)
+ax.set_ylabel(r"Current ($\mu A$)", fontsize=28)
+ax.legend(fontsize=25);
+```
+
+Using a semicircular spectral density
+
+```{code-cell} ipython3
+def semicircle(w,W):
+    result=np.zeros(len(w))
+    mask=W>=w
+    mask2=-W<=w
+    result[mask&mask2]=np.sqrt(1-(w[mask&mask2]/W)**2)
+    return result
+    
+```
+
+```{code-cell} ipython3
+w_list=np.linspace(-300,300,10_000)
+plt.plot(w_list,semicircle(w_list,30))
+```
+
+```{code-cell} ipython3
+J=semicircle(w_list,30)
+```
+
+```{code-cell} ipython3
+g=0.01
+w_list=np.linspace(-300,300,10_000)
+fenvL=FermionicEnvironment.from_spectral_density(g*semicircle(w_list,10),w_list,T=0.025851991,mu=1)
+fenvR=FermionicEnvironment.from_spectral_density(g*semicircle(w_list,10),w_list,T=0.025851991,mu=-1)
+```
+
+```{code-cell} ipython3
+# Times to solve for and initial system state:
+tlist = np.linspace(0, 30, 1000)
+rho0 = basis(2, 0) * basis(2, 0).dag()
+tk=np.linspace(0,200,1_000)
+#tk=np.linspace(0,50,1_000)
+#tk=np.linspace(0,5000,6000)
+k=12
+# fpenvL=fenvL._approx_by_prony(method="esprit",tlist=tk,Np=k,Nm=6,tag="L")
+fpenvL=fenvL._approx_by_prony(method="espira-I",tlist=tk,Np=k,Nm=k,tag="L")
+fpenvR=fenvR._approx_by_prony(method="espira-I",tlist=tk,Np=k,Nm=k,tag="R")
+
+with timer("RHS construction time"):
+    solver_fit = HEOMSolver(H, [(fpenvL,d1), (fpenvR,d1)], max_depth=2, options=options)
+
+# with timer("ODE solver time"):
+#     result_fit = solver_fit.run(rho0, tlist)
 
+with timer("Steady state solver time"):
+    rho_ss_fit, ado_ss_fit = solver_fit.steady_state()
+```
+
+```{code-cell} ipython3
+rho_ss_fit
+```
+
+```{code-cell} ipython3
+
+curr_ss_p_L = state_current(ado_ss_fit, "L")
+curr_ss_p_R = state_current(ado_ss_fit, "R")
+print(f"Fit (R): {curr_ss_p_R}")
+```
+
+```{code-cell} ipython3
+curr_ss_analytic = analytical_steady_state_current(fenvL, fenvR, e1)
+
+print(f"Analytical steady state current: {curr_ss_analytic}")
+```
+
+```{code-cell} ipython3
+plt.plot(tlist,fenvL.correlation_function_plus(tlist).real)
+plt.plot(tlist,fpenvL.correlation_function_plus(tlist).real,"--")
+plt.show()
+plt.plot(tlist,fenvL.correlation_function_plus(tlist).imag)
+plt.plot(tlist,fpenvL.correlation_function_plus(tlist).imag,"--")
+```
+
+```{code-cell} ipython3
+plt.plot(tlist,fenvL.correlation_function_minus(tlist).real)
+plt.plot(tlist,fpenvL.correlation_function_minus(tlist).real,"--")
+plt.show()
+plt.plot(tlist,fenvL.correlation_function_minus(tlist).imag)
+plt.plot(tlist,fpenvL.correlation_function_minus(tlist).imag,"--")
+```
+
+```{code-cell} ipython3
+w_list=np.linspace(-20,20,10_000)
+```
+
+```{code-cell} ipython3
+plt.plot(w_list,fenvL.power_spectrum_plus(w_list).real)
+plt.plot(w_list,fpenvL.power_spectrum_plus(w_list).real,"--")
+```
+
+```{code-cell} ipython3
+plt.plot(w_list,fenvL.power_spectrum_minus(w_list).real)
+plt.plot(w_list,fpenvL.power_spectrum_minus(w_list).real,"--")
+```
+
+```{code-cell} ipython3
+w_list=np.linspace(-100,100,10_000)
+plt.plot(w_list,g*semicircle(w_list,10))
+plt.plot(w_list,fenvL.spectral_density(w_list),"--")
+plt.plot(w_list,fpenvL.spectral_density(w_list),"--")
+# plt.xlim(50,70)
+# plt.ylim(-1e-12,1e-12)
+```
+
+```{code-cell} ipython3
+curr_ss_analytic = analytical_steady_state_current(fenvL, fenvR, e1)
+
+print(f"Analytical steady state current: {curr_ss_analytic}")
+```
+
+```{code-cell} ipython3
+def current_analytic_for_theta(e1, theta):
+    """ Return the analytic current for a given theta. """
+    fenvL=FermionicEnvironment.from_spectral_density(g*semicircle(w_list,3),w_list,T=0.025851991,mu=theta/2)
+    fenvR=FermionicEnvironment.from_spectral_density(g*semicircle(w_list,3),w_list,T=0.025851991,mu=-theta/2)
+    current = analytical_steady_state_current(
+        fenvL,
+        fenvR,
+        e1,
+        theta
+    )
+    progress.value += 1
+    return np.real(current)
+
+
+def current_pade_for_theta(H, theta, Nk,method):
+    """ Return the steady state current using the Pade approximation. """
+    fenvL=FermionicEnvironment.from_spectral_density(g*semicircle(w_list,3),w_list,T=0.025851991,mu=theta/2)
+    fenvR=FermionicEnvironment.from_spectral_density(g*semicircle(w_list,3),w_list,T=0.025851991,mu=-theta/2)
+    if method=="pade":
+        bathL=fenvL.approx_by_pade(Nk=Nk,tag="L")
+        bathR=fenvR.approx_by_pade(Nk=Nk,tag="R")
+    else:
+        bathL=fenvL._approx_by_prony(method="espira-I",tlist=tk,Np=Nk,Nm=Nk,tag="L")
+        bathR=fenvR._approx_by_prony(method="espira-I",tlist=tk,Np=Nk,Nm=Nk,tag="R")
+    solver_pade = HEOMSolver(H, [(bathL,d1), (bathR,d1)], max_depth=2, options=options)
+    rho_ss_pade, ado_ss_pade = solver_pade.steady_state()
+    current = state_current(ado_ss_pade, bath_tag="R")
+
+    progress.value += 1
+    return np.real(current)
+
+
+curr_ss_analytic_thetas = [
+    current_analytic_for_theta(e1, theta)
+    for theta in thetas
+]
+```
+
+```{code-cell} ipython3
+from tqdm import tqdm
+```
+
+```{code-cell} ipython3
+curr_ss_fit_theta = [
+    current_pade_for_theta(H, theta, Nk=5,method="pde")
+    for theta in tqdm(thetas)
+]
+```
+
+```{code-cell} ipython3
+fig, ax = plt.subplots(figsize=(12, 7))
+
+
+ax.plot(
+    thetas, 2.434e-4 * 1e6 * np.array(curr_ss_analytic_thetas),
+    color="black", linewidth=3,
+    label=r"Analytical",
+)
+ax.plot(
+    thetas, 2.434e-4 * 1e6 * np.array(curr_ss_fit_theta),
+    'b--', linewidth=3,
+    label=rf"HEOM Fit $N_k={6}$,"+r"$n_{\mathrm{max}}=2$",
+)
+# ax.plot(
+#     thetas, 2.434e-4 * 1e6 * np.array(curr_ss_pade_theta),
+#     'r--', linewidth=3,
+#     label=rf"HEOM Fit $N_k={6}$,"+r"$n_{\mathrm{max}}=2$",
+# )
 ax.locator_params(axis='y', nbins=4)
 ax.locator_params(axis='x', nbins=4)
 
diff --git a/tutorials-v5/heom/heom-5b-fermions-discrete-boson-model.md b/tutorials-v5/heom/heom-5b-fermions-discrete-boson-model.md
index bc2e23e5..de466f2a 100644
--- a/tutorials-v5/heom/heom-5b-fermions-discrete-boson-model.md
+++ b/tutorials-v5/heom/heom-5b-fermions-discrete-boson-model.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.16.1
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -95,7 +95,7 @@ The complete setup now consists of four parts:
 
 ## Setup
 
-```{code-cell} ipython3
+```{code-cell}
 import contextlib
 import dataclasses
 import time
@@ -122,7 +122,7 @@ from IPython.display import display
 
 ## Helpers
 
-```{code-cell} ipython3
+```{code-cell}
 @contextlib.contextmanager
 def timer(label):
     """ Simple utility for timing functions:
@@ -136,7 +136,7 @@ def timer(label):
     print(f"{label}: {end - start}")
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def state_current(ado_state, bath_tag):
     """ Determine current from the given bath (either "R" or "L") to
         the system in the given ADO state.
@@ -158,7 +158,7 @@ def state_current(ado_state, bath_tag):
     )
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Solver options:
 
 # We set store_ados to True so that we can
@@ -181,7 +181,7 @@ options = {
 
 Let us set up the system Hamiltonian and specify the properties of the two reservoirs.
 
-```{code-cell} ipython3
+```{code-cell}
 # Define the system Hamiltonian:
 
 @dataclasses.dataclass
@@ -208,7 +208,7 @@ class SystemParameters:
 sys_p = SystemParameters()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Define parameters for left and right fermionic baths.
 # Each bath is a lead (i.e. a wire held at a potential)
 # with temperature T and chemical potential mu.
@@ -262,7 +262,7 @@ bath_R = LorentzianBathParameters(W=10**4, lead="R")
 
 Next let's plot the emission and absorption by the leads.
 
-```{code-cell} ipython3
+```{code-cell}
 w_list = np.linspace(-2, 2, 100)
 
 fig, ax = plt.subplots(figsize=(12, 7))
@@ -310,7 +310,7 @@ Here we just give one example of the current as a function of bias voltage, but
 
 One note:  for very large problems, this can be slow.
 
-```{code-cell} ipython3
+```{code-cell}
 def steady_state_pade_for_theta(sys_p, bath_L, bath_R, theta, Nk, Nc, Nbos):
     """ Return the steady state current using the Pade approximation. """
 
@@ -336,7 +336,7 @@ def steady_state_pade_for_theta(sys_p, bath_L, bath_R, theta, Nk, Nc, Nbos):
     return np.real(2.434e-4 * 1e6 * current)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Parameters:
 
 Nk = 6
@@ -360,7 +360,7 @@ for theta in thetas:
     progress.value += 1
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 fig, ax = plt.subplots(figsize=(12, 10))
 
 ax.plot(
@@ -382,7 +382,7 @@ ax.legend(loc=4);
 
 ## About
 
-```{code-cell} ipython3
+```{code-cell}
 qutip.about()
 ```
 
@@ -390,6 +390,6 @@ qutip.about()
 
 This section can include some tests to verify that the expected outputs are generated within the notebook. We put this section at the end of the notebook, so it's not interfering with the user experience. Please, define the tests using assert, so that the cell execution fails if a wrong output is generated.
 
-```{code-cell} ipython3
+```{code-cell}
 assert 1 == 1
 ```
diff --git a/tutorials-v5/heom/heom-6-parity.md b/tutorials-v5/heom/heom-6-parity.md
new file mode 100644
index 00000000..f7677742
--- /dev/null
+++ b/tutorials-v5/heom/heom-6-parity.md
@@ -0,0 +1,244 @@
+---
+jupyter:
+  jupytext:
+    text_representation:
+      extension: .md
+      format_name: markdown
+      format_version: '1.3'
+      jupytext_version: 1.16.1
+  kernelspec:
+    display_name: Python 3 (ipykernel)
+    language: python
+    name: python3
+---
+
+```python
+import numpy as np
+import matplotlib.pyplot as plt
+from qutip import tensor,destroy,qeye,spre,spost,operator_to_vector,expect,qeye
+from qutip.solver.heom import FermionicBath,LorentzianPadeBath,HEOMSolver
+import scipy.sparse as sp
+from scipy.sparse.linalg import lgmres
+```
+
+<!-- #region -->
+# Kondo  Peak
+
+In this tutorial we will find the kondo peak for a system of two spins connected to two fermionic reservoirs 
+
+The Hamiltonian of this setup is given by:
+
+$H_{T}=H_{S}+H_{f}+H_{ef}$
+
+Where $H_{S}$ is the system Hamiltonian, which is divided into three contributions. The electron system hamiltonian ($H_{e}$), the single-mode cavity ($H_{c}$) and the light matter coupling ($H_{ef}$):
+
+$H_{S} = H_{e} + H_{ec}$
+
+
+$H_{e}=  \sum_{n=g,e}\sum_{\sigma=\uparrow,\downarrow} \epsilon_{n} \hat{n}^{n \sigma} + U_{n} \hat{n}_{n \uparrow} \hat{n}_{n \downarrow} $
+
+$H_{ec} = \sum_{\sigma=\uparrow,\downarrow} g_{ct} (d_{g\sigma}^{\dagger} d_{e\sigma} + d_{e\sigma}^{\dagger} d_{g\sigma} )(a^{\dagger}+a)$
+
+
+The other two terms in the total Hamiltonian are given by
+
+$H_{f}= \sum_{\alpha} \sum_{k} \epsilon_{\alpha,k} c^{\dagger}_{\alpha,k}c_{\alpha,k}$
+
+and
+
+$H_{ef}= \sum_{k} \sum_{\alpha=L,R} \sum_{\sigma=\uparrow,\downarrow} g_{\alpha,k} (c_{\alpha,k}^{\dagger} d_{g\sigma}+c_{\alpha,k} d_{g\sigma}^{\dagger})$
+
+
+The interaction between the electronic systm and fermionic leads can be fully characterized by the Lorentzian spectral density
+
+$J_{f_{\alpha}}(\omega)=\frac{1}{2 \pi} \frac{\Gamma_{\alpha} W_{f}^{2}}{(\omega-\mu_{\alpha})^{2}+W_{f}^{2}}$
+
+
+<!-- #endregion -->
+
+We first define some utility functions to construct the Hailtonian and  the system's Liouvillian supper operator
+
+$\mathcal{L}(A)= - i [ H , A]$
+
+```python
+def _Is(i): return [qeye(2) for j in range(0, i)]
+def _sm(N, i): return tensor(_Is(i) + [destroy(2)] + _Is(N - i - 1))
+def _sz(N, i): return tensor(_Is(i) + [qeye(2)-2*destroy(2).dag() * destroy(2)] + _Is(N - i - 1))
+def _oprd(lst): return np.prod(np.array(lst, dtype=object))
+def _f(N, n): return _oprd([_sz(N, j) for j in range(n)])*_sm(N,n)
+def liouvillian(H):
+    return - 1j * (spre(H) - spost(H))
+```
+
+We define the system and bath parameters and  construct our fermionic baths using the `LorentzianPadeBath`  class from `Qutip`
+
+```python
+Ncc = 5
+Nk=2
+#########################
+### system: two fermions
+N = 2
+d_1 = _f(N,0)
+d_2 = _f(N,1)
+
+#bath params:
+mu = 0.  #chemical potential
+Gamma  = 1  #coupling strenght
+W = 2.5 #bath width
+
+#system params:
+#coulomb repulsion
+U = 3 * np.pi * Gamma
+#impurity energy
+w0 = - U / 2.
+
+beta = 1 / (0.2 * Gamma) # Inverse Temperature
+
+Qops = [d_1.dag(),d_1,d_2.dag(),d_2] # coupling operators 
+
+He = w0 *(d_1.dag() * d_1 + d_2.dag() * d_2) + U * d_1.dag() * d_1 * d_2.dag() * d_2
+
+L = liouvillian(He) 
+
+times = np.linspace(0,10,20000)
+PadeBath=LorentzianPadeBath(Q=sum(Qops),gamma=2*Gamma,w=W,mu=mu,T=1/beta,Nk=Nk)
+```
+
+Since we have two leads, we now extract the coefficients of the bath to contruct a bath for each lead
+
+```python
+coeffsplus=PadeBath._corr(2*Gamma,W,mu,1/beta,Nk,sigma=1) #absorption coefficients
+coeffsminus=PadeBath._corr(2*Gamma,W,mu,1/beta,Nk,sigma=-1) #emission coefficients
+```
+
+We then construct the baths from the lead coefficients
+
+```python
+bath1 = FermionicBath(d_1, coeffsplus[0], coeffsplus[1], coeffsminus[0], coeffsminus[1], tag ="Lead 1")
+bath2 = FermionicBath(d_2, coeffsplus[0], coeffsplus[1], coeffsminus[0], coeffsminus[1], tag ="Lead 2")
+```
+
+At this point we obtain the steady state of the system and obtain the expectation value of $d_{1}^{\dagger} d_{1}$ which should be around $\frac{1}{2}$ for standard parameters
+
+```python
+resultHEOMPade = HEOMSolver(L, [bath1,bath2], Ncc)  #<---- normal parity HEOM to get normal steadystate
+rhoss, fullss= resultHEOMPade.steady_state()
+expect(rhoss, d_1.dag()*d_1)
+```
+
+We now use the steady state to construct the density of states, using the Generator of the dynamics and the quantum regression theorem one has:
+
+$\langle d^{\dagger}(\tau) d(0) \rangle = Tr[d^{\dagger}(0) e^{\mathcal{L}\tau} \big(d(0) \rho_{ss}\big)]$
+
+If we define 
+
+$C_{R}(\tau)=\theta(\tau)(\langle d^{\dagger}(\tau) d(0) \rangle+\langle  d(0) d^{\dagger}(-\tau) \rangle)$
+
+$C_{A}(\tau)=-\theta(\tau)(\langle d^{\dagger}(\tau) d(0) \rangle+\langle  d(0) d^{\dagger}(-\tau) \rangle)$
+
+One can then write the density of states as 
+
+$A(\omega)=\frac{1}{2 \pi} \int_{-\infty}^{\infty} dt e^{i \omega t} (C_{R}(t)+C_{A}(t))$
+
+below there are some auxiliary functions that compute this quantities from the steady state of the system and the generator (the right hand side of HEOM). 
+
+```python
+
+wlist = np.linspace(-15,15,100)
+def prepare_matrices(result,fullss):
+    M_1 = result.rhs
+    sup_dim = result._sup_shape
+    N_he  = result._n_ados
+    rhoptemp = fullss._ado_state.reshape((N_he*sup_dim, 1))
+    unit_helems = sp.identity(N_he, format='csr')
+    M2=sp.csr_matrix(M_1.to_list()[0].full())
+    ###################################
+    d_1_big = sp.kron(unit_helems, sp.csr_matrix(spre(d_1).full()))
+    rho0d1 = np.array(d_1_big @ rhoptemp, dtype=complex)
+    # ######################################
+    d_1_bigdag = sp.kron(unit_helems, sp.csr_matrix(spre(d_1.dag()).full()))
+    rho0d1dag = np.array(d_1_bigdag @  rhoptemp, dtype=complex)
+    I = tensor([qeye(n) for n in d_1.dims[0]])
+    I_vec3 = sp.csr_matrix(operator_to_vector(I).full().T)
+    Nfull = N_he*sup_dim
+    c_I = sp.identity(Nfull)
+    return  I_vec3,c_I,rho0d1dag,rho0d1,sup_dim,M2, d_1_bigdag, d_1_big
+def D2(w,M2,c_I):
+    return (M2-1.0j*w*c_I)
+
+def D3(w,M2,c_I):
+    return (1.0j*w*c_I+M2)
+
+def density_of_states(wlist,result,fullss):
+    ddagd = []
+    dddag = []
+    I_vec3,c_I,rho0d1dag,rho0d1,sup_dim,M2, d_1_bigdag, d_1_big=prepare_matrices(result,fullss)
+    for idx in range(len(wlist)):
+        w=wlist[idx]
+        x2, _= lgmres(D2(w,M2,c_I), rho0d1,atol=1e-8)
+        Cw21 = d_1_bigdag  @ x2
+        Cw22 =  (I_vec3 @ Cw21[:sup_dim])
+        ddagd.append(Cw22)
+        
+        x3, _= lgmres(D3(w,M2,c_I),  rho0d1dag,atol=1e-8)
+        Cw31 = d_1_big @ x3
+        Cw32 =  (I_vec3 @ Cw31[:sup_dim])
+        dddag.append(Cw32)
+        
+    return -2*(np.array(ddagd).flatten()+np.array(dddag).flatten())
+```
+
+```python
+ddos=density_of_states(wlist,resultHEOMPade,fullss)
+```
+
+### We can now visualize our density of states
+
+```python
+
+plt.rcParams["figure.figsize"] = (12,10)
+plt.plot(wlist,np.real(ddos),linestyle='--',label=r"HEOM Padé $N_k = 2$",linewidth=4)
+plt.legend(fontsize=20,loc = (0.23,0.7))
+
+plt.xlim(-10,15)
+plt.yticks([0.,1,2],[0,1,2])
+plt.xlabel(r"$\omega/\Gamma$",fontsize=28,labelpad=-10)
+plt.ylabel(r"$2\pi \Gamma A(\omega)$ ",fontsize=28)
+
+         
+plt.show()
+```
+
+Notice we don't see any peak as the RHS used for evolution without taking the parity into account, to clarify why this is wrong consider that the generator is now acting as
+
+$Tr[d^{\dagger}(0) e^{\mathcal{L}\tau} \big(d(0) \rho_{ss}\big)]$
+
+rather than the usual 
+
+$ e^{\mathcal{L}\tau} \rho$
+
+as the state of the system has even parity, when we apply creation or anihilation operatos on it we make it odd, and as such we need to evolve it with an odd parity generator
+
+To take the parity into account we just set the odd_parity argument to True on the HEOMSolver, we now repeat the calculations
+
+```python
+
+resultHEOMPade2 = HEOMSolver(L, [bath1,bath2], Ncc, odd_parity=True) #<------ use ODD parity 
+ddos=density_of_states(wlist,resultHEOMPade2,fullss) # recalculate density of states with odd parity
+```
+
+```python
+
+
+plt.rcParams["figure.figsize"] = (12,10)
+plt.plot(wlist,np.real(ddos),linestyle='--',label=r"HEOM Padé $N_k = 2$",linewidth=4)
+plt.legend(fontsize=20,loc = (0.23,0.7))
+
+plt.xlim(-10,15)
+plt.yticks([0.,1,2],[0,1,2])
+plt.xlabel(r"$\omega/\Gamma$",fontsize=28,labelpad=-10)
+plt.ylabel(r"$2\pi \Gamma A(\omega)$ ",fontsize=28)
+
+         
+plt.show()
+```
diff --git a/tutorials-v5/heom/heom-index.md b/tutorials-v5/heom/heom-index.md
index dcb1fa89..4bc08d9d 100644
--- a/tutorials-v5/heom/heom-index.md
+++ b/tutorials-v5/heom/heom-index.md
@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.4
+    jupytext_version: 1.16.1
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python