Accurate Lindblad-form master equation for weakly damped quantum systems across all regimes

Abstract

Realistic models of quantum systems must include dissipative interactions with a thermal environment. For weakly-damped systems, while the Lindblad-form Markovian master equation is invaluable for this task, it applies only when the frequencies of any subset of the system’s transitions are degenerate, or their differences are much greater than the transitions’ linewidths. Outside of these regimes the only available efficient description has been the Bloch–Redfield master equation, the efficacy of which has long been controversial due to its failure to guarantee the positivity of the density matrix. The ability to efficiently simulate weakly-damped systems across all regimes is becoming increasingly important, especially in quantum technologies. Here we solve this long-standing problem by deriving a Lindblad-form master equation for weakly-damped systems that is accurate for all regimes. We further show that when this master equation breaks down, so do all time-independent Markovian equations, including the B-R equation. We thus obtain a replacement for the B-R equation for thermal damping that is no less accurate, simpler in structure, completely positive, allows simulation by efficient quantum trajectory methods, and unifies the previous Lindblad master equations. We also show via exact simulations that the new master equation can describe systems in which slowly-varying transition frequencies cross each other during the evolution. System identification tools, developed in systems engineering, play an important role in our analysis. We expect these tools to prove useful in other areas of physics involving complex systems.

Introduction

Weakly damped open systems are important across a wide range of areas in both physics and chemistry, from quantum thermodynamics^1,2,3 to the control of chemical reactions⁴, to quantum technologies^{5,6,7,8,9,10,11,12}. So long as the thermal environment that induces the weak damping has a high cut-off frequency (something we assume throughout), Lindblad-form Markovian master equations are tremendously useful for modeling these systems as they avoid the computationally expensive, and often prohibitive, task of simulating the thermal environment^13,14. However, recent developments have made it clear that the regime that existing Lindblad master equations cannot describe—the "near degenerate” regime in which non-degenerate transition frequencies are close together—while long ignored, is crucial for investigating important questions in a range of topics, including reservoir engineering and cascaded systems^{15,16,17,18,19,20}, adiabatic computation^21,22, super and sub-radiance^23,24,25,26, and "weak lasing”^27,28, with the possibility that this regime will also reveal new tools for controlling quantum systems.

Weakly-damped quantum systems can be divided into three regimes depending on the frequency difference between pairs of transitions. These regimes are degenerate (the frequency difference is zero), non-degenerate (the frequency difference is much greater than the transitions’ linewidths), and near-degenerate (everything else). The degenerate and non-degenerate regimes are described, respectively, by two quite different Lindblad master equations^29,30 (for ease of reference we present these master equations in the Supplementary Discussion). The difference between them is exemplified by the fact that degenerate transitions exhibit super and sub-radiance, whereas non-degenerate transitions do not. These two Lindblad master equations are obtained from the Bloch–Redfield master equation by making the secular (rotating-wave) approximation. However, no Lindblad master equation has been obtained for the near-degenerate regime^31,32. Thus to simulate systems in which two or more distinct transitions are separated by less than a few linewidths, one must resort to the Bloch–Redfield (B–R) master equation^33,34. This equation has long been the subject of debate because it is not guaranteed to preserve the positivity of the density matrix^35,36. In some subfields (e.g., photo-chemistry^37,38), the B–R master equation is used as the standard vehicle for treating weakly-damped systems. Practitioners in other fields, for example quantum optics and many areas of quantum technologies, do not use it because its failure to ensure such a fundamental property as positivity is seen as an indication that it cannot be trusted.

There have been a number of papers, some quite recent, arguing that the B–R equation is a valid and effective model so long as the system is close to Markovian^28,39,40. In particular, Eastham et al.²⁸ considered a model of two coupled linear oscillators that can be solved exactly, and examined how well the B–R equation describes the near-degenerate regime (since this is the regime in which it is needed). They found that the B–R equation was both very accurate and preserved positivity to a very good approximation. They attributed this to the fact that the dynamics of the coupled oscillators stays close to Markovian, which was in turn due to the relatively slow variation of the spectral density. Jeske et al.⁴⁰ also noted that when transitions are close enough that they share the same value of the spectral density, the B–R equation reduces to the degenerate master equation, which is a key element in our analysis here. These recent works raise an interesting question: other authors have assumed that the near-degenerate regime is non-Markovian, due to the apparent lack of a Lindblad-form master equation in that regime^31,32. If the dynamics of weakly-damped systems is indeed Markovian in the near-degenerate regime, then it is not unreasonable to suggest that there may be a completely positive Markovian master equation that accurately describes it.

Here we show that there is a single, Lindblad-form master equation that describes weakly-damped systems across all three regimes. This master equation, which can be found in Eqs. (37) and (50), applies to the Ohmic spectrum, and any spectrum that varies sufficiently slowly on the scale of the Lamb shifts and linewidths. For baths with spectra that satisfy this "slow variation" condition, our master equation agrees to very high accuracy with the B–R equation, something that follows from our derivation and is confirmed by numerical simulations. We further show, using exact simulations, that when our "slow variation" condition is broken not only does the B–R equation for thermal damping break down, but so do all time-independent Markovian master equations. We thus show that this B–R equation cannot be trusted outside the regime in which our master equation is valid, and in this sense our master equation is a complete replacement for it. The Lindblad master equation has a simpler form that the Bloch–Redfield equation, and thus provides new insight into the behavior of the near-degenerate regime.

Numerical simulations reveal that our master equation describes non-degenerate transitions more accurately than the existing Lindblad master equation for non-degenerate transitions. Not only does our master equation finally provide a non-controversial method for simulating all weakly damped systems; being in the Lindblad form it can also be simulated using efficient Monte Carlo methods^41,42,43,44, and provides a formulation of the action of a thermal bath as a continuous measurement on the system. This formulation quantifies the way in which information flows from the system to the bath.

We expect that many important problems involving the near-degenerate regime will also involve transition frequencies that change with time, and possibly cross during the evolution. Examples of this are the Landau-Zener transition^45,46 and the control of super- and sub-radiance by shifting energy levels. We show that the adiabatic extension of our master equation is able to accurately describe such time-dependent problems, so long as the rate of change of the transition frequencies is not too fast.

Our derivation of the Lindblad master equation provides the following corollaries:

• 1. It shows that the secular approximation is unnecessary: weak damping, a high bath cut-off frequency, and sufficient flatness of the spectral density suffice to guarantee positivity and Markovianity. We show that the flatness of the spectral density is also a necessary condition for Markovianity for the oscillator bath model of thermal damping.

• 2. It resolves the controversy regarding the Bloch–Redfield master equation for thermal damping: it shows that when the spectral density is sufficiently flat, this equation is very close to a Lindblad master equation, and will thus approximately preserve positivity. Conversely, outside this flatness condition this Bloch–Redfield equation is, in general, no longer valid, confirming the conjecture in ref. ²⁸).

We obtain the new master equation in two steps. First, we use exact simulations of a V system coupled to an Ohmic bath, together with the method of system identification, developed in systems engineering, to show that weakly-damped quantum systems with an Ohmic spectrum are not only Markovian but also time-independent across all three regimes. This method also allows us to directly back-out the Lindblad-form equation of motion for this V system. Second, aided by the form obtained in step one, we show how to derive the new Lindblad master equation from the Bloch–Redfield equation valid for all regimes and all temperatures. We provide additional confirmation of its accuracy for the Ohmic bath by comparing its predictions to those of exact simulations for two further systems, a trident system and two co-located qubits.

We note that here we are specifically concerned with the standard oscillator-bath model of thermal damping. Thus the B–R equation we are concerned with is the one that results from coupling the system to quadrature operators of the bath oscillators. Naturally, other kinds of system-bath coupling may lead to different B–R equations. It is an open question as to whether our method for obtaining a Lindblad equation that replaces the B–R equation extends to models in which the system couples to other kinds of bath operators.

The results section is laid out as follows. First we perform simulations of a V system with an Ohmic bath, and obtain an accurate master equation for this system for all regimes using system identification. Next we use the information obtained from the V system to derive a master equation for all weakly damped systems and all temperatures given a constraint on the derivative of the spectral density. We then perform numerical simulations to provide further confirmation of the accuracy of the master equation, and show how both it and the B–R equation break down as the derivative of the spectral density is increased. These simulations also show that weakly damped system become non-Markovian as these equations break down. Finally, we use numerical simulations to show that the master equation is able to describe time-dependent systems whose energy levels cross during the evolution.

Results

System identification in the near-degenerate regime

The methods of system identification provide us with a way to determine, from the time series of a linear time-invariant system, the minimal number of variables required to generate this time-series (that is, the dimension of the system), as well as its equations of motion. System identification (SID) methods are typically concerned with input/output systems. SID involves obtaining the outputs of a system for a large enough set of distinct inputs that the equations of motion can be determined. While our system does not have inputs, SID methods are easily adapted to replace the set of inputs with a set of initial states. (The SID method that we use is given in Methods.) Since the evolution of the V system is non-trivial only when the upper levels are populated, and the evolution does not generate coherence with the lower level, the two upper populations together with the real and imaginary parts of their coherence form a closed four-dimensional system. Performing exact simulations of the V system, SID provides us with the dynamics of a fictitious (and possibly larger) system that generates the four-dimensional dynamics. Specifically, if we denote the state of the fictitious system at time t by v(t), then SID provides us with a matrix M(τ) where v(τ) = M(τ)v(0) for a specified time τ. The number of appreciable eigenvalues of M is the effective size of the fictitious system.

Since it is only the ratios between the various rate parameters that determine the dynamical behavior (up to a scaling of time) we specify all frequencies in terms of an arbitrary frequency, $\tilde{\nu }$. We perform SID on the V system depicted in Fig. 1a with bath cut-off frequency $\Omega =80\pi \tilde{\nu }$ (we give the details of the bath model when we derive the master equation below), fix the mean transition frequency $\bar{\omega }\equiv ({\omega }_{1}+{\omega }_{2})/2=3\pi \tilde{\nu }$, and choose the coupling constants g₁ and g₂ (defined in Eq. (10)) so as to give the decay rates ${\gamma }_{1}=0.1\tilde{\nu }$ and ${\gamma }_{2}=0.05\tilde{\nu }$. Note that since we are simulating the Hamiltonian in Eq. (8), it is the coupling constants g_j, rather than the damping rates γ_j, that define the simulation, with γ_j ∝ ∣g_j∣². Because we have fixed $\bar{\omega }$ while changing the detuning, and since γ_j ∝ ω_j for the Ohmic spectral density, when the detuning changes we must change g_j to keep the damping rates fixed. We wish to examine the evolution when the detuning, Δω ≡ ω₂ − ω₁, is not large compared to the damping rates, so we simulate the evolution for the following four values of Δω: 0, 0.28πγ₁, 2πγ₁, and 4.8πγ₁.

Fig. 1: The four open systems we use as examples.

Here we depict four systems with transitions that decay due to a coupling with a thermal bath at zero temperature. The red bars are the energy eigenstates of the system and the blue wiggly lines indicate the transitions. These system are a V system, b trident system, c four-level system, and d two co-located qubits. Given that the relative energy of each level in the diagram is indicated by its vertical position, Δ denotes the detuning between the transitions in a, c, and d. In system (b) there are three transitions and thus two independent detunings denoted by Δ₁ and Δ₂. The decay rate of the j^th transition is denoted by γ_j. The transition operators and frequencies for each of these systems are given in the supplementary discussion.

To perform the exact simulations we use the method detailed in^13,14 which employs the matrix-product-state method of Vidal^47,48. This in turn requires a split operator method, for which we use a second-order method valid for time-dependent systems, and choose a time-step small enough to obtain an accuracy of about six digits of precision.

Obtaining the matrix M for each value of the detuning, Δω, we find that the largest four eigenvalues of M account for almost all of the dynamical behavior for all four values: the magnitudes of all the rest of the eigenvalues contribute a fraction of less than 3 × 10⁻⁴ to the 1-norm of M. This result implies that the dynamics of the system in the near-degenerate regime is both time-independent and Markovian to very good approximation.

A 4-dimensional dynamical model for the four independent variables of the V system can now be obtained merely by taking the log of the matrix M(t) for some appropriate value of t. We note that the time index t must be smaller than the smallest period of the dynamics to avoid multiple branches of the complex logarithm. Writing the four variables as the vector x, the approximate model is $\dot{{\bf{x}}}=D{\bf{x}}$, with $D=\mathrm{ln}\,[M(t)]/t$. To determine the Lindblad-form master equation specified by this model, we need to translate from the elements of D to the familiar terms used to express such master equations. The simplest way to do this is to take a general degenerate master equation for a V system and derive its D matrix. The degenerate master equation for a V system (Fig. 1a), in which both transitions have frequency ω, is given by

$$\dot{\rho }=-\frac{\,\text{i}\,}{\hslash }\left[{H}_{0}+{H}_{{\rm{L}}},\rho \right]+{\mathcal{D}}[\Sigma ]\rho .$$

(1)

Here

$${H}_{0}=\hslash {\omega }_{0}(\left|1\right\rangle \left\langle 1\right|+\left|2\right\rangle \left\langle 2\right|),$$

(2)

$$\Sigma =\sqrt{{\gamma }_{1}}{\sigma }_{1}+{e}^{{\rm{i}}\phi }\sqrt{{\gamma }_{2}}{\sigma }_{2},$$

(3)

with ω₀ the frequency of both transitions, ${\mathcal{D}}$ is a superoperator defined by

$${\mathcal{D}}[c]\rho \equiv \frac{1}{2}\left(2c\rho {c}^{\dagger }-{c}^{\dagger }c\rho -\rho {c}^{\dagger }c\right)$$

(4)

for an arbitrary operator c, and H_L is the Lamb shift Hamiltonian, given by

$${H}_{{\rm{L}}}=-\hslash \left[\sum _{j}{\Delta }_{j}{\sigma }_{j}^{\dagger }{\sigma }_{j}+\sqrt{{\Delta }_{1}{\Delta }_{2}}({e}^{{\rm{i}}\phi }{\sigma }_{1}^{\dagger }{\sigma }_{2}+{\rm{H}}.{\rm{c}}.)\right]$$

(5)

with ${\sigma }_{1}=\left|0\right\rangle \left\langle 1\right|$, ${\sigma }_{2}=\left|0\right\rangle \left\langle 2\right|$. The phase ϕ is determined by the phases of the interactions between the transitions and the bath (see Eq. (10)). We also note that H_L can be factored as H_L = − ħf(ω)D^†D in which $D=\sqrt{{\Delta }_{1}}{\sigma }_{1}+{e}^{{\rm{i}}\phi }\sqrt{{\Delta }_{2}}{\sigma }_{2}$.

From the derivation of the degenerate master equation²⁹ (also see below) we know that the decay rates γ_j depend on the spectral density of the bath evaluated at their corresponding transition frequencies ω_j. The Lamb shifts depend both on the damping rates (to which they are proportional) as well as a factor that is an integral of the entire spectral density. Thus if the frequencies of the transitions are changed while leaving the spectral density the same, the decay rates and Lamb shifts also change.

We find that the backed-out model for Δω ≠ 0, in which the Hamiltonian is now

$${H}_{0}=\hslash {\omega }_{0}\left|1\right\rangle \left\langle 1\right|+\hslash ({\omega }_{0}+\Delta \omega )\left|2\right\rangle \left\langle 2\right|,$$

(6)

has exactly the same form as the degenerate master equation. That is, it can be written as Eq. (1) with ${H}_{{\rm{L}}}={\sum }_{j}{\zeta }_{j}{\sigma }_{j}^{\dagger }{\sigma }_{j}+(\xi {\sigma }_{1}^{\dagger }{\sigma }_{2}+\,\text{H.c.})$ and σ = ∑_jη_jσ_j for some set of {ζ_j, ξ, η_j}. This may be considered a little surprising, given that the non-degenerate master equation has no terms in the Lamb shift Hamiltonian proportional to σ₁σ₂. A simple guess for the parameters ζ_j, ξ, and η_j as functions of the Lamb shifts and the decay rates is to take exactly the expression for the degenerate master equation, but to replace Δ₂(ω₀) and γ₂(ω₀) by new values implied by the new value of ω₂, namely Δ₂(ω₀ + Δω) and γ₂(ω₀ + Δω).

We find that this trial master equation does indeed match the model backed out using SID for all three values of Δω. We compare further the evolution predicted by this master equation to the exact evolution in Fig. 2, for a range of values of Δω. For these simulations we use ${\omega }_{0}=10\pi \tilde{\nu }$ and ${\gamma }_{1}({\omega }_{0})=2{\gamma }_{2}({\omega }_{0})=0.1\tilde{\nu }$. In Fig. 2c we show a measure of the difference between the evolution given by the master equation and exact simulations for a range of values of the detuning. This measure is an average of the absolute values of the differences between the populations and coherences of the density matrix averaged over time. The measure is below 10⁻³ for all values of the detuning shown.

Fig. 2: Accuracy of the master equation inferred from system identification.

Here we compare the exact evolution of the weakly-damped V system with that predicted by the master equation in Eq. (1) with the replacements Δ₂(ω₀) → Δ₂(ω₀ + Δω) and γ₂(ω₀) → γ₂(ω₀ + Δω), which for the V-system is equivalent to the new Lindblad master equation we derive here. In (a) and (b) we show the populations of levels $\left|1\right\rangle$ (red) and $\left|2\right\rangle$ (blue) as a function of time with the initial state $\left|{\psi }_{0}\right\rangle =(\left|1\right\rangle +\left|2\right\rangle )/\sqrt{2}$. The damping rates are ${\gamma }_{1}=2{\gamma }_{2}=0.1\tilde{\nu }$, in which $\tilde{\nu }$ is an arbitrary frequency specifying our frequency units. The solid curves are the exact evolution and the dashed curves are that of the master equation. The detuning is (a) Δω = 4γ₁, (b) Δω = 100γ₁. In (c), for the values of the detuning shown in the legend, and with the parameters above, we plot a measure of the deviation of the master equation from the exact dynamics as a function of time. This measure is an average taken over the deviations of the four relevant elements of the density matrix: the populations of levels $\left|1\right\rangle$ and $\left|2\right\rangle$, and the real and imaginary part of the coherence between them. Notably, this average error remains less than 8 × 10⁻⁴ for all the values of detuning we explored.

In Fig. 2b we plot the evolution of the populations of the upper levels for both the master equation and the exact simulation for a value of Δω that might be considered well into the non-degenerate regime (Δω = 100γ). We see that the evolution contains "wiggles” that are correctly predicted by our trial master equation, but are not predicted by the non-degenerate master equation. These wiggles vanish as Δω → ∞.

In showing that, for the V-system, the degenerate, near-degenerate, and far-detuned regimes are all described by a Lindblad master equation that is essentially the degenerate master equation, SID has provided us with the insight we need to derive this Lindblad equation from the B–R equation. An inspection of the B–R equation, Eq. (25) in the next section, shows that the only way that it can take this Lindblad form is if the frequency-dependent parameters R_j and I_j that appear in this equation are effectively equal for transitions whose frequency differences are on the order of the damping rates and Lamb shifts, respectively. Combining this with the fact that certain terms containing these parameters are effectively removed by the rotating-wave approximation when the same frequency differences are much larger than the damping rates and Lamb shifts allows us to derive the Lindblad master equation. We perform this derivation in the next section.

Derivation of the master equation

Having shown numerically that there is a Markovian master equation describing thermal damping of two arbitrarily detuned transitions, at least for the Ohmic bath, as well as obtaining the form that this equation takes, a close examination of the usual derivation of the existing Markovian master equations reveals how this more general master equation can be derived for an arbitrary number of levels. For simplicity we present this derivation first for a bath at zero temperature. We then outline the derivation for non-zero temperature, since it is essentially the same.

The Hamiltonian for the standard model of thermal damping, in which a system is coupled to a continuum of independent harmonic oscillators, is given by

$$\begin{array}{*{20}{l}}H ={H}_{{\rm{sys}}}+\hslash (A+{A}^{\dagger }+D){\int\nolimits_{0}^{\Omega}} \,\sqrt{J(\omega )}\left[b(\omega )+{b}^{\dagger }(\omega )\right]d\omega \\+\,{\int \nolimits_{0}^{\Omega}}\, \hslash \omega b{(\omega )}^{\dagger }b(\omega )\ d\omega .\end{array}$$

(7)

Here H_sys is the Hamiltonian of the system which has a discrete set of energy levels. The operators b(ω) are the annihilation operators for a continuum of harmonic oscillators indexed by their frequencies ω (equivalently the modes of a quantum field). The function J(ω) is the density of oscillators per unit frequency, usually referred to as the spectral density of the bath. The maximum frequency of the bath oscillators is Ω and is called the "cut-off” frequency. Instead of having a sharp cut-off at frequency Ω one can instead arrange J(ω) to exhibit a smooth drop-off above some frequency. We use a sharp cut-off purely for simplicity. As will be clear in what follows, in the weak-damping regime the addition of a smooth cut-off merely modifies the values of the damping rates and Lamb shifts.

We have written the Hermitian system operator with which the bath couples as A + A^† + D. In the eigenbasis of H_sys, D is defined as being diagonal with real elements and A is upper-triangular. Thus A has all the matrix elements that transform higher energy levels to lower ones. The diagonal component D generates only dephasing for the transitions. In Methods we derive the terms in the master equation that result from a non-zero D, and show that they play no role in near-degenerate transitions. Further, the terms that arise from D are already in the Lindblad form when deriving the B-R equation. Thus we can set D to zero in our analysis here.

If we denote the energy levels of the system by $\left|n\right\rangle$, thus writing

$${H}_{{\rm{sys}}}=\sum _{n}{E}_{n}\left|n\right\rangle \left\langle n\right|,$$

(8)

and define transition operators (or "lowering operators")

$${\sigma }_{j}=\left|{n}_{j}\right\rangle \left\langle {m}_{j}\right|,\ \ {E}_{{m}_{j}}\;>\;{E}_{{n}_{j}},\ \ j=1,\ldots ,J,$$

(9)

then A can be written as

$$A=\mathop{\sum }\limits_{j=1}^{J}{g}_{j}{\sigma }_{j},$$

(10)

where ${g}_{j}=\left|{g}_{j}\right|{e}^{{\rm{i}}{\phi }_{j}}$ are complex numbers giving the magnitude and phase of the coupling of transition j to the bath. The frequency of transition j is

$${\omega }_{j}=\frac{{E}_{{m}_{j}}-{E}_{{n}_{j}}}{\hslash },$$

(11)

and the evolution of the transition operators σ_j and bath operators b(ω) in the interaction picture is

$${\sigma }_{j}^{{\rm{I}}}={\sigma }_{j}{e}^{-{\rm{i}}{\omega }_{j}t},$$

(12)

$${b}^{{\rm{I}}}(\omega )=b(\omega ){e}^{-{\rm{i}}\omega t}.$$

(13)

To derive Markovian master equations one applies a rotating-wave approximation (RWA) to the Hamiltonian above. (Without the RWA one obtains non-Makovian evolution, as exemplified by the Caldeira-Leggett master equation⁴⁹.) The RWA that is made at this point should not be confused with a second rotating-wave approximation which is the final step that turns the Bloch-Redfeild master equation into the degenerate and non-degenerate Lindblad master equations. To apply this first RWA we move into the interaction picture. The terms in the interaction Hamiltonian that contain the products Ab(ω) and A^†b^†(ω) (the so-called "off-resonant” terms) become

$${H}_{{\rm{OR}}}^{{\rm{I}}}=\hslash {\int \nolimits_{0}^{\Omega}} \, \sqrt{J(\omega )}\left[\sum _{j}{g}_{j}{\sigma }_{j}b(\omega ){e}^{-{\rm{i}}({\omega }_{j}+\omega )}+{\rm{H}}.{\rm{c}}\right]d\omega .$$

(14)

Since the minimum frequency at which each of the terms in the sum over j oscillates is ω_j, when the damping rates and Lamb shifts (to be derived below) are much less than all the ω_j, these terms will average to zero on the timescale of the dynamics induced by the bath, and can be discarded. With this approximation the Hamiltonian of the system and bath becomes

$$\begin{array}{rl}{H}_{{\rm{RWA}}}= \!\!\!\! & {H}_{{\rm{sys}}}+\hslash {\int \nolimits_{0}^{\Omega}}\, \sqrt{J(\omega )}\left[{A}^{\dagger }b(\omega )+A{b}^{\dagger }(\omega )\right]d\omega \\ & +\,{\int \nolimits_{0}^{\Omega }}\, \hslash \omega b{(\omega )}^{\dagger }b(\omega )\ d\omega .\end{array}$$

(15)

The regime of weak damping is defined as the regime in which we are close enough to the limit in which $\mathop{\min }\nolimits_{j}({\omega }_{j})/\mathop{\max }\nolimits_{j}({\gamma }_{j})\to \infty$ so that this approximation is a good one.

To proceed one now applies what are known as the Born-Markov approximations to the evolution generated by H_RWA. For the details of these approximations we refer the reader to^29,36,50). The result is the following expression for the evolution of the density matrix of the system in the interaction picture:

$$\frac{d{\rho }^{{\rm{I}}}}{dt}=\frac{1}{{\hslash }^{2}}{\int \nolimits_{0}^{\infty}}\,{\text{Tr}}_{{\rm{B}}}\left[\left[{H}_{{\rm{R}}}^{{\rm{I}}}(t-s),{\rho }^{{\rm{I}}}(t)\otimes {\rho }_{{\rm{B}}}(0)\right],{H}_{{\rm{R}}}^{{\rm{I}}}(t)\right]ds,$$

(16)

where

$${H}_{{\rm{R}}}^{{\rm{I}}}=\hslash {\int \nolimits_{0}^{\Omega}}\, \sqrt{J(\omega )}\left[\sum _{j}{g}_{j}{\sigma }_{j}^{{\rm{I}}\dagger }{b}^{{\rm{I}}}(\omega )+{\rm{H}}.{\rm{c}}\right]d\omega$$

(17)

is the interaction between the system and the bath that appears in H_RWA, in the interaction picture. The operator ρ^I(t) is the density matrix of the system in the interaction picture, ρ_B(0) is the initial density matrix of the bath, and ${{\rm{Tr}}}_{{\rm{B}}}[\cdot ]$ denotes the trace over the bath.

To proceed now we will examine a single term from the expression above, since all the terms are similar and each is processed in the same way. Substituting ${H}_{{\rm{R}}}^{{\rm{I}}}(t)$ and ${H}_{{\rm{R}}}^{{\rm{I}}}(s)$ into the expression above, one of the terms we obtain is

$$K={\int \nolimits_{0}^{\infty}}\left[{\int \nolimits_{0}^{\Omega}}\, G(\omega )\, d\omega \right]{\sigma }_{k}^{{\rm{I}}\dagger }(t){\sigma }_{j}^{{\rm{I}}}(t-s){\rho }^{{\rm{I}}}(t)\ ds,$$

(18)

where

$$\begin{array}{*{20}{l}}G(\omega) ={e}^{-{\rm{i}}\omega s}{\int \nolimits_{0}^{\Omega}}\ \left\langle b(\omega ){b}^{\dagger }(\omega ^{\prime} )\right\rangle \sqrt{J(\omega )J(\omega ^{\prime} )}\ d\omega ^{\prime} \\=\,J(\omega ){e}^{-{\rm{i}}\omega s}.\end{array}$$

(19)

Here we have used the relation $[b(\omega ),{b}^{\dagger }(\omega ^{\prime} )]=\delta (\omega -\omega ^{\prime} )$ and chosen the field to be at zero temperature so that $\left\langle {b}^{\dagger }(\omega ^{\prime} )b(\omega )\right\rangle =0$. Note that all terms in Eq. (16) that contain two lowering or two raising operators are eliminated by the rotating-wave approximation.

Now substituting G(ω) into K and rearranging we obtain

$$\begin{array}{*{20}{l}}K={\int \nolimits_{0}^{\infty}}\left[{\int \nolimits_{0}^{\Omega}}\, J(\omega ){\text{e}}^{-{\rm{i}}\omega s}\, d\omega \right]{\sigma }_{k}^{{\rm{I}}\dagger }(t){\sigma }_{j}^{{\rm{I}}}(t-s){\rho }^{{\rm{I}}}(t) \, ds\\=\, \left[{\int \nolimits_{0}^{\Omega}}\,{\int \nolimits_{0}^{\infty}}\, J(\omega ){\text{e}}^{-{\rm{i}}(\omega -{\omega }_{j})s}\ ds\ d\omega \right]{\sigma }_{k}^{{\rm{I}}\dagger }(t){\sigma }_{j}^{{\rm{I}}}(t){\rho }^{{\rm{I}}}(t).\end{array}$$

(20)

It is useful to define

$${\Gamma }_{j}\equiv \Gamma ({\omega }_{j})\equiv {\int \nolimits_{0}^{\Omega}}\, {\int \nolimits_{0}^{\infty}}\, J(\omega ){\text{e}}^{-{\rm{i}}(\omega -{\omega }_{j})s}\ ds\ d\omega ,$$

(21)

$${R}_{j}=\,\text{Re}\,[{\Gamma }_{j}],$$

(22)

$${I}_{j}=\,\text{Im}\,[{\Gamma }_{j}].$$

(23)

Moving back into the Schrödinger picture, and writing down all the terms, we obtain the Bloch–Redfield equation for thermal damping, which is

$$\begin{array}{*{20}{l}}\dot{\rho }=-\frac{{\rm{i}}}{{\hslash}}[{H}_{{\rm{sys}}},\rho ]-\,{\rm{i}}\,{\mathop{\sum}\limits_{j}^{\,}}{I}_{j}[{{\Sigma}_j^\dagger}{{\Sigma}_j},\rho ]+ {\mathop{\sum}\limits_{j}}2{R}_{j}{\mathcal{D}}[{{\Sigma}_j}]\rho \\-\,{\rm{i}}\,{\mathop{\sum}\limits_{k\ne j}^{\,}}{I}_{k}\left[{{\Sigma}_j^\dagger }{{\Sigma}_k}\rho -\rho {{\Sigma}_k^\dagger}{{\Sigma}_j}+{{\Sigma}_j}\rho {{\Sigma}_k^\dagger }-{{\Sigma}_k}\rho {{\Sigma}_j^\dagger }\right]\\-{\sum \limits _{k\ne j}^{\,}}{R}_{k}\left[{{\Sigma}_j^\dagger }{{\Sigma}_k}\rho +\rho {{\Sigma} _k^\dagger }{{\Sigma}_j}-{{\Sigma}_j}\rho {{\Sigma}_k^\dagger }-{{\Sigma}_k}\rho {{\Sigma}_j^\dagger }\right].\end{array}$$

(24)

Here, for compactness, we have defined

$${\Sigma }_{j}\equiv {g}_{j}{\sigma }_{j},$$

(25)

and ${\mathcal{D}}$ is the superoperator defined in Eq. (4).

We note that the decay rates γ_j will be

$$\begin{array}{*{20}{l}}{\gamma }_{j}= 2| {g}_{j}{| }^{2}{R}_{j}\\= 2| {g}_{j}{| }^{2} {\int \nolimits_{0}^{\Omega}}\ J(\omega )\left[{\int \nolimits_{0}^{\infty}}\, \cos ([\omega -{\omega }_{j}]s)ds\right]d\omega \\= 2| {g}_{j}{| }^{2} {\int \nolimits_{0}^{\Omega}}\, J(\omega )\pi \delta (\omega -{\omega }_{j})d\omega \\= 2\pi | {g}_{j}{| }^{2}J({\omega }_{j})\\\equiv \gamma ({\omega }_{j})\end{array}$$

(26)

and the Lamb shifts will be

$$\begin{array}{*{20}{l}}{\Delta }_{j}= | {g}_{j}{| }^{2}\ {I}_{j}\\= | {g}_{j}{| }^{2}{\int \nolimits_{0}^{\Omega}}\ J(\omega )\left[{\int \nolimits_{0}^{\infty}}\, \sin ([\omega -{\omega }_{j}]s)ds\right]d\omega \\= | {g}_{j}{| }^{2}\ {\mathbb{P}}\left[{\int\nolimits_{0}^{\Omega}}\ J(\omega )\left(\frac{1}{\omega -{\omega }_{j}}\right)d\omega \right]\\= | {g}_{j}{| }^{2}\ {\mathbb{P}}\left[{\int \nolimits_{-{\omega }_{j}}^{\Omega -{\omega }_{j}}}\frac{J(\omega +{\omega }_{j})}{\omega }d\omega \right] \\\equiv \Delta ({\omega }_{j}).\end{array}$$

(27)

Here ${\mathbb{P}}[\cdot ]$ denotes the principle value of an integral. For readers not familiar with this quantity we give the definition and an example in Methods. So long as the spectral density does not decrease with ω, and Ω > ω_j, the Lamb shift Δ_j can be expected to be greater than the damping rate γ_j (this is true for the Ohmic bath, see below).

The master equation we have derived in Eq. (25), the Bloch–Redfield equation for thermal damping, includes arbitrary detuning between the levels, but it is not in the Lindblad form, and does not guarantee that the density matrix will remain positive. We wish to obtain a master equation in the Lindblad form that is still valid for all detunings between the transitions.

Note first that when transitions j and k are degenerate Γ_j = Γ_k. In this case the last two lines of Eq. (25) combine respectively with the last two terms on the first line to give the degenerate master equation for these transitions. In this case the Lamb-shift Hamiltonian is H_L = ${\hbar}$D^†D with $D={\sum }_{j}\sqrt{{I}_{j}}{\Sigma }_{j}$ and the Lindblad damping term is ${\mathcal{D}}[\Sigma ]$ with $\Sigma ={\sum }_{j}\sqrt{2{R}_{j}}{\Sigma }_{j}$.

There is another situation in which the B–R equation will reduce immediately to a Lindblad master equation: if the spectrum is flat, meaning that it is the same for all ω. (The assumption of a flat spectrum is often called the "white noise approximation", as it is useful for deriving quantum Langevin equations for open systems in the non-degenerate regime^51,52.) If the spectrum is flat then Γ_j = Γ_k for all values of ω_j and ω_k, and so in this case the resulting master equation, which has the same form as the degenerate master equation, is valid for all regimes. Unfortunately, physically relevant spectra are not flat.

What we do now is to show that a flat spectrum is not required to derive a master equation valid for all regimes; the rate of change of the spectrum with respect to ω does not need to be zero, it merely needs to be small enough. We determine the necessary condition on this rate of change, and use it to derive the Lindblad master equation. This master equation has the same form as that for a flat spectrum, but of course the damping rates and Lamb shifts will not be the same as those for a flat spectrum, since the spectral density will now in general be different for different transitions.

To begin we observe that when two transition frequencies, ω_j and ω_k, are different, every term in the last two lines of Eq. (25) will oscillate at the difference frequency Δω_jk = ∣ω_j − ω_k∣. If this difference frequency is sufficiently high then these terms will average to zero and we will be left only with the first line of Eq. (25), which is the non-degenerate master equation. How large does Δω_jk need to be to eliminate the last two lines of Eq. (25)? It needs to be much larger than the magnitudes of the rest of the dynamical terms in the master equation (excluding the Hamiltonian of the system, since this does not change the populations of the system’s eigenstates). The magnitudes of the second and third terms on the top line are Δ_j and γ_j, respectively, and those on the last two lines are

$${M}_{jk}\equiv | {g}_{j}{g}_{k}| {R}_{j} \sim \sqrt{{\gamma }_{j}{\gamma }_{k}},$$

(28)

$${O}_{jk}\equiv | {g}_{j}{g}_{k}| {I}_{j} \sim \sqrt{{\Delta }_{j}{\Delta }_{k}}.$$

(29)

So the terms on the last two lines are eliminated when $\mathop{\min }\nolimits_{jk}\Delta {\omega }_{jk}\gg \mathop{\max }\nolimits_{lm}{M}_{lm}$ and $\mathop{\min }\nolimits_{jk}\Delta {\omega }_{jk}\gg \mathop{\max }\nolimits_{lm}{O}_{lm}$. Without loss of generality we will assume that the Lamb shifts are greater than the damping rates (in the opposite case one merely switches the roles of O_j and M_j). To find a set of terms that are in the Lindblad form, and that are an excellent approximation to the last two lines of Eq. (25), we only need concern ourselves with the regime

$$\Delta {\omega }_{jk}\lesssim \mathop{\min }\nolimits_{jk}{O}_{jk}$$

(30)

(since we have assumed O_jk ≥ M_jk, the regime $\Delta {\omega }_{jk}\lesssim \mathop{\min }\nolimits_{jk}{O}_{jk}$ automatically includes the regime $\Delta {\omega }_{jk}\lesssim \mathop{\min }\nolimits_{jk}{M}_{jk}$). Outside of this regime, the last two lines will be eliminated by the oscillations at the detuning frequency Δω_jk.

We now recall that γ_j and Δ_j, and therefore M_jk and O_jk, must be much smaller than both transition frequencies in order for the master equation to be valid. This is a requirement of the initial rotating wave approximation discussed above. Combining this with Eq. (30) we need only consider the regime in which

$$\Delta {\omega }_{jk}\lesssim {O}_{jk} \sim \sqrt{{\Delta }_{j}{\Delta }_{k}}\ll \sqrt{{\omega }_{j}{\omega }_{k}}.$$

(31)

Now if Γ_j (and thus R_j and I_j) does not vary rapidly on the scale of ${O}_{jk} \sim \sqrt{{\Delta }_{j}{\Delta }_{k}}$ (which is the scale of the Lamb shifts), then in the regime we need to consider we have Γ_j ≈ Γ_k. More specifically, we consider systems for which the spectral density satisfies

$$J({\omega }_{j}+{\Delta }_{k})\approx J({\omega }_{j}),\ \forall j,k,$$

(32)

for which the more precise statement is

$$| J({\omega }_{j}+{\Delta }_{k})-J({\omega }_{j})| \ll J({\omega }_{j}),\ \forall j,k,$$

(33)

since this implies that Γ(ω) also satisfies the same "slow variation” conditions. Under the condition in Eq. (33) we have

$$| {\omega }_{j}-{\omega }_{k}| \lesssim {O}_{jk}\ \ \Rightarrow \ \ \left\{\begin{array}{l}{R}_{j}\approx {R}_{k}\\ {I}_{j}\approx {I}_{k}\end{array}\right.$$

(34)

while at the same time allowing

$${R}_{j}\;\ne\;{R}_{k},\ {I}_{j}\;\ne\;{I}_{k}\ \ \,\text{when}\,\ \ | {\omega }_{j}-{\omega }_{k}| \gg {O}_{jk},$$

(35)

The relation (34) allows us to make the replacements ${R}_{j}\approx \sqrt{{R}_{j}{R}_{k}}$ and ${I}_{j}\approx \sqrt{{I}_{j}{I}_{k}}$ in Eq. (25) because i) when ∣ω_j − ω_k∣ ≲ O_jk these replacements are well justified, and ii) when ∣ω_j − ω_k∣ is larger the terms containing $\sqrt{{R}_{j}{R}_{k}}$ and $\sqrt{{I}_{j}{I}_{k}}$ are eliminated by the rotating wave approximation. The result is

$$\begin{array}{*{20}{l}}\dot{\rho }=-\frac{{\rm{i}}}{{\hslash}}[{H}_{{\rm{sys}}},\rho ]-\,{\rm{i}}\,{\mathop{\sum} \limits_{j}}{I}_{j}[{\Sigma}_{j}^{\dagger }{\Sigma }_{j},\rho ]+{\mathop{\sum}\limits_{j}}2{R}_{j}{\mathcal{D}}[{\Sigma }_{j}]\rho \\-\,{\rm{i}}{\mathop{\sum}\limits_{k\ne j}}\sqrt{{I}_{j}{I}_{k}}\left[{\Sigma }_{j}^{\dagger}{\Sigma }_{k}+{{\Sigma}_{k}^{\dagger}}{{\Sigma}_{j}},\rho \right]\\ - \, 2{\mathop{\sum}\limits_{k\ne j}}\sqrt{{R}_{j}{R}_{k}}\left[{\Sigma }_{j}^{\dagger }{\Sigma }_{k}\rho +\rho {\Sigma }_{k}^{\dagger }{\Sigma }_{j}-2{\Sigma }_{j}\rho {\Sigma }_{k}^{\dagger }\right].\end{array}$$

(36)

The terms in this equation can be re-factored so as to write it in a much neater form, namely

$$\dot{\rho }=-\frac{\,\text{i}\,}{\hslash }\left[{H}_{0}-\hslash {D}^{\dagger }D,\rho \right]+{\mathcal{D}}[\Sigma ]\rho$$

(37)

This is the zero-temperature Lindblad-form master equation for all regimes. The operators Σ and D are

$$\Sigma =\mathop{\sum }\limits_{j=1}^{N}\sqrt{{\gamma }_{j}}{\text{e}}^{{\rm{i}}{\phi }_{j}}{\sigma }_{j},$$

(38)

$$D=\mathop{\sum }\limits_{j=1}^{N}\sqrt{{\Delta }_{j}}{\text{e}}^{{\rm{i}}{\phi }_{j}}{\sigma }_{j},$$

(39)

in which ${\phi }_{j}=\arg [{g}_{j}]$ as defined below Eq. (10). The term H_L ≡ −ℏD^†D is the Lamb-shift Hamiltonian. If we expand it out we see that, so long as ${\sigma }_{j}^{\dagger }{\sigma }_{k}\;\ne\;0$, the upper levels of the transitions j and k are coupled together via the bath:

$$\begin{array}{l}{H}_{{\rm{L}}}=-\hslash \mathop{\sum }\limits_{j=1}^{N}{\Delta }_{j}{\sigma }_{j}^{\dagger }{\sigma }_{j} -\ \hslash \mathop{\sum }\limits_{j=1}^{N}\mathop{\sum }\limits_{k=j+1}^{N}\sqrt{{\Delta }_{j}{\Delta }_{k}}\left({\text{e}}^{{\rm{i}}\Delta {\phi }_{jk}}{\sigma }_{j}^{\dagger }{\sigma }_{k}+{\rm{H}}.{\rm{c}}.\right),\end{array}$$

(40)

where the phases are given by

$$\Delta {\phi }_{jk}={\phi }_{k}-{\phi }_{j}.$$

(41)

The decay rates γ_j are given in Eq. (28) and are determined by the value of the spectral density J(ω) only at transition frequency ω_j. The Lamb shifts, on the other hand, depend on the whole spectral density, and in particular on the cut-off frequency. As an example, for the Ohmic spectrum with a sharp cut-off at Ω, in which J(ω) ∝ ω (we choose to define J(ω) = ω/Ω²) the Lamb shifts are

$$\begin{array}{*{20}{l}}{\Delta }_{j}= | {g}_{j}{| }^{2}\ {\mathbb{P}}\left[{\int\nolimits_{-{\omega }_{j}}^{\Omega -{\omega }_{j}}}\frac{J(\omega\,+\,{\omega }_{j})}{\omega }d\omega \right]\\= \frac{| {g}_{j}{| }^{2}}{{\Omega }^{2}}\ {\mathbb{P}}\left[{\int\nolimits_{-{\omega }_{j}}^{\Omega -{\omega }_{j}}}\frac{(\omega\,+\,{\omega }_{j})}{\omega }d\omega \right]\\= \frac{{\gamma }_{j}}{2\pi }\left[\frac{\Omega }{{\omega }_{j}}+\mathrm{ln}\!\left(\frac{\Omega }{{\omega }_{j}}-1\right)\right].\end{array}$$

(42)

We see that so long as Ω is larger than 2ω_j the Lamb shift is larger than the damping rate by at least a factor of Ω/ω_j. Recall that the Lamb shifts are required to be much less that the transition frequencies. Using the expression for Δ_j above, we have

$$\begin{array}{*{20}{l}}\frac{{\Delta }_{j}}{{\omega }_{j}}= \frac{{\gamma }_{j}}{2\pi }\left[\frac{\Omega }{{\omega }_{j}}+\mathrm{ln}\!\left(\frac{\Omega }{{\omega }_{j}}-1\right)\right]\\\approx \frac{1}{2\pi }\left(\frac{{\gamma }_{j}}{{\omega }_{j}}\right)\left(\frac{\Omega }{{\omega }_{j}}\right).\end{array}$$

(43)

Thus to satisfy the condition Δ_j ≪ ω_j requires that the cut-off frequency is not too large. In particular $\Omega \ll 2\pi {\omega }_{j}^{2}/{\gamma }_{j}$.

We can now evaluate the fidelity of our approximation explicitly for the Ohmic spectrum. Recall that we require Δ_k/Δ_j ≈ 1, when Δω = ∣ω_j − ω_k∣ ≲ Δ_j. Denoting Δ_j by Δ(ω_j), and writing ω_k = ω_j + Δ_j, we have

$$\begin{array}{*{20}{l}}\frac{\Delta ({\omega }_{j}\,+\,{\Delta }_{j})}{\Delta ({\omega }_{j})}-1 = \frac{{\Delta }_{j}}{\Omega }{\mathrm{ln}}\!\left(\frac{\Omega }{{\omega }_{j}}-1\right)+{\mathcal{O}}(\Delta {\omega }^{2})\\\approx \left[\frac{{\Delta }_{j}}{{\omega }_{j}}\right]\left[\frac{{\omega }_{j}}{\Omega }{\mathrm{ln}}\!\left(\frac{\Omega }{{\omega }_{j}}\right)\right]\ll 1.\end{array}$$

(44)

Since the master equation is already derived under the conditions that Ω ≫ ω_j, and ω_j ≫ Δ_j, the expressions in both of the square brackets are individually much less than unity. Thus the slowly varying spectrum approximation that we have introduced is automatically a very good approximation for the Ohmic bath.

Master equation for arbitrary temperature

To derive the master equation for non-zero temperature we merely replace the zero-temperature expectation values of the bath operators with their expectation values at non-zero temperature, which are

$$\left\langle {b}^{\dagger }(\omega ^{\prime} )b(\omega )\right\rangle ={n}_{T}(\omega )\delta (\omega -\omega ^{\prime} ),$$

(45)

$$\left\langle b(\omega ^{\prime} ){b}^{\dagger }(\omega )\right\rangle =[1+{n}_{T}(\omega )]\delta (\omega -\omega ^{\prime} ),$$

(46)

in which

$${n}_{T}(\omega )=\frac{1}{\exp \left[\hslash \omega /\left({k}_{{\rm{B}}}T\right)\right]-1}.$$

(47)

Here T is the temperature of the bath and k_B is Boltzmann’s constant. With these new expectation values we now obtain more terms in the master equation. For the new terms the spectral density J(ω) is multiplied by n_T(ω), so the new terms give a new integral for which we have to calculate the principle value. For the Ohmic bath this integral is

$$\begin{array}{*{20}{l}}{\rm{Im}}\,[{\tilde{\Gamma }}_{j}^{{\rm{T}}}]= {\mathbb{P}}\left[{\int \nolimits_{-{\omega }_{j}}^{\Omega -{\omega }_{j}}}\frac{J(\omega +{\omega }_{j}){n}_{T}(\omega +{\omega }_{j})}{\omega }d\omega \right]\\= \frac{1}{{\Omega }^{2}}{\mathbb{P}}\left[{\int \nolimits_{-{\omega }_{j}}^{\Omega -{\omega }_{j}}} \frac{1+({\omega }_{j}/\omega )}{\exp [\hslash (\omega +{\omega }_{j})/({k}_{{\rm{B}}}T)]-1}d\omega \right].\end{array}$$

(48)

Unfortunately this integral does not have an analytic solution, so we leave it as an integral and define a new set of Lamb shifts

$${\Delta }_{j}^{{\rm{T}}}\equiv | {g}_{j}{| }^{2}\,{\rm{Im}}\,[{\tilde{\Gamma }}_{j}^{{\rm{T}}}].$$

(49)

The approximations we used for the zero temperature part of the master equation can be applied in exactly the same way to the new terms that appear at non-zero temperature. The resulting master equation for arbitrary temperatures is

$$\dot{\rho }=-\frac{\,\text{i}\,}{\hslash }[{H}_{0}+{H}_{{\rm{L}}},\rho ]-{\mathcal{D}}[\Theta (T)]\rho -{\mathcal{D}}[\Upsilon (T)]\rho ,$$

(50)

where

$$\Theta (T)=\mathop{\sum }\limits_{j=1}^{N}\sqrt{{\gamma }_{j}\left[1+{n}_{T}({\omega }_{j})\right]}\ {\text{e}}^{{\rm{i}}{\phi }_{j}}{\sigma }_{j},$$

(51)

$$\Upsilon (T)=\mathop{\sum }\limits_{j=1}^{N}\sqrt{{\gamma }_{j}{n}_{T}({\omega }_{j})}\ {\text{e}}^{-{\rm{i}}{\phi }_{j}}{\sigma }_{j}^{\dagger },$$

(52)

and

$${H}_{{\rm{L}}}=-\hslash \left[{B}^{\dagger }B-C{C}^{\dagger }\right],$$

(53)

with

$$B=\mathop{\sum }\limits_{j=1}^{N}\sqrt{{\Delta }_{j}+{\Delta }_{j}^{{\rm{T}}}}\ {\text{e}}^{{\rm{i}}{\phi }_{j}}{\sigma }_{j},$$

(54)

$$C=\mathop{\sum }\limits_{j=1}^{N}\sqrt{{\Delta }_{j}^{{\rm{T}}}}\ {\text{e}}^{{\rm{i}}{\phi }_{j}}{\sigma }_{j}.$$

(55)

Note that when T > 0 the bath induces a Hamiltonian coupling not only between the upper levels of the different transitions, but also the lower levels.

The non-degenerate master equation and numerical efficiency

In deriving the master equation valid for all regimes, we did not make the secular approximation, which involves dropping terms that oscillate at the frequency difference between different transitions. However, when the difference between the frequencies of two transitions is much larger than the Lamb shifts and linewidths, keeping the resulting rapidly oscillating terms greatly increases the numerical overhead while contributing little to the evolution. In this case one should drop these terms for numerical efficiency. Doing so transforms the master equation into the non-degenerate master equation, but only for those pairs of transitions for which the detuning is very large. If we write the master equation in the form given in Eq. (37), then dropping the rapidly oscillating terms means merely dropping terms in the second and third lines for the pairs of values of j and k whose transitions are detuned by much more that their Lamb shifts and linewidths.

For readers very familiar with the degenerate and non-degenerate master equations, the result of applying the secular approximation to a subset of pairs of transitions will likely be clear. For readers without this familiarity, we give an explicit example. Let us say that we can divide our transitions into two sets, where the frequencies of those in the first set differ from the frequencies of those in the second set by at least $1{0}^{3}{\Delta }_{\max }$ in which ${\Delta }_{\max }$ is the maximum Lamb shift among all the transitions. If we denote the transition operators in the first set by ${\sigma }_{j}^{(1)}$, with j = 1, …, N₁, and those in the second by ${\sigma }_{j}^{(2)}$, with j = 1, …, N₂, then the result of making the secular approximation on the master equation in Eq. (37) is

$$\dot{\rho }=-\frac{\,\text{i}\,}{\hslash }\left[{H}_{0}-\hslash \mathop{\sum }\limits_{m=1}^{2}{D}_{m}^{\dagger }{D}_{m},\rho \right]+\mathop{\sum }\limits_{m=1}^{2}{\mathcal{D}}[{\Sigma }_{m}]\rho ,$$

(56)

with

$${\Sigma }_{m}=\mathop{\sum }\limits_{j=1}^{{N}_{m}}\sqrt{{\gamma }_{j}^{(m)}}{\text{e}}^{{\rm{i}}{\phi }_{j}^{(m)}}{\sigma }_{j}^{(m)},$$

(57)

$${D}_{m}=\mathop{\sum }\limits_{j=1}^{{N}_{m}}\sqrt{{\Delta }_{j}^{(m)}}{\text{e}}^{{\rm{i}}{\phi }_{j}^{(m)}}{\sigma }_{j}^{(m)}.$$

(58)

Additional confirmation of accuracy for the Ohmic bath

It is clear from the derivation in the previous section that the master equation we have obtained, given in Eqs. (37) and (50), will be valid so long as the variation of the Lamb shifts Δ(ω_j) and damping rates γ(ω_j) on the scale of these same Lamb shifts and damping rates is sufficiently small. This variation of the Lamb shifts and damping rates will be small if the variation of the spectral density J(ω) on the scale of the Lamb shifts and damping rates is sufficiently small. Below we will verify quantitatively the accuracy of the new master equation for the Ohmic spectrum at zero temperature, using exact simulations of two example systems. The Ohmic spectrum is appropriate for systems such as atoms, color centers, or superconducting qubits coupled to one-dimensional wave-guides or transmission lines.

Since the efficacy of the approximation used to derive the master equation does not depend on the temperature of the bath or the specific functional form of the spectral density (it depends only on the local variation of the resulting Lamb shifts and damping rates around their respective transition frequencies), simulations for the Ohmic bath at zero temperature provide a high level of confidence in the accuracy of the master equation more generally. Further, since we have exact analytic expressions for the Lamb shifts and damping rates in this case, if desired the variations of these quantities at the transition frequencies can be related directly to the accuracy determined in our simulations.

We have already compared the evolution of the master equation to exact simulations for the V system in Fig. 2. We now consider two further systems. The first is the "trident” system depicted in Fig. 1b. This system has three transitions, and thus also three pairs of transitions which we can place simultaneously in the near-degenerate regime. Again using $\tilde{\nu }$ as our arbitrary frequency reference, we choose parameters ${\omega }_{1}=10\pi \tilde{\nu }$, ω_j = ω₁ + γ₂/(j − 1), for j = 2, 3, and ${\gamma }_{j}=[(5-j)/40]\tilde{\nu }$, for j = 1, 2, 3, with the cut-off frequency $\Omega =80\pi \tilde{\nu }$. Choosing the initial state $\left|{\psi }_{0}\right\rangle =(7i\left|1\right\rangle +3\left|2\right\rangle )/\sqrt{58}$, we plot the evolution of the populations predicted by the master equation along with the exact evolution in Fig. 3a. The maximum error in the evolution of the master equation over the time period plotted in Fig. 3a is less than 2 × 10⁻³.

Fig. 3: Comfirming the accuracy of the Lindblad master equation for two examples.

A comparison of the master equation given in Eq. (37) with the exact evolution for two open systems. The evolution of the master equation is shown as dashed lines and the exact evolution as solid lines. a The populations of the three upper levels of the trident system depicted in Fig. 1b with initial state $\left|{\psi }_{0}\right\rangle =(7\,\text{i}\,\left|1\right\rangle +3\left|2\right\rangle )/\sqrt{58}$. b The populations of the two upper levels of the two-qubit system depicted in Fig. 1d with the initial state $\left|{\psi }_{0}\right\rangle =(\,\text{i}\,{\left|1\right\rangle }_{1}{\left|0\right\rangle }_{2}+{\left|0\right\rangle }_{1}{\left|1\right\rangle }_{2})/\sqrt{2}$ and the parameters ${\gamma }_{1}={\gamma }_{2}=0.1\tilde{\nu }$, ${\omega }_{1}=10\pi \tilde{\nu }$, and ω₂ = ω₁ + 2γ.

We now perform simulations for two co-located qubits, whose level structure is depicted in Fig. 1d. We choose the parameters ${\omega }_{1}=10\pi \tilde{\nu }$, ω₂ = ω₁ + 2γ₁, and ${\gamma }_{1}={\gamma }_{2}=0.1\tilde{\nu }$, with the same cut-off frequency as before. We find that this system requires significantly larger values of the weak damping parameters ("quality factors”), Q_j ≡ ω_j/γ_j, in order for the master equation to accurately model the dynamics. Since available numerical resources place restrictions on the sizes of the Q_j’s that we can practically simulate, for this system we apply the first rotating-wave approximation to our model Hamiltonian prior to performing the exact simulations. That is, we simulate the Hamiltonian H_RWA (Eq. 15) instead of the full model in Eq. 7. These simulations thus show us how well the master equation will perform so long as the Q_j’s are large enough to satisfy the first rotating-wave approximation. We stress that the values of the Q_j’s we actually simulate here are not large enough to satisfy this approximation for this system. This fact is interesting in itself, because it shows that different systems, even with only a few levels, can require quite different quality factors to reach the weak damping regime. We believe this is due to the availability of more channels via which the off-resonant terms in the system bath interaction can excite the two-qubit system over the V and trident systems.

Choosing the initial state $\left|{\psi }_{0}\right\rangle =(i{\left|1\right\rangle }_{1}{\left|0\right\rangle }_{2}+{\left|0\right\rangle }_{1}{\left|1\right\rangle }_{2})/\sqrt{2}$, we show the evolution of the populations for both the master equation and the exact simulations of H_RWA in Fig. 3b. The error in the evolution of the master equation over the duration shown in Fig. 3a is less than 5 × 10⁻³.

Finally, we compare the evolution of the Bloch–Redfield equation (Eq. 24) both to our Lindblad master equation and the exact simulations. These comparisons, plots of which are given in the Supplementary Figs. 1 and 2, confirm that the B-R equation and our master equation have essentially the same accuracy. This result is implied, of course, by the derivation of the Lindblad master equation.

Regime of validity

We have shown that when the spectral density is sufficiently flat the B–R equation for thermal damping can be replaced by a Lindblad equation, and that this is an excellent approximation for the Ohmic bath. The question we need to answer now is whether the slowly varying spectrum (SVS) approximation remains an excellent approximation over the entire domain in which the B–R equation itself is valid, or whether there is a regime in which the B–R equation is valid but the Lindblad equation is not. Since the SVS approximation depends solely on the slope of the spectral density (strictly, the difference between the values of the spectral density at the frequencies of nearby transitions) the question is, as we increase this slope, whether the B–R equation deviates from the exact evolution before or after the Lindblad and B–R equations deviate from each other?

If we use an Ohmic spectrum, then we cannot increase the slope of the spectral density without similarly increasing the decay rate(s). We already know that the master equations break down (deviate from the exact evolution) at large damping. To explore how the slope of the spectral density affects the master equations we thus require a new spectral density, and we use the one depicted in Fig. 4. This density, which we will denote by J_r(ω), is divided into three segments. It has a constant slope in each segment, with all but the middle segment having the same slope as the Ohmic spectrum (dJ_r(ω)/dω = 1/Ω²). In the middle segment, which includes the transition frequency(ies) of the system, this slope is increased by a factor of r.

Fig. 4: The piece-wise linear spectrum, J_r(ω).

We use this piece-wise spectrum to examine the effect of the slope of the spectral density on the accuracy of the master equations. On the two outer segments this spectrum has the slope of the Ohmic spectrum, while on the middle segment the slope is increased by a factor r.

It has already been established that the spectral density must be sufficiently flat in order that the damping induced by the bath be exponential^32,53. We can use system identification, which we used above to find the master equation for the V-system, to determine the number of dynamical variables required to reproduce the evolution of the open system as the slope of the spectral density is increased. When this number is greater than that possessed by the system, the evolution is non-Markovian and all time-independent Markovian master equations will break down. This allows us to determine not only how the B–R equation performs, but whether any time-independent Markovian master equation is able to model baths with steeply varying spectra.

We examine the exact evolution of both a two-level system and the three-level V-system of Fig. 1a. The frequency of the two-level transition is ${\omega }_{0}=10\pi \tilde{\nu }$, with a damping rate of $\gamma =0.296\tilde{\nu }$ when r = 1, while the V-system has ω₁ = ω₀ − γ/2, ω₂ = ω₀ + γ/2, and ${\gamma }_{1}={\gamma }_{2}=0.2\tilde{\nu }$. Since we have chosen a higher damping rate for the two-level system the "baseline" error of the master equations for this system will be a little higher than that for the V-system.

As we increase the slope of J_r, all the damping rates do increase (at first only a little) because the values of J_r(ω) at the transition frequencies also increase. The increase in the error of the master equations will thus have two sources, the increasing slope and the increasing damping rate(s). By determining the error of the master equations as we increase the damping rates while keeping the spectrum Ohmic, we can largely distinguish the relative contributions of the two sources of error.

In Fig. 5a we plot the population of the excited state of the two-level system as it decays into the bath, for r = 2ⁿ, with n = 0 (the Ohmic spectrum) and n = 3, 4, …, 12. We see that as the slope increases the evolution distorts away from exponential decay. Since the system starts in the excited state, the coherences remain zero and there is only one independent variable in the evolution. Since a time-independent linear equation can only generate real non-exponential evolution if it has more than one dynamical variable, any deviation from exponential behavior necessarily implies non-Markovian evolution. In Fig. 5c we plot the error in the evolution of the master equations (note that both the Lindblad and B–R master equations are identical for a two-level system) as a function of n (equivalently ${\mathrm{log}\,}_{2}r$). Note that the evolution of the master equation is determined solely by the damping rate γ(r) = 2πg²J_r(ω₀). Rather than using the damping rate as given by the master equation, for each value of r we can alternatively choose the damping rate that minimizes the error. We also plot this "best fit" error in Fig. 5c, as well as the error of the master equations for an Ohmic spectrum with the damping rate γ(r). We see that the effect of the increased slope on the error is much greater than that of the increased damping rate alone. We also see that the "best-fit" error, while smaller than that given by the value of γ specified by the master equation, increases just as rapidly as the latter.

Fig. 5: Breakdown of the Lindblad and Bloch–Redfield (B–R) master equations.

Here we show how the Lindblad and Bloch–Redfield (B–R) master equations break down as the condition on the flatness of the spectral density is relaxed. Plots a and c show this break down for a two-level system and plots b and d for the three-level V-system of Fig. 1a. a The population of the excited state as it decays, for a range of values of the slope of the spectral density at the transition frequency, $J^{\prime} ({\omega }_{0})=r/{\Omega }^{2}$. The blue curve is the Ohmic spectrum, and rest of the curves, from green through violet, show the evolution for r = 2ⁿ, with n = 3, 4, …, 12. b The population of level $\left|2\right\rangle$ of the V-system with r = 2⁷ (blue) along with the predictions of the Lindblad and B–R master equations (green and turquoise, respectively). The orange and pink curves show the "best-fit" to the exact dynamics for both master equations, obtained by fitting the parameters R_j and I_j. c The blue curve is error of the master equation in describing the decay of the excited state when the spectrum is J_r, as r is increased. The red curve is the error of the master equation when the spectrum is Ohmic and the damping rate is the equivalent of that set by J_r, namely γ(r) = 2πg²J_r(ω₀). The orange curve shows the error of the master equation for the spectrum J_r when we fit the damping rate. d The blue and magenta curves give, respectively, the errors of the Linblad and B-R master equations in describing the evolution of the population of the initial state $\left|2\right\rangle$, when the spectrum is J_r. (Note that these errors are sufficiently similar that the magenta curve is almost entirely obscured by the blue curve.) The red and dark red curves, respectively, give the same errors when the spectrum is Ohmic and the damping rates are the equivalents of those set by J_r. The orange and purple curves give the same errors when we choose the values of R_j and I_j so as to give the best fit to the exact dynamics for each master equation.

We use system identification to determine the dimension of the dynamics that generates the exact evolution for each value of r. As above, we define this dimension as the number of dynamical variables required to account for 0.999 of the combined magnitude of the dynamical eigenvalues. The results are presented in Table 1. We see that the dimension, and thus the non-Markovianity, starts increasing as soon as we increase the slope, showing that all time-independent Markovian master equations break down.

Table 1 Effective dimension of the two-level system and V-system.

We now examine whether there is difference between the accuracy of the Lindblad and B–R master equations as the slope, r, increases. For this we need to explore the evolution of the V-system. We choose the initial state of this system to be $\left|2\right\rangle$, and plot in Fig. 5d the error of both master equations as a function of n, as well as the error when R_j and I_j are chosen so as to give the best fit to the exact dynamics. As Fig. 5d shows, the behavior of the V-system is very similar to that of the two-level system. Even when both master equations have deviated significantly from the exact dynamics, they remain so close to each other that the difference in their respective errors is hardly distinguishable on the plot. Thus the master equations break down well before they deviate from each other, and thus well before the SVS approximation breaks down. We also note that the B-R equation continues to maintain positivity to good approximation up through r = 10; the magnitude of the most negative eigenvalue of the density matrix remains below 1 × 10⁻¹².

In Fig. 5b We show how the master equations compare to the exact evolution for the V-system when n = 7 (r = 128). Notable is how different the master equations are from the exact evolution, while being close to eachother. The fact that the master equations are able to approximate the exact evolution considerably better given the optimal choices for R_j and I_j is rather interesting. It suggests that there might be some way to develop improved formulae for R_j and I_j.

Time-dependent problems: accuracy of the adiabatic extension

Many important problems involve open systems whose Hamiltonians change with time. Our master equation can be used to describe these systems if the time-dependence is not too fast. To do so one takes the master equation and changes the parameters and operators that appear in it, namely γ_j(ω_j), Δ_j(ω_j), and σ_j (which depend on the system eigenstates and thus on the system Hamiltonian), so that at each time they take the values determined by the Hamiltonian of the system at that time⁵⁴. The resulting time-dependent "adiabatic” master equation will be effective for sufficiently slow changes in the Hamiltonian.

Here we examine some examples to confirm that the adiabatic version of the master equation is accurate even when two levels of an open system cross each other, or move from degenerate to near-degenerate, during the evolution. We consider first the V system in which both damping rates are equal and the detuning changes with time. We start the system in the state $\left|{\psi }_{-}\right\rangle \equiv (\left|1\right\rangle -\left|2\right\rangle )/\sqrt{2}$, which for Δω = 0 will not decay since it is the (sub-radiant) dark state^15,25. We then change the detuning with time as determined by the following two functions:

$${F}_{1}(t)=\left\{\begin{array}{*{20}{l}}0,&\ \ 0< t <\frac{T}{4},\\ \frac{\pi }{64}\left(t-\frac{T}{4}\right),&\ \ \frac{T}{4}< t < T,\end{array}\right.$$

(59)

$${F}_{2}(t)=\left\{\begin{array}{*{20}{l}}0,&\ \ 0< t < \frac{T}{4},\\ \frac{\pi }{2},&\ \ \frac{T}{4}< t <\frac{3T}{4},\\ 0,&\ \ \frac{3T}{4}< t < T.\end{array}\right.$$

(60)

The function F₁ is chosen so that the detuning increases gradually, while F₂ involves rapid changes. In Fig. 6 we compare the adiabatic version of the master equation with the exact evolution for the two cases. For Δω(t) = F₁(t) the maximum error of the adiabatic master equation is 3.4 × 10⁻³, and for Δω(t) = F₂(t) the maximum error is 1.5 × 10⁻².

Fig. 6: Accuracy of the Lindblad master equation for two time-dependent examples.

A comparison of the evolution predicted by the adiabatic extension of the new master equation, Eq. (37), with an exact simulation of a V system coupled to an Ohmic bath (Fig. 1a), in which the detuning between the transitions, Δω, has the time-dependence given in Eqs. (59) and (60). The initial state is $\left|{\psi }_{-}\right\rangle =(\left|1\right\rangle -\left|2\right\rangle )/\sqrt{2}$ which is the dark state for degenerate transitions with equal damping rates. The transition frequency ${\omega }_{1}=3\pi \tilde{\nu }$, the damping rates are ${\gamma }_{1}={\gamma }_{2}=\gamma =0.1\tilde{\nu }$, and the simulation time is T = 8. (a,c,e) Δω(t) = F₁(t) (Eq. (59)). a The populations of levels $\left|1\right\rangle$ (blue) and $\left|2\right\rangle$ (red), with the evolution of the adiabatic master equation denoted by dashed lines and that of the exact simulation by solid lines. c The absolute value of the difference between the populations predicted by the adiabatic master equation and the exact evolution. e The detuning as a function of time. b, d, f The same set of results but with the detuning given by F₂(t) (Eq. (60)).

As our final example we consider a generalized version of the Landau-Zener transition^45,46, in which the energies of two coupled levels cross each other. In particular we use the 4-level system depicted in Fig. 1c, in which we add a coupling between the upper two levels. In the original Landua-Zener transition, for which there is an analytic solution, the energy of one of the levels is fixed and the other increases linearly with time. We generalize this by choosing the following sinusoidal time-dependence for ω₁:

$${\omega }_{1}(t)={\omega }_{2}-{\Delta }_{0}\cos (\nu t),\ 0\le t\le \frac{\pi }{\nu }.$$

(61)

The detuning, Δω ≡ ω₁ − ω₂, starts at − Δ₀ and increases as a sinusoid through zero to end at Δ₀. We plot the evolution of the populations of the two levels in Fig. 7. The maximum error of the adiabatic extension of the master equation is less than 5.4 × 10⁻³.

Fig. 7: Accuracy of the Linblad master equation for an example Landau-Zener transition.

Here we plot the evolution resulting from a generalized Landau-Zener transition in which the energies of two coupled levels cross. The levels are the upper two levels of the four-level system depicted in Fig. 1c. The energy of level $\left|2\right\rangle$ is fixed so that ${\omega }_{2}=2\pi \tilde{\nu }$, and that of level $\left|1\right\rangle$, in which all the population starts, increases with time. The detuning between the levels is shown as a function of time in the inset. We plot the populations of the two levels, both the exact evolution (solid) and that predicted by the master equation (dashed). The damping rates of the two levels are ${\gamma }_{2}=2{\gamma }_{1}=0.05\tilde{\nu }$, the coupling between them is $c=0.2\tilde{\nu }$, the initial detuning is ${\Delta }_{0}=(\pi /2)\tilde{\nu }$, and the period of the sinusoid is $2\pi /\nu =64/\tilde{\nu }$.

Discussion

We have shown that a slow variation condition on the spectral density is necessary for Markovianity in the oscillator-bath model of thermalization, and along with weak damping and a high cut-off frequency is a sufficient condition for the existence of an accurate Lindblad master equation for all regimes of detuning. In doing so we have shown that the Bloch–Redfield equation for thermal damping can be replaced with this Lindblad equation. This resolves the long-running controversy with the Bloch–Redfield equation, and confirms the conjectures of Eastham et al.²⁸.

The new master equation unifies the existing Lindblad master equations for degenerate and non-degenerate systems and in doing so provides insight into the dynamics of the near-degenerate regime. It also allows both the use of efficient Monte Carlo methods and a measurement description of the action of a thermal bath for all regimes. Further, its adiabatic extension provides a powerful tool for simulating systems in which transitions are time-dependent and cross during the evolution, so long as this time-dependence is not too fast. This suggests that further exploration of the accuracy of the adiabatic extension as a function of the speed of the time-dependence may be a worthwhile endeavor. Such an exploration would help to delineate the class of controlled systems for which it is effective.

The technique of system identification played an important role in obtaining the master equation, as well as determining when open systems are Markovian. As far as we are aware, system identification has not been used before as a tool to understand the dynamics of open quantum systems, or emergent phenomena in many-body systems more generally. We expect that it will prove to be powerful for exploring a wide range of problems in open systems and many-body physics.

Methods

System identification for linear systems

Here we present the method we used to determine the minimal model of a linear system given a knowledge of the evolution of a subset of the system’s state-space. In our case the linear system consists of a low-dimensional quantum system interacting with a high-dimensional bath, and the subset of the state space that we can observe is the density matrix of the low-dimensional system. The following method is one of a family of elegant methods referred to as subspace identification methods, adapted so as to use a set of initial conditions rather than a set of inputs. Further information on subspace identification methods can be found in^55,56,57,58.

Let us say we have a high-dimensional system with dimension J (in our case the open system and the bath), and we have the ability to observe N < J variables of the system, as well as to evolve the system with any choice of initial conditions for the N variables we can observe. We would like to find an accurate model (another linear system) that generates the evolution of the N variables but is only M dimensional with N ≤ M < J. Let us denote a state of the full J-dimensional system by the vector v, and the subset of N variables in which we are interested by the N-dimensional vector x. The map that gives the state of the total J-dimensional system at time τ given an initial state v(0) we will call Z(τ) so that v(τ) = Z(τ)v(0). Defining Z ≡ Z(τ) we note that Z(nτ) = Zⁿ. We also define the non-square projector P that projects onto the N variables so that x = Pv.

Given the ability to evolve the total system with any choice of initial conditions for the N-dimensional subsystem, along with a single choice for the initial values of the rest of the variables (of which there are N − J), we can obtain the matrices Y_n = PZ(nτ)P^T that maps the N variables at time 0 to their values at time τ, for any time τ.

We now construct the following two symmetric "block Hankel” matrices:

$${H}_{0}\equiv \left(\begin{array}{*{20}{l}}{Y}_{0}&{Y}_{1}&\cdots \ &{Y}_{n}\\ {Y}_{1}&{Y}_{2}&&{Y}_{n+1}\\ \vdots &&\ddots &\vdots \\ {Y}_{n}&{Y}_{n+1}&\cdots \ &{Y}_{2n}\end{array}\right),$$

(62)

$${H}_{1}\equiv \left(\begin{array}{*{20}{l}}{Y}_{1}&{Y}_{2}&\cdots \ &{Y}_{n+1}\\ {Y}_{2}&{Y}_{3}&&{Y}_{n+2}\\ \vdots &&\ddots &\vdots \\ {Y}_{n+1}&{Y}_{n+2}&\cdots \ &{Y}_{2n+1}\end{array}\right).$$

(63)

We now note that H₀ can be written as an outer product H₀ = CD^T of (non-square) matrices given by

$$C=\left(\begin{array}{*{20}{l}}P\\ PZ\\ \vdots \\ P{Z}^{m}\end{array}\right)=P\left(\begin{array}{*{20}{l}}I\\ Z\\ \vdots \\ {Z}^{m}\end{array}\right),$$

(64)

$${D}^{{\rm{T}}}=\left(\ Z\ {Z}^{2}\ \cdots \ {Z}^{m}\right){P}^{{\rm{T}}}.$$

(65)

We can now determine C and D by doing a singular value decomposition of H₀ to give H₀ = USV^T. Note that since the smaller dimension of the matrices C and D is smaller than that of H₀ we expect many of the columns of U and the rows of V^T will be zero, as will many of the eigenvalues of H₀ which are given in the diagonal matrix S. Note that we can now view P and Z as defining a linear model that generates evolution of the N-dimensional subsystem. The number of eigenvalues that are appreciably non-zero tells us the dimension of the model.

To distinguish the model from the original total system we started with, we can write the matrices C and D as

$${C}^{{\rm{T}}}=\left(\ {M}^{{\rm{T}}}\ {[{M}^{{\rm{T}}}]}^{2}\ \cdots \ {[{M}^{{\rm{T}}}]}^{m}\right){Q}^{{\rm{T}}},$$

(66)

$${D}^{{\rm{T}}}=\left(\ M\ {M}^{2}\ \cdots \ {M}^{m}\right){Q}^{{\rm{T}}},$$

(67)

where M is the evolutionary map for the model and Q is the projector onto the subsystem. Let us now decompose H₀ = USV^T into the outer product of two vectors $\tilde{C}\equiv U\sqrt{S}{P}^{{\rm{T}}}$ and $\tilde{D}={(P\sqrt{S}{V}^{{\rm{T}}})}^{{\rm{T}}}$. Noting that

$${H}_{1}=CZ{D}^{{\rm{T}}}$$

(68)

we can obtain the evolutionary map for the model, M, from H₁ using

$$M={(\tilde{C}{\tilde{C}}^{{\rm{T}}})}^{-1}{\tilde{C}}^{{\rm{T}}}{H}_{1}\tilde{D}{({\tilde{D}}^{{\rm{T}}}\tilde{D})}^{-1}.$$

(69)

We note that the above method determines the equations of motion of the linear system, but does not itself give the change of basis between the original system variables and those that appear in the obtained equations. Obtaining this change of basis requires additional methods.

Including a diagonal system-bath coupling operator: dephasing

In our derivation of the master equation above we set the diagonal component of the interaction operator, D, to zero. To include D we note that it can be expressed as

$$D={\kappa }_{0}I+\mathop{\sum }\limits_{j=1}^{N}{\kappa }_{j}{z}_{j}$$

(70)

for some real numbers κ₀, κ_j. Here we have defined ${z}_{j}\equiv {\sigma }_{z}^{(j)}=\left|{m}_{j}\right\rangle \left\langle {m}_{j}\right|-\left|{n}_{j}\right\rangle \left\langle {n}_{j}\right|$ as the Pauli z-operator for transition j. Note that the term κ₀I does not contribute to the dynamics of the system.

In driving the master equation we now have additional terms of the form of K in Eq. (18), but with the operators ${\sigma }_{j}^{{\rm{I}}}$ replaced by ${z}_{j}^{{\rm{I}}}={z}_{j}$ (The z_j operators have no evolution in the interaction picture because they commute with H_sys). Note that there are no terms of the form of K in which there is one z_j and one ${\sigma }_{k}^{{\rm{I}}}$ because the rotating-wave approximation eliminates them. Since the ${z}_{j}^{{\rm{I}}}$ do not oscillate, the terms they generate are "at zero frequency" and thus cannot be near-degenerate with any of the damping terms. As such, in proceeding to process the new terms to derive the B–R equation, one immediately obtains a Lindblad form for these terms and the new approximation we have introduced is not required. The contribution to the master equation is

$$\dot{\rho }=-\,\text{i}\,\left[{Z}^{2},\rho \right]-\left[Y,\left[Y,\rho \right]\right]$$

(71)

where

$$Z=\sum _{j}\sqrt{2\pi {\kappa }_{j}^{2}J(0)}\ {z}_{j}$$

(72)

$$Y=\sum _{j}{\kappa }_{j}{\left[\mathop{\int}\nolimits_{0}^{\Omega }\frac{J(\omega )}{\omega }d\omega \right]}^{1/2}{z}_{j}$$

(73)

Note that the integral in Eq. (73) will only exist so long as J(ω) goes to zero sufficiently fast as ω → 0. But in this case J(0) = 0 so the Hamiltonian contribution vanishes. The double commutator in Eq. (71) gives collective dephasing of the transitions.

Principle value of an integral

In deriving our Lindblad master equation we used the fact that

$$\int\ F(x)\left[{\int \nolimits_{0}^{\infty }}\, \sin ([\omega -{\omega }_{0}]s)ds\right]dx={\mathbb{P}}\left[\int\frac{F(x)}{x}dx\right]$$

(74)

for any smooth function F(x), in which ${\mathbb{P}}\left[\cdots \ \right]$ denotes the principle value of a divergent integral. The principle value of an integral that diverges at a point a (where a ∈ (b, c)), is defined by

$${\mathbb{P}}\left[{\int \nolimits_{b}^{c}}f(x)dx\right]\equiv {\mathop {\lim }\limits_{\varepsilon \to 0}}\left[{\int \nolimits_{b}^{a-\varepsilon}} f(x)dx\,+\,{\int \nolimits_{a+\varepsilon }^{c}}\ f(x)dx\right].$$

(75)

Since the divergent function f(x) = 1/x is anti-symmetric and diverges at a = 0, it is simple to evaluate the principle value of $\mathop{\int}\nolimits_{-b}^{c}(1/x)dx$. Assuming that b and c are positive and c > b we have

$$\begin{array}{*{20}{l}}{\mathbb{P}}\left[{\int \nolimits_{-b}^{c}}\frac{dx}{x}\right]\equiv \mathop {\lim }\limits_{\varepsilon \to 0}\left[{\int \nolimits_{-b}^{-\varepsilon }}\frac{dx}{x}+{\int \nolimits_{\varepsilon }^{c}}\frac{dx}{x}\right]\\ = \mathop {\lim }\limits_{\varepsilon \to 0}\left[{\int \nolimits_{-b}^{-\varepsilon}}\frac{dx}{x}+ {\int \nolimits_{\varepsilon }^{b}}\frac{dx}{x}\right]+ {\int \nolimits_{b}^{c}}\frac{dx}{x}\\= {\int \nolimits_{b}^{c}}\frac{dx}{x}=\mathrm{ln}\,(c/b).\end{array}$$

(76)