Environments of Open Quantum Systems

written by Paul Menczel and Gerardo Suarez

QuTiP can describe environments of open quantum systems. They can be passed to various solvers, where their influence is taken into account exactly or approximately. In the following, we will discuss bosonic and fermionic thermal environments. In our definitions, we follow [BoFiN23].

Note that currently, we only support a single coupling term per environment. If a more generalized coupling would be useful to you, please let us know on GitHub.

Bosonic Environments

A bosonic environment is described by a continuum of harmonic oscillators. The open quantum system and its environment are thus described by the Hamiltonian

\[H = H_{\rm s} + Q \otimes X + \sum_k \hbar\omega_k\, a_k^\dagger a_k , \qquad X = \sum_k g_k (a_k + a_k^\dagger)\]

(in the case of a single bosonic environment). Here, \(H_{\rm s}\) is the Hamiltonian of the open system and \(Q\) a system coupling operator. The sums enumerate the environment modes with frequencies \(\omega_k\) and couplings \(g_k\). The couplings are described by the spectral density

\[J(\omega) = \pi \sum_k g_k^2 \delta(\omega - \omega_k) .\]

Equivalently, a bosonic environment can be described by its auto-correlation function

\[C(t) = \langle X(t) X \rangle = \sum_{k} g_{k}^{2} \Big( \cos(\omega_{k} t) \underbrace{( 2 n_{k}+1)}_{\coth(\frac{\beta \omega_{k}}{2})} - i \sin(\omega_{k} t) \Big)\]

(where \(X(t)\) is the interaction picture operator) or its power spectrum

\[S(\omega) = \int_{-\infty}^\infty \mathrm dt\, C(t)\, \mathrm e^{\mathrm i\omega t} ,\]

which is the inverse Fourier transform of the correlation function. The correlation function satisfies the symmetry relation \(C(-t) = C(t)^\ast\).

Assuming that the initial state is thermal with inverse temperature \(\beta\), in the continuum limit the correlation function and the power spectrum can be calculated from the spectral density via

(1)\[\begin{split}\begin{aligned} C(t) &= \int_0^\infty \frac{\mathrm d\omega}{\pi}\, J(\omega) \Bigl[ \coth\Bigl( \frac{\beta\omega}{2} \Bigr) \cos\bigl( \omega t \bigr) - \mathrm i \sin\bigl( \omega t \bigr) \Bigr] , \\ S(\omega) &= \operatorname{sign}(\omega)\, J(|\omega|) \Bigl[ \coth\Bigl( \frac{\beta\omega}{2} \Bigr) + 1 \Bigr] . \end{aligned}\end{split}\]

Here, \(\operatorname{sign}(\omega) = \pm 1\) depending on the sign of \(\omega\). At zero temperature, these equations become \(C(t) = \int_0^\infty \frac{\mathrm d\omega}{\pi} J(\omega) \mathrm e^{-\mathrm i\omega t}\) and \(S(\omega) = 2 \Theta(\omega) J(|\omega|)\), where \(\Theta\) is the Heaviside function.

If the environment is coupled weakly to the open system, the environment induces quantum jumps with transition rates \(\gamma \propto S(-\Delta\omega)\), where \(\Delta\omega\) is the energy change in the system corresponding to the quantum jump. In the strong coupling case, QuTiP provides exact integration methods based on multi-exponential decompositions of the correlation function, see below.

Note

We generally assume that the frequencies \(\omega_k\) are all positive and hence \(J(\omega) = 0\) for \(\omega \leq 0\). To handle a spectral density \(J(\omega)\) with support at negative frequencies, one can use an effective spectral density \(J'(\omega) = J(\omega) - J(-\omega)\). This process might result in the desired correlation function, because

\[\int_0^\infty \frac{\mathrm d\omega}{\pi}\, J'(\omega) \bigl[ \cdots \bigr] = \int_{-\infty}^\infty \frac{\mathrm d\omega}{\pi}\, J(\omega) \bigl[ \cdots \bigr] ,\]

where \([\cdots]\) stands for the square brackets in Eq. (1). Note however that the derivation of Eq. (1) is not valid in this situation, since the environment does not have a thermal state.

Note

Note that the expressions provided above for \(S(\omega)\) are ill-defined at \(\omega=0\). For zero temperature, we simply have \(S(0) = 0\).

For non-zero temperature, one obtains

\[S(0) = 2\beta^{-1}\, J'(0) .\]

Hence, \(S(0)\) diverges if the spectral density is sub-ohmic.

Pre-defined Environments

Ohmic Environment

Ohmic environments can be constructed in QuTiP using the class OhmicEnvironment. They are characterized by spectral densities of the form

(2)\[J(\omega) = \alpha \frac{\omega^s}{\omega_c^{s-1}} e^{-\omega / \omega_c} ,\]

where \(\alpha\) is a dimensionless parameter that indicates the coupling strength, \(\omega_{c}\) is the cutoff frequency, and \(s\) is a parameter that determines the low-frequency behaviour. Ohmic environments are usually classified according to this parameter as

  • Sub-Ohmic (\(s<1\))

  • Ohmic (\(s=1\))

  • Super-Ohmic (\(s>1\)).

Note

In the literature, the Ohmic spectral density can often be found as \(J(\omega) = \alpha \frac{\omega^s}{\omega_c^{s-1}} f(\omega)\), where \(f(\omega)\) with \(\lim\limits_{\omega \to \infty} f(\omega) = 0\) is known as the cutoff function. The cutoff function ensures that the spectral density and its integrals (for example (1)) do not diverge. Sometimes, with sub-Ohmic spectral densities, an infrared cutoff is used as well so that \(\lim\limits_{\omega \to 0} J(\omega) = 0\). This pre-defined Ohmic environment class is restricted to an exponential cutoff function, which is one of the most commonly used in the literature. Other cutoff functions can be used in QuTiP with user-defined environments as explained below.

Substituting the Ohmic spectral density (2) into (1), the correlation function can be computed analytically:

\[C(t)= \frac{\alpha}{\pi} w_{c}^{1-s} \beta^{-(s+1)} \Gamma(s+1) \left[ \zeta\left(s+1,\frac{1+\beta w_{c} -i w_{c} t}{\beta w_{c}} \right) +\zeta\left(s+1,\frac{1+ i w_{c} t}{\beta w_{c}}\right) \right] ,\]

where \(\beta\) is the inverse temperature, \(\Gamma\) the Gamma function, and \(\zeta\) the Hurwitz zeta function. The zero temperature case can be obtained by taking the limit \(\beta \to \infty\), which results in

\[C(t) = \frac{\alpha}{\pi} \omega_c^2\, \Gamma(s+1) (1+ i \omega_{c} t)^{-(s+1)} .\]

The evaluation of the zeta function for complex arguments requires mpmath, so certain features of the Ohmic enviroment are only available if mpmath is installed.

Multi-exponential approximations to Ohmic environments can be obtained through the fitting procedures approx_by_cf_fit and approx_by_sd_fit. The following example shows how to create a sub-Ohmic environment, and how to use approx_by_cf_fit to fit the real and imaginary parts of the correlation function with two exponential terms each.

import numpy as np
import qutip as qt
import matplotlib.pyplot as plt

# Define a sub-Ohmic environment with the given temperature, coupling strength and cutoff
env = qt.OhmicEnvironment(T=0.1, alpha=1, wc=3, s=0.7)

# Fit the correlation function with three exponential terms
tlist = np.linspace(0, 3, 250)
approx_env, info = env.approx_by_cf_fit(tlist, target_rsme=None, Nr_max=3, Ni_max=3, maxfev=1e8)

The environment approx_env created here could be used, for example, with the HEOM solver. The variable info contains info about the convergence of the fit; here, we will just plot the fit together with the analytical correlation function. Note that a larger number of exponential terms would have yielded a better result.

plt.plot(tlist, np.real(env.correlation_function(tlist)), label='Real part (analytic)')
plt.plot(tlist, np.real(approx_env.correlation_function(tlist)), '--', label='Real part (fit)')

plt.plot(tlist, np.imag(env.correlation_function(tlist)), label='Imag part (analytic)')
plt.plot(tlist, np.imag(approx_env.correlation_function(tlist)), '--', label='Imag part (fit)')

plt.xlabel('Time')
plt.ylabel('Correlation function')
plt.tight_layout()
plt.legend()
../_images/guide-environments-2.png

Drude-Lorentz Environment

Drude-Lorentz environments, also known as overdamped environments, can be constructed in QuTiP using the class DrudeLorentzEnvironment. They are characterized by spectral densities of the form

\[J(\omega) = \frac{2 \lambda \gamma \omega}{\gamma^{2}+\omega^{2}} ,\]

where \(\lambda\) is a coupling strength (with the dimension of energy) and \(\gamma\) the cutoff frequency.

To compute the corresponding correlation function, one can apply the Matsubara expansion:

\[C(t) = \sum_{k=0}^{\infty} c_k e^{- \nu_k t}\]

The coefficients of this expansion are

\[\begin{split}\nu_{k} = \begin{cases} \gamma & k = 0\\ {2 \pi k} / {\beta} & k \geq 1\\ \end{cases} \;, \qquad c_k = \begin{cases} \lambda \gamma [\cot(\beta \gamma / 2) - i] & k = 0\\ \frac{4 \lambda \gamma \nu_k }{ (\nu_k^2 - \gamma^2)\beta} & k \geq 1\\ \end{cases} \;.\end{split}\]

The function approx_by_matsubara creates a multi-exponential approximation to the Drude-Lorentz environment by truncating this series at a finite index \(N_k\). This approximation can then be used with the HEOM solver, for example. The HEOM section of this guide contains further examples using the Drude-Lorentz enviroment.

Similarly, the function approx_by_pade can be used to apply and truncate the numerically more efficient Padé expansion.

Underdamped Environment

Underdamped environments can be constructed in QuTiP using the class UnderDampedEnvironment. They are characterized by spectral densities of the form

\[J(\omega) = \frac{\lambda^{2} \Gamma \omega}{(\omega_0^{2}- \omega^{2})^{2}+ \Gamma^{2} \omega^{2}} ,\]

where \(\lambda\), \(\Gamma\) and \(\omega_0\) are the coupling strength (with dimension \((\text{energy})^{3/2}\)), the cutoff frequency and the resonance frequency.

Similar to the Drude-Lorentz environment, the correlation function can be approximated by a Matsubara expansion. This functionality is available with the approx_by_matsubara function.

For small temperatures, the Matsubara expansion converges slowly. It is recommended to instead use a fitting procedure for the Matsubara contribution as described in [Lambert19].

User-Defined Environments

As stated in the introduction, a bosonic environment is fully characterized by its temperature and spectral density (SD), or alternatively by its correlation function (CF) or its power spectrum (PS). QuTiP allows for the creation of an user-defined environment by specifying either the spectral density, the correlation function, or the power spectrum.

QuTiP then computes the other two functions based on the provided one. To do so, it converts between the SD and the PS using the formula \(S(\omega) = \operatorname{sign}(\omega)\, J(|\omega|) \bigl[ \coth( \beta\omega / 2 ) + 1 \bigr]\) introduced earlier, and between the PS and the CF using the fast Fourier transform. The former conversion requires the bath temperature to be specified; the latter requires a cutoff frequency (or cutoff time) to be provided together with the specified function (SD, CF or PS). In this way, all characteristic functions can be computed from the specified one.

The following example manually creates an environment with an underdamped spectral density. It then compares the correlation function obtained via fast Fourier transformation with the Matsubara expansion. The slow convergence of the Matsubara expansion is visible around \(t=0\).

# Define underdamped environment parameters
T = 0.1
lam = 1
gamma = 2
w0 = 5

# User-defined environment based on SD
def underdamped_sd(w):
    return lam**2 * gamma * w / ((w**2 - w0**2)**2 + (gamma*w)**2)
env = qt.BosonicEnvironment.from_spectral_density(underdamped_sd, wMax=50, T=T)

tlist = np.linspace(-2, 2, 250)
plt.plot(tlist, np.real(env.correlation_function(tlist)), label='FFT')

# Pre-defined environment and Matsubara approximations
env2 = qt.UnderDampedEnvironment(T, lam, gamma, w0)
for Nk in range(0, 11, 2):
    approx_env = env2.approx_by_matsubara(Nk)
    plt.plot(tlist, np.real(approx_env.correlation_function(tlist)), label=f'Nk={Nk}')

plt.xlabel('Time')
plt.ylabel('Correlation function (real part)')
plt.tight_layout()
plt.legend()
../_images/guide-environments-3.png

Multi-Exponential Approximations

Many approaches to simulating the dynamics of an open quantum system strongly coupled to an environment assume that the environment correlation function can be approximated by a multi-exponential expansion like

\[C(t) = C_R(t) + \mathrm i C_I(t) , \qquad C_{R,I}(t) = \sum_{k=1}^{N_{R,I}} c^{R,I}_k \exp[-\nu^{R,I}_k t]\]

with small numbers \(N_{R,I}\) of exponents. Note that \(C_R(t)\) and \(C_I(t)\) are the real and imaginary parts of the correlation function, but the coefficients \(c^{R,I}_k\) and exponents \(\nu^{R,I}_k\) are not required to be real in general.

In the previous sections, various methods of obtaining multi-exponential approximations were introduced. The output of these approximation functions are ExponentialBosonicEnvironment objects. An ExponentialBosonicEnvironment is basically a collection of CFExponent s, which store (in the bosonic case) the coefficient, the exponent, and whether the exponent contributes to the real part, the imaginary part, or both. As we have already seen above, one can then compute the spectral density, correlation function and power spectrum corresponding to the exponents, in order to compare them to the original, exact environment.

Let \(c_k \mathrm e^{-\nu_k t}\) be a term in the correlation function (i.e., \(c_k = c^R_k\) or \(c_k = \mathrm i c^I_k\)). The corresponding term in the power spectrum is

\[S_k(\omega) = 2\Re\Bigr[ \frac{c_k}{\nu_k - \mathrm i\omega} \Bigr]\]

and, if a temperature has been specified, the corresponding term in the spectral density can be computed as described above.

The following example shows how to manually create an ExponentialBosonicEnvironment for the simple example \(C(t) = c \mathrm e^{-\nu t}\) with real \(c\), \(\nu\). The power spectrum then is a Lorentzian, \(S(\omega) = 2c\nu / (\nu^2 + \omega^2)\).

c = 1
nu = 2
wlist = np.linspace(-3, 3, 250)

env = qt.ExponentialBosonicEnvironment([c], [nu], [], [])

plt.figure(figsize=(4, 3))
plt.plot(wlist, env.power_spectrum(wlist))
plt.plot(wlist, 2 * c * nu / (nu**2 + wlist**2), '--')
plt.xlabel('Frequency')
plt.ylabel('Power spectrum')
plt.tight_layout()
../_images/guide-environments-4.png

Fermionic environments

The implementation of fermionic environments in QuTiP is not yet as advanced as the bosonic environments. Currently, user-defined fermionic environments and fitting are not implemented.

However, the overall structure of fermionic environments in QuTiP is analogous to the bosonic environments. There is one pre-defined fermionic environment, the Lorentzian environment, and multi-exponential fermionic environments. Lorentzian environments can be approximated by multi-exponential Matsubara or Padé expansions.

In the fermionic case, we consider the number-conserving Hamiltonian

\[H = H_s + (B^\dagger c + c^\dagger B) + \sum_k \hbar\omega_k\, b^\dagger_k b_k , \qquad B = \sum_k f_k b_k ,\]

where the bath operators \(b_k\) and the system operator \(c\) obey fermionic anti-commutation relations. In analogy to the bosonic case, we define the spectral density

\[J(\omega) = 2\pi \sum_k f_k^2\, \delta[\omega - \omega_k] ,\]

which may however now be defined for all (including negative) frequencies, since the spectrum of each mode is bounded.

The fermionic environment is characterized either by its spectral density, inverse temperature \(\beta\) and chemical potential \(\mu\), or equivalently by two correlation functions or by two power spectra. The correlation functions are

\[C^\sigma(t) = \langle B^\sigma(t) B^{-\sigma} \rangle = \int_{-\infty}^\infty \frac{\mathrm d\omega}{2\pi}\, J(\omega)\, \mathrm e^{\sigma \mathrm i\omega t}\, f_F(\sigma \beta[\omega - \mu]) ,\]

where \(\sigma = \pm 1\), \(B^+ = B^\dagger\) and \(B^- = B\). Further, \(f_F(x) = (\mathrm e^x + 1)^{-1}\) is the Fermi-Dirac function. Note that we still have \(C^\sigma(-t) = C^\sigma(t)^\ast\). The power spectra are the Fourier transformed correlation functions,

\[S^\sigma(\omega) = \int_{-\infty}^\infty \mathrm dt\, C^\sigma(t)\, \mathrm e^{-\sigma \mathrm i\omega t} = J(\omega) f_F(\sigma\beta[\omega - \mu]) .\]

Since \(f_F(x) + f_F(-x) = 1\), we have \(S^+(\omega) + S^-(\omega) = J(\omega)\).

Note

The relationship between the spectral density and the two power spectra (or the two correlation functions) is not one-to-one. A pair of functions \(S^\pm(\omega)\) is physical if they satisfy the condition

\[S^-(\omega) = \mathrm e^{\beta(\omega - \mu)}\, S^+(\omega) .\]

For the correlation functions, the condition becomes \(C^-(t) = \mathrm e^{-\beta\mu}\, C^+(t - \mathrm i\beta)^\ast\). For flexibility, we do not enforce the power spectra / correlation functions to be physical in this sense.

Lorentzian Environment

Fermionic Lorentzian environments are represented by the class LorentzianEnvironment. They are characterized by spectral densities of the form

\[J(\omega) = \frac{\gamma W^2}{(\omega - \omega_0)^2 + W^2} ,\]

where \(\gamma\) is the coupling strength, \(W\) the spectral width and \(\omega_0\) the resonance frequency. Often, the resonance frequency is taken to be equal to the chemical potential of the environment.

As with the bosonic Drude-Lorentz environments, multi-exponential approximations of the correlation functions,

\[C^\sigma(t) \approx \sum_{k=0}^{N_k} c^\sigma_k e^{- \nu^\sigma_k t} ,\]

can be obtained using the Matsubara or Padé expansions. The functions approx_by_matsubara and approx_by_pade implement these approximations in QuTiP, yielding approximated environments that can be used, for example, with the HEOM solver. Note that for this type of environment, the Matsubara expansion is very inefficient, converging much more slowly than the Padé expansion. Typically, at least \(N_k \geq 20\) is required for good convergence.

For reference, we tabulate the values of the coefficients and exponents in the following. For the Matsubara expansion, they are

\[\begin{split}\nu^\sigma_{k} = \begin{cases} W - \mathrm i \sigma\, \omega_0 & k = 0\\ \frac{(2k - 1) \pi}{\beta} - \mathrm i \sigma\, \mu & k \geq 1\\ \end{cases} \;, \qquad c^\sigma_k = \begin{cases} \frac{\gamma W}{2} f_F[\sigma\beta(\omega_0 - \mu) + \mathrm i\, \beta W] & k = 0\\ \frac{\mathrm i \gamma W^2}{\beta} \frac{1}{[\mathrm i \sigma (\omega_0 - \mu) + (2k - 1) \pi / \beta]^2 - W^2} & k \geq 1\\ \end{cases} \;.\end{split}\]

The Padé decomposition approximates the Fermi distribution as:

\[f_F(x) \approx f_F^{\mathrm{approx}}(x) = \frac{1}{2} - \sum_{k=0}^{N_k} \frac{2\kappa_k x}{x^2 + \epsilon_k^2}\]

where \(\kappa_k\) and \(\epsilon_k\) are coefficients that depend on \(N_k\) and are defined in J. Chem Phys 133, “Efficient on the fly calculation of time correlation functions in computer simulations”. This approach results in the exponents

\[\begin{split}\nu^\sigma_{k} = \begin{cases} W - \mathrm i \sigma\, \omega_0 & k = 0\\ \frac{\epsilon_k}{\beta} - \mathrm i \sigma\, \mu & k \geq 1\\ \end{cases} \;, \qquad c^\sigma_k = \begin{cases} \frac{\gamma W}{2} f_F^{\mathrm{approx}}[\sigma\beta(\omega_0 - \mu) + \mathrm i\, \beta W] & k = 0\\ \frac{\mathrm i\, \kappa_k \gamma W^2}{\beta} \frac{1}{[\mathrm i\sigma(\omega_0 - \mu) + \epsilon_k / \beta]^2 - W^2} & k \geq 1\\ \end{cases} \;.\end{split}\]

Multi-Exponential Fermionic Environment

Analogous to the ExponentialBosonicEnvironment in the bosonic case, the ExponentialFermionicEnvironment describes fermionic environments where the correlation functions are given by multi-exponential decompositions,

\[C^\sigma(t) \approx \sum_{k=0}^{N_k^\sigma} c^\sigma_k e^{-\nu^\sigma_k t} .\]

Like in the bosonic case, the class allows us to automatically compute the spectral density and power spectra that correspond to the multi-exponential correlation functions. In this case, they are

\[S^\sigma(\omega) = \sum_{k=0}^{N_k^\sigma} 2\Re\Bigr[ \frac{c_k^\sigma}{\nu_k^\sigma + \mathrm i \sigma\, \omega} \Bigr]\]

and \(J(\omega) = S^+(\omega) + S^-(\omega)\).