diff --git a/tutorials-v4/heom/heom-5a-fermions-single-impurity-model.md b/tutorials-v4/heom/heom-5a-fermions-single-impurity-model.md
index 623af4db..d56b4219 100644
--- a/tutorials-v4/heom/heom-5a-fermions-single-impurity-model.md
+++ b/tutorials-v4/heom/heom-5a-fermions-single-impurity-model.md
@@ -20,7 +20,7 @@ kernelspec:
 
 ## Introduction
 
-Here we model a single fermion coupled to two electronic leads or reservoirs (e.g.,  this can describe a single quantum dot, a molecular transistor, etc).  Note that in this implementation we primarily follow the definitions used by Christian Schinabeck in his dissertation https://opus4.kobv.de/opus4-fau/files/10984/DissertationChristianSchinabeck.pdf and related publications.
+Here we model a single fermion coupled to two electronic leads or reservoirs (e.g.,  this can describe a single quantum dot, a molecular transistor, etc).  Note that in this implementation we primarily follow the definitions used by Christian Schinabeck in his dissertation https://open.fau.de/items/36fdd708-a467-4b59-bf4e-4a2110fbc431 and related publications.
 
 Notation:
 
diff --git a/tutorials-v4/heom/heom-5b-fermions-discrete-boson-model.md b/tutorials-v4/heom/heom-5b-fermions-discrete-boson-model.md
index 19d63126..f2afc4f8 100644
--- a/tutorials-v4/heom/heom-5b-fermions-discrete-boson-model.md
+++ b/tutorials-v4/heom/heom-5b-fermions-discrete-boson-model.md
@@ -20,7 +20,7 @@ kernelspec:
 
 Here we model a single fermion coupled to two electronic leads or reservoirs (e.g.,  this can describe a single quantum dot, a molecular transistor, etc), also coupled to a discrete bosonic (vibronic) mode.
 
-Note that in this implementation we primarily follow the definitions used by Christian Schinabeck in his Dissertation https://opus4.kobv.de/opus4-fau/files/10984/DissertationChristianSchinabeck.pdf and related publications. In particular this example reproduces some results from https://journals.aps.org/prb/abstract/10.1103/PhysRevB.94.201407
+Note that in this implementation we primarily follow the definitions used by Christian Schinabeck in his Dissertation https://open.fau.de/items/36fdd708-a467-4b59-bf4e-4a2110fbc431 and related publications. In particular this example reproduces some results from https://journals.aps.org/prb/abstract/10.1103/PhysRevB.94.201407
 
 Notation:
 
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..292ff8aa 100644
--- a/tutorials-v5/heom/heom-1a-spin-bath-model-basic.md
+++ b/tutorials-v5/heom/heom-1a-spin-bath-model-basic.md
@@ -1,19 +1,19 @@
 ---
-jupyter:
-  jupytext:
-    text_representation:
-      extension: .md
-      format_name: markdown
-      format_version: '1.3'
-      jupytext_version: 1.14.5
-  kernelspec:
-    display_name: Python 3 (ipykernel)
-    language: python
-    name: python3
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.17.0
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
 ---
 
 # HEOM 1a: Spin-Bath model (introduction)
 
++++
 
 ## Introduction
 
@@ -66,77 +66,57 @@ density is given by:
     c_k = \begin{cases}
                \lambda \gamma (\cot(\beta \gamma / 2) - i) \
                & k = 0\\
-               4 \lambda \gamma \nu_k / \{(nu_k^2 - \gamma^2)\beta \} \
+               4 \lambda \gamma \nu_k / \{(\nu_k^2 - \gamma^2)\beta \} \
                & k \geq 1\\
            \end{cases}
 \end{equation*}
 
 Note that in the above, and the following, we set $\hbar = k_\mathrm{B} = 1$.
 
++++
 
 ## Setup
 
-```python
+```{code-cell} ipython3
 import contextlib
 import time
 
 import numpy as np
-from matplotlib import pyplot as plt
 from scipy.optimize import curve_fit
+from matplotlib import pyplot as plt
 
-import qutip
 from qutip import (
-    basis,
-    brmesolve,
-    destroy,
-    expect,
-    liouvillian,
-    qeye,
-    sigmax,
-    sigmaz,
-    spost,
-    spre,
-    tensor,
+    about, basis, brmesolve, destroy, expect, liouvillian,
+    qeye, sigmax, sigmaz, spost, spre, tensor
 )
-
-from qutip.solver.heom import (
-    BosonicBath,
-    DrudeLorentzBath,
-    DrudeLorentzPadeBath,
-    HEOMSolver,
-    HSolverDL,
+from qutip.core.environment import (
+    DrudeLorentzEnvironment, ExponentialBosonicEnvironment, system_terminator
 )
+from qutip.solver.heom import HEOMSolver, HSolverDL
 
 %matplotlib inline
 ```
 
-
 ## Helper functions
 
 Let's define some helper functions for calculating correlation function expansions, plotting results and timing how long operations take:
 
-```python
+```{code-cell} ipython3
 def cot(x):
     """Vectorized cotangent of x."""
     return 1.0 / np.tan(x)
 ```
 
-```python
+```{code-cell} ipython3
 def dl_matsubara_params(lam, gamma, T, nk):
     """Calculation of the real and imaginary expansions of the Drude-Lorenz
     correlation functions.
     """
     ckAR = [lam * gamma * cot(gamma / (2 * T))]
     ckAR.extend(
-        4
-        * lam
-        * gamma
-        * T
-        * 2
-        * np.pi
-        * k
-        * T
-        / ((2 * np.pi * k * T) ** 2 - gamma**2)
+        8 * lam * gamma * T * np.pi * k * T / (
+            (2 * np.pi * k * T) ** 2 - gamma**2
+        )
         for k in range(1, nk + 1)
     )
     vkAR = [gamma]
@@ -148,22 +128,7 @@ def dl_matsubara_params(lam, gamma, T, nk):
     return ckAR, vkAR, ckAI, vkAI
 ```
 
-```python
-def dl_corr_approx(t, nk):
-    """Drude-Lorenz correlation function approximation.
-
-    Approximates the correlation function at each time t to nk exponents.
-    """
-    c = lam * gamma * (-1.0j + cot(gamma / (2 * T))) * np.exp(-gamma * t)
-    for k in range(1, nk):
-        vk = 2 * np.pi * k * T
-        c += (4 * lam * gamma * T * vk / (vk**2 - gamma**2)) * np.exp(
-            -vk * t
-        )
-    return c
-```
-
-```python
+```{code-cell} ipython3
 def plot_result_expectations(plots, axes=None):
     """Plot the expectation values of operators as functions of time.
 
@@ -191,7 +156,7 @@ def plot_result_expectations(plots, axes=None):
     return fig
 ```
 
-```python
+```{code-cell} ipython3
 @contextlib.contextmanager
 def timer(label):
     """Simple utility for timing functions:
@@ -205,7 +170,7 @@ def timer(label):
     print(f"{label}: {end - start}")
 ```
 
-```python
+```{code-cell} ipython3
 # Default solver options:
 
 default_options = {
@@ -222,19 +187,19 @@ default_options = {
 
 And let us set up the system Hamiltonian, bath and system measurement operators:
 
-```python
+```{code-cell} ipython3
 # Defining the system Hamiltonian
 eps = 0.5  # Energy of the 2-level system.
 Del = 1.0  # Tunnelling term
 Hsys = 0.5 * eps * sigmaz() + 0.5 * Del * sigmax()
 ```
 
-```python
+```{code-cell} ipython3
 # Initial state of the system.
 rho0 = basis(2, 0) * basis(2, 0).dag()
 ```
 
-```python
+```{code-cell} ipython3
 # System-bath coupling (Drude-Lorentz spectral density)
 Q = sigmaz()  # coupling operator
 
@@ -252,7 +217,7 @@ Nk = 2  # terms in the Matsubara expansion of the correlation function
 tlist = np.linspace(0, 50, 1000)
 ```
 
-```python
+```{code-cell} ipython3
 # 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()
@@ -265,7 +230,7 @@ P12p = basis(2, 0) * basis(2, 1).dag()
 
 Now we are ready to begin. Let's look at the shape of the spectral density given the bath parameters we defined above:
 
-```python
+```{code-cell} ipython3
 def plot_spectral_density():
     """Plot the Drude-Lorentz spectral density"""
     w = np.linspace(0, 5, 1000)
@@ -287,12 +252,12 @@ 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.
 
-```python
+```{code-cell} ipython3
 ckAR, vkAR, ckAI, vkAI = dl_matsubara_params(nk=Nk, lam=lam, gamma=gamma, T=T)
 ```
 
 Having created the lists which specify the bath correlation functions, we
-create a `BosonicBath` from them and pass the bath to the `HEOMSolver` class.
+create an `ExponentialBosonicEnvironment` from them and pass the environment to the `HEOMSolver` class.
 
 The solver constructs the "right hand side" (RHS) determinining how the
 system and auxiliary density operators evolve in time. This can then be used
@@ -301,18 +266,18 @@ to solve for dynamics or steady-state.
 Below we create the bath and solver and then solve for the dynamics by
 calling `.run(rho0, tlist)`.
 
-```python
+```{code-cell} ipython3
 options = {**default_options}
 
 with timer("RHS construction time"):
-    bath = BosonicBath(Q, ckAR, vkAR, ckAI, vkAI)
-    HEOMMats = HEOMSolver(Hsys, bath, NC, options=options)
+    env = ExponentialBosonicEnvironment(ckAR, vkAR, ckAI, vkAI)
+    HEOMMats = HEOMSolver(Hsys, (env, Q), NC, options=options)
 
 with timer("ODE solver time"):
     resultMats = HEOMMats.run(rho0, tlist)
 ```
 
-```python
+```{code-cell} ipython3
 plot_result_expectations(
     [
         (resultMats, P11p, "b", "P11 Mats"),
@@ -323,37 +288,67 @@ plot_result_expectations(
 
 In practice, one would not perform this laborious expansion for the
 Drude-Lorentz correlation function, because QuTiP already has a class,
-`DrudeLorentzBath`, that can construct this bath for you. Nevertheless,
+`DrudeLorentzEnvironment`, that can construct this bath for you. Nevertheless,
 knowing how to perform this expansion will allow you to construct your own
 baths for other spectral densities.
 
 Below we show how to use this built-in functionality:
 
-```python
+```{code-cell} ipython3
 # Compare to built-in Drude-Lorentz bath:
 
 with timer("RHS construction time"):
-    bath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
-    HEOM_dlbath = HEOMSolver(Hsys, bath, NC, options=options)
+    # Abstract representation of D-L Environment
+    dlenv = DrudeLorentzEnvironment(lam=lam, gamma=gamma, T=T)
+    # Matsubara approximation of D-L Environment
+    dlenv_approx = dlenv.approximate(method="matsubara", Nk=Nk)
+    HEOM_dlbath = HEOMSolver(Hsys, (dlenv_approx, Q), NC, options=options)
 
 with timer("ODE solver time"):
-    result_dlbath = HEOM_dlbath.run(rho0, tlist)  # normal  115
+    result_dlbath = HEOM_dlbath.run(rho0, tlist)
 ```
 
-```python
+```{code-cell} ipython3
 plot_result_expectations(
     [
-        (result_dlbath, P11p, "b", "P11 (DrudeLorentzBath)"),
-        (result_dlbath, P12p, "r", "P12 (DrudeLorentzBath)"),
+        (result_dlbath, P11p, "b", "P11 (DrudeLorentzEnvironment)"),
+        (result_dlbath, P12p, "r", "P12 (DrudeLorentzEnvironment)"),
     ]
 );
 ```
 
+The `DrudeLorentzEnvironment` class also allows us to easily obtain the power spectrum, correlation function, and spectral density. The approximated Environment is a `BosonicEnvironment` where the approximated correlation function, and the corresponding effective power spectrum and spectral density are also accessible. In the following plots, the solid lines are the exact expressions, and the dashed lines are based on our Matsubara approximation of the correlation function with a finite number of exponents.
+
+The `DrudeLorentzEnvironment` computes the exact correlation function using the Pade approximation. The number of terms to use defaults to $10$, but when simulating low temperatures, $10$ Pade exponents may noy be enough. More details about the Pade approximation are provided below. Next we show how to use this built-in functionality:
+
+```{code-cell} ipython3
+w = np.linspace(-10, 20, 1000)
+w2 = np.linspace(0, 20, 1000)
+
+fig, axs = plt.subplots(2, 2)
+
+axs[0, 0].plot(w, dlenv.power_spectrum(w))
+axs[0, 0].plot(w, dlenv_approx.power_spectrum(w), "--")
+axs[0, 0].set(xlabel=r"$\omega$", ylabel=r"$S(\omega)$")
+axs[0, 1].plot(w2, dlenv.spectral_density(w2))
+axs[0, 1].plot(w2, dlenv_approx.spectral_density(w2), "--")
+axs[0, 1].set(xlabel=r"$\omega$", ylabel=r"$J(\omega)$")
+axs[1, 0].plot(w2, np.real(dlenv.correlation_function(w2, Nk=100)))  # 100 Pade
+axs[1, 0].plot(w2, np.real(dlenv_approx.correlation_function(w2)), "--")
+axs[1, 0].set(xlabel=r"$t$", ylabel=r"$C_{R}(t)$")
+axs[1, 1].plot(w2, np.imag(dlenv.correlation_function(w2, Nk=100)))
+axs[1, 1].plot(w2, np.imag(dlenv_approx.correlation_function(w2)), "--")
+axs[1, 1].set(xlabel=r"$t$", ylabel=r"$C_{I}(t)$")
+
+fig.tight_layout()
+plt.show()
+```
+
 We also provide a legacy class, `HSolverDL`, which calculates the
 Drude-Lorentz correlation functions automatically, to be backwards
 compatible with the previous HEOM solver in QuTiP:
 
-```python
+```{code-cell} ipython3
 # Compare to legacy class:
 
 # The legacy class performs the above collation of coefficients automatically,
@@ -366,7 +361,7 @@ with timer("ODE solver time"):
     resultLegacy = HEOMlegacy.run(rho0, tlist)  # normal  115
 ```
 
-```python
+```{code-cell} ipython3
 plot_result_expectations(
     [
         (resultLegacy, P11p, "b", "P11 Legacy"),
@@ -375,6 +370,20 @@ plot_result_expectations(
 );
 ```
 
+Another legacy class kept for convenience is the `DrudeLorentzBath`. The code
+```python
+dlenv = DrudeLorentzEnvironment(lam=lam, gamma=gamma, T=T)
+dlenv_approx = dlenv.approximate(method="matsubara", Nk=Nk)  # Computes Matsubara exponents
+HEOM_dlbath = HEOMSolver(Hsys, (dlenv_approx, Q), NC, options=options)
+```
+that we used above is equivalent to the following code:
+```python
+dlbath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
+HEOM_dlbath = HEOMSolver(Hsys, dlbath, NC, options=options)
+```
+
++++
+
 ## Ishizaki-Tanimura Terminator
 
 To speed up convergence (in terms of the number of exponents kept in the
@@ -394,38 +403,39 @@ $ e^{-\nu_k t} \approx \delta(t)/\nu_k$, and $C(t)=\sum_{k=N_k}^{\infty}
 \frac{c_k}{\nu_k} \delta(t)$
 
 It is convenient to calculate the whole sum
-$ C(t)=\sum_{k=0}^{\infty} \frac{c_k}{\nu_k} =  2 \lambda / (\beta \gamma)
-- i\lambda $
+$ C(t)=\sum_{k=0}^{\infty} \frac{c_k}{\nu_k} =  2 \lambda / (\beta \gamma)- i\lambda $
 , and subtract off the contribution from the finite number of Matsubara terms
-that are kept in the hierarchy, and treat the residual as a Lindblad.
+that are kept in the hierarchy, and treat the residual as a contribution in 
+Lindblad form.
 
 This is clearer if we plot the correlation function with a large number of
-Matsubara terms:
+Matsubara terms. To create the plot, we use the utility function of the
+`DrudeLorentzEnvironment` class mentioned above.
 
-```python
+```{code-cell} ipython3
 def plot_correlation_expansion_divergence():
     """We plot the correlation function with a large number of Matsubara terms
     to show that the real part is slowly diverging at t = 0.
     """
     t = np.linspace(0, 2, 100)
 
-    # correlation coefficients with 15k and 2 terms
-    corr_15k = dl_corr_approx(t, 15_000)
-    corr_2 = dl_corr_approx(t, 2)
+    # correlation coefficients with 100 pade and with 2 matsubara terms
+    corr_100 = dlenv.correlation_function(t, Nk=100)
+    corr_2 = dlenv_approx.correlation_function(t)
 
     fig, ax1 = plt.subplots(figsize=(12, 7))
 
     ax1.plot(
-        t, np.real(corr_2), color="b", linewidth=3, label=r"Mats = 2 real"
+        t, np.real(corr_2), color="b", linewidth=3, label=rf"Mats = {Nk} real"
     )
     ax1.plot(
-        t, np.imag(corr_2), color="r", linewidth=3, label=r"Mats = 2 imag"
+        t, np.imag(corr_2), color="r", linewidth=3, label=rf"Mats = {Nk} imag"
     )
     ax1.plot(
-        t, np.real(corr_15k), "b--", linewidth=3, label=r"Mats = 15000 real"
+        t, np.real(corr_100), "b--", linewidth=3, label=r"Pade = 100 real"
     )
     ax1.plot(
-        t, np.imag(corr_15k), "r--", linewidth=3, label=r"Mats = 15000 imag"
+        t, np.imag(corr_100), "r--", linewidth=3, label=r"Pade = 100 imag"
     )
 
     ax1.set_xlabel("t")
@@ -433,18 +443,21 @@ def plot_correlation_expansion_divergence():
     ax1.legend()
 
 
-plot_correlation_expansion_divergence();
+plot_correlation_expansion_divergence()
 ```
 
 Let us evaluate the result including this Ishizaki-Tanimura terminator:
 
-```python
-# Run HEOM solver and include the Ishizaki-Tanimura terminator
+```{code-cell} ipython3
+# Run HEOM solver including the Ishizaki-Tanimura terminator
 
 # Notes:
 #
-# * when using the built-in DrudeLorentzBath, the terminator (L_bnd) is
-#   available from bath.terminator().
+# * here, we will first show how to compute the terminator manually
+#
+# * when using the built-in DrudeLorentzEnvironment the terminator (L_bnd) is
+#   available from by setting the parameter compute_delta to True in the
+#   approximate method
 #
 # * in the legacy HSolverDL function the terminator is included automatically
 #   if the parameter bnd_cut_approx=True is used.
@@ -467,14 +480,14 @@ Ltot = liouvillian(Hsys) + L_bnd
 options = {**default_options, "rtol": 1e-14, "atol": 1e-14}
 
 with timer("RHS construction time"):
-    bath = BosonicBath(Q, ckAR, vkAR, ckAI, vkAI)
-    HEOMMatsT = HEOMSolver(Ltot, bath, NC, options=options)
+    env = ExponentialBosonicEnvironment(ckAR, vkAR, ckAI, vkAI)
+    HEOMMatsT = HEOMSolver(Ltot, (env, Q), NC, options=options)
 
 with timer("ODE solver time"):
     resultMatsT = HEOMMatsT.run(rho0, tlist)
 ```
 
-```python
+```{code-cell} ipython3
 plot_result_expectations(
     [
         (resultMatsT, P11p, "b", "P11 Mats + Term"),
@@ -483,47 +496,41 @@ plot_result_expectations(
 );
 ```
 
-Or using the built-in Drude-Lorentz bath we can write simply:
+Or using the built-in Drude-Lorentz environment we can write simply:
 
-```python
+```{code-cell} ipython3
 options = {**default_options, "rtol": 1e-14, "atol": 1e-14}
 
 with timer("RHS construction time"):
-    bath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
-    _, terminator = bath.terminator()
-    Ltot = liouvillian(Hsys) + terminator
-    HEOM_dlbath_T = HEOMSolver(Ltot, bath, NC, options=options)
+    dlenv_approx, delta = dlenv.approximate(
+        "matsubara", Nk=Nk, compute_delta=True
+    )
+    Ltot = liouvillian(Hsys) + system_terminator(Q, delta)
+    HEOM_dlbath_T = HEOMSolver(Ltot, (dlenv_approx, Q), NC, options=options)
 
 with timer("ODE solver time"):
     result_dlbath_T = HEOM_dlbath_T.run(rho0, tlist)
 ```
 
-```python
-plot_result_expectations(
-    [
-        (result_dlbath_T, P11p, "b", "P11 Mats (DrudeLorentzBath + Term)"),
-        (result_dlbath_T, P12p, "r", "P12 Mats (DrudeLorentzBath + Term)"),
-    ]
-);
+```{code-cell} ipython3
+plot_result_expectations([
+    (result_dlbath_T, P11p, "b", "P11 Mats (DrudeLorentzEnvironment + Term)"),
+    (result_dlbath_T, P12p, "r", "P12 Mats (DrudeLorentzEnvironment + Term)"),
+]);
 ```
 
 We can compare the solution obtained from the QuTiP Bloch-Redfield solver:
 
-```python
-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))) "
-    f"* ((1/(exp((w) * {beta})-1))+1)"
-)
+```{code-cell} ipython3
 options = {**default_options}
 
 with timer("ODE solver time"):
     resultBR = brmesolve(
-        Hsys, rho0, tlist, a_ops=[[sigmaz(), DL]], options=options
+        Hsys, rho0, tlist, a_ops=[[sigmaz(), dlenv]], options=options
     )
 ```
 
-```python
+```{code-cell} ipython3
 plot_result_expectations(
     [
         (resultMats, P11p, "b", "P11 Mats"),
@@ -538,11 +545,12 @@ plot_result_expectations(
 
 ## Padé decomposition
 
++++
 
 The Matsubara decomposition is not the only option.  We can also use the
 faster-converging Pade decomposition.
 
-```python
+```{code-cell} ipython3
 def deltafun(j, k):
     if j == k:
         return 1.0
@@ -580,7 +588,7 @@ def pade_chi(lmax):
             )
 
     eigvalsAP = np.linalg.eigvalsh(AlphaP)
-    chi = [-2 / val for val in eigvalsAP[0:lmax - 1]]
+    chi = [-2 / val for val in eigvalsAP[0:(lmax - 1)]]
     return chi
 
 
@@ -641,8 +649,8 @@ def pade_corr(tlist, lmax):
 
 tlist_corr = np.linspace(0, 2, 100)
 cppLP, etapLP, gampLP = pade_corr(tlist_corr, 2)
-corr_15k = dl_corr_approx(tlist_corr, 15_000)
-corr_2k = dl_corr_approx(tlist_corr, 2)
+corr_100 = dlenv.correlation_function(tlist_corr, Nk=100)
+corr_2 = dlenv.approximate("matsubara", Nk=2).correlation_function(tlist_corr)
 
 fig, ax1 = plt.subplots(figsize=(12, 7))
 ax1.plot(
@@ -654,35 +662,35 @@ ax1.plot(
 )
 ax1.plot(
     tlist_corr,
-    np.real(corr_15k),
+    np.real(corr_100),
     "r--",
     linewidth=3,
-    label=r"real mats 15000 terms",
+    label=r"real pade 100 terms",
 )
 ax1.plot(
     tlist_corr,
-    np.real(corr_2k),
+    np.real(corr_2),
     "g--",
     linewidth=3,
     label=r"real mats 2 terms",
 )
 
 ax1.set_xlabel("t")
-ax1.set_ylabel(r"$C$")
+ax1.set_ylabel(r"$C_{R}(t)$")
 ax1.legend()
 
 fig, ax1 = plt.subplots(figsize=(12, 7))
 
 ax1.plot(
     tlist_corr,
-    np.real(cppLP) - np.real(corr_15k),
+    np.real(cppLP) - np.real(corr_100),
     color="b",
     linewidth=3,
     label=r"pade error",
 )
 ax1.plot(
     tlist_corr,
-    np.real(corr_2k) - np.real(corr_15k),
+    np.real(corr_2) - np.real(corr_100),
     "r--",
     linewidth=3,
     label=r"mats error",
@@ -693,7 +701,7 @@ ax1.set_ylabel(r"Error")
 ax1.legend();
 ```
 
-```python
+```{code-cell} ipython3
 # put pade parameters in lists for heom solver
 ckAR = [np.real(eta) + 0j for eta in etapLP]
 ckAI = [np.imag(etapLP[0]) + 0j]
@@ -703,14 +711,14 @@ vkAI = [gampLP[0] + 0j]
 options = {**default_options, "rtol": 1e-14, "atol": 1e-14}
 
 with timer("RHS construction time"):
-    bath = BosonicBath(Q, ckAR, vkAR, ckAI, vkAI)
-    HEOMPade = HEOMSolver(Hsys, bath, NC, options=options)
+    bath = ExponentialBosonicEnvironment(ckAR, vkAR, ckAI, vkAI)
+    HEOMPade = HEOMSolver(Hsys, (bath, Q), NC, options=options)
 
 with timer("ODE solver time"):
     resultPade = HEOMPade.run(rho0, tlist)
 ```
 
-```python
+```{code-cell} ipython3
 plot_result_expectations(
     [
         (resultMats, P11p, "b", "P11 Mats"),
@@ -723,61 +731,67 @@ plot_result_expectations(
 );
 ```
 
-The Padé decomposition of the Drude-Lorentz bath is also available via a
-built-in class, `DrudeLorentzPadeBath` bath. Like `DrudeLorentzBath`, one
-can obtain the terminator by calling `bath.terminator()`.
+As mentioned previously, the Padé decomposition of the Drude-Lorentz bath is also available via the
+built-in `DrudeLorentzEnvironment`. Similarly to the terminator
+section when approximating by Padé one can calculate the terminator easily by
+requesting the approximation function to compute delta
 
-Below we show how to use the built-in Padé Drude-Lorentz bath and its
-terminator (although the terminator does not provide much improvement here,
-because the Padé expansion already fits the correlation function well):
+Below we show how to use the built-in Drude-Lorentz Environment to obtain a
+Padé decomposition approximation and its terminator (although the terminator 
+does not provide much improvement here,because the Padé expansion already fits 
+the correlation function well):
 
-```python
+```{code-cell} ipython3
 options = {**default_options, "rtol": 1e-14, "atol": 1e-14}
 
 with timer("RHS construction time"):
-    bath = DrudeLorentzPadeBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
-    _, terminator = bath.terminator()
-    Ltot = liouvillian(Hsys) + terminator
-    HEOM_dlpbath_T = HEOMSolver(Ltot, bath, NC, options=options)
+    env_approx, delta = dlenv.approximate("pade", Nk=2, compute_delta=True)
+    Ltot = liouvillian(Hsys) + system_terminator(Q, delta)
+    HEOM_dlpbath_T = HEOMSolver(Ltot, (env_approx, Q), NC, options=options)
 
 with timer("ODE solver time"):
     result_dlpbath_T = HEOM_dlpbath_T.run(rho0, tlist)
 ```
 
-```python
+```{code-cell} ipython3
 plot_result_expectations(
     [
-        (result_dlpbath_T, P11p, "b", "P11 Padé (DrudeLorentzBath + Term)"),
-        (result_dlpbath_T, P12p, "r", "P12 Padé (DrudeLorentzBath + Term)"),
+        (result_dlpbath_T, P11p, "b", "P11 Padé + Term"),
+        (result_dlpbath_T, P12p, "r", "P12 Padé + Term"),
     ]
 );
 ```
 
 ### Next we compare the Matsubara and Pade correlation function fits
 
-This is not efficient for this example, but can be extremely useful in
-situations where large number of exponents are needed (e.g., near zero
-temperature).
+Fitting the correlation function is not efficient for this example, but
+can be extremely useful in situations where large number of exponents
+are needed (e.g., near zero temperature). We will perform the fitting
+manually below, and then show how to do it with the built-in tools.
 
-First we collect a large sum of Matsubara terms for many time steps:
+For the manual fit we first we collect a large sum of Pade terms for 
+many time steps:
 
-```python
+```{code-cell} ipython3
 tlist2 = np.linspace(0, 2, 10000)
 
-corr_15k_t10k = dl_corr_approx(tlist2, 15_000)
+corr_100_t10k = dlenv.correlation_function(tlist2, Nk=100)
+# Nk specifies the number of pade terms to be used for the correlation function
+
+corrRana = np.real(corr_100_t10k)
+corrIana = np.imag(corr_100_t10k)
 
-corrRana = np.real(corr_15k_t10k)
-corrIana = np.imag(corr_15k_t10k)
+corrRMats = np.real(dlenv_approx.correlation_function(tlist2))
 ```
 
 We then fit this sum with standard least-squares approach:
 
-```python
+```{code-cell} ipython3
 def wrapper_fit_func(x, N, args):
     """ Fit function wrapper that unpacks its arguments. """
     x = np.array(x)
     a = np.array(args[:N])
-    b = np.array(args[N:2 * N])
+    b = np.array(args[N:(2 * N)])
     return fit_func(x, a, b)
 
 
@@ -824,15 +838,13 @@ def fitter(ans, tlist, k):
     return (a, b)
 ```
 
-```python
+```{code-cell} ipython3
 kR = 4  # number of exponents to use for real part
 poptR = []
 with timer("Correlation (real) fitting time"):
     for i in range(kR):
         poptR.append(fitter(corrRana, tlist2, i + 1))
 
-corrRMats = np.real(dl_corr_approx(tlist2, Nk))
-
 kI = 1  # number of exponents for imaginary part
 poptI = []
 with timer("Correlation (imaginary) fitting time"):
@@ -842,32 +854,55 @@ with timer("Correlation (imaginary) fitting time"):
 
 And plot the results of the fits:
 
-```python
-plt.plot(tlist2, corrRana, label="Analytic")
-plt.plot(tlist2, corrRMats, label="Matsubara")
+```{code-cell} ipython3
+# Define line styles and colors
+linestyles = ["-", "--", "-.", ":", (0, (3, 1, 1, 1)), (0, (5, 1))]
+colors = ["blue", "green", "purple", "orange", "red", "brown"]
+
+# Define a larger linewidth
+linewidth = 2.5
+
+# Create a single figure with two subplots
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
+
+# Plot the real part on the first subplot (ax1)
+ax1.plot(tlist2, corrRana, label="Analytic", color=colors[0],
+         linestyle=linestyles[0], linewidth=linewidth)
+ax1.plot(tlist2, corrRMats, label="Matsubara", color=colors[1],
+         linestyle=linestyles[1], linewidth=linewidth)
 
 for i in range(kR):
     y = fit_func(tlist2, *poptR[i])
-    plt.plot(tlist2, y, label=f"Fit with {i} terms")
+    ax1.plot(tlist2, y, label=f"Fit with {i+1} terms", color=colors[i + 2],
+             linestyle=linestyles[i + 2], linewidth=linewidth)
 
-plt.title("Fit to correlations (real part)")
-plt.legend()
-plt.show()
-```
+ax1.set_ylabel(r"$C_{R}(t)$")
+ax1.set_xlabel(r"$t$")
+ax1.legend()
 
-```python
-plt.plot(tlist2, corrIana, label="Analytic")
+# Plot the imaginary part on the second subplot (ax2)
+ax2.plot(tlist2, corrIana, label="Analytic", color=colors[0],
+         linestyle=linestyles[0], linewidth=linewidth)
 
 for i in range(kI):
     y = fit_func(tlist2, *poptI[i])
-    plt.plot(tlist2, y, label=f"Fit with {i} terms")
-
-plt.title("Fit to correlations (imaginary part)")
-plt.legend()
+    ax2.plot(tlist2, y, label=f"Fit with {i+1} terms", color=colors[i + 3],
+             linestyle=linestyles[i + 3], linewidth=linewidth)
+
+ax2.set_ylabel(r"$C_{I}(t)$")
+ax2.set_xlabel(r"$t$")
+ax2.legend()
+
+# Add overall plot title and show the figure
+fig.suptitle(
+    "Comparison of Analytic and Fit to Correlations"
+    " (Real and Imaginary Parts)",
+    fontsize=16,
+)
 plt.show()
 ```
 
-```python
+```{code-cell} ipython3
 # Set the exponential coefficients from the fit parameters
 
 ckAR1 = poptR[-1][0]
@@ -883,7 +918,7 @@ vkAI1 = poptI[-1][1]
 vkAI = [-x + 0j for x in vkAI1]
 ```
 
-```python
+```{code-cell} ipython3
 # overwrite imaginary fit with analytical value (not much reason to use the
 # fit for this)
 
@@ -891,20 +926,20 @@ ckAI = [lam * gamma * (-1.0) + 0.0j]
 vkAI = [gamma + 0.0j]
 ```
 
-```python
+```{code-cell} ipython3
 options = {**default_options}
 
 NC = 4
 
 with timer("RHS construction time"):
-    bath = BosonicBath(Q, ckAR, vkAR, ckAI, vkAI)
-    HEOMFit = HEOMSolver(Hsys, bath, NC, options=options)
+    bath = ExponentialBosonicEnvironment(ckAR, vkAR, ckAI, vkAI)
+    HEOMFit = HEOMSolver(Hsys, (bath, Q), NC, options=options)
 
 with timer("ODE solver time"):
     resultFit = HEOMFit.run(rho0, tlist)
 ```
 
-```python
+```{code-cell} ipython3
 plot_result_expectations(
     [
         (resultFit, P11p, "b", "P11 Fit"),
@@ -913,15 +948,90 @@ plot_result_expectations(
 );
 ```
 
+Now we use the built-in fitting functions. The `BosonicEnvironment` class includes a 
+method that performs this fit automatically. More information on how the
+built-in functios work can be found in `HEOM 1d: Spin-Bath model, fitting of spectrum and correlation functions`
+
+```{code-cell} ipython3
+max_val = dlenv.correlation_function(0).real
+guess = [max_val / 3, 0, 0, 0]
+lower = [-max_val, -np.inf, -np.inf, -np.inf]
+upper = [max_val, 0, 0, 0]
+envfit, fitinfo = dlenv.approximate(
+    "cf", tlist=tlist2, full_ansatz=True, Ni_max=1, Nr_max=3,
+    upper=upper, lower=lower, guess=guess)
+```
+
+The `approximate("cf", ...)` method outputs an `ExponentialBosonicEnvironment` object,
+which contains a decaying exponential representation of the original 
+environment , and a dictionary containing all information related to the fit.
+The dictionary contains a summary of the fir information and the normalized 
+root mean squared error, which assesses how good the fit is.
+
+```{code-cell} ipython3
+print(fitinfo["summary"])
+```
+
+We can then compare the result of the built-in fit with the manual fit
+
+```{code-cell} ipython3
+# Create a single figure with two subplots
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
+
+# Plot the real part on the first subplot (ax1)
+ax1.plot(tlist2, corrRana, label="Original", marker="o", markevery=500)
+ax1.plot(tlist2, fit_func(tlist2, *poptR[-1]), color="r", label="Manual Fit")
+ax1.plot(tlist2, np.real(envfit.correlation_function(tlist2)), "k--",
+         label="Built-in fit")
+ax1.set_ylabel(r"$C_{R}(t)$")
+ax1.set_xlabel(r"$t$")
+ax1.legend()
+
+# Plot the imaginary part on the second subplot (ax2)
+ax2.plot(tlist2, corrIana, label="Original", marker="o", markevery=500)
+ax2.plot(tlist2, fit_func(tlist2, *poptI[-1]), color="r", label="Manual Fit")
+ax2.plot(tlist2, np.imag(envfit.correlation_function(tlist2)), "k--",
+         label="Built-in fit")
+ax2.set_ylabel(r"$C_{I}(t)$")
+ax2.set_xlabel(r"$t$")
+ax2.legend()
+
+# Adjust layout
+plt.tight_layout(rect=[0, 0.03, 1, 0.95])
+plt.show()
+```
+
+```{code-cell} ipython3
+options = {**default_options}
+
+with timer("RHS construction time"):
+    HEOMFit_2 = HEOMSolver(Hsys, (envfit, Q), NC, options=options)
+
+with timer("ODE solver time"):
+    resultFit_2 = HEOMFit_2.run(rho0, tlist)
+```
+
+```{code-cell} ipython3
+plot_result_expectations(
+    [
+        (resultFit, P11p, "b", "P11 Fit"),
+        (resultFit, P12p, "r", "P12 Fit"),
+        (resultFit_2, P11p, "r--", "P11 Built-in-Fit"),
+        (resultFit_2, P12p, "b--", "P12 Built-in-Fit"),
+    ]
+);
+```
+
 ## A reaction coordinate approach
 
++++
 
 Here we construct a reaction coordinate inspired model to capture the
 steady-state behavior, and compare to the HEOM prediction. This result is
 more accurate for narrow spectral densities.  We will use the population and
 coherence from this cell in our final plot below.
 
-```python
+```{code-cell} ipython3
 dot_energy, dot_state = Hsys.eigenstates()
 deltaE = dot_energy[1] - dot_energy[0]
 
@@ -930,9 +1040,7 @@ wa = 2 * np.pi * gamma2 * gamma  # reaction coordinate frequency
 g = np.sqrt(np.pi * wa * lam / 2.0)  # reaction coordinate coupling
 # reaction coordinate coupling factor over 2 because of diff in J(w)
 # (it is 2 lam now):
-g = np.sqrt(
-    np.pi * wa * lam / 4.0
-)  #
+g = np.sqrt(np.pi * wa * lam / 4.0)  #
 
 NRC = 10
 
@@ -970,7 +1078,7 @@ P11RC = tensor(qeye(NRC), basis(2, 0) * basis(2, 0).dag())
 
 Finally, let's plot all of our different results to see how they shape up against each other.
 
-```python
+```{code-cell} ipython3
 rcParams = {
     "axes.titlesize": 25,
     "axes.labelsize": 30,
@@ -987,7 +1095,7 @@ rcParams = {
 }
 ```
 
-```python
+```{code-cell} ipython3
 fig, axes = plt.subplots(2, 1, sharex=False, figsize=(12, 15))
 
 with plt.rc_context(rcParams):
@@ -1009,13 +1117,14 @@ with plt.rc_context(rcParams):
                 resultFit,
                 P11p,
                 "r",
-                r"Fit $N_f = 4$, $N_k=15\times 10^3$",
+                r"Fit $N_f = 4$, Pade $N_k=100$",
                 {"dashes": [3, 2]},
             ),
             (
                 resultRC,
                 P11RC,
-                "--", "Thermal",
+                "--",
+                "Thermal",
                 {"linewidth": 2, "color": "black"},
             ),
         ],
@@ -1046,7 +1155,7 @@ with plt.rc_context(rcParams):
                 resultFit,
                 P12p,
                 "r",
-                r"Fit $N_f = 4$, $N_k=15\times 10^3$",
+                r"Fit $N_f = 4$, Pade $N_k=100$",
                 {"dashes": [3, 2]},
             ),
             (
@@ -1067,18 +1176,19 @@ with plt.rc_context(rcParams):
 
 And that's the end of a detailed first dive into modeling bosonic environments with the HEOM.
 
++++
 
 ## About
 
-```python
-qutip.about()
+```{code-cell} ipython3
+about()
 ```
 
 ## Testing
 
 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.
 
-```python
+```{code-cell} ipython3
 # Check P11p
 assert np.allclose(
     expect(P11p, resultMatsT.states),
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..57a29e32 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
@@ -1,15 +1,14 @@
 ---
 jupytext:
-  formats: ipynb,md:myst
   text_representation:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.17.0
 kernelspec:
+  name: python3
   display_name: Python 3 (ipykernel)
   language: python
-  name: python3
 ---
 
 # HEOM 1b: Spin-Bath model (very strong coupling)
@@ -86,25 +85,11 @@ import contextlib
 import time
 
 import numpy as np
-from scipy.optimize import curve_fit
 import matplotlib.pyplot as plt
 
-import qutip
-from qutip import (
-    basis,
-    brmesolve,
-    expect,
-    liouvillian,
-    sigmax,
-    sigmaz,
-)
-from qutip.solver.heom import (
-    HEOMSolver,
-    BosonicBath,
-    DrudeLorentzBath,
-    DrudeLorentzPadeBath,
-    BathExponent,
-)
+from qutip import about, basis, brmesolve, expect, liouvillian, sigmax, sigmaz
+from qutip.core.environment import DrudeLorentzEnvironment, system_terminator
+from qutip.solver.heom import HEOMSolver
 
 %matplotlib inline
 ```
@@ -116,7 +101,7 @@ Let's define some helper functions for calculating correlation function expansio
 ```{code-cell} ipython3
 def cot(x):
     """ Vectorized cotangent of x. """
-    return 1. / np.tan(x)
+    return 1.0 / np.tan(x)
 ```
 
 ```{code-cell} ipython3
@@ -152,8 +137,8 @@ And let us set up the system Hamiltonian, bath and system measurement operators:
 
 ```{code-cell} ipython3
 # Defining the system Hamiltonian
-eps = .0     # Energy of the 2-level system.
-Del = .2     # Tunnelling term
+eps = 0.0  # Energy of the 2-level system.
+Del = 0.2  # Tunnelling term
 Hsys = 0.5 * eps * sigmaz() + 0.5 * Del * sigmax()
 ```
 
@@ -167,10 +152,10 @@ rho0 = basis(2, 0) * basis(2, 0).dag()
 Q = sigmaz()  # coupling operator
 
 # Bath properties (see Shi et al., J. Chem. Phys. 130, 084105 (2009)):
-gamma = 1.  # cut off frequency
-lam = 2.5   # coupling strength
-T = 1.      # in units where Boltzmann factor is 1
-beta = 1. / T
+gamma = 1.0  # cut off frequency
+lam = 2.5    # coupling strength
+T = 1.0      # in units where Boltzmann factor is 1
+beta = 1.0 / T
 
 # HEOM parameters:
 
@@ -199,8 +184,9 @@ P12p = basis(2, 0) * basis(2, 1).dag()
 Let us briefly inspect the spectral density.
 
 ```{code-cell} ipython3
+env = DrudeLorentzEnvironment(lam=lam, gamma=gamma, T=T, Nk=500)
 w = np.linspace(0, 5, 1000)
-J = w * 2 * lam * gamma / ((gamma**2 + w**2))
+J = env.spectral_density(w)
 
 # Plot the results
 fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 8))
@@ -213,8 +199,8 @@ axes.set_ylabel(r'J', fontsize=28);
 
 ```{code-cell} ipython3
 with timer("RHS construction time"):
-    bath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
-    HEOMMats = HEOMSolver(Hsys, bath, NC, options=options)
+    matsEnv = env.approximate(method="matsubara", Nk=Nk)
+    HEOMMats = HEOMSolver(Hsys, (matsEnv, Q), NC, options=options)
 
 with timer("ODE solver time"):
     resultMats = HEOMMats.run(rho0, tlist)
@@ -224,10 +210,12 @@ with timer("ODE solver time"):
 
 ```{code-cell} ipython3
 with timer("RHS construction time"):
-    bath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
-    _, terminator = bath.terminator()
+    matsEnv, delta = env.approximate(
+        method="matsubara", Nk=Nk, compute_delta=True
+    )
+    terminator = system_terminator(Q, delta)
     Ltot = liouvillian(Hsys) + terminator
-    HEOMMatsT = HEOMSolver(Ltot, bath, NC, options=options)
+    HEOMMatsT = HEOMSolver(Ltot, (matsEnv, Q), NC, options=options)
 
 with timer("ODE solver time"):
     resultMatsT = HEOMMatsT.run(rho0, tlist)
@@ -258,51 +246,21 @@ axes.legend(loc=0, fontsize=12);
 
 ```{code-cell} ipython3
 # 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)
-
-# We will compare against a summation of {lmaxmats} Matsubara terms
-lmaxmats = 15000
-exactBath = DrudeLorentzBath(
-    Q, lam=lam, gamma=gamma, T=T, Nk=lmaxmats, combine=False,
-)
-
-
-def CR(bath, t):
-    """ C_R, the real part of the correlation function. """
-    result = 0
-    for exp in bath.exponents:
-        if (
-            exp.type == BathExponent.types['R'] or
-            exp.type == BathExponent.types['RI']
-        ):
-            result += exp.ck * np.exp(-exp.vk * t)
-    return result
-
-
-def CI(bath, t):
-    """ C_I, the imaginary part of the correlation function. """
-    result = 0
-    for exp in bath.exponents:
-        if exp.type == BathExponent.types['I']:
-            result += exp.ck * np.exp(exp.vk * t)
-        if exp.type == BathExponent.types['RI']:
-            result += exp.ck2 * np.exp(exp.vk * t)
-    return result
+padeEnv = env.approximate("pade", Nk=Nk)
 
 
 fig, (ax1, ax2) = plt.subplots(ncols=2, sharey=True, figsize=(16, 8))
 
 ax1.plot(
-    tlist, CR(exactBath, tlist),
-    "r", linewidth=2, label=f"Mats (Nk={lmaxmats})",
+    tlist, np.real(env.correlation_function(tlist)),
+    "r", linewidth=2, label="Exact",
 )
 ax1.plot(
-    tlist, CR(matsBath, tlist),
+    tlist, np.real(matsEnv.correlation_function(tlist)),
     "g--", linewidth=2, label=f"Mats (Nk={Nk})",
 )
 ax1.plot(
-    tlist, CR(padeBath, tlist),
+    tlist, np.real(padeEnv.correlation_function(tlist)),
     "b--", linewidth=2, label=f"Pade (Nk={Nk})",
 )
 
@@ -312,11 +270,13 @@ ax1.legend(loc=0, fontsize=12)
 
 tlist2 = tlist[0:50]
 ax2.plot(
-    tlist2, np.abs(CR(matsBath, tlist2) - CR(exactBath, tlist2)),
+    tlist2, np.abs(matsEnv.correlation_function(tlist2)
+                   - env.correlation_function(tlist2)),
     "g", linewidth=2, label="Mats Error",
 )
 ax2.plot(
-    tlist2, np.abs(CR(padeBath, tlist2) - CR(exactBath, tlist2)),
+    tlist2, np.abs(padeEnv.correlation_function(tlist2)
+                   - env.correlation_function(tlist2)),
     "b--", linewidth=2, label="Pade Error",
 )
 
@@ -326,8 +286,7 @@ ax2.legend(loc=0, fontsize=12);
 
 ```{code-cell} ipython3
 with timer("RHS construction time"):
-    bath = DrudeLorentzPadeBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
-    HEOMPade = HEOMSolver(Hsys, bath, NC, options=options)
+    HEOMPade = HEOMSolver(Hsys, (padeEnv, Q), NC, options=options)
 
 with timer("ODE solver time"):
     resultPade = HEOMPade.run(rho0, tlist)
@@ -358,107 +317,77 @@ axes.legend(loc=0, fontsize=12);
 
 ## Simulation 4: Fitting approach
 
-```{code-cell} ipython3
-def wrapper_fit_func(x, N, args):
-    """ Fit function wrapper that unpacks its arguments. """
-    x = np.array(x)
-    a = np.array(args[:N])
-    b = np.array(args[N:2 * N])
-    return fit_func(x, a, b)
-
-
-def fit_func(x, a, b):
-    """ Fit function. Calculates the value of the
-        correlation function at each x, given the
-        fit parameters in a and b.
-    """
-    return np.sum(
-        a[:, None] * np.exp(np.multiply.outer(b, x)),
-        axis=0,
-    )
-
-
-def fitter(ans, tlist, k):
-    """ Compute fit with k exponents. """
-    upper_a = abs(max(ans, key=abs)) * 10
-    # sets initial guesses:
-    guess = (
-        [upper_a / k] * k +  # guesses for a
-        [0] * k  # guesses for b
-    )
-    # sets lower bounds:
-    b_lower = (
-        [-upper_a] * k +  # lower bounds for a
-        [-np.inf] * k  # lower bounds for b
-    )
-    # sets higher bounds:
-    b_higher = (
-        [upper_a] * k +  # upper bounds for a
-        [0] * k  # upper bounds for b
-    )
-    param_bounds = (b_lower, b_higher)
-    p1, p2 = curve_fit(
-        lambda x, *params_0: wrapper_fit_func(x, k, params_0),
-        tlist,
-        ans,
-        p0=guess,
-        sigma=[0.01 for t in tlist],
-        bounds=param_bounds,
-        maxfev=1e8,
-    )
-    a, b = p1[:k], p1[k:]
-    return (a, b)
-```
+In `HEOM 1a: Spin-Bath model (introduction)` there is an example of a manually performed fit, here
+we will only use the built-in tools. More details about them can be seen in 
+`HEOM 1d: Spin-Bath model, fitting of spectrum and correlation functions`
 
 ```{code-cell} ipython3
-# Fitting the real part of the correlation function:
-
-# Correlation function values to fit:
-tlist_fit = np.linspace(0, 6, 10000)
-corrRana = CR(exactBath, tlist_fit)
-
-# Perform the fit:
-kR = 3  # number of exponents to use for real part
-poptR = []
-with timer("Correlation (real) fitting time"):
-    for i in range(kR):
-        poptR.append(fitter(corrRana, tlist_fit, i + 1))
+tfit = np.linspace(0, 10, 10000)
+lower = [0, -np.inf, -1e-6, -3]
+guess = [np.real(env.correlation_function(0)) / 10, -10, 0, 0]
+upper = [5, 0, 1e-6, 0]
+# for better fits increase the first element in upper, or change approximate
+# method that makes the simulation much slower (Larger C(t) as C(0) is fit
+# better)
+
+envfit, fitinfo = env.approximate(
+    "cf", tlist=tfit, Nr_max=2, Ni_max=1, full_ansatz=True,
+    sigma=0.1, maxfev=1e6, target_rmse=None,
+    lower=lower, upper=upper, guess=guess,
+)
 ```
 
 ```{code-cell} ipython3
-plt.plot(tlist_fit, corrRana, label="Analytic")
-
-for i in range(kR):
-    y = fit_func(tlist_fit, *poptR[i])
-    plt.plot(tlist_fit, y, label=f"Fit with {i} terms")
-
-plt.title("Fit to correlations (real part)")
-plt.legend()
-plt.show()
+print(fitinfo["summary"])
 ```
 
+We can quickly compare the result of the fit with the Pade expansion
+
 ```{code-cell} ipython3
-# Set the exponential coefficients from the fit parameters
+fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(16, 8))
 
-ckAR1 = poptR[-1][0]
-ckAR = [x + 0j for x in ckAR1]
+ax1.plot(
+    tlist, np.real(env.correlation_function(tlist)),
+    "r", linewidth=2, label="Exact",
+)
+ax1.plot(
+    tlist, np.real(envfit.correlation_function(tlist)),
+    "g--", linewidth=2, label="Fit", marker="o", markevery=50,
+)
+ax1.plot(
+    tlist, np.real(padeEnv.correlation_function(tlist)),
+    "b--", linewidth=2, label=f"Pade (Nk={Nk})",
+)
 
-vkAR1 = poptR[-1][1]
-vkAR = [-x + 0j for x in vkAR1]
+ax1.set_xlabel(r"t", fontsize=28)
+ax1.set_ylabel(r"$C_R(t)$", fontsize=28)
+ax1.legend(loc=0, fontsize=12)
 
-# Imaginary part: use analytical value
+ax2.plot(
+    tlist, np.imag(env.correlation_function(tlist)),
+    "r", linewidth=2, label="Exact",
+)
+ax2.plot(
+    tlist, np.imag(envfit.correlation_function(tlist)),
+    "g--", linewidth=2, label="Fit", marker="o", markevery=50,
+)
+ax2.plot(
+    tlist, np.imag(padeEnv.correlation_function(tlist)),
+    "b--", linewidth=2, label=f"Pade (Nk={Nk})",
+)
 
-ckAI = [lam * gamma * (-1.0) + 0j]
-vkAI = [gamma + 0j]
+ax2.set_xlabel(r"t", fontsize=28)
+ax2.set_ylabel(r"$C_I(t)$", fontsize=28)
+ax2.legend(loc=0, fontsize=12)
+plt.show()
 ```
 
 ```{code-cell} ipython3
 with timer("RHS construction time"):
-    bath = BosonicBath(Q, ckAR, vkAR, ckAI, vkAI)
     # We reduce NC slightly here for speed of execution because we retain
     # 3 exponents in ckAR instead of 1. Please restore full NC for
     # convergence though:
-    HEOMFit = HEOMSolver(Hsys, bath, int(NC * 0.7), options=options)
+    HEOMFit = HEOMSolver(Hsys, (envfit, Q), int(NC * 0.7), options=options)
 
 with timer("ODE solver time"):
     resultFit = HEOMFit.run(rho0, tlist)
@@ -467,16 +396,10 @@ with timer("ODE solver time"):
 ## Simulation 5: Bloch-Redfield
 
 ```{code-cell} ipython3
-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))) "
-    "* ((1 / (exp(w * {beta}) - 1)) + 1)"
-).format(gamma=gamma, beta=beta, lam=lam)
-
 with timer("ODE solver time"):
     resultBR = brmesolve(
         Hsys, rho0, tlist,
-        a_ops=[[sigmaz(), DL]], sec_cutoff=0, options=options,
+        a_ops=[[sigmaz(), env]], sec_cutoff=0, options=options,
     )
 ```
 
@@ -551,7 +474,7 @@ with plt.rc_context(rcParams):
 ## About
 
 ```{code-cell} ipython3
-qutip.about()
+about()
 ```
 
 ## Testing
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..af48a007 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
@@ -1,15 +1,14 @@
 ---
 jupytext:
-  formats: ipynb,md:myst
   text_representation:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.17.0
 kernelspec:
+  name: python3
   display_name: Python 3 (ipykernel)
   language: python
-  name: python3
 ---
 
 # HEOM 1c: Spin-Bath model (Underdamped Case)
@@ -83,22 +82,12 @@ import time
 import numpy as np
 from matplotlib import pyplot as plt
 
-import qutip
-from qutip import (
-    basis,
-    brmesolve,
-    destroy,
-    expect,
-    qeye,
-    sigmax,
-    sigmaz,
-    tensor,
-)
-from qutip.solver.heom import (
-    HEOMSolver,
-    BosonicBath,
-    UnderDampedBath,
+from qutip import (about, basis, brmesolve, destroy, expect, qeye,
+                   sigmax, sigmaz, tensor)
+from qutip.core.environment import (
+    ExponentialBosonicEnvironment, UnderDampedEnvironment
 )
+from qutip.solver.heom import HEOMSolver
 
 %matplotlib inline
 ```
@@ -110,13 +99,13 @@ Let's define some helper functions for calculating correlation function expansio
 ```{code-cell} ipython3
 def cot(x):
     """ Vectorized cotangent of x. """
-    return 1. / np.tan(x)
+    return 1.0 / np.tan(x)
 ```
 
 ```{code-cell} ipython3
 def coth(x):
     """ Vectorized hyperbolic cotangent of x. """
-    return 1. / np.tanh(x)
+    return 1.0 / np.tanh(x)
 ```
 
 ```{code-cell} ipython3
@@ -125,8 +114,8 @@ def underdamped_matsubara_params(lam, gamma, T, nk):
         underdamped correlation functions.
     """
     Om = np.sqrt(w0**2 - (gamma / 2)**2)
-    Gamma = gamma / 2.
-    beta = 1. / T
+    Gamma = gamma / 2.0
+    beta = 1.0 / T
 
     ckAR = [
         (lam**2 / (4*Om)) * coth(beta * (Om + 1.0j * Gamma) / 2),
@@ -142,12 +131,9 @@ def underdamped_matsubara_params(lam, gamma, T, nk):
         -1.0j * Om + Gamma,
         1.0j * Om + Gamma,
     ]
-    vkAR.extend(
-        2 * np.pi * k * T + 0.j
-        for k in range(1, nk + 1)
-    )
+    vkAR.extend(2 * np.pi * k * T + 0.0j for k in range(1, nk + 1))
 
-    factor = 1. / 4
+    factor = 1.0 / 4
 
     ckAI = [
         -factor * lam**2 * 1.0j / Om,
@@ -222,8 +208,8 @@ And let us set up the system Hamiltonian, bath and system measurement operators:
 
 ```{code-cell} ipython3
 # Defining the system Hamiltonian
-eps = .5     # Energy of the 2-level system.
-Del = 1.0    # Tunnelling term
+eps = 0.5  # Energy of the 2-level system.
+Del = 1.0  # Tunnelling term
 Hsys = 0.5 * eps * sigmaz() + 0.5 * Del * sigmax()
 ```
 
@@ -237,11 +223,11 @@ rho0 = basis(2, 0) * basis(2, 0).dag()
 Q = sigmaz()  # coupling operator
 
 # Bath properties:
-gamma = .1  # cut off frequency
-lam = .5  # coupling strength
-w0 = 1.  # resonance frequency
-T = 1.
-beta = 1. / T
+gamma = 0.1  # cut off frequency
+lam = 0.5    # coupling strength
+w0 = 1.0     # resonance frequency
+T = 1.0
+beta = 1.0 / T
 
 # HEOM parameters:
 
@@ -292,7 +278,7 @@ The correlation functions are now very oscillatory, because of the Lorentzian pe
 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)
-    Gamma = gamma / 2.
+    Gamma = gamma / 2.0
     ek = 2 * np.pi * k / beta
 
     return (
@@ -304,7 +290,7 @@ def Mk(t, k, gamma, w0, beta):
 def c(t, Nk, lam, gamma, w0, beta):
     """ Calculate the correlation function for a vector of times, t. """
     Om = np.sqrt(w0**2 - (gamma / 2)**2)
-    Gamma = gamma / 2.
+    Gamma = gamma / 2.0
 
     Cr = (
         coth(beta * (Om + 1.0j * Gamma) / 2) * np.exp(1.0j * Om * t)
@@ -345,15 +331,13 @@ def plot_matsubara_correlation_function_contributions():
     """ Plot the underdamped correlation function. """
     t = np.linspace(0, 20, 1000)
 
-    M_Nk2 = np.sum([
-        Mk(t, k, gamma=gamma, w0=w0, beta=beta)
-        for k in range(1, 2 + 1)
-    ], 0)
+    M_Nk2 = np.sum(
+        [Mk(t, k, gamma=gamma, w0=w0, beta=beta) for k in range(1, 2 + 1)], 0
+    )
 
-    M_Nk100 = np.sum([
-        Mk(t, k, gamma=gamma, w0=w0, beta=beta)
-        for k in range(1, 100 + 1)
-    ], 0)
+    M_Nk100 = np.sum(
+        [Mk(t, k, gamma=gamma, w0=w0, beta=beta) for k in range(1, 100 + 1)], 0
+    )
 
     fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 8))
     axes.plot(t, np.real(M_Nk2), '-', color="black", label="Re[M(t)] Nk=2")
@@ -380,7 +364,7 @@ ckAR, vkAR, ckAI, vkAI = underdamped_matsubara_params(
 )
 ```
 
-Having created the lists which specify the bath correlation functions, we create a `BosonicBath` from them and pass the bath to the `HEOMSolver` class.
+Having created the lists which specify the bath correlation functions, we create an `ExponentialBosonicEnvironment` from them and pass the bath to the `HEOMSolver` class.
 
 The solver constructs the "right hand side" (RHS) determinining how the system and auxiliary density operators evolve in time. This can then be used to solve for dynamics or steady-state.
 
@@ -388,8 +372,8 @@ Below we create the bath and solver and then solve for the dynamics by calling `
 
 ```{code-cell} ipython3
 with timer("RHS construction time"):
-    bath = BosonicBath(Q, ckAR, vkAR, ckAI, vkAI)
-    HEOMMats = HEOMSolver(Hsys, bath, NC, options=options)
+    env = ExponentialBosonicEnvironment(ckAR, vkAR, ckAI, vkAI)
+    HEOMMats = HEOMSolver(Hsys, (env, Q), NC, options=options)
 
 with timer("ODE solver time"):
     resultMats = HEOMMats.run(rho0, tlist)
@@ -403,8 +387,8 @@ plot_result_expectations([
 ```
 
 In practice, one would not perform this laborious expansion for the underdamped correlation function, because
-QuTiP already has a class, `UnderDampedBath`, that can construct this bath for you. Nevertheless, knowing how
-to perform this expansion will allow you to construct your own baths for other spectral densities.
+QuTiP already has a class, `UnderDampedEnvironment`, that can construct this bath for you. Nevertheless, knowing how
+to perform this expansion is an useful skill.
 
 Below we show how to use this built-in functionality:
 
@@ -412,8 +396,9 @@ Below we show how to use this built-in functionality:
 # Compare to built-in under-damped bath:
 
 with timer("RHS construction time"):
-    bath = UnderDampedBath(Q, lam=lam, gamma=gamma, w0=w0, T=T, Nk=Nk)
-    HEOM_udbath = HEOMSolver(Hsys, bath, NC, options=options)
+    env = UnderDampedEnvironment(lam=lam, gamma=gamma, w0=w0, T=T)
+    env_approx = env.approximate("matsubara", Nk=Nk)
+    HEOM_udbath = HEOMSolver(Hsys, (env_approx, Q), NC, options=options)
 
 with timer("ODE solver time"):
     result_udbath = HEOM_udbath.run(rho0, tlist)
@@ -421,11 +406,38 @@ with timer("ODE solver time"):
 
 ```{code-cell} ipython3
 plot_result_expectations([
-    (result_udbath, P11p, 'b', "P11 (UnderDampedBath)"),
-    (result_udbath, P12p, 'r', "P12 (UnderDampedBath)"),
+    (result_udbath, P11p, 'b', "P11 (UnderDampedEnvironment)"),
+    (result_udbath, P12p, 'r', "P12 (UnderDampedEnvironment)"),
 ]);
 ```
 
+The `UnderDampedEnvironment` class also allows us to easily evaluate analytical expressions for the power spectrum, correlation function, and spectral density. In the following plots, the solid lines are the exact expressions, and the dashed lines are based on our approximation of the correlation function with a finite number of exponents. In this case, there is an excellent agreement.
+
+```{code-cell} ipython3
+w = np.linspace(-3, 3, 1000)
+w2 = np.linspace(0, 3, 1000)
+t = np.linspace(0, 10, 1000)
+env_cf = env.correlation_function(t)
+
+fig, axs = plt.subplots(2, 2)
+
+axs[0, 0].plot(w, env.power_spectrum(w))
+axs[0, 0].plot(w, env_approx.power_spectrum(w), "--")
+axs[0, 0].set(xlabel=r"$\omega$", ylabel=r"$S(\omega)$")
+axs[0, 1].plot(w2, env.spectral_density(w2))
+axs[0, 1].plot(w2, env_approx.spectral_density(w2), "--")
+axs[0, 1].set(xlabel=r"$\omega$", ylabel=r"$J(\omega)$")
+axs[1, 0].plot(t, np.real(env_cf))
+axs[1, 0].plot(t, np.real(env_approx.correlation_function(t)), "--")
+axs[1, 0].set(xlabel=r"$t$", ylabel=r"$C_{R}(t)$")
+axs[1, 1].plot(t, np.imag(env_cf))
+axs[1, 1].plot(t, np.imag(env_approx.correlation_function(t)), "--")
+axs[1, 1].set(xlabel=r"$t$", ylabel=r"$C_{I}(t)$")
+
+fig.tight_layout()
+plt.show()
+```
+
 ## Compare the results
 
 +++
@@ -433,17 +445,10 @@ plot_result_expectations([
 ### We can compare these results to those of the Bloch-Redfield solver in QuTiP:
 
 ```{code-cell} ipython3
-UD = (
-    f"2 * {lam}**2 * {gamma} / ( {w0}**4 * {beta}) if (w==0)"
-    " else "
-    f"2 * ({lam}**2 * {gamma} * w / (({w0}**2 - w**2)**2 + {gamma}**2 * w**2))"
-    f" * ((1 / (exp(w * {beta}) - 1)) + 1)"
-)
-
 with timer("ODE solver time"):
     resultBR = brmesolve(
         Hsys, rho0, tlist,
-        a_ops=[[sigmaz(), UD]], options=options,
+        a_ops=[[sigmaz(), env]], options=options,
     )
 ```
 
@@ -478,14 +483,14 @@ a = tensor(destroy(NRC), qeye(2))
 
 H0 = wa * a.dag() * a + Hsys_exp
 # interaction
-H1 = (g * (a.dag() + a) * Q_exp)
+H1 = g * (a.dag() + a) * Q_exp
 
 H = H0 + H1
 
 energies, states = H.eigenstates()
 rhoss = 0 * states[0] * states[0].dag()
 for kk, energ in enumerate(energies):
-    rhoss += (states[kk] * states[kk].dag() * np.exp(-beta * energies[kk]))
+    rhoss += states[kk] * states[kk].dag() * np.exp(-beta * energies[kk])
 rhoss = rhoss / rhoss.norm()
 
 P12RC = tensor(qeye(NRC), basis(2, 0) * basis(2, 1).dag())
@@ -542,7 +547,7 @@ with plt.rc_context(rcParams):
 ## About
 
 ```{code-cell} ipython3
-qutip.about()
+about()
 ```
 
 ## Testing
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..bf06c0f8 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
@@ -1,15 +1,14 @@
 ---
 jupytext:
-  formats: ipynb,md:myst
   text_representation:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.17.0
 kernelspec:
+  name: python3
   display_name: Python 3 (ipykernel)
   language: python
-  name: python3
 ---
 
 # HEOM 1d: Spin-Bath model, fitting of spectrum and correlation functions
@@ -27,127 +26,40 @@ The properties of the system are encoded in Hamiltonian, and a coupling operator
 
 The bosonic environment is implicitly assumed to obey a particular Hamiltonian ([see paper](https://arxiv.org/abs/2010.10806)), the parameters of which are encoded in the spectral density, and subsequently the free-bath correlation functions.
 
-In the example below we show how to model an Ohmic environment with exponential cut-off in two ways:
+In the example below we show how to model an Ohmic environment with exponential cut-off in three ways:
 
 * First we fit the spectral density with a set of underdamped brownian oscillator functions.
-
 * Second, we evaluate the correlation functions, and fit those with a certain choice of exponential functions.
+* Third, we use the built-in OhmicBath class, and explore the other approximation methods QuTiP offers
 
-In each case we will use the fit parameters to determine the correlation function expansion co-efficients needed to construct a description of the bath (i.e. a `BosonicBath` object) to supply to the `HEOMSolver` so that we can solve for the system dynamics.
+In each case we will use the fit parameters to determine the correlation function expansion co-efficients needed to construct a description of the bath (i.e. a `BosonicEnvironment` object) to supply to the `HEOMSolver` so that we can solve for the system dynamics.
 
 +++
 
 ## Setup
 
 ```{code-cell} ipython3
-import contextlib
-import dataclasses
-import time
+# mpmath is required for this tutorial,
+# for the evaluation of gamma and zeta
+# functions in the expression for the correlation:
+from mpmath import mp
 
 import numpy as np
 from matplotlib import pyplot as plt
-from scipy.optimize import curve_fit
-
-import qutip
-from qutip import (
-    basis,
-    expect,
-    liouvillian,
-    sigmax,
-    sigmaz,
-    spost,
-    spre,
-)
-from qutip.solver.heom import (
-    HEOMSolver,
-    BosonicBath,
-)
 
-# Import mpmath functions for evaluation of gamma and zeta
-# functions in the expression for the correlation:
+from qutip import about, basis, expect, sigmax, sigmaz
+from qutip.core.environment import BosonicEnvironment, OhmicEnvironment
+from qutip.solver.heom import HEOMSolver
 
-from mpmath import mp
+%matplotlib inline
 
 mp.dps = 15
 mp.pretty = True
-
-%matplotlib inline
-```
-
-## Helper functions
-
-Let's define some helper functions for plotting results and timing how long operations take:
-
-```{code-cell} ipython3
-def coth(x):
-    """ Vectorized hyperbolic cotangent of x. """
-    return 1. / np.tanh(x)
-```
-
-```{code-cell} ipython3
-def plot_result_expectations(plots, axes=None):
-    """ Plot the expectation values of operators as functions of time.
-
-        Each plot in plots consists of (solver_result,
-        measurement_operation, color, label).
-    """
-    if axes is None:
-        fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 8))
-        fig_created = True
-    else:
-        fig = None
-        fig_created = False
-
-    # add kw arguments to each plot if missing
-    plots = [p if len(p) == 5 else p + ({},) for p in plots]
-    for result, m_op, color, label, kw in plots:
-        exp = np.real(expect(result.states, m_op))
-        kw.setdefault("linewidth", 2)
-        if color == 'rand':
-            axes.plot(
-                result.times, exp,
-                c=np.random.rand(3,), label=label, **kw,
-            )
-        else:
-            axes.plot(result.times, exp, color, label=label, **kw)
-
-    if fig_created:
-        axes.legend(loc=0, fontsize=12)
-        axes.set_xlabel("t", fontsize=28)
-
-    return fig
-```
-
-```{code-cell} ipython3
-@contextlib.contextmanager
-def timer(label):
-    """ Simple utility for timing functions:
-
-        with timer("name"):
-            ... code to time ...
-    """
-    start = time.time()
-    yield
-    end = time.time()
-    print(f"{label}: {end - start}")
-```
-
-```{code-cell} ipython3
-# Solver options:
-
-options = {
-    "nsteps": 15000,
-    "store_states": True,
-    "rtol": 1e-14,
-    "atol": 1e-14,
-    "method": "vern9",
-    "progress_bar": "enhanced",
-}
 ```
 
 ## System and bath definition
 
-And let us set up the system Hamiltonian, bath and system measurement operators:
+Let us set up the system Hamiltonian, bath and system measurement operators:
 
 +++
 
@@ -155,12 +67,10 @@ And let us set up the system Hamiltonian, bath and system measurement operators:
 
 ```{code-cell} ipython3
 # Defining the system Hamiltonian
-eps = 0.0    # Energy of the 2-level system.
-Del = 0.2    # Tunnelling term
+eps = 0    # Energy of the 2-level system.
+Del = 0.2  # Tunnelling term
 Hsys = 0.5 * eps * sigmaz() + 0.5 * Del * sigmax()
-```
 
-```{code-cell} ipython3
 # Initial state of the system.
 rho0 = basis(2, 0) * basis(2, 0).dag()
 ```
@@ -208,17 +118,19 @@ def ohmic_correlation(t, alpha, wc, beta, s=1):
     """ The Ohmic bath correlation function as a function of t
         (and the bath parameters).
     """
-    corr = (
-        (1 / np.pi) * alpha * wc**(1 - s) * beta**(-(s + 1)) * mp.gamma(s + 1)
-    )
+    corr = (1 / np.pi) * alpha * wc ** (1 - s)
+    corr *= beta ** (-(s + 1)) * mp.gamma(s + 1)
     z1_u = (1 + beta * wc - 1.0j * wc * t) / (beta * wc)
     z2_u = (1 + 1.0j * wc * t) / (beta * wc)
     # Note: the arguments to zeta should be in as high precision as possible.
     # See http://mpmath.org/doc/current/basics.html#providing-correct-input
-    return np.array([
-        complex(corr * (mp.zeta(s + 1, u1) + mp.zeta(s + 1, u2)))
-        for u1, u2 in zip(z1_u, z2_u)
-    ], dtype=np.complex128)
+    return np.array(
+        [
+            complex(corr * (mp.zeta(s + 1, u1) + mp.zeta(s + 1, u2)))
+            for u1, u2 in zip(z1_u, z2_u)
+        ],
+        dtype=np.complex128,
+    )
 ```
 
 ```{code-cell} ipython3
@@ -226,49 +138,35 @@ def ohmic_spectral_density(w, alpha, wc):
     """ The Ohmic bath spectral density as a function of w
         (and the bath parameters).
     """
-    return w * alpha * np.e**(-w / wc)
+    return w * alpha * np.e ** (-w / wc)
 ```
 
 ```{code-cell} ipython3
 def ohmic_power_spectrum(w, alpha, wc, beta):
     """ The Ohmic bath power spectrum as a function of w
         (and the bath parameters).
+        We here obtain it naively using the Fluctuation-Dissipation Theorem,
+        but this fails at w=0 where the limit should be taken properly
     """
-    return (
-        w * alpha * np.e**(-abs(w) / wc) *
-        ((1 / (np.e**(w * beta) - 1)) + 1) * 2
-    )
+    bose = (1 / (np.e ** (w * beta) - 1)) + 1
+    return w * alpha * np.e ** (-abs(w) / wc) * 2 * bose
 ```
 
 ### Bath and HEOM parameters
 
 +++
 
-Finally, let's set the bath parameters we will work with and write down some measurement operators:
+Finally, let's set the bath parameters we will work with
 
 ```{code-cell} ipython3
-# Bath parameters:
-
-@dataclasses.dataclass
-class OhmicBathParameters:
-    """ Ohmic bath parameters. """
-    Q: object = dataclasses.field(default_factory=sigmaz, repr=False)
-    alpha: float = 3.25
-    T: float = 0.5
-    wc: float = 1.0
-    s: float = 1
-
-    def __post_init__(self):
-        self.beta = 1 / self.T
-
-    def replace(self, **kw):
-        return dataclasses.replace(self, **kw)
-
-
-obp = OhmicBathParameters()
+Q = sigmaz()
+alpha = 3.25
+T = 0.5
+wc = 1.0
+s = 1
 ```
 
-And set the cut-off for the HEOM hierarchy:
+and set the cut-off for the HEOM hierarchy:
 
 ```{code-cell} ipython3
 # HEOM parameters:
@@ -276,823 +174,670 @@ And set the cut-off for the HEOM hierarchy:
 # The max_depth defaults to 5 so that the notebook executes more
 # quickly. Change it to 11 to wait longer for more accurate results.
 max_depth = 5
+
+# options used for the differential equation solver
+# "bdf" integration method is faster here
+options = {
+    "nsteps": 15000,
+    "store_states": True,
+    "rtol": 1e-12,
+    "atol": 1e-12,
+    "method": "bdf",
+}
 ```
 
-## Building the HEOM bath by fitting the spectral density
+# Obtaining a decaying Exponential description of the environment
 
-+++
+In order to carry out our HEOM simulation, we need to express the correlation 
+function as a sum of decaying exponentials, that is we need to express it as 
 
-We begin by fitting the spectral density, using a series of $k$ underdamped harmonic oscillators case with the Meier-Tannor form (J. Chem. Phys. 111, 3365 (1999); https://doi.org/10.1063/1.479669):
+$$C(\tau)= \sum_{k=0}^{N-1}c_{k}e^{-\nu_{k}t}$$
 
-\begin{equation}
-J_{\mathrm approx}(\omega; a, b, c) = \sum_{i=0}^{k-1} \frac{2 a_i b_i w}{((w + c_i)^2 + b_i^2) ((w - c_i)^2 + b_i^2)}
-\end{equation}
+As the correlation function of the environment is tied to it's power spectrum via 
+a Fourier transform, such a representation of the correlation function implies a 
+power spectrum of the form
 
-where $a, b$ and $c$ are the fit parameters and each is a vector of length $k$.
+$$S(\omega)= \sum_{k}2 Re\left( \frac{c_{k}}{\nu_{k}- i \omega}\right)$$
+
+There are several ways one can obtain such a decomposition, in this tutorial we 
+will cover the following approaches:
+
+- Non-Linear Least Squares:
+    - On the Spectral Density (`sd`)
+    - On the Correlation function (`cf`)
+    - On the Power spectrum (`ps`)
+- Methods based on the Prony Polynomial
+    - Prony on the correlation function(`prony`)
+    - ESPRIT on the correlation function(`esprit`)
+- Methods based on rational Approximations
+    - The AAA algorithm on the Power Spectrum (`aaa`)
+    - ESPIRA-I on the correlation function and its FFT (`espira-I`)
+    - ESPIRA-II on the correlation function and its FFT (`espira-II`)
+
++++
+
+## Building User defined Bosonic Environments
+
+Before obtaining exponential approximations, we first need to construct a 
+`BosonicEnviroment` describing the exact environment.
+Here, we will briefly explain how to create a user-defined 
+`BosonicEnviroment` by specifying the spectral density. The same can be done 
+using either the correlation function or the power spectrum. For this example, we will
+use the Ohmic Spectral density we defined above:
 
 ```{code-cell} ipython3
-# Helper functions for packing the paramters a, b and c into a single numpy
-# array as required by SciPy's curve_fit:
+w = np.linspace(0, 25, 20000)
+J = ohmic_spectral_density(w, alpha, wc)
+```
 
-def pack(a, b, c):
-    """ Pack parameter lists for fitting. """
-    return np.concatenate((a, b, c))
+The `BosonicEnvironment` class has special constructors that can be used to 
+create environments from arbitrary spectral densities, correlation functions, or
+power spectrums. For example:
 
+```{code-cell} ipython3
+# From an array
+sd_env = BosonicEnvironment.from_spectral_density(J=J, wlist=w, T=T)
+```
 
-def unpack(params):
-    """ Unpack parameter lists for fitting. """
-    N = len(params) // 3
-    a = np.array(params[:N])
-    b = np.array(params[N:2 * N])
-    c = np.array(params[2 * N:])
-    return a, b, c
+Specifying the temperature is optional, but it allows us to automatically compute the corresponding power spectrum and correlation function. For example, the automatically computed power spectrum matches the analytically defined `ohmic_power_spectrum` function from above:
+
+```{code-cell} ipython3
+# Here we avoid w=0
+np.allclose(
+    sd_env.power_spectrum(w[1:]), ohmic_power_spectrum(w[1:], alpha, wc, 1 / T)
+)
 ```
 
+Specifying the Temperature also allows QuTiP to automatically compute the correlation function by fast Fourier transformation:
+
 ```{code-cell} ipython3
-# The approximate spectral density and a helper for fitting the approximate
-# spectral density to values calculated from the analytical formula:
+tlist = np.linspace(0, 10, 500)
+plt.plot(
+    tlist,
+    sd_env.correlation_function(tlist).real,
+    label="BosonicEnvironment FFT (Real Part)",
+)
+plt.plot(
+    tlist,
+    ohmic_correlation(tlist, alpha, wc, 1 / T).real,
+    "--",
+    label="Analytical (Real Part)",
+)
+plt.plot(
+    tlist,
+    np.imag(sd_env.correlation_function(tlist)),
+    label="BosonicEnvironment FFT (Imaginary Part)",
+)
+plt.plot(
+    tlist,
+    np.imag(ohmic_correlation(tlist, alpha, wc, 1 / T)),
+    "--",
+    label="Analytical (Imaginary Part)",
+)
+plt.ylabel("C(t)")
+plt.xlabel("t")
+plt.legend()
+plt.show()
+```
 
-def spectral_density_approx(w, a, b, c):
-    """ Calculate the fitted value of the function for the given
-        parameters.
-    """
-    return np.sum(
-        2 * a[:, None] * np.multiply.outer(b, w) / (
-            ((w + c[:, None])**2 + b[:, None]**2) *
-            ((w - c[:, None])**2 + b[:, None]**2)
-        ),
-        axis=0,
-    )
+Note that above, we constructed the `BosonicEnvironment` from the arrays `w` and `J`.
+Instead, one can also use a pure Python function.
+In that case, it is important to specify the parameter `wMax`, which is the cutoff frequency where the 
+spectral density or power spectrum has effectively decayed to zero. That is, for $\omega > \omega_{max}$, the function can be 
+considered to be essentially zero. The following is therefore equivalent to the environment that we used above:
 
+```{code-cell} ipython3
+# From a function
+sd_env2 = BosonicEnvironment.from_spectral_density(
+    ohmic_spectral_density, T=T, wMax=25 * wc, args={"alpha": alpha, "wc": wc}
+)
+```
 
-def fit_spectral_density(J, w, alpha, wc, N):
-    """ Fit the spectral density with N underdamped oscillators. """
-    sigma = [0.0001] * len(w)
+## Building the Exponential environment by fitting the spectral density
 
-    J_max = abs(max(J, key=abs))
+Once our `BosonicEnvironment` has been constructed, we can obtain a decaying
+exponential representation of the environment, via fitting either the spectral
+density, power spectrum or the correlation function. 
 
-    guesses = pack([J_max] * N, [wc] * N, [wc] * N)
-    lower_bounds = pack([-100 * J_max] * N, [0.1 * wc] * N, [0.1 * wc] * N)
-    upper_bounds = pack([100 * J_max] * N, [100 * wc] * N, [100 * wc] * N)
+We begin with a nonlinear-least-squares fit of the spectral density, using a series of $k$ underdamped harmonic oscillators case with the Meier-Tannor form (J. Chem. Phys. 111, 3365 (1999); https://doi.org/10.1063/1.479669):
 
-    params, _ = curve_fit(
-        lambda x, *params: spectral_density_approx(w, *unpack(params)),
-        w, J,
-        p0=guesses,
-        bounds=(lower_bounds, upper_bounds),
-        sigma=sigma,
-        maxfev=1000000000,
-    )
+\begin{equation}
+J_{\mathrm approx}(\omega; a, b, c) = \sum_{i=0}^{k-1} \frac{2 a_i b_i w}{((w + c_i)^2 + b_i^2) ((w - c_i)^2 + b_i^2)}
+\end{equation}
 
-    return unpack(params)
-```
+where $a, b$ and $c$ are the fit parameters and each is a vector of length $k$.
+The idea here is that we express our arbitrary spectral density as a sum of 
+underdamped spectral densities with different coefficients, for which the
+Matsubara decomposition is available.
 
-With the spectral density approximation $J_{\mathrm approx}(w; a, b, c)$ implemented above, we can now perform the fit and examine the results.
+The fit can be done easily using the `approximate` method. Its output is a tuple containing an `ExponentialBosonicEnvironment`
+and a dictionary that has all the relevant information about the fit performed.
+The goodness of the fit is measured via the normalized root mean squared error.
 
-```{code-cell} ipython3
-w = np.linspace(0, 25, 20000)
-J = ohmic_spectral_density(w, alpha=obp.alpha, wc=obp.wc)
+By default, the number of terms in the fit is increased automatically until the target accuracy 
+is reached or the maximum number allowed terms `Nmax` is reached. (The target accuracy can be set to None,
+then the fit is performed only with the specified number `Nmax` of exponents.)
 
-params_k = [
-    fit_spectral_density(J, w, alpha=obp.alpha, wc=obp.wc, N=i+1)
-    for i in range(4)
-]
+```{code-cell} ipython3
+# adding a small uncertainty "sigma" helps the fit routine
+approx_env, fitinfo = sd_env.approximate("sd", w, Nmax=4, sigma=0.0001)
 ```
 
-Let's plot the fit for each $k$ and examine how it improves with an increasing number of terms:
+To obtain an overview of the results of the fit we may take a look at the summary from the ``fitinfo``
 
 ```{code-cell} ipython3
-for k, params in enumerate(params_k):
-    lam, gamma, w0 = params
-    y = spectral_density_approx(w, lam, gamma, w0)
-    print(f"Parameters [k={k}]: lam={lam}; gamma={gamma}; w0={w0}")
-    plt.plot(w, J, w, y)
-    plt.show()
+print(fitinfo["summary"])
 ```
 
-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:
+Since the effective spectral density and power spectrum corresponding to the approximated correlation function are available through the `approx_env` object, we can compare them to the original:
 
 ```{code-cell} ipython3
-# The parameters for the fit with four terms:
+fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 5))
 
-lam, gamma, w0 = params_k[-1]
+ax1.plot(w, J, label="Original spectral density")
+ax1.plot(w, approx_env.spectral_density(w), "--", label="Effective fitted SD")
+ax1.set_xlabel(r"$\omega$")
+ax1.set_ylabel(r"$J$")
+ax1.legend()
 
-print(f"Parameters [k={len(params_k) - 1}]: lam={lam}; gamma={gamma}; w0={w0}")
+ax2.plot(w, np.abs(J - approx_env.spectral_density(w)), label="Error")
+ax2.set_xlabel(r"$\omega$")
+ax2.set_ylabel(r"$|J-J_{approx}|$")
+ax2.legend()
+
+fig.tight_layout()
+plt.show()
 ```
 
+Here we see a surprisingly large discrepancy in our approximated or effective spectral density. This happens because we are not using enough exponentials (i.e., not enough Matsubara terms) from each of the underdamped modes to have an appropiate fit. All modes use the same number of Matsubara terms; when not specified, the number defaults to $1$, which is not enough to model a bath with the temperature considered here. Let us repeat this with a larger number of exponents.
+
 ```{code-cell} ipython3
-# Plot the components of the fit separately:
+# 3 Matsubara terms per mode instead of one (default)
+approx_env, fitinfo = sd_env.approximate("sd", w, Nmax=4, Nk=3, sigma=0.0001)
 
-def spectral_density_ith_component(w, i, lam, gamma, w0):
-    """ Return the i'th term of the approximation for the spectral density. """
-    return (
-        2 * lam[i] * gamma[i] * w /
-        (((w + w0[i])**2 + gamma[i]**2) * ((w - w0[i])**2 + gamma[i]**2))
-    )
+fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10, 5))
+
+ax1.plot(w, J, label="Original spectral density")
+ax1.plot(w, approx_env.spectral_density(w), "--", label="Effective fitted SD")
+ax1.set_xlabel(r"$\omega$")
+ax1.set_ylabel(r"$J$")
+ax1.legend()
 
+ax2.plot(w, np.abs(J - approx_env.spectral_density(w)), label="Error")
+ax2.set_xlabel(r"$\omega$")
+ax2.set_ylabel(r"$|J-J_{approx}|$")
+ax2.legend()
 
-def plot_spectral_density_fit_components(J, w, lam, gamma, w0):
-    """ Plot the individual components of a fit to the spectral density. """
-    fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 8))
-    axes.plot(w, J, 'r--', linewidth=2, label="original")
+fig.tight_layout()
+plt.show()
+```
+
+Since the number of exponents increases simulation time, one should go with the least amount of exponents that correctly describe the bath properties (Power spectrum, Spectral density and the correlation function).
+
+Let's take a closer look at our last fit by plotting the contribution of each term of the fit:
+
+```{code-cell} ipython3
+def plot_fit_components(func, J, w, lam, gamma, w0):
+    """ Plot the individual components of a fit to the spectral density
+        and how they contribute to the full fit"""
+    plt.figure(figsize=(10, 5))
+    plt.plot(w, J, "r--", linewidth=2, label="original")
     for i in range(len(lam)):
-        axes.plot(
-            w, spectral_density_ith_component(w, i, lam, gamma, w0),
-            linewidth=2,
-            label=f"fit component {i}",
-        )
-
-    axes.set_xlabel(r'$w$', fontsize=28)
-    axes.set_ylabel(r'J', fontsize=28)
-    axes.legend()
+        component = func(w, lam[i], gamma[i], w0[i])
+        plt.plot(w, component, label=rf"$k={i+1}$")
+    plt.xlabel(r"$\omega$", fontsize=20)
+    plt.ylabel(r"$J(\omega)$", fontsize=20)
+    plt.legend(bbox_to_anchor=(1.04, 1))
+    plt.show()
 
-    return fig
+lam = fitinfo["params"][:, 0]
+gamma = fitinfo["params"][:, 1]
+w0 = fitinfo["params"][:, 2]
 
+def _sd_fit_model(wlist, a, b, c):
+    return (
+        2 * a * b * wlist
+        / (((wlist + c) ** 2 + b**2) * ((wlist - c) ** 2 + b**2))
+    )
 
-plot_spectral_density_fit_components(J, w, lam, gamma, w0);
+plot_fit_components(_sd_fit_model, J, w, lam, gamma, w0)
 ```
 
 And let's also compare the power spectrum of the fit and the analytical spectral density:
 
 ```{code-cell} ipython3
-def plot_power_spectrum(alpha, wc, beta, lam, gamma, w0, save=True):
+def plot_power_spectrum(alpha, wc, beta):
     """ Plot the power spectrum of a fit against the actual power spectrum. """
     w = np.linspace(-10, 10, 50000)
-
     s_orig = ohmic_power_spectrum(w, alpha=alpha, wc=wc, beta=beta)
-    s_fit = (
-        spectral_density_approx(w, lam, gamma, w0) *
-        ((1 / (np.e**(w * beta) - 1)) + 1) * 2
-    )
-
-    fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 8))
-    axes.plot(w, s_orig, 'r', linewidth=2, label="original")
-    axes.plot(w, s_fit, 'b', linewidth=2, label="fit")
+    s_fit = approx_env.power_spectrum(w)
+    fig, axes = plt.subplots(1, 1, sharex=True)
+    axes.plot(w, s_orig, "r", linewidth=2, label="original")
+    axes.plot(w, np.real(s_fit), "b", linewidth=2, label="fit")
 
-    axes.set_xlabel(r'$\omega$', fontsize=28)
-    axes.set_ylabel(r'$S(\omega)$', fontsize=28)
+    axes.set_xlabel(r"$\omega$", fontsize=20)
+    axes.set_ylabel(r"$S(\omega)$", fontsize=20)
     axes.legend()
 
-    if save:
-        fig.savefig('powerspectrum.eps')
-
-
-plot_power_spectrum(obp.alpha, obp.wc, obp.beta, lam, gamma, w0, save=False)
+plot_power_spectrum(alpha, wc, 1 / T)
 ```
 
-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.
+Now if we want to see the systems's behaviour as we change the parameters of the fit and the simulation, we may use this auxiliary function.
 
 ```{code-cell} ipython3
-def matsubara_coefficients_from_spectral_fit(lam, gamma, w0, beta, Q, Nk):
-    """ Calculate the Matsubara co-efficients for a fit to the spectral
-        density.
-    """
-    # initial 0 value with the correct dimensions:
-    terminator = 0. * spre(Q)
-    # the number of matsubara expansion terms to include in the terminator:
-    terminator_max_k = 1000
-
-    ckAR = []
-    vkAR = []
-    ckAI = []
-    vkAI = []
-
-    for lamt, Gamma, Om in zip(lam, gamma, w0):
-
-        ckAR.extend([
-            (lamt / (4 * Om)) * coth(beta * (Om + 1.0j * Gamma) / 2),
-            (lamt / (4 * Om)) * coth(beta * (Om - 1.0j * Gamma) / 2),
-        ])
-        for k in range(1, Nk + 1):
-            ek = 2 * np.pi * k / beta
-            ckAR.append(
-                (-2 * lamt * 2 * Gamma / beta) * ek /
-                (
-                    ((Om + 1.0j * Gamma)**2 + ek**2) *
-                    ((Om - 1.0j * Gamma)**2 + ek**2)
-                )
-            )
-
-        terminator_factor = 0
-        for k in range(Nk + 1, terminator_max_k):
-            ek = 2 * np.pi * k / beta
-            ck = (
-                (-2 * lamt * 2 * Gamma / beta) * ek /
-                (
-                    ((Om + 1.0j * Gamma)**2 + ek**2) *
-                    ((Om - 1.0j * Gamma)**2 + ek**2)
-                )
-            )
-            terminator_factor += ck / ek
-        terminator += terminator_factor * (
-            2 * spre(Q) * spost(Q.dag())
-            - spre(Q.dag() * Q)
-            - spost(Q.dag() * Q)
-        )
-
-        vkAR.extend([
-            -1.0j * Om + Gamma,
-            1.0j * Om + Gamma,
-        ])
-        vkAR.extend([
-            2 * np.pi * k * obp.T + 0.j
-            for k in range(1, Nk + 1)
-        ])
-
-        ckAI.extend([
-            -0.25 * lamt * 1.0j / Om,
-            0.25 * lamt * 1.0j / Om,
-        ])
-        vkAI.extend([
-            -(-1.0j * Om - Gamma),
-            -(1.0j * Om - Gamma),
-        ])
-
-    return ckAR, vkAR, ckAI, vkAI, terminator
-```
-
-```{code-cell} ipython3
-def generate_spectrum_results(obp, params, Nk, max_depth):
+def generate_spectrum_results(N, Nk, max_depth):
     """ Run the HEOM with the given bath parameters and
         and return the results of the evolution.
     """
-    lam, gamma, w0 = params
-    ckAR, vkAR, ckAI, vkAI, terminator = (
-        matsubara_coefficients_from_spectral_fit(
-            lam, gamma, w0, beta=obp.beta, Q=obp.Q, Nk=Nk,
-        )
+    approx_env, fitinfo = sd_env.approximate(
+        "sd", w, Nmax=N, Nk=Nk, sigma=0.0001, target_rmse=None
     )
-    Ltot = liouvillian(Hsys) + terminator
     tlist = np.linspace(0, 30 * np.pi / Del, 600)
 
-    with timer("RHS construction time"):
-        bath = BosonicBath(obp.Q, ckAR, vkAR, ckAI, vkAI)
-        HEOM_spectral_fit = HEOMSolver(
-            Ltot, bath, max_depth=max_depth, options=options,
-        )
-
-    with timer("ODE solver time"):
-        results_spectral_fit = (HEOM_spectral_fit.run(rho0, tlist))
+    print(f"Starting calculations for N={N}, Nk={Nk}"
+          f" and max_depth={max_depth} ... ")
 
+    HEOM_spectral_fit = HEOMSolver(
+        Hsys, (approx_env, Q), max_depth=max_depth,
+        options={**options, 'progress_bar': False},
+    )
+    results_spectral_fit = HEOM_spectral_fit.run(rho0, tlist)
     return results_spectral_fit
 ```
 
+```{code-cell} ipython3
+def plot_result_expectations(plots, axes=None):
+    """ Plot the expectation values of operators as functions of time.
+
+        Each plot in plots consists of (solver_result,
+        measurement_operation, color, label).
+    """
+    if axes is None:
+        fig, axes = plt.subplots(1, 1, sharex=True)
+        fig_created = True
+    else:
+        fig = None
+        fig_created = False
+
+    # add kw arguments to each plot if missing
+    plots = [p if len(p) == 5 else p + ({},) for p in plots]
+    for result, m_op, color, label, kw in plots:
+        exp = np.real(expect(result.states, m_op))
+        kw.setdefault("linewidth", 2)
+        if color == "rand":
+            axes.plot(result.times, exp, color=np.random.rand(3),
+                      label=label, **kw)
+        else:
+            axes.plot(result.times, exp, color, label=label, **kw)
+
+    if fig_created:
+        axes.legend(loc=0, fontsize=12)
+        axes.set_xlabel("t", fontsize=20)
+
+    return fig
+```
+
 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
 # Generate results for different number of lorentzians in fit:
 
 results_spectral_fit_pk = [
-    generate_spectrum_results(obp, params, Nk=1, max_depth=max_depth)
-    for params in params_k
+    generate_spectrum_results(n, Nk=1, max_depth=max_depth)
+    for n in range(1, 5)
 ]
 
 plot_result_expectations([
-    (
-        result, P11p, 'rand',
-        f"P11 (spectral fit) $k_J$={pk + 1}",
-    )
+    (result, P11p, "rand", f"P11 (spectral fit) $k$={pk + 1}")
     for pk, result in enumerate(results_spectral_fit_pk)
 ]);
 ```
 
 ```{code-cell} ipython3
-# generate results for different number of Matsubara terms per Lorentzian
-# for max number of Lorentzians:
+# generate results for different number of Matsubara terms per Lorentzian:
 
-Nk_list = range(2, 4)
+Nk_list = range(1, 4)
 results_spectral_fit_nk = [
-    generate_spectrum_results(obp, params_k[-1], Nk=Nk, max_depth=max_depth)
+    generate_spectrum_results(4, Nk=Nk, max_depth=max_depth)
     for Nk in Nk_list
 ]
 
 plot_result_expectations([
-    (
-        result, P11p, 'rand',
-        f"P11 (spectral fit) K={nk}",
-    )
+    (result, P11p, "rand", f"P11 (spectral fit) N_k={nk}")
     for nk, result in zip(Nk_list, results_spectral_fit_nk)
 ]);
 ```
 
 ```{code-cell} ipython3
-# Generate results for different depths:
+# generate results for different hierarchy depths:
 
-Nc_list = range(2, max_depth)
+Nc_list = range(3, max_depth+1)
 results_spectral_fit_nc = [
-    generate_spectrum_results(obp, params_k[-1], Nk=1, max_depth=Nc)
+    generate_spectrum_results(4, Nk=1, max_depth=Nc)
     for Nc in Nc_list
 ]
 
 plot_result_expectations([
-    (
-        result, P11p, 'rand',
-        f"P11 (spectral fit) $N_C={nc}$",
-    )
+    (result, P11p, "rand", f"P11 (spectral fit) $N_C={nc}$")
     for nc, result in zip(Nc_list, results_spectral_fit_nc)
 ]);
 ```
 
-We now combine the fitting and correlation function data into one large plot.
-
-```{code-cell} ipython3
-def correlation_approx_matsubara(t, ck, vk):
-    """ Calculate the approximate real or imaginary part of the
-        correlation function from the matsubara expansion co-efficients.
-    """
-    ck = np.array(ck)
-    vk = np.array(vk)
-    return np.sum(ck[:, None] * np.exp(-vk[:, None] * t), axis=0)
-```
-
-```{code-cell} ipython3
-def plot_cr_fit_vs_actual(t, ckAR, vkAR, C, axes):
-    """ Plot the C_R(t) fit. """
-    yR = correlation_approx_matsubara(t, ckAR, vkAR)
+## Obtaining a decaying exponential description by fitting the correlation function
 
-    axes.plot(
-        t, np.real(C),
-        "r", linewidth=3, label="Original",
-    )
-    axes.plot(
-        t, np.real(yR),
-        "g", dashes=[3, 3], linewidth=2, label="Reconstructed",
-    )
++++
 
-    axes.legend(loc=0)
-    axes.set_ylabel(r'$C_R(t)$', fontsize=28)
-    axes.set_xlabel(r'$t\;\omega_c$', fontsize=28)
-    axes.locator_params(axis='y', nbins=4)
-    axes.locator_params(axis='x', nbins=4)
-    axes.text(0.15, 0.85, "(a)", fontsize=28, transform=axes.transAxes)
+Having successfully fitted the spectral density and used the result to calculate the Matsubara expansion and terminator for the HEOM bosonic bath, we now proceed to the second case of fitting the correlation function itself instead.
 
+Here we fit the real and imaginary parts separately, using the following ansatz 
 
-def plot_ci_fit_vs_actual(t, ckAI, vkAI, C, axes):
-    """ Plot the C_I(t) fit. """
-    yI = correlation_approx_matsubara(t, ckAI, vkAI)
+$$C_R^F(t) = \sum_{i=1}^{k_R} c_R^ie^{-\gamma_R^i t}\cos(\omega_R^i t)$$
 
-    axes.plot(
-        t, np.imag(C),
-        "r", linewidth=3, label="Original",
-    )
-    axes.plot(
-        t, np.real(yI),
-        "g", dashes=[3, 3], linewidth=2, label="Reconstructed",
-    )
+$$C_I^F(t) = \sum_{i=1}^{k_I} c_I^ie^{-\gamma_I^i t}\sin(\omega_I^i t)$$
 
-    axes.legend(loc=0)
-    axes.set_ylabel(r'$C_I(t)$', fontsize=28)
-    axes.set_xlabel(r'$t\;\omega_c$', fontsize=28)
-    axes.locator_params(axis='y', nbins=4)
-    axes.locator_params(axis='x', nbins=4)
-    axes.text(0.80, 0.80, "(b)", fontsize=28, transform=axes.transAxes)
+Also this fit can easily be performed using the `approximate` method. The main difference with respect to the spectral density fit is that now we are perfoming two fits, one for the real part and another one for the imaginary part.
 
++++
 
-def plot_jw_fit_vs_actual(bath_fit, obp, axes):
-    """ Plot the J(w) fit. """
-    [lam, gamma, w0] = bath_fit
-    [alpha, wc] = [obp.alpha, obp.wc]
+Note that the ansatz is not good if $C_I^F(0) \neq 0$. In this case, the option `full_ansatz=True` allows for the usage of a 
+more general ansatz. The fit however tends to be significantly slower. We refer to the documentation for details.
 
-    w = np.linspace(0, 25, 20000)
+```{code-cell} ipython3
+def generate_corr_results(N, max_depth):
+    tlist = np.linspace(0, 30 * np.pi / Del, 600)
+    approx_env, fitinfo = sd_env.approximate(
+        "cf", tlist=tlist, Ni_max=N, Nr_max=N, maxfev=1e8, target_rmse=None
+    )
 
-    J_orig = ohmic_spectral_density(w, alpha=alpha, wc=wc)
-    J_fit = spectral_density_approx(w, lam, gamma, w0)
+    print(f"Starting calculations for N={N}"
+          f" and max_depth={max_depth} ... ")
 
-    axes.plot(
-        w, J_orig,
-        "r", linewidth=3, label=r"$J(\omega)$ original",
-    )
-    axes.plot(
-        w, J_fit,
-        "g", dashes=[3, 3], linewidth=2, label=r"$J(\omega)$ Fit $k_J = 4$",
+    HEOM_corr_fit = HEOMSolver(
+        Hsys, (approx_env, Q), max_depth=max_depth,
+        options={**options, 'progress_bar': False},
     )
+    results_corr_fit = HEOM_corr_fit.run(rho0, tlist)
+    return results_corr_fit
+```
 
-    axes.legend(loc=0)
-    axes.set_ylabel(r'$J(\omega)$', fontsize=28)
-    axes.set_xlabel(r'$\omega/\omega_c$', fontsize=28)
-    axes.locator_params(axis='y', nbins=4)
-    axes.locator_params(axis='x', nbins=4)
-    axes.text(0.15, 0.85, "(c)", fontsize=28, transform=axes.transAxes)
-
+```{code-cell} ipython3
+# Generate results for different number of exponentials in fit:
+results_corr_fit_pk = [
+    generate_corr_results(i, max_depth=max_depth)
+    for i in range(1, 4)
+]
 
-def plot_sw_fit_vs_actual(bath_fit, obp, axes):
-    """ Plot the S(w) fit. """
-    [lam, gamma, w0] = bath_fit
-    [alpha, wc, beta] = [obp.alpha, obp.wc, obp.beta]
+plot_result_expectations([
+    (result, P11p, "rand", f"P11 (correlation fit) k_R=k_I={pk + 1}")
+    for pk, result in enumerate(results_corr_fit_pk)
+]);
+```
 
-    # avoid the pole in the fit around zero:
-    w = np.concatenate(
-        [np.linspace(-10, -0.1, 5000),
-         np.linspace(0.1, 10, 5000)],
-    )
+```{code-cell} ipython3
+# Comparison plot
 
-    s_orig = ohmic_power_spectrum(w, alpha=alpha, wc=wc, beta=beta)
-    s_fit = (
-        spectral_density_approx(w, lam, gamma, w0) *
-        ((1 / (np.e**(w * beta) - 1)) + 1) * 2
-    )
+fig, axes = plt.subplots(1, 1, sharex=True, figsize=(10, 6))
 
-    axes.plot(w, s_orig, "r", linewidth=3, label="Original")
-    axes.plot(w, s_fit, "g", dashes=[3, 3], linewidth=2, label="Reconstructed")
+plot_result_expectations([
+    (results_corr_fit_pk[0], P11p, "y", "Correlation Fct. Fit $k_R=k_I=1$"),
+    (results_corr_fit_pk[2], P11p, "k", "Correlation Fct. Fit $k_R=k_I=3$"),
+    (results_spectral_fit_pk[0], P11p, "b", "Spectral Density Fit $k_J=1$"),
+    (results_spectral_fit_pk[3], P11p, "r-.", "Spectral Density Fit $k_J=4$"),
+], axes=axes)
 
-    axes.legend()
-    axes.set_ylabel(r'$S(\omega)$', fontsize=28)
-    axes.set_xlabel(r'$\omega/\omega_c$', fontsize=28)
-    axes.locator_params(axis='y', nbins=4)
-    axes.locator_params(axis='x', nbins=4)
-    axes.text(0.15, 0.85, "(d)", fontsize=28, transform=axes.transAxes)
+axes.set_yticks([0.6, 0.8, 1])
+axes.set_ylabel(r"$\rho_{11}$", fontsize=20)
+axes.set_xlabel(r"$t\;\omega_c$", fontsize=20)
+axes.legend(loc=0, fontsize=15);
+```
 
+# Using the Ohmic Environment class
 
-def plot_matsubara_spectrum_fit_vs_actual(
-    t, C, matsubara_fit, bath_fit, obp,
-):
-    """ Plot the Matsubara fit of the spectrum . """
-    fig = plt.figure(figsize=(12, 10))
-    grid = plt.GridSpec(2, 2, wspace=0.4, hspace=0.3)
+As the ohmic spectrum is popular in the modeling of open quantum systems, it has its own dedicated class. The results above can be reproduced quickly by using the `OhmicEnvironment` class. This allows for rapid implementation of fitted Ohmic baths.
 
-    [ckAR, vkAR, ckAI, vkAI] = matsubara_fit
+```{code-cell} ipython3
+obs = OhmicEnvironment(T, alpha, wc, s=1)
+tlist = np.linspace(0, 30 * np.pi / Del, 600)
+```
 
-    plot_cr_fit_vs_actual(
-        t, ckAR, vkAR, C,
-        axes=fig.add_subplot(grid[0, 0]),
-    )
-    plot_ci_fit_vs_actual(
-        t, ckAI, vkAI, C,
-        axes=fig.add_subplot(grid[0, 1]),
-    )
-    plot_jw_fit_vs_actual(
-        bath_fit, obp,
-        axes=fig.add_subplot(grid[1, 0]),
-    )
-    plot_sw_fit_vs_actual(
-        bath_fit, obp,
-        axes=fig.add_subplot(grid[1, 1]),
-    )
+Just like the other `BosonicEnvironment` we can obtain a decaying exponential 
+representation of the environment via the `approximate` function. Let us first do the same 
+methods we explored before:
 
-    return fig
+```{code-cell} ipython3
+sd_approx_env, fitinfo = obs.approximate(
+    method="sd", wlist=w, Nmax=4, Nk=3, sigma=0.0001, target_rmse=None
+)
+print(fitinfo["summary"])
+HEOM_ohmic_sd_fit = HEOMSolver(
+    Hsys, (sd_approx_env, Q), max_depth=max_depth, options=options
+)
+results_ohmic_sd_fit = HEOM_ohmic_sd_fit.run(rho0, tlist)
 ```
 
 ```{code-cell} ipython3
-t = np.linspace(0, 15, 100)
-C = ohmic_correlation(t, alpha=obp.alpha, wc=obp.wc, beta=obp.beta)
-
-ckAR, vkAR, ckAI, vkAI, terminator = (
-    matsubara_coefficients_from_spectral_fit(
-        lam, gamma, w0, beta=obp.beta, Q=obp.Q, Nk=1,
-    )
+cf_approx_env, fitinfo = obs.approximate(
+    method="cf", tlist=tlist, Nr_max=4, Ni_max=4, maxfev=1e8, target_rmse=None
 )
-
-matsubara_fit = [ckAR, vkAR, ckAI, vkAI]
-bath_fit = [lam, gamma, w0]
-
-plot_matsubara_spectrum_fit_vs_actual(
-    t, C, matsubara_fit,
-    bath_fit, obp,
-);
+print(fitinfo["summary"])
+HEOM_ohmic_corr_fit = HEOMSolver(
+    Hsys, (cf_approx_env, Q), max_depth=max_depth, options=options
+)
+results_ohmic_corr_fit = HEOM_ohmic_corr_fit.run(rho0, tlist)
 ```
 
-## Building the HEOM bath by fitting the correlation function
+## Other Approximation methods
+### Methods based on the Prony Polynomial 
 
-+++
+The Prony polynomial forms the mathematical foundation for many spectral analysis techniques that estimate frequencies, damping factors, and amplitudes of signals. These methods work by interpreting a given signal as a sum of complex exponentials and deriving a polynomial whose roots correspond to the frequencies or poles of the system.
 
-Having successfully fitted the spectral density and used the result to calculate the Matsubara expansion and terminator for the HEOM bosonic bath, we now proceed to the second case of fitting the correlation function itself instead.
+The methods consider a signal 
 
-Here we fit the real and imaginary parts seperately, using the following ansatz
+$$f(t)=\sum_{k=0}^{N-1} c_{k} e^{-\nu_{k} t} =\sum_{k=0}^{N-1} c_{k} z_{k}^{t}  $$
 
-$$C_R^F(t) = \sum_{i=1}^{k_R} c_R^ie^{-\gamma_R^i t}\cos(\omega_R^i t)$$
+The $z_{k}$ can be seen as the generalized eigenvalues of the matrix pencil
 
-$$C_I^F(t) = \sum_{i=1}^{k_I} c_I^ie^{-\gamma_I^i t}\sin(\omega_I^i t)$$
-
-```{code-cell} ipython3
-# The approximate correlation functions and a helper for fitting
-# the approximate correlation function to values calculated from
-# the analytical formula:
-
-def correlation_approx_real(t, a, b, c):
-    """ Calculate the fitted value of the function for the given parameters.
-    """
-    a = np.array(a)
-    b = np.array(b)
-    c = np.array(c)
-    return np.sum(
-        a[:, None] * np.exp(b[:, None] * t) * np.cos(c[:, None] * t),
-        axis=0,
-    )
-
-
-def correlation_approx_imag(t, a, b, c):
-    """ Calculate the fitted value of the function for the given parameters.
-    """
-    a = np.array(a)
-    b = np.array(b)
-    c = np.array(c)
-    return np.sum(
-        a[:, None] * np.exp(b[:, None] * t) * np.sin(c[:, None] * t),
-        axis=0,
-    )
+\begin{align}
+z_{j}  {\mathbf H}_{2N-L,L}(0) - {\mathbf H}_{2N-L,L}(1) = {\mathbf V}_{2N-L,M}
+({\mathbf z})   \mathrm{diag}  \Big( \left( (z_{j} - z_{k})\gamma_{k} 
+\right)_{k=1}^{M} \Big) {\mathbf V}_{L,M}({\mathbf z})^{T}
+\end{align}
 
 
-def fit_correlation_real(C, t, wc, N):
-    """ Fit the spectral density with N underdamped oscillators. """
-    sigma = [0.1] * len(t)
 
-    C_max = abs(max(C, key=abs))
+The amplitudes ($c_{k}$) can later be obtained by solving the least-squares Vandermonde system given by
 
-    guesses = pack([C_max] * N, [-wc] * N, [wc] * N)
-    lower_bounds = pack([-20 * C_max] * N, [-np.inf] * N, [0.] * N)
-    upper_bounds = pack([20 * C_max] * N, [0.1] * N, [np.inf] * N)
+$$ V_{N,M}(z)c = f $$
 
-    params, _ = curve_fit(
-        lambda x, *params: correlation_approx_real(t, *unpack(params)),
-        t, C,
-        p0=guesses,
-        bounds=(lower_bounds, upper_bounds),
-        sigma=sigma,
-        maxfev=1000000000,
-    )
+where $V_{N,M}(z)$ is the Vandermonde matrix given by
 
-    return unpack(params)
 
+$$V_{M,N}(z)=\begin{pmatrix} 
+1 &1 &\dots &1 \\
+z_{1} & z_{2} &\dots & z_{N} \\
+z_{1}^{2} & z_{2}^{2} &\dots & z_{N}^{2} \\
+\vdots & \vdots & \ddots & \vdots  \\
+z_{1}^{M} & z_{2}^{M} &\dots & z_{N}^{M} \\
+\end{pmatrix}$$
 
-def fit_correlation_imag(C, t, wc, N):
-    """ Fit the spectral density with N underdamped oscillators. """
-    sigma = [0.0001] * len(t)
+and $M$ is the length of the signal, and $N$ the number of exponents, and $f=f(t_{sample})$ is the signal evaluated in the sampling points,is a vector $c = (c_{1}, \dots, c_{N})$.
 
-    C_max = abs(max(C, key=abs))
+The main difference between the methods is the way one obtains the roots of the polynomial, typically whether this system is solved or a low rank approximation is found for the polynomial. [This article](https://academic.oup.com/imajna/article-abstract/43/2/789/6525860?redirectedFrom=fulltext) is a good reference, the QuTiP implementations are based on it and on the matlab implementations made available by the authors.
 
-    guesses = pack([-C_max] * N, [-2] * N, [1] * N)
-    lower_bounds = pack([-5 * C_max] * N, [-100] * N, [0.] * N)
-    upper_bounds = pack([5 * C_max] * N, [0.01] * N, [100] * N)
+The prony like methods include:
 
-    params, _ = curve_fit(
-        lambda x, *params: correlation_approx_imag(t, *unpack(params)),
-        t, C,
-        p0=guesses,
-        bounds=(lower_bounds, upper_bounds),
-        sigma=sigma,
-        maxfev=1000000000,
-    )
+- Prony
+- ESPRIT
+- ESPIRA
 
-    return unpack(params)
-```
+Though ESPIRA is prony like, since it is based on rational polynomial approximations
+we group it with other methods.
 
-```{code-cell} ipython3
-t = np.linspace(0, 15, 15000)
-C = ohmic_correlation(t, alpha=obp.alpha, wc=obp.wc, beta=obp.beta)
++++
 
-params_k_real = [
-    fit_correlation_real(np.real(C), t, wc=obp.wc, N=i+1)
-    for i in range(3)
-]
+##### Using the Original Prony Method on the Correlation Function
 
-params_k_imag = [
-    fit_correlation_imag(np.imag(C), t, wc=obp.wc, N=i+1)
-    for i in range(3)
-]
-```
+The method is available via `approximate` passing "prony" as method. Compared to the other approaches showed so far. The Prony based methods, shine on their simplicity no information needs to be known about the function, and one just needs to provide the sampling points, and the Number of Exponents one desires
 
 ```{code-cell} ipython3
-for k, params in enumerate(params_k_real):
-    lam, gamma, w0 = params
-    y = correlation_approx_real(t, lam, gamma, w0)
-    print(f"Parameters [k={k}]: lam={lam}; gamma={gamma}; w0={w0}")
-    plt.plot(t, np.real(C), label="C_R(t) analytic")
-    plt.plot(t, y, label=f"C_R(t) k={k + 1}")
-    plt.legend()
-    plt.show()
+tlist2 = np.linspace(0, 40, 100)
 ```
 
 ```{code-cell} ipython3
-for k, params in enumerate(params_k_imag):
-    lam, gamma, w0 = params
-    y = correlation_approx_imag(t, lam, gamma, w0)
-    print(f"Parameters [k={k}]: lam={lam}; gamma={gamma}; w0={w0}")
-    plt.plot(t, np.imag(C), label="C_I(t) analytic")
-    plt.plot(t, y, label=f"C_I(t) k={k + 1}")
-    plt.legend()
-    plt.show()
+prony_approx_env, fitinfo = obs.approximate("prony", tlist2, Nr=4)
+print(fitinfo["summary"])
+HEOM_ohmic_prony_fit = HEOMSolver(
+    Hsys, (prony_approx_env, Q), max_depth=max_depth, options=options
+)
+results_ohmic_prony_fit = HEOM_ohmic_prony_fit.run(rho0, tlist)
 ```
 
-Now we construct the `BosonicBath` co-efficients and frequencies from the fit to the correlation function:
+Similar to how we approximated via prony we can use ESPRIT, the main difference
+between both methods lies in the construction of the pencil matrix
 
 ```{code-cell} ipython3
-def matsubara_coefficients_from_corr_fit_real(lam, gamma, w0):
-    """ Return the matsubara coefficients for the imaginary part
-        of the correlation function.
-    """
-    ckAR = [0.5 * x + 0j for x in lam]  # the 0.5 is from the cosine
-    # extend the list with the complex conjugates:
-    ckAR.extend(np.conjugate(ckAR))
-
-    vkAR = [-x - 1.0j * y for x, y in zip(gamma, w0)]
-    vkAR.extend([-x + 1.0j * y for x, y in zip(gamma, w0)])
-
-    return ckAR, vkAR
-
-
-def matsubara_coefficients_from_corr_fit_imag(lam, gamma, w0):
-    """ Return the matsubara coefficients for the imaginary part
-        of the correlation function.
-    """
-    ckAI = [-0.5j * x for x in lam]  # the 0.5 is from the sine
-    # extend the list with the complex conjugates:
-    ckAI.extend(np.conjugate(ckAI))
-
-    vkAI = [-x - 1.0j * y for x, y in zip(gamma, w0)]
-    vkAI.extend([-x + 1.0j * y for x, y in zip(gamma, w0)])
-
-    return ckAI, vkAI
-```
-
-```{code-cell} ipython3
-ckAR, vkAR = matsubara_coefficients_from_corr_fit_real(*params_k_real[-1])
-ckAI, vkAI = matsubara_coefficients_from_corr_fit_imag(*params_k_imag[-1])
+esprit_approx_env, fitinfo = obs.approximate(
+    "esprit", tlist2, Nr=4, separate=False
+)
+print(fitinfo["summary"])
+HEOM_ohmic_es_fit = HEOMSolver(
+    Hsys, (esprit_approx_env, Q), max_depth=max_depth, options=options
+)
+results_ohmic_es_fit = HEOM_ohmic_es_fit.run(rho0, tlist)
 ```
 
-```{code-cell} ipython3
-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)
-    for ck, vk in zip(ckAR, vkAR):
-        S += (
-            2 * ck * np.real(vk) /
-            ((w - np.imag(vk))**2 + (np.real(vk)**2))
-        )
-    for ck, vk in zip(ckAI, vkAI):
-        S += (
-            2 * 1.0j * ck * np.real(vk) /
-            ((w - np.imag(vk))**2 + (np.real(vk)**2))
-        )
-    return S
-```
+## Fitting the power spectrum 
 
-```{code-cell} ipython3
-def plot_jw_correlation_fit_vs_actual(matsubara_fit, obp, axes):
-    """ Plot J(w) from the correlation fit. """
-    [ckAR, vkAR, ckAI, vkAI] = matsubara_fit
-    [alpha, wc] = [obp.alpha, obp.wc]
+So far all the methods covered fitted either the spectral density or the
+correlation function. Here we will fit the power spectrum.
 
-    w = np.linspace(0.001, 25, 20000)
+### AAA algorithm 
 
-    J_orig = ohmic_spectral_density(w, alpha=alpha, wc=wc)
-    J_fit = np.real(
-        corr_spectrum_approx(w, ckAR, vkAR, ckAI, vkAI) /
-        (((1 / (np.e**(w * obp.beta) - 1)) + 1) * 2)
-    )
+The Adaptive Antoulas–Anderson algorithm (AAA) is a method for the 
+approximation of a function in terms of a quotient of polynomials
 
-    axes.plot(
-        w, J_orig,
-        "r", linewidth=3, label=r"$J(\omega)$ original",
-    )
-    axes.plot(
-        w, J_fit,
-        "g", dashes=[3, 3], linewidth=2, label=r"$J(\omega)$ fit",
-    )
+\begin{align}
+    f(z) =\frac{q(z)}{p(z)} \approx \sum_{j=1}^{m} \frac{residues}{z-poles}
+\end{align}
 
-    axes.legend(loc=0)
-    axes.set_ylabel(r'$J(\omega)$', fontsize=28)
-    axes.set_xlabel(r'$\omega/\omega_c$', fontsize=28)
-    axes.locator_params(axis='y', nbins=4)
-    axes.locator_params(axis='x', nbins=4)
-    axes.text(3, 1.1, "(c)", fontsize=28)
+We don't use this method on the correlation function directly, but on the power spectrum .  After obtaining this rational polynomial form of the power spectrum one can recover the correlation function by noticing that 
 
+\begin{align}
+    s(\omega) = \int_{-\infty}^{\infty} dt e^{i \omega t} C(t)  = 2 \Re \left(\sum_{k} \frac{c_{k}}{\nu_{k}-i \omega} \right)
+\end{align}
 
-def plot_sw_correlation_fit_vs_actual(matsubara_fit, obp, axes):
-    """ Plot S(W) from the correlation fit. """
-    [ckAR, vkAR, ckAI, vkAI] = matsubara_fit
-    [alpha, wc, beta] = [obp.alpha, obp.wc, obp.beta]
+Which allows us to identify 
 
-    # avoid the pole in the fit around zero:
-    w = np.concatenate([
-        np.linspace(-10, -0.1, 5000),
-        np.linspace(0.1, 10, 5000),
-    ])
+\begin{align}
+    \nu_{k}= i \times poles \\
+    c_{k} = -i \times residues
+\end{align}
 
-    s_orig = ohmic_power_spectrum(w, alpha=alpha, wc=wc, beta=beta)
-    s_fit = corr_spectrum_approx(w, ckAR, vkAR, ckAI, vkAI)
+This method works best when the sampling points provided are in the logarithmic scale:
 
-    axes.plot(
-        w, s_orig,
-        "r", linewidth=3, label="Original",
-    )
-    axes.plot(
-        w, s_fit,
-        "g", dashes=[3, 3], linewidth=2, label="Reconstructed",
-    )
+```{code-cell} ipython3
+wlist = np.concatenate((-np.logspace(3, -8, 3500), np.logspace(-8, 3, 3500)))
 
-    axes.legend()
-    axes.set_ylabel(r'$S(\omega)$', fontsize=28)
-    axes.set_xlabel(r'$\omega/\omega_c$', fontsize=28)
-    axes.locator_params(axis='y', nbins=4)
-    axes.locator_params(axis='x', nbins=4)
-    axes.text(0.15, 0.85, "(d)", fontsize=28, transform=axes.transAxes)
+aaa_aprox_env, fitinfo = obs.approximate("aaa", wlist, Nmax=8, tol=1e-15)
+print(fitinfo["summary"])
+HEOM_ohmic_aaa_fit = HEOMSolver(
+    Hsys, (aaa_aprox_env, Q), max_depth=max_depth, options=options
+)
+results_ohmic_aaa_fit = HEOM_ohmic_aaa_fit.run(rho0, tlist)
+```
 
+### NLSQ on the power spectrum
 
-def plot_matsubara_correlation_fit_vs_actual(t, C, matsubara_fit, obp):
-    fig = plt.figure(figsize=(12, 10))
-    grid = plt.GridSpec(2, 2, wspace=0.4, hspace=0.3)
+On the first part of the tutorials we dealt with methods based on non-linear
+least squares. This is another one of those methods, but applied on the power 
+spectrum, compared to fitting the spectral density this is advantageous since
+we don't need to specify $N_k$, while compared to the correlation fit, it is 
+easier to obtain approximations that hold the KMS relation.
 
-    ckAR, vkAR, ckAI, vkAI = matsubara_fit
+we fit the power spectrum to a function of the form
 
-    plot_cr_fit_vs_actual(
-        t, ckAR, vkAR, C,
-        axes=fig.add_subplot(grid[0, 0]),
-    )
-    plot_ci_fit_vs_actual(
-        t, ckAI, vkAI, C,
-        axes=fig.add_subplot(grid[0, 1]),
-    )
-    plot_jw_correlation_fit_vs_actual(
-        matsubara_fit, obp,
-        axes=fig.add_subplot(grid[1, 0]),
-    )
-    plot_sw_correlation_fit_vs_actual(
-        matsubara_fit, obp,
-        axes=fig.add_subplot(grid[1, 1]),
-    )
-```
+$$S(\omega) = \sum_{k=1}^{N}\frac{2(a_k c_k + b_k (d_k - \omega))}
+{(\omega - d_k)^2 + c_k^2}= 2 \Re \left(\sum_{k} \frac{c_{k}}{\nu_{k}-i \omega} \right)$$
 
 ```{code-cell} ipython3
-t = np.linspace(0, 15, 100)
-C = ohmic_correlation(t, alpha=obp.alpha, wc=obp.wc, beta=obp.beta)
+w2 = np.concatenate((-np.linspace(10, 1e-2, 100), np.linspace(1e-2, 10, 100)))
 
-matsubara_fit = [ckAR, vkAR, ckAI, vkAI]
-
-plot_matsubara_correlation_fit_vs_actual(
-    t, C, matsubara_fit, obp,
+ps_approx_env, fitinfo = obs.approximate("ps", w2, Nmax=4)
+print(fitinfo["summary"])
+HEOM_ohmic_ps_fit = HEOMSolver(
+    Hsys, (ps_approx_env, Q), max_depth=max_depth, options=options
 )
+results_ohmic_ps_fit = HEOM_ohmic_ps_fit.run(rho0, tlist)
 ```
 
-```{code-cell} ipython3
-def generate_corr_results(params_real, params_imag, max_depth):
-    ckAR, vkAR = matsubara_coefficients_from_corr_fit_real(
-        *params_real
-    )
-    ckAI, vkAI = matsubara_coefficients_from_corr_fit_imag(
-        *params_imag
-    )
+### ESPIRA
 
-    tlist = np.linspace(0, 30 * np.pi / Del, 600)
+ESPIRA is a Prony-like method. While it takes a correlation function as 
+input, it exploits the relationship between parameter estimation (what we do 
+in Prony) and rational approximations. The rational approximation is done on the 
+DFT via the AAA algorithm, effectively using both information about the 
+power spectrum and the correlation function in the same fit.  
 
-    with timer("RHS construction time"):
-        bath = BosonicBath(obp.Q, ckAR, vkAR, ckAI, vkAI)
-        HEOM_corr_fit = HEOMSolver(
-            Hsys, bath, max_depth=max_depth, options=options,
-        )
+We have two implementations of ESPIRA, ESPIRA-I and ESPIRA-II. ESPIRA-I 
+is typically better, but in many cases especially when `separate=True`, 
+ESPIRA-II will yield better results. ESPIRA-II is recommended if 
+extremely slowly decaying exponents are expected. Otherwise ESPIRA-I is 
+recommended.
 
-    with timer("ODE solver time"):
-        results_corr_fit = (HEOM_corr_fit.run(rho0, tlist))
++++
 
-    return results_corr_fit
+##### ESPIRA I
 
+```{code-cell} ipython3
+tlist4 = np.linspace(0, 20, 1000)
 
-# Generate results for different number of lorentzians in fit:
-results_corr_fit_pk = [
-    print(f"{pk + 1}") or generate_corr_results(
-        params_real, params_imag, max_depth=max_depth,
-    )
-    for pk, (params_real, params_imag)
-    in enumerate(zip(params_k_real, params_k_imag))
-]
+espi_approx_env, fitinfo = obs.approximate("espira-I", tlist4, Nr=4)
+print(fitinfo["summary"])
+HEOM_ohmic_espira_fit = HEOMSolver(
+    Hsys, (espi_approx_env, Q), max_depth=max_depth, options=options
+)
+results_ohmic_espira_fit = HEOM_ohmic_espira_fit.run(rho0, tlist)
 ```
 
+##### ESPIRA-II
+
 ```{code-cell} ipython3
-plot_result_expectations([
-    (
-        result, P11p, 'rand',
-        f"P11 (correlation fit) k_R=k_I={pk + 1}",
-    )
-    for pk, result in enumerate(results_corr_fit_pk)
-]);
+espi2_approx_env, fitinfo = obs.approximate(
+    "espira-II", tlist4, Nr=4, Ni=4, separate=True
+)
+print(fitinfo["summary"])
+HEOM_ohmic_espira_fit2 = HEOMSolver(
+    Hsys, (espi2_approx_env, Q), max_depth=max_depth, options=options
+)
+results_ohmic_espira2_fit = HEOM_ohmic_espira_fit2.run(rho0, tlist)
 ```
 
+Finally we plot the dynamics obtained by the different methods
+
 ```{code-cell} ipython3
-fig, axes = plt.subplots(1, 1, sharex=True, figsize=(12, 7))
+fig, axes = plt.subplots(1, 1, sharex=True, figsize=(10, 6))
 
 plot_result_expectations([
-    (
-        results_corr_fit_pk[0], P11p,
-        'y', "Correlation Function Fit $k_R=k_I=1$",
-    ),
-    (
-        results_corr_fit_pk[2], P11p,
-        'y-.', "Correlation Function Fit $k_R=k_I=3$",
-    ),
-    (results_spectral_fit_pk[0], P11p, 'b', "Spectral Density Fit $k_J=1$"),
-    (results_spectral_fit_pk[2], P11p, 'g--', "Spectral Density Fit $k_J=3$"),
-    (results_spectral_fit_pk[3], P11p, 'r-.', "Spectral Density Fit $k_J=4$"),
+    (results_ohmic_corr_fit, P11p, "r", "Correlation Fct. Fit"),
+    (results_ohmic_sd_fit, P11p, "g--", "Spectral Density Fit"),
+    (results_ohmic_ps_fit, P11p, "g--", "Power Spectrum Fit Ohmic Bath"),
+    (results_ohmic_prony_fit, P11p, "k", " Prony Fit"),
+    (results_ohmic_es_fit, P11p, "b-.", "ESPRIT Fit"),
+    (results_ohmic_aaa_fit, P11p, "r-.", "Matrix AAA Fit"),
+    (results_ohmic_espira_fit, P11p, "k", "ESPIRA I Fit"),
+    (results_ohmic_espira2_fit, P11p, "--", "ESPIRA II Fit"),
 ], axes=axes)
 
-axes.set_yticks([0.6, 0.8, 1])
-axes.set_ylabel(r'$\rho_{11}$', fontsize=30)
-axes.set_xlabel(r'$t\;\omega_c$', fontsize=30)
-axes.legend(loc=0, fontsize=20);
+axes.set_ylabel(r"$\rho_{11}$", fontsize=20)
+axes.set_xlabel(r"$t\;\omega_c$", fontsize=20)
+axes.legend(loc=0, fontsize=15)
+axes.set_yscale("log")
 ```
 
 ## About
 
 ```{code-cell} ipython3
-qutip.about()
+about()
 ```
 
 ## Testing
@@ -1105,4 +850,30 @@ assert np.allclose(
     expect(P11p, results_spectral_fit_pk[3].states),
     rtol=1e-2,
 )
+assert np.allclose(
+    expect(P11p, results_ohmic_aaa_fit.states),
+    expect(P11p, results_spectral_fit_pk[3].states),
+    rtol=1e-2,
+)
+assert np.allclose(
+    expect(P11p, results_ohmic_prony_fit.states),
+    expect(P11p, results_spectral_fit_pk[3].states),
+    rtol=1e-2,
+)
+
+assert np.allclose(
+    expect(P11p, results_ohmic_es_fit.states),
+    expect(P11p, results_spectral_fit_pk[3].states),
+    rtol=1e-2,
+)
+assert np.allclose(
+    expect(P11p, results_ohmic_espira_fit.states),
+    expect(P11p, results_spectral_fit_pk[3].states),
+    rtol=1e-2,
+)
+assert np.allclose(
+    expect(P11p, results_ohmic_espira2_fit.states),
+    expect(P11p, results_spectral_fit_pk[3].states),
+    rtol=1e-2,
+)
 ```
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..af1fb51d 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
@@ -1,15 +1,14 @@
 ---
 jupytext:
-  formats: ipynb,md:myst
   text_representation:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.17.0
 kernelspec:
+  name: python3
   display_name: Python 3 (ipykernel)
   language: python
-  name: python3
 ---
 
 # HEOM 1e: Spin-Bath model (pure dephasing)
@@ -24,7 +23,7 @@ In this example we show the evolution of a single two-level system in contact wi
 
 The Bosonic environment is implicitly assumed to obey a particular Hamiltonian (see paper), the parameters of which are encoded in the spectral density, and subsequently the free-bath correlation functions.
 
-In the example below we show how to model the overdamped Drude-Lorentz Spectral Density, commonly used with the HEOM. We show how to do the Matsubara and Pade analytical decompositions, as well as how to fit the latter with a finite set of approximate exponentials.  This differs from examble 1a in that we assume that the system and coupling parts of the Hamiltonian commute, hence giving an analytically solvable ''pure dephasing'' model. This is a useful example to look at when introducing other approximations  (e.g., fitting of correlation functions) to check for validity/convergence against the analytical results.  (Note that, generally, for the fitting examples, the pure dephasing model is the 'worst possible case'.  
+In the example below we show how to model the overdamped Drude-Lorentz Spectral Density, commonly used with the HEOM. We show how to do the Matsubara and Pade analytical decompositions, as well as how to fit the latter with a finite set of approximate exponentials.  This differs from example 1a in that we assume that the system and coupling parts of the Hamiltonian commute, hence giving an analytically solvable ''pure dephasing'' model. This is a useful example to look at when introducing other approximations  (e.g., fitting of correlation functions) to check for validity/convergence against the analytical results.  (Note that, generally, for the fitting examples, the pure dephasing model is the 'worst possible case'.)
 
 ### Drude-Lorentz spectral density
 
@@ -71,24 +70,12 @@ import contextlib
 import time
 
 import numpy as np
-from matplotlib import pyplot as plt
 import scipy
-from scipy.optimize import curve_fit
-
-import qutip
-from qutip import (
-    basis,
-    expect,
-    liouvillian,
-    sigmax,
-    sigmaz,
-)
-from qutip.solver.heom import (
-    HEOMSolver,
-    BosonicBath,
-    DrudeLorentzBath,
-    DrudeLorentzPadeBath,
-)
+from matplotlib import pyplot as plt
+
+from qutip import about, basis, expect, liouvillian, sigmax, sigmaz
+from qutip.core.environment import DrudeLorentzEnvironment, system_terminator
+from qutip.solver.heom import HEOMSolver
 
 %matplotlib inline
 ```
@@ -100,12 +87,12 @@ Let's define some helper functions for calculating correlation function expansio
 ```{code-cell} ipython3
 def cot(x):
     """ Vectorized cotangent of x. """
-    return 1. / np.tan(x)
+    return 1.0 / np.tan(x)
 
 
 def coth(x):
     """ Vectorized hyperbolic cotangent of x. """
-    return 1. / np.tanh(x)
+    return 1.0 / np.tanh(x)
 ```
 
 ```{code-cell} ipython3
@@ -190,7 +177,7 @@ Q = sigmaz()  # coupling operator
 gamma = 0.5  # cut off frequency
 lam = 0.1  # coupling strength
 T = 0.5
-beta = 1. / T
+beta = 1.0 / T
 
 # HEOM parameters:
 # cut off parameter for the bath:
@@ -205,10 +192,10 @@ tlist = np.linspace(0, 50, 1000)
 
 ```{code-cell} ipython3
 # Define some operators with which we will measure the system
-# 1,1 element of density matrix - corresonding to groundstate
+# 1,1 element of density matrix - corresponding to groundstate
 P11p = basis(2, 0) * basis(2, 0).dag()
 P22p = basis(2, 1) * basis(2, 1).dag()
-# 1,2 element of density matrix  - corresonding to coherence
+# 1,2 element of density matrix  - corresponding to coherence
 P12p = basis(2, 0) * basis(2, 1).dag()
 ```
 
@@ -220,12 +207,19 @@ psi = (basis(2, 0) + basis(2, 1)).unit()
 rho0 = psi * psi.dag()
 ```
 
+We then define our environment, from which all the different simulations will 
+be obtained
+
+```{code-cell} ipython3
+env = DrudeLorentzEnvironment(lam=lam, gamma=gamma, T=T, Nk=Nk)
+```
+
 ## Simulation 1: Matsubara decomposition, not using Ishizaki-Tanimura terminator
 
 ```{code-cell} ipython3
 with timer("RHS construction time"):
-    bath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
-    HEOMMats = HEOMSolver(Hsys, bath, NC, options=options)
+    env_mats = env.approximate(method="matsubara", Nk=Nk)
+    HEOMMats = HEOMSolver(Hsys, (env_mats, Q), NC, options=options)
 
 with timer("ODE solver time"):
     resultMats = HEOMMats.run(rho0, tlist)
@@ -243,10 +237,11 @@ plot_result_expectations([
 
 ```{code-cell} ipython3
 with timer("RHS construction time"):
-    bath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
-    _, terminator = bath.terminator()
-    Ltot = liouvillian(Hsys) + terminator
-    HEOMMatsT = HEOMSolver(Ltot, bath, NC, options=options)
+    env_mats, delta = env.approximate(
+        method="matsubara", Nk=Nk, compute_delta=True
+    )
+    Ltot = liouvillian(Hsys) + system_terminator(Q, delta)
+    HEOMMatsT = HEOMSolver(Ltot, (env_mats, Q), NC, options=options)
 
 with timer("ODE solver time"):
     resultMatsT = HEOMMatsT.run(rho0, tlist)
@@ -268,8 +263,8 @@ As in example 1a, we can compare to Pade and Fitting approaches.
 
 ```{code-cell} ipython3
 with timer("RHS construction time"):
-    bath = DrudeLorentzPadeBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
-    HEOMPade = HEOMSolver(Hsys, bath, NC, options=options)
+    env_pade = env.approximate(method="pade", Nk=Nk)
+    HEOMPade = HEOMSolver(Hsys, (env_pade, Q), NC, options=options)
 
 with timer("ODE solver time"):
     resultPade = HEOMPade.run(rho0, tlist)
@@ -288,121 +283,12 @@ plot_result_expectations([
 ## Simulation 4: Fitting approach
 
 ```{code-cell} ipython3
-def c(t, Nk):
-    """ Calculates real and imag. parts of the correlation function
-        using Nk Matsubara terms.
-    """
-    vk = 2 * np.pi * T * np.arange(1, Nk)
-
-    result = (
-        lam * gamma * (-1.0j + cot(gamma * beta / 2.)) *
-        np.exp(-gamma * t[None, :])
-    )
-    result += np.sum(
-        (4 * lam * gamma * T * vk[:, None] / (vk[:, None]**2 - gamma**2)) *
-        np.exp(-vk[:, None] * t[None, :]),
-        axis=0,
-    )
-    result = result.squeeze(axis=0)
-
-    return result
-
-
-tlist_fit = np.linspace(0, 2, 10000)
-lmaxmats = 15000
-
-corr_ana = c(tlist_fit, lmaxmats)
-corrRana, corrIana = np.real(corr_ana), np.imag(corr_ana)
-```
-
-```{code-cell} ipython3
-def wrapper_fit_func(x, N, *args):
-    """ Wrapper for fitting function. """
-    a, b = args[0][:N], args[0][N:2*N]
-    return fit_func(x, a, b)
-
-
-def fit_func(x, a, b):
-    """ Fitting function. """
-    a = np.array(a)
-    b = np.array(b)
-    x = np.atleast_1d(np.array(x))
-    return np.sum(
-        a[:, None] * np.exp(b[:, None] * x[None, :]),
-        axis=0,
-    )
-
-
-def fitter(ans, tlist, i):
-    """ Compute the fit. """
-    upper_a = abs(max(ans, key=np.abs)) * 10
-    # set initial guess:
-    guess = [upper_a] * i + [0] * i
-    # set bounds: a's = anything, b's = negative
-    # sets lower bound
-    b_lower = [-upper_a] * i + [-np.inf] * i
-    # sets higher bound
-    b_higher = [upper_a] * i + [0] * i
-    param_bounds = (b_lower, b_higher)
-    p1, p2 = curve_fit(
-        lambda x, *params: wrapper_fit_func(x, i, params),
-        tlist,
-        ans,
-        p0=guess,
-        sigma=[0.01] * len(tlist),
-        bounds=param_bounds,
-        maxfev=1e8,
-    )
-    return p1[:i], p1[i:]
-
-
-# Fits of the real part with up to 4 exponents
-popt1 = []
-for i in range(4):
-    a, b = fitter(corrRana, tlist_fit, i + 1)
-    popt1.append((a, b))
-    y = fit_func(tlist_fit, a, b)
-    plt.plot(tlist_fit, corrRana, label="C_R(t)")
-    plt.plot(tlist_fit, y, label=f"Fit with k={i + 1}")
-    plt.xlabel("t")
-    plt.ylabel("C_R(t)")
-    plt.legend()
-    plt.show()
-
-# Fit of the imaginary part with 1 exponent
-popt2 = []
-for i in range(1):
-    a, b = fitter(corrIana, tlist_fit, i + 1)
-    popt2.append((a, b))
-    y = fit_func(tlist_fit, a, b)
-    plt.plot(tlist_fit, corrIana, label="C_I(t)")
-    plt.plot(tlist_fit, y, label=f"Fit with k={i + 1}")
-    plt.xlabel("t")
-    plt.ylabel("C_I(t)")
-    plt.legend()
-    plt.show()
-```
-
-```{code-cell} ipython3
-# Set the exponential coefficients from the fit parameters
-
-ckAR = popt1[-1][0]
-vkAR = -1 * popt1[-1][1]
-
-ckAI = popt2[-1][0]
-vkAI = -1 * popt2[-1][1]
-
-# The imaginary fit can also be determined analytically and is
-# a single term:
-#
-# ckAI = [complex(lam * gamma * (-1.0))]
-# vkAI = [complex(gamma)]
-```
-
-```{code-cell} ipython3
+tfit = np.linspace(0, 10, 1000)
 with timer("RHS construction time"):
-    bath = BosonicBath(Q, ckAR, vkAR, ckAI, vkAI)
-    HEOMFit = HEOMSolver(Hsys, bath, NC, options=options)
+    env_fit, _ = env.approximate(
+        method="cf", tlist=tfit, Ni_max=1, Nr_max=3, target_rmse=None
+    )
+    HEOMFit = HEOMSolver(Hsys, (env_fit, Q), NC, options=options)
 
 with timer("ODE solver time"):
     resultFit = HEOMFit.run(rho0, tlist)
@@ -434,10 +320,9 @@ def pure_dephasing_evolution_analytical(tlist, wq, ck, vk):
     integral: float
         The value of the integral function at time t.
     """
-    evolution = np.array([
-        np.exp(-1j * wq * t - correlation_integral(t, ck, vk))
-        for t in tlist
-    ])
+    evolution = np.array(
+        [np.exp(-1j * wq * t - correlation_integral(t, ck, vk)) for t in tlist]
+    )
     return evolution
 
 
@@ -471,17 +356,9 @@ def correlation_integral(t, ck, vk):
     integral: float
         The value of the integral function at time t.
     """
-    t1 = np.sum(
-        (ck / vk**2) *
-        (np.exp(vk * t) - 1)
-    )
-    t2 = np.sum(
-        (ck.conj() / vk.conj()**2) *
-        (np.exp(vk.conj() * t) - 1)
-    )
-    t3 = np.sum(
-        (ck / vk + ck.conj() / vk.conj()) * t
-    )
+    t1 = np.sum((ck / vk**2) * (np.exp(vk * t) - 1))
+    t2 = np.sum((ck.conj() / vk.conj()**2) * (np.exp(vk.conj() * t) - 1))
+    t3 = np.sum((ck / vk + ck.conj() / vk.conj()) * t)
     return 2 * (t1 + t2 - t3)
 ```
 
@@ -491,16 +368,12 @@ For the pure dephasing analytics, we just sum up as many matsubara terms as we c
 lmaxmats2 = 15000
 
 vk = [complex(-gamma)]
-vk.extend([
-    complex(-2. * np.pi * k * T)
-    for k in range(1, lmaxmats2)
-])
+vk.extend([complex(-2.0 * np.pi * k * T) for k in range(1, lmaxmats2)])
 
-ck = [complex(lam * gamma * (-1.0j + cot(gamma * beta / 2.)))]
-ck.extend([
-    complex(4 * lam * gamma * T * (-v) / (v**2 - gamma**2))
-    for v in vk[1:]
-])
+ck = [complex(lam * gamma * (-1.0j + cot(gamma * beta / 2.0)))]
+ck.extend(
+    [complex(4 * lam * gamma * T * (-v) / (v**2 - gamma**2)) for v in vk[1:]]
+)
 
 P12_ana = 0.5 * pure_dephasing_evolution_analytical(
     tlist, 0, np.asarray(ck), np.asarray(vk)
@@ -511,13 +384,13 @@ Alternatively, we can just do the integral of the propagator directly, without u
 
 ```{code-cell} ipython3
 def JDL(omega, lamc, omega_c):
-    return 2. * lamc * omega * omega_c / (omega_c**2 + omega**2)
+    return 2.0 * lamc * omega * omega_c / (omega_c**2 + omega**2)
 
 
 def integrand(omega, lamc, omega_c, Temp, t):
     return (
-        (-4. * JDL(omega, lamc, omega_c) / omega**2) *
-        (1. - np.cos(omega*t)) * (coth(omega/(2.*Temp)))
+        (-4.0 * JDL(omega, lamc, omega_c) / omega**2) *
+        (1.0 - np.cos(omega*t)) * (coth(omega/(2*Temp)))
         / np.pi
     )
 
@@ -564,7 +437,7 @@ axes.legend(loc=0, fontsize=12);
 ## About
 
 ```{code-cell} ipython3
-qutip.about()
+about()
 ```
 
 ## Testing
diff --git a/tutorials-v5/heom/heom-2-fmo-example.md b/tutorials-v5/heom/heom-2-fmo-example.md
index 7a301108..185da1ab 100644
--- a/tutorials-v5/heom/heom-2-fmo-example.md
+++ b/tutorials-v5/heom/heom-2-fmo-example.md
@@ -1,15 +1,14 @@
 ---
 jupytext:
-  formats: ipynb,md:myst
   text_representation:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.17.0
 kernelspec:
+  name: python3
   display_name: Python 3 (ipykernel)
   language: python
-  name: python3
 ---
 
 # HEOM 2: Dynamics in Fenna-Mathews-Olsen complex (FMO)
@@ -38,67 +37,15 @@ import time
 import numpy as np
 from matplotlib import pyplot as plt
 
-import qutip
-from qutip import (
-    Qobj,
-    basis,
-    brmesolve,
-    expect,
-    liouvillian,
-    mesolve,
-)
-from qutip.solver.heom import (
-    HEOMSolver,
-    DrudeLorentzBath,
-)
+from qutip import Qobj, about, basis, brmesolve, expect, liouvillian, mesolve
+from qutip.core.environment import DrudeLorentzEnvironment, system_terminator
+from qutip.solver.heom import HEOMSolver
 
 %matplotlib inline
 ```
 
 ## Helper functions
 
-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
-def cot(x):
-    """ Vectorized cotangent of x. """
-    return 1 / np.tan(x)
-```
-
-```{code-cell} ipython3
-def J0(energy):
-    """ Under-damped brownian oscillator spectral density. """
-    return 2 * lam * gamma * energy / (energy**2 + gamma**2)
-```
-
-```{code-cell} ipython3
-def J0_dephasing():
-    """ Under-damped brownian oscillator dephasing probability.
-
-        This returns the limit as w -> 0 of J0(w) * n_th(w, T) / T.
-    """
-    return 2 * lam * gamma / gamma**2
-```
-
-```{code-cell} ipython3
-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
-def dl_corr_approx(t, nk):
-    """ Drude-Lorenz correlation function approximation.
-
-        Approximates the correlation function at each time t to nk exponents.
-    """
-    c = lam * gamma * (-1.0j + cot(gamma / (2 * T))) * np.exp(-gamma * t)
-    for k in range(1, nk):
-        vk = 2 * np.pi * k * T
-        c += (4 * lam * gamma * T * vk / (vk**2 - gamma**2)) * np.exp(-vk * t)
-    return c
-```
-
 ```{code-cell} ipython3
 @contextlib.contextmanager
 def timer(label):
@@ -162,11 +109,15 @@ 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
+env = DrudeLorentzEnvironment(T=T, lam=lam, gamma=gamma)
+```
+
 ```{code-cell} ipython3
 wlist = np.linspace(0, 200 * 3e10 * 2 * np.pi, 100)
 tlist = np.linspace(0, 1e-12, 1000)
 
-J = J0(wlist) / (3e10*2*np.pi)
+J = env.spectral_density(wlist) / (3e10 * 2 * np.pi)
 
 fig, axes = plt.subplots(1, 2, sharex=False, figsize=(10, 3))
 
@@ -182,11 +133,11 @@ axes[0].legend()
 # Correlation plot:
 
 axes[1].plot(
-    tlist, np.real(dl_corr_approx(tlist, 10)),
+    tlist, np.real(env.correlation_function(tlist, 10)),
     color='r', ls='--', label="C(t) real",
 )
 axes[1].plot(
-    tlist, np.imag(dl_corr_approx(tlist, 10)),
+    tlist, np.imag(env.correlation_function(tlist, 10)),
     color='g', ls='--', label="C(t) imaginary",
 )
 axes[1].set_xlabel(r'$t$', fontsize=20)
@@ -211,24 +162,20 @@ NC = 4  # Use NC=8 for more precise results
 Nk = 0
 
 Q_list = []
-baths = []
 Ltot = liouvillian(Hsys)
+env_approx, delta = env.approximate(
+    method="matsubara", Nk=Nk, compute_delta=True
+)
 for m in range(7):
     Q = basis(7, m) * basis(7, m).dag()
     Q_list.append(Q)
-    baths.append(
-        DrudeLorentzBath(
-            Q, lam=lam, gamma=gamma, T=T, Nk=Nk,
-            tag=str(m)
-        )
-    )
-    _, terminator = baths[-1].terminator()
-    Ltot += terminator
+    Ltot += system_terminator(Q, delta)
 ```
 
 ```{code-cell} ipython3
 with timer("RHS construction time"):
-    HEOMMats = HEOMSolver(Hsys, baths, NC, options=options)
+    HEOMMats = HEOMSolver(Hsys, [(env_approx, Q) for Q in Q_list],
+                          NC, options=options)
 
 with timer("ODE solver time"):
     outputFMO_HEOM = HEOMMats.run(rho0, tlist)
@@ -260,8 +207,8 @@ for m in range(7):
 axes.set_title('HEOM solution', fontsize=24)
 axes.legend(loc=0)
 axes.set_xlim(0, 1000)
-plt.yticks([0., 0.5, 1], [0, 0.5, 1])
-plt.xticks([0., 500, 1000], [0, 500, 1000]);
+plt.yticks([0.0, 0.5, 1], [0, 0.5, 1])
+plt.xticks([0.0, 500, 1000], [0, 500, 1000]);
 ```
 
 ## Comparison with Bloch-Redfield solver
@@ -271,16 +218,10 @@ 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
-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))) * "
-    f"((1 / (exp((w) * {beta}) - 1)) + 1)"
-)
-
 with timer("BR ODE solver time"):
     outputFMO_BR = brmesolve(
         Hsys, rho0, tlist,
-        a_ops=[[Q, DL] for Q in Q_list],
+        a_ops=[[Q, env] for Q in Q_list],
         options=options,
     )
 ```
@@ -315,6 +256,15 @@ 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
+def J0_dephasing():
+    """ Under-damped brownian oscillator dephasing probability.
+
+        This returns the limit as w -> 0 of J0(w) * n_th(w, T) / T.
+    """
+    return 2 * lam * gamma / gamma**2
+```
+
 ```{code-cell} ipython3
 def get_collapse(H, T, dephasing=1):
     """ Calculate collapse operators for a given system H and
@@ -323,10 +273,7 @@ def get_collapse(H, T, dephasing=1):
     all_energy, all_state = H.eigenstates(sort="low")
     Nmax = len(all_energy)
 
-    Q_list = [
-        basis(Nmax, n) * basis(Nmax, n).dag()
-        for n in range(Nmax)
-    ]
+    Q_list = [basis(Nmax, n) * basis(Nmax, n).dag() for n in range(Nmax)]
 
     collapse_list = []
 
@@ -335,24 +282,18 @@ def get_collapse(H, T, dephasing=1):
             for k in range(j + 1, Nmax):
                 Deltajk = abs(all_energy[k] - all_energy[j])
                 if abs(Deltajk) > 0:
-                    rate = (
-                        np.abs(Q.matrix_element(
-                            all_state[j].dag(), all_state[k]
-                        ))**2 *
-                        2 * J0(Deltajk) * (n_th(Deltajk, T) + 1)
-                    )
+                    rate = np.abs(
+                        Q.matrix_element(all_state[j].dag(), all_state[k])
+                    ) ** 2 * env.power_spectrum(Deltajk)
                     if rate > 0.0:
                         # emission:
                         collapse_list.append(
                             np.sqrt(rate) * all_state[j] * all_state[k].dag()
                         )
 
-                    rate = (
-                        np.abs(Q.matrix_element(
-                            all_state[k].dag(), all_state[j]
-                        ))**2 *
-                        2 * J0(Deltajk) * n_th(Deltajk, T)
-                    )
+                    rate = np.abs(
+                        Q.matrix_element(all_state[k].dag(), all_state[j])
+                    ) ** 2 * env.power_spectrum(-Deltajk)
                     if rate > 0.0:
                         # absorption:
                         collapse_list.append(
@@ -362,10 +303,9 @@ def get_collapse(H, T, dephasing=1):
         if dephasing:
             for j in range(Nmax):
                 rate = (
-                    np.abs(Q.matrix_element(
-                        all_state[j].dag(), all_state[j])
-                    )**2 *
-                    J0_dephasing() * T
+                    np.abs(
+                        Q.matrix_element(all_state[j].dag(), all_state[j])
+                    ) ** 2 * env.power_spectrum(0) / 2
                 )
                 if rate > 0.0:
                     # emission:
@@ -376,7 +316,7 @@ def get_collapse(H, T, dephasing=1):
     return collapse_list
 ```
 
-Now we are able to switch the pure dephasing tersms on and off.
+Now we are able to switch the pure dephasing terms 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.
 
@@ -437,14 +377,14 @@ plt.xticks([0, 500, 1000], [0, 500, 1000])
 axes.legend(fontsize=18);
 ```
 
-And now we see that without the dephasing, the oscillations reappear. The full dynamics capture by the HEOM are still not capture by this simpler model, however.
+And now we see that without the dephasing, the oscillations reappear. The full dynamics captured by the HEOM are still not capture by this simpler model, however.
 
 +++
 
 ## About
 
 ```{code-cell} ipython3
-qutip.about()
+about()
 ```
 
 ## Testing
diff --git a/tutorials-v5/heom/heom-3-quantum-heat-transport.md b/tutorials-v5/heom/heom-3-quantum-heat-transport.md
index 8e956d6c..95611961 100644
--- a/tutorials-v5/heom/heom-3-quantum-heat-transport.md
+++ b/tutorials-v5/heom/heom-3-quantum-heat-transport.md
@@ -1,15 +1,14 @@
 ---
 jupytext:
-  formats: ipynb,md:myst
   text_representation:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.17.0
 kernelspec:
+  name: python3
   display_name: Python 3 (ipykernel)
   language: python
-  name: python3
 ---
 
 # HEOM 3: Quantum Heat Transport
@@ -57,11 +56,9 @@ import numpy as np
 import matplotlib.pyplot as plt
 
 import qutip as qt
-from qutip.solver.heom import (
-    DrudeLorentzPadeBath,
-    BathExponent,
-    HEOMSolver,
-)
+from qutip.core.environment import (CFExponent, DrudeLorentzEnvironment,
+                                    system_terminator)
+from qutip.solver.heom import HEOMSolver
 
 from ipywidgets import IntProgress
 from IPython.display import display
@@ -91,6 +88,7 @@ options = {
 @dataclasses.dataclass
 class SystemParams:
     """ System parameters and Hamiltonian. """
+
     epsilon: float = 1.0
     J12: float = 0.1
 
@@ -143,8 +141,16 @@ class BathParams:
         return qt.tensor(Q)
 
     def bath(self, Nk, tag=None):
-        return DrudeLorentzPadeBath(
-            self.Q(), self.lam, self.gamma, self.T, Nk, tag=tag
+        env = DrudeLorentzEnvironment(
+            lam=self.lam, gamma=self.gamma, T=self.T, tag=tag
+        )
+        env_approx, delta = env.approximate(
+            "pade", Nk=Nk, compute_delta=True, tag=tag
+        )
+        return (
+            (env_approx, self.Q()),
+            system_terminator(self.Q(), delta),
+            delta,
         )
 
     def replace(self, **kw):
@@ -201,9 +207,9 @@ def bath_heat_current(bath_tag, ado_state, hamiltonian, coupling_op, delta=0):
         [exp] = ado_state.exps(label)
         result += exp.vk * (coupling_op * ado_state.extract(label)).tr()
 
-        if exp.type == BathExponent.types['I']:
+        if exp.type == CFExponent.types["I"]:
             cI0 += exp.ck
-        elif exp.type == BathExponent.types['RI']:
+        elif exp.type == CFExponent.types["RI"]:
             cI0 += exp.ck2
 
     result -= 2 * cI0 * (coupling_op * coupling_op * ado_state.rho).tr()
@@ -290,14 +296,12 @@ tlist = np.linspace(0, 50, 250)
 ```{code-cell} ipython3
 H = sys.H()
 
-bath1 = bath_p1.bath(Nk, tag='bath 1')
+bath1, b1term, b1delta = bath_p1.bath(Nk, tag="bath 1")
 Q1 = bath_p1.Q()
 
-bath2 = bath_p2.bath(Nk, tag='bath 2')
+bath2, b2term, b2delta = bath_p2.bath(Nk, tag="bath 2")
 Q2 = bath_p2.Q()
 
-b1delta, b1term = bath1.terminator()
-b2delta, b2term = bath2.terminator()
 solver = HEOMSolver(
     qt.liouvillian(H) + b1term + b2term,
     [bath1, bath2],
@@ -318,7 +322,7 @@ We first plot $\langle \sigma_z^1 \rangle$ to see the time evolution of the syst
 
 ```{code-cell} ipython3
 fig, axes = plt.subplots(figsize=(8, 8))
-axes.plot(tlist, result.expect[0], 'r', linewidth=2)
+axes.plot(tlist, np.real(result.expect[0]), 'r', linewidth=2)
 axes.set_xlabel('t', fontsize=28)
 axes.set_ylabel(r"$\langle \sigma_z^1 \rangle$", fontsize=28);
 ```
@@ -382,15 +386,13 @@ def heat_currents(sys, bath_p1, bath_p2, Nk, NC, options):
     """ Calculate the steady sate heat currents for the given system and
         bath.
     """
-    bath1 = bath_p1.bath(Nk, tag="bath 1")
+
+    bath1, b1term, b1delta = bath_p1.bath(Nk, tag="bath 1")
     Q1 = bath_p1.Q()
 
-    bath2 = bath_p2.bath(Nk, tag="bath 2")
+    bath2, b2term, b2delta = bath_p2.bath(Nk, tag="bath 2")
     Q2 = bath_p2.Q()
 
-    b1delta, b1term = bath1.terminator()
-    b2delta, b2term = bath2.terminator()
-
     solver = HEOMSolver(
         qt.liouvillian(sys.H()) + b1term + b2term,
         [bath1, bath2],
@@ -445,18 +447,9 @@ def calculate_heat_current(J12, zb, Nk, progress=progress):
 
 # Calculate steady state currents for range of zeta_bars
 # for J12 = 0.01, 0.1 and 0.5:
-j1s = [
-    calculate_heat_current(0.01, zb, Nk)
-    for zb in zeta_bars
-]
-j2s = [
-    calculate_heat_current(0.1, zb, Nk)
-    for zb in zeta_bars
-]
-j3s = [
-    calculate_heat_current(0.5, zb, Nk)
-    for zb in zeta_bars
-]
+j1s = [calculate_heat_current(0.01, zb, Nk) for zb in zeta_bars]
+j2s = [calculate_heat_current(0.1, zb, Nk) for zb in zeta_bars]
+j3s = [calculate_heat_current(0.5, zb, Nk) for zb in zeta_bars]
 ```
 
 ## Create Plot
@@ -495,11 +488,3 @@ axes.legend(loc=0);
 ```{code-cell} ipython3
 qt.about()
 ```
-
-## Testing
-
-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
-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..538c5b4d 100644
--- a/tutorials-v5/heom/heom-4-dynamical-decoupling.md
+++ b/tutorials-v5/heom/heom-4-dynamical-decoupling.md
@@ -1,15 +1,14 @@
 ---
 jupytext:
-  formats: ipynb,md:myst
   text_representation:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.17.0
 kernelspec:
+  name: python3
   display_name: Python 3 (ipykernel)
   language: python
-  name: python3
 ---
 
 # HEOM 4: Dynamical decoupling of a non-Markovian environment
@@ -31,36 +30,17 @@ We first show the standard example of equally spaced pulses, and then consider t
 import numpy as np
 import matplotlib.pyplot as plt
 
-import qutip
-from qutip import (
-    QobjEvo,
-    basis,
-    expect,
-    ket2dm,
-    sigmax,
-    sigmaz,
-)
-from qutip.solver.heom import (
-    HEOMSolver,
-    DrudeLorentzPadeBath,
-)
+from qutip import (DrudeLorentzEnvironment, QobjEvo, about,
+                   basis, expect, ket2dm, sigmax, sigmaz)
+from qutip.solver.heom import HEOMSolver
 
-from ipywidgets import IntProgress
 from IPython.display import display
+from ipywidgets import IntProgress
 
 %matplotlib inline
 ```
 
-## Helper functions
-
-Let's define some helper functions for calculating the spectral density:
-
-```{code-cell} ipython3
-def dl_spectrum(w, lam, gamma):
-    """ Return the Drude-Lorentz spectral density. """
-    J = w * 2 * lam * gamma / (gamma**2 + w**2)
-    return J
-```
+## Solver options
 
 ```{code-cell} ipython3
 # Solver options:
@@ -95,10 +75,10 @@ H_sys = 0 * sigmaz()
 
 ```{code-cell} ipython3
 # Define some operators with which we will measure the system
-# 1,1 element of density matrix - corresonding to groundstate
+# 1,1 element of density matrix - corresponding to groundstate
 P11p = basis(2, 0) * basis(2, 0).dag()
 P22p = basis(2, 1) * basis(2, 1).dag()
-# 1,2 element of density matrix  - corresonding to coherence
+# 1,2 element of density matrix  - corresponding to coherence
 P12p = basis(2, 0) * basis(2, 1).dag()
 ```
 
@@ -115,7 +95,8 @@ Q = sigmaz()
 # number of terms to keep in the expansion of the bath correlation function:
 Nk = 3
 
-bath = DrudeLorentzPadeBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
+env = DrudeLorentzEnvironment(lam=lam, gamma=gamma, T=T)
+env_approx = env.approximate(method="pade", Nk=Nk)
 ```
 
 ```{code-cell} ipython3
@@ -168,12 +149,15 @@ Let's start by plotting the spectral density of our Drude-Lorentz bath:
 
 ```{code-cell} ipython3
 wlist = np.linspace(0, 0.5, 1000)
-J = dl_spectrum(wlist, lam, gamma)
+J = env.spectral_density(wlist)
+J_approx = env_approx.spectral_density(wlist)
 
 fig, axes = plt.subplots(1, 1, figsize=(8, 8))
-axes.plot(wlist, J, 'r', linewidth=2)
-axes.set_xlabel(r'$\omega$', fontsize=28)
-axes.set_ylabel(r'J', fontsize=28);
+axes.plot(wlist, J, "r", linewidth=2)
+axes.plot(wlist, J_approx, "b--", linewidth=2)
+
+axes.set_xlabel(r"$\omega$", fontsize=28)
+axes.set_ylabel(r"J", fontsize=28);
 ```
 
 ## Dynamic decoupling with fast and slow pulses
@@ -194,14 +178,14 @@ rho0 = (basis(2, 1) + basis(2, 0)).unit()
 rho0 = ket2dm(rho0)
 
 # without pulses
-hsolver = HEOMSolver(H_sys, bath, NC, options=options)
+hsolver = HEOMSolver(H_sys, (env_approx, Q), NC, options=options)
 outputnoDD = hsolver.run(rho0, tlist)
 
 # with pulses
 drive_fast = drive(amplitude=0.5, delay=20, integral=np.pi / 2)
-H_d = qutip.QobjEvo([H_sys, [H_drive, drive_fast]])
+H_d = QobjEvo([H_sys, [H_drive, drive_fast]])
 
-hsolver = HEOMSolver(H_d, bath, NC, options=options)
+hsolver = HEOMSolver(H_d, (env_approx, Q), NC, options=options)
 outputDD = hsolver.run(rho0, tlist)
 ```
 
@@ -211,14 +195,14 @@ And now the longer slower pulses:
 # Slow driving (longer, small amplitude pulses)
 
 # without pulses
-hsolver = HEOMSolver(H_sys, bath, NC, options=options)
+hsolver = HEOMSolver(H_sys, (env_approx, Q), NC, options=options)
 outputnoDDslow = hsolver.run(rho0, tlist)
 
 # with pulses
-drive_slow = drive(amplitude=0.01, delay=20, integral=np.pi/2)
+drive_slow = drive(amplitude=0.01, delay=20, integral=np.pi / 2)
 H_d = QobjEvo([H_sys, [H_drive, drive_slow]])
 
-hsolver = HEOMSolver(H_d, bath, NC, options=options)
+hsolver = HEOMSolver(H_d, (env_approx, Q), NC, options=options)
 outputDDslow = hsolver.run(rho0, tlist)
 ```
 
@@ -233,9 +217,9 @@ def plot_dd_results(outputnoDD, outputDD, outputDDslow):
     tlist = outputDD.times
 
     P12 = basis(2, 1) * basis(2, 0).dag()
-    P12DD = qutip.expect(outputDD.states, P12)
-    P12noDD = qutip.expect(outputnoDD.states, P12)
-    P12DDslow = qutip.expect(outputDDslow.states, P12)
+    P12DD = expect(outputDD.states, P12)
+    P12noDD = expect(outputnoDD.states, P12)
+    P12DDslow = expect(outputDDslow.states, P12)
 
     plt.sca(axes[0])
     plt.yticks([0, 0.25, 0.5], [0, 0.25, 0.5])
@@ -267,7 +251,7 @@ def plot_dd_results(outputnoDD, outputDD, outputDDslow):
     pulseslow = [drive_slow(t) for t in tlist]
 
     plt.sca(axes[1])
-    plt.yticks([0., 0.25, 0.5], [0, 0.25, 0.5])
+    plt.yticks([0, 0.25, 0.5], [0, 0.25, 0.5])
 
     axes[1].plot(
         tlist, pulse,
@@ -331,10 +315,9 @@ def cummulative_delay_fractions(N):
 
         as the cummulative delay after the j'th pulse.
     """
-    return np.array([
-        np.sin((np.pi / 2) * (j / (N + 1)))**2
-        for j in range(0, N + 1)
-    ])
+    return np.array(
+        [np.sin((np.pi / 2) * (j / (N + 1)))**2 for j in range(0, N + 1)]
+    )
 
 
 def drive_opt(amplitude, avg_delay, integral, N):
@@ -368,7 +351,7 @@ On the same axes we plot the individual $j^{th}$ delays as a fraction of the ave
 def plot_cummulative_delay_fractions(N):
     cummulative = cummulative_delay_fractions(N)
     individual = (cummulative[1:] - cummulative[:-1]) * N
-    plt.plot(np.arange(0, N + 1), cummulative, label="Cummulative delay")
+    plt.plot(np.arange(0, N + 1), cummulative, label="Cumulative delay")
     plt.plot(np.arange(0, N), individual, label="j'th delay")
     plt.xlabel("j")
     plt.ylabel("Fraction of delay")
@@ -448,14 +431,14 @@ def simulate_100_pulses(lam, gamma, T, NC, Nk):
     duration = integral / amplitude
     delay = avg_cycle_time - duration
 
-    bath = DrudeLorentzPadeBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
-
+    env = DrudeLorentzEnvironment(lam=lam, gamma=gamma, T=T)
+    env_approx = env.approximate("pade", Nk=Nk)
     # Equally spaced pulses:
 
     pulse_eq = drive(amplitude, delay, integral)
     H_d = QobjEvo([H_sys, [H_drive, pulse_eq]])
 
-    hsolver = HEOMSolver(H_d, bath, NC, options=options)
+    hsolver = HEOMSolver(H_d, (env_approx, Q), NC, options=options)
     result = hsolver.run(rho0, tlist)
 
     P12_eq = expect(result.states[-1], P12p)
@@ -466,7 +449,7 @@ def simulate_100_pulses(lam, gamma, T, NC, Nk):
     pulse_opt = drive_opt(amplitude, delay, integral, N)
     H_d = QobjEvo([H_sys, [H_drive, pulse_opt]])
 
-    hsolver = HEOMSolver(H_d, bath, NC, options=options)
+    hsolver = HEOMSolver(H_d, (env_approx, Q), NC, options=options)
     result = hsolver.run(rho0, tlist)
 
     P12_opt = expect(result.states[-1], P12p)
@@ -519,13 +502,5 @@ And now you know about dynamically decoupling a qubit from its environment!
 ## About
 
 ```{code-cell} ipython3
-qutip.about()
-```
-
-## Testing
-
-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
-assert 1 == 1
+about()
 ```
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..7746cb15 100644
--- a/tutorials-v5/heom/heom-5a-fermions-single-impurity-model.md
+++ b/tutorials-v5/heom/heom-5a-fermions-single-impurity-model.md
@@ -1,15 +1,14 @@
 ---
 jupytext:
-  formats: ipynb,md:myst
   text_representation:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.17.0
 kernelspec:
+  name: python3
   display_name: Python 3 (ipykernel)
   language: python
-  name: python3
 ---
 
 # HEOM 5a: Fermionic single impurity model
@@ -18,7 +17,7 @@ kernelspec:
 
 ## Introduction
 
-Here we model a single fermion coupled to two electronic leads or reservoirs (e.g.,  this can describe a single quantum dot, a molecular transistor, etc).  Note that in this implementation we primarily follow the definitions used by Christian Schinabeck in his dissertation https://opus4.kobv.de/opus4-fau/files/10984/DissertationChristianSchinabeck.pdf and related publications.
+Here we model a single fermion coupled to two electronic leads or reservoirs (e.g.,  this can describe a single quantum dot, a molecular transistor, etc).  Note that in this implementation we primarily follow the definitions used by Christian Schinabeck in his dissertation https://open.fau.de/items/36fdd708-a467-4b59-bf4e-4a2110fbc431 and related publications.
 
 Notation:
 
@@ -75,23 +74,15 @@ import dataclasses
 import time
 
 import numpy as np
-import matplotlib.pyplot as plt
 from scipy.integrate import quad
+import matplotlib.pyplot as plt
 
-import qutip
-from qutip import (
-    basis,
-    destroy,
-    expect,
-)
-from qutip.solver.heom import (
-    HEOMSolver,
-    LorentzianBath,
-    LorentzianPadeBath,
-)
+from qutip import about, basis, destroy, expect
+from qutip.core.environment import LorentzianEnvironment
+from qutip.solver.heom import HEOMSolver
 
-from ipywidgets import IntProgress
 from IPython.display import display
+from ipywidgets import IntProgress
 
 %matplotlib inline
 ```
@@ -164,7 +155,7 @@ class LorentzianBathParameters:
         if self.lead == "L":
             self.mu = self.theta / 2.0
         else:
-            self.mu = - self.theta / 2.0
+            self.mu = -self.theta / 2.0
 
     def J(self, w):
         """ Spectral density. """
@@ -176,7 +167,7 @@ class LorentzianBathParameters:
         return fF(x)
 
     def lamshift(self, w):
-        """ Return the lamshift. """
+        """ Return the lamb shift. """
         return 0.5 * (w - self.mu) * self.J(w) / self.W
 
     def replace(self, **kw):
@@ -280,17 +271,22 @@ rho0 = basis(2, 0) * basis(2, 0).dag()
 
 Nk = 10  # Number of exponents to retain in the expansion of each bath
 
-bathL = LorentzianPadeBath(
-    bath_L.Q, bath_L.gamma, bath_L.W, bath_L.mu, bath_L.T,
-    Nk, tag="L",
+envL = LorentzianEnvironment(
+    bath_L.T, bath_L.mu, bath_L.gamma, bath_L.W,
 )
-bathR = LorentzianPadeBath(
-    bath_R.Q, bath_R.gamma, bath_R.W, bath_R.mu, bath_R.T,
-    Nk, tag="R",
+envL_pade = envL.approx_by_pade(Nk=Nk, tag="L")
+envR = LorentzianEnvironment(
+    bath_R.T, bath_R.mu, bath_R.gamma, bath_R.W,
 )
+envR_pade = envR.approx_by_pade(Nk=Nk, tag="R")
 
 with timer("RHS construction time"):
-    solver_pade = HEOMSolver(H, [bathL, bathR], max_depth=2, options=options)
+    solver_pade = HEOMSolver(
+        H,
+        [(envL_pade, bath_L.Q), (envR_pade, bath_R.Q)],
+        max_depth=2,
+        options=options,
+    )
 
 with timer("ODE solver time"):
     result_pade = solver_pade.run(rho0, tlist)
@@ -325,17 +321,17 @@ Now let us do the same for the Matsubara expansion:
 ```{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",
-)
+envL_mats = envL.approx_by_matsubara(Nk=Nk, tag="L")
+envR_mats = envR.approx_by_matsubara(Nk=Nk, tag="R")
+
 
 with timer("RHS construction time"):
-    solver_mats = HEOMSolver(H, [bathL, bathR], max_depth=2, options=options)
+    solver_mats = HEOMSolver(
+        H,
+        [(envL_mats, bath_L.Q), (envR_mats, bath_R.Q)],
+        max_depth=2,
+        options=options,
+    )
 
 with timer("ODE solver time"):
     result_mats = solver_mats.run(rho0, tlist)
@@ -391,7 +387,7 @@ def analytical_steady_state_current(bath_L, bath_R, e1):
             bath_L.J(w) * bath_R.J(w) * (bath_L.fF(w) - bath_R.fF(w)) /
             (
                 (bath_L.J(w) + bath_R.J(w))**2 +
-                4*(w - e1 - bath_L.lamshift(w) - bath_R.lamshift(w))**2
+                4 * (w - e1 - bath_L.lamshift(w) - bath_R.lamshift(w))**2
             )
         )
 
@@ -437,8 +433,7 @@ def state_current(ado_state, bath_tag):
         return exp.Q if exp.type == exp.types["+"] else exp.Q.dag()
 
     return -1.0j * sum(
-        exp_sign(exp) * (exp_op(exp) * aux).tr()
-        for aux, exp in level_1_aux
+        exp_sign(exp) * (exp_op(exp) * aux).tr() for aux, exp in level_1_aux
     )
 ```
 
@@ -513,16 +508,15 @@ def current_pade_for_theta(H, bath_L, bath_R, theta, Nk):
     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",
-    )
+    envL = LorentzianEnvironment(bath_L.T, bath_L.mu, bath_L.gamma, bath_L.W)
+    bathL = envL.approx_by_pade(Nk=Nk)
+    envR = LorentzianEnvironment(bath_R.T, bath_R.mu, bath_R.gamma, bath_R.W)
 
-    solver_pade = HEOMSolver(H, [bathL, bathR], max_depth=2, options=options)
+    bathR = envR.approx_by_pade(Nk=Nk, tag="R")
+
+    solver_pade = HEOMSolver(
+        H, [(bathL, bath_L.Q), (bathR, bath_R.Q)], 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,15 +525,13 @@ 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)
-    for theta in thetas
+    current_analytic_for_theta(e1, bath_L, bath_R, 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, bath_L, bath_R, theta, Nk=6)
-    for theta in thetas
+    current_pade_for_theta(H, bath_L, bath_R, theta, Nk=6) for theta in thetas
 ]
 ```
 
@@ -573,7 +565,7 @@ ax.legend(fontsize=25);
 ## About
 
 ```{code-cell} ipython3
-qutip.about()
+about()
 ```
 
 ## Testing
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 7e139d8d..372a78b5 100644
--- a/tutorials-v5/heom/heom-5b-fermions-discrete-boson-model.md
+++ b/tutorials-v5/heom/heom-5b-fermions-discrete-boson-model.md
@@ -1,15 +1,14 @@
 ---
 jupytext:
-  formats: ipynb,md:myst
   text_representation:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.17.0
 kernelspec:
+  name: python3
   display_name: Python 3 (ipykernel)
   language: python
-  name: python3
 ---
 
 # HEOM 5b: Discrete boson coupled to an impurity and fermionic leads
@@ -100,22 +99,15 @@ import contextlib
 import dataclasses
 import time
 
-import matplotlib.pyplot as plt
 import numpy as np
+import matplotlib.pyplot as plt
 
-import qutip
-from qutip import (
-    destroy,
-    qeye,
-    tensor,
-)
-from qutip.solver.heom import (
-    HEOMSolver,
-    LorentzianPadeBath,
-)
+from qutip import about, destroy, qeye, tensor
+from qutip.core.environment import LorentzianEnvironment
+from qutip.solver.heom import HEOMSolver
 
-from ipywidgets import IntProgress
 from IPython.display import display
+from ipywidgets import IntProgress
 
 %matplotlib inline
 ```
@@ -153,8 +145,7 @@ def state_current(ado_state, bath_tag):
         return exp.Q if exp.type == exp.types["+"] else exp.Q.dag()
 
     return -1.0j * sum(
-        exp_sign(exp) * (exp_op(exp) * aux).tr()
-        for aux, exp in level_1_aux
+        exp_sign(exp) * (exp_op(exp) * aux).tr() for aux, exp in level_1_aux
     )
 ```
 
@@ -227,7 +218,7 @@ class LorentzianBathParameters:
         if self.lead == "L":
             self.mu = self.theta / 2.0
         else:
-            self.mu = - self.theta / 2.0
+            self.mu = -self.theta / 2.0
 
     def J(self, w):
         """ Spectral density. """
@@ -239,7 +230,7 @@ class LorentzianBathParameters:
         return fF(x)
 
     def lamshift(self, w):
-        """ Return the lamshift. """
+        """ Return the lamb shift. """
         return 0.5 * (w - self.mu) * self.J(w) / self.W
 
     def replace(self, **kw):
@@ -318,17 +309,17 @@ def steady_state_pade_for_theta(sys_p, bath_L, bath_R, theta, Nk, Nc, Nbos):
     bath_L = bath_L.replace(theta=theta)
     bath_R = bath_R.replace(theta=theta)
 
-    bathL = LorentzianPadeBath(
-        sys_p.Q, bath_L.gamma, bath_L.W, bath_L.mu, bath_L.T,
-        Nk, tag="L",
-    )
-    bathR = LorentzianPadeBath(
-        sys_p.Q, bath_R.gamma, bath_R.W, bath_R.mu, bath_R.T,
-        Nk, tag="R",
-    )
+    envL = LorentzianEnvironment(bath_L.T, bath_L.mu, bath_L.gamma, bath_L.W)
+    envR = LorentzianEnvironment(bath_R.T, bath_R.mu, bath_R.gamma, bath_R.W)
+
+    bathL = envL.approx_by_matsubara(Nk, tag="L")
+    bathR = envR.approx_by_matsubara(Nk, tag="R")
 
     solver_pade = HEOMSolver(
-        sys_p.H, [bathL, bathR], max_depth=2, options=options,
+        sys_p.H,
+        [(bathL, sys_p.Q), (bathR, sys_p.Q)],
+        max_depth=2,
+        options=options,
     )
     rho_ss_pade, ado_ss_pade = solver_pade.steady_state()
     current = state_current(ado_ss_pade, bath_tag="R")
@@ -383,13 +374,5 @@ ax.legend(loc=4);
 ## About
 
 ```{code-cell} ipython3
-qutip.about()
-```
-
-## Testing
-
-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
-assert 1 == 1
+about()
 ```