Accurate Lindblad-Form Master Equation for Weakly Damped Quantum Systems Across All Regimes

TL;DR

This paper introduces a new Lindblad-form master equation that accurately describes weakly-damped quantum systems across all regimes, overcoming limitations of previous models and ensuring positivity and efficiency.

Contribution

The authors derive a universal Lindblad-form master equation valid for all regimes of weak damping, unifying previous approaches and ensuring complete positivity and computational efficiency.

Findings

The new master equation is accurate across all regimes of damping.

It guarantees the positivity of the density matrix.

It enables efficient simulation of complex quantum systems.

Abstract

Realistic models of quantum systems must include dissipative interactions with an environment. For weakly-damped systems the Lindblad-form Markovian master equation is invaluable for this task due to its tractability and efficiency. This equation only applies, however, when the frequencies of any subset of the system's transitions are either equal (degenerate), or their differences are much greater than the transitions' linewidths (far-detuned). Outside of these two regimes the only available efficient description has been the Bloch-Redfield (B-R) 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 the area of quantum technologies. Here we solve this long-standing problem. We discover that a condition on the slope of the spectral density is sufficient to derive a Lindblad form master equation that is accurate for all regimes. We further show that this condition is necessary for weakly-damped systems to be described by the B-R equation or indeed any Markovian master equation. We thus obtain a replacement for the B-R equation over its entire domain of applicability 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.

Tables

Table 1. Table 1: The effective dimension of the damped two-level system, D 2 subscript 𝐷 2 D_{2} , and that of the damped three-level V-system, D V subscript 𝐷 𝑉 D_{V} , as the slope of the spectral density, J ′ ​ ( ω ) = r / Ω 2 superscript 𝐽 ′ 𝜔 𝑟 superscript Ω 2 J^{\prime}(\omega)=r/\Omega^{2} , is increased for a fixed cut-off frequency, Ω Ω \Omega

$r$	$2^0$	$2^1$	$2^2$	$2^3$	$2^4$	$2^5$	$2^6$	$2^7$	$2^8$	$2^9$	$2^{10}$	$2^{11}$	$2^{12}$	$2^{13}$
$D_2$	1	2	2	2	3	4	4	4	4	5	5	5	5	5
$D_V$	4	5	6	7	8	9	10	10	11	11	12	12	12	12

Equations

$$\dot{\rho }$$

$$H_0$$

$$\Sigma$$

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

$$H_{\text{mbox}L}$$

$$H_0$$

$$H$$

$$+{\int }_0^{\Omega }\hslash \omega b(\omega {)}^{†}b(\omega )d\omega .$$

$$H_{\text{mbox}sys}=\sum _n{E}_n|n⟩⟨n|,$$

$${\sigma }_j=|n_j⟩⟨m_j|,E_{m_j}>E_{n_j},j=1,\dots ,J,$$

$$A$$

$${\omega }_j=\frac{E_{m_j}-E_{n_j}}{\hslash },$$

$${\sigma }_j^{\text{mbox}I}$$

$$b^{\text{mbox}I}(\omega )$$

$$H_{\text{mbox}OR}^{\text{mbox}I}$$

$$H_{\text{mbox}RWA}$$

$$+{\int }_0^{\Omega }\hslash \omega b(\omega {)}^{†}b(\omega )d\omega .$$

$$\frac{d{\rho }^{\text{mbox}I}}{dt}$$

$$\displaystyle H_{\mbox{\scriptsize R}}^{\mbox{\scriptsize I}}=\hbar\int_{0}^{\Omega}\!\!\!\!\sqrt{J(\omega)}\Biggl{[}\sum_{j}g_{j}\sigma_{j}^{\mbox{\scriptsize I}\dagger}b^{\mbox{\scriptsize I}}(\omega)+\mbox{H.c}\Biggr{]}d\omega$$

$$K$$

$$G(\omega )$$

$$=J(\omega )e^{-i\omega (t-s)}.$$

$$K$$

$$=\left[{\int }_0^{\infty }{\int }_0^{\Omega }G(\omega -{\omega }_j)d\omega ds\right]e^{i({\omega }_k-{\omega }_j)t}{\sigma }_k^{†}{\sigma }_j{\rho }^{\text{mbox}I}(t)$$

$$=\left[{\int }_0^{\Omega }{\int }_0^{\infty }G(\omega -{\omega }_j)dsd\omega \right]{\sigma }_k^{\text{mbox}I†}(t){\sigma }_j^{\text{mbox}I}(t){\rho }^{\text{mbox}I}(t).$$

$${\Gamma }_j$$

$$R_j$$

$$I_j$$

$$\dot{\rho }$$

$$-i\sum _{k\ne j}{I}_k\left[{\Sigma }_j^{†}{\Sigma }_k\rho -\rho {\Sigma }_k^{†}{\Sigma }_j+{\Sigma }_j\rho {\Sigma }_k^{†}-{\Sigma }_k\rho {\Sigma }_j^{†}\right]$$

$$-\sum _{k\ne j}{R}_k\left[{\Sigma }_j^{†}{\Sigma }_k\rho +\rho {\Sigma }_k^{†}{\Sigma }_j-{\Sigma }_j\rho {\Sigma }_k^{†}-{\Sigma }_k\rho {\Sigma }_j^{†}\right].$$

$${\Sigma }_j\equiv g_j{\sigma }_j,$$

$${\gamma }_j$$

$$=2|g_j{|}^2{\int }_0^{\Omega }J(\omega )\left[{\int }_0^{\infty }\cos ([\omega -{\omega }_j]s)ds\right]d\omega$$

$$=2|g_j{|}^2{\int }_0^{\Omega }J(\omega )\pi \delta (\omega -{\omega }_j)d\omega$$

$$=2\pi |g_j{|}^{2}J({\omega }_j)$$

$$\equiv \gamma ({\omega }_j)$$

$${\Delta }_j$$

Full text

I Introduction

Weakly damped open systems are important across a wide range of areas in both physics and chemistry, from quantum thermodynamics Vinjanampathy and Anders (2016); Horowitz and Jacobs (2015); Shizume (1995) to the control of chemical reactions Brif et al. (2010), to quantum technologies MacFarlane et al. (2003); Taylor and Bowen (2016); Degen et al. (2017); Dunjko and Briegel (2017); Braun et al. (2018); Crosson and Harrow (2016); Farhi et al. (2011); Reiserer and Rempe (2015). 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 modelling these systems as they avoid the computationally expensive, and often prohibitive, task of simulating the thermal environment Prior et al. (2010); Chin et al. (2010). 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 Stannigel et al. (2012); Tomadin et al. (2012); Chang et al. (2012); Cook et al. (2018); Goerz and Jacobs (2018); Vuglar et al. (2018), adiabatic computation Albash et al. (2012); Amin et al. (2008), super and sub-radiance Dantan et al. (2006); Clark et al. (2003); Gross and Haroche (1982); Le Kien et al. (2005), and “weak lasing” Aleiner et al. (2012); Eastham et al. (2016), 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 Breuer and Petruccione (2007); Jacobs (2014) (for ease of reference we present these master equations in the supplement). 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 Majenz et al. (2013); Santra et al. (2017). 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 Redfield (1957); Bloch (1957). This equation has long been the subject of debate because it is not guaranteed to preserve the positivity of the density matrix Dümcke and Spohn (1979); Gardiner and Zoller (2010). In some subfields (e.g. photo-chemistry Zhang et al. (1998); Yang and Fleming (2002)), 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 Whitney (2008); Jeske et al. (2015); Eastham et al. (2016). In particular, Eastham et al. Eastham et al. (2016) 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. In Jeske et al. (2015) 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 Majenz et al. (2013); Santra et al. (2017). 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 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 break down, but so do all Markovian master equations. We thus show that the 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 the B-R equation. The Linblad 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 Wiseman and Milburn (1993); Molmer et al. (1993); Diósi (1986); Jacobs (2010), and provides a formulation of the action of a thermal bath as a continuous measurement on the system. This 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 Wubs et al. (2006); Wittig (2005) 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:

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 flatness of the spectral density is also a necessary condition for Markovianity 2. 2.

It resolves the controversy regarding the Bloch-Redfield master equation: 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 the Bloch-Redfield equation is, in general, no longer valid, confirming the conjecture in Eastham et al. (2016)).

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.

In the next section 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. In Section III we use the information obtained in Section II to derive a master equation for all weakly damped systems and all temperatures given a constraint on the derivative of the spectral density. In Section IV we further confirm the accuracy of the master equation with numerical simulations. In Section VI we use numerical simulations to show that the master equation is able to describe time-dependent systems whose levels cross. Section Discussion concludes with a discussion of some open questions.

II 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. (We give the details of the SID method that we use in Appendix A.) Since 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 their complex 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 $\mathbf{v}(t)$, then SID provides us with a matrix $M(\tau)$ where $\mathbf{v}(\tau)=M(\tau)\mathbf{v}(0)$ for a specified time $\tau$. The number of appreciable eigenvalues of $M$ is the effective size of the fictitious system.

Since it is only the ratios between the 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$ (details of the bath model are given below in Section III), fix the mean transition frequency $\bar{\omega}\equiv(\omega_{1}+\omega_{2})/2=3\pi\tilde{\nu}$, and choose the coupling constants $g_{1}$ and $g_{2}$ (defined in Eq.(10)) so as to give the decay rates $\gamma_{1}=0.1\tilde{\nu}$ and $\gamma_{2}=0.05\tilde{\nu}$ 111Since we are simulating the Hamiltonian in Eq.(7), it is the coupling constants $g_{j}$, rather than the damping rates $\gamma_{j}$, that define the simulation, with $\gamma_{j}\propto|g_{j}|^{2}$. Because we have fixed $\bar{\omega}$ while changing the detuning, and since $\gamma_{j}\propto\omega_{j}$ for the Ohmic spectral density, when the detuning changes we must change $g_{j}$ to keep the damping rates fixed. Since we wish to examine the evolution when the detuning, $\Delta\omega\equiv\omega_{2}-\omega_{1}$, is not large compared to the damping rates, we simulate the evolution for the following four values of $\Delta\omega$: [math], $0.28\pi\gamma_{1}$, $2\pi\gamma_{1}$, and $4.8\pi\gamma_{1}$.

To perform the exact simulations we use the method detailed in Prior et al. (2010); Chin et al. (2010) which employs the matrix-product-state method of Vidal Vidal (2004, 2003). 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, $\Delta\omega$, 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\times 10^{-4}$ 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$ 222The 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 $\mathbf{x}$, the approximate model is $\dot{\mathbf{x}}=D\mathbf{x}$, with $D=\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 $\omega$, is given by

[TABLE]

Here

[TABLE]

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

[TABLE]

for an arbitrary operator $c$, and $H_{\mbox{\scriptsize L}}$ is the Lamb shift Hamiltonian, given by

[TABLE]

with $\sigma_{1}=|0\rangle\langle 1|$, $\sigma_{2}=|0\rangle\langle 2|$. The phase $\phi$ is determined by the phases of the interactions between the transitions and the bath (see Eq.(10)). We also note that $H_{\mbox{\scriptsize L}}$ can be factored as $H_{L}=-\hbar f(\omega)D^{\dagger}D$ in which $D=\sqrt{\Delta_{1}}\sigma_{1}+e^{i\phi}\sqrt{\Delta_{2}}\sigma_{2}$.

From the derivation of the degenerate master equation Breuer and Petruccione (2007) (see also Section III) we know that the decay rates $\gamma_{j}$ depend on the spectral density of the bath evaluated at their corresponding transition frequencies $\omega_{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 $\Delta\omega\not=0$, in which the Hamiltonian is now

[TABLE]

has exactly the same form as the degenerate master equation. That is, it can be written as Eq.(1) with $H_{\mbox{\scriptsize L}}=\sum_{j}\zeta_{j}\sigma_{j}^{\dagger}\sigma_{j}+(\xi\sigma_{1}^{\dagger}\sigma_{2}+\mbox{H.c.})$ and $\sigma=\sum_{j}\eta_{j}\sigma_{j}$ for some set of $\{\zeta_{j},\xi,\eta_{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 $\sigma_{1}\sigma_{2}$. A simple guess for the parameters $\zeta_{j}$, $\xi$, and $\eta_{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 $\Delta_{2}(\omega_{0})$ and $\gamma_{2}(\omega_{0})$ by new values implied by the new value of $\omega_{2}$, namely $\Delta_{2}(\omega_{0}+\Delta\omega)$ and $\gamma_{2}(\omega_{0}+\Delta\omega)$.

We find that this trial master equation does indeed match the model backed out using SID for all three values of $\Delta\omega$. We compare further the evolution predicted by this master equation to the exact evolution in Fig. 2, for a range of values of $\Delta\omega$. For these simulations we use $\omega_{0}=10\pi$ and $\gamma_{1}(\omega_{0})=2\gamma_{2}(\omega_{0})=0.1$. 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^{-3}$ for all values of the detuning shown.

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 $\Delta\omega$ that might be considered well into the non-degenerate regime ($\Delta\omega=100\gamma$). 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 $\Delta\omega\rightarrow\infty$.

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.(24) 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.

III Derivation of the master equation

Having shown numerically that there is a Markovian master equation describing 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 system and the bath, in which the latter consists of a continuum of independent harmonic oscillators, is given by

[TABLE]

Here $H_{\mbox{\scriptsize sys}}$ is the Hamiltonian of the system which has a discrete set of energy levels. The operators $b(\omega)$ are the annihilation operators for a continuum of harmonic oscillators indexed by their frequencies $\omega$ (equivalently the modes of a quantum field). The function $J(\omega)$ 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 $\Omega$ and is called the “cut-off” frequency. Instead of having sharp cut-off at frequency $\Omega$ one can instead arrange $J(\omega)$ 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 operator of the system that couples to the bath as $A+A^{\dagger}$. Here $A$ is defined as containing all the matrix elements of this Hermitian operator that transform higher energy levels to lower ones. If we denote the energy levels of the system by $|n\rangle$, thus writing

[TABLE]

and define transition (or “decay”) operators

[TABLE]

then $A$ can be written as

[TABLE]

where $g_{j}=|g_{j}|e^{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

[TABLE]

and the evolution of the transition operators $\sigma_{j}$ and bath operators $b(\omega)$ in the interaction picture is

[TABLE]

To derive Markovian master equations one applies a rotating-wave approximation (RWA) to the Hamiltonian above. This 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 the first RWA we move into the interaction picture. The terms in the interaction Hamiltonian that contain the products $Ab(\omega)$ and $A^{\dagger}b^{\dagger}(\omega)$ (the so-called “off-resonant” terms) become

[TABLE]

Since the minimum frequency at which each of the terms in the sum over $j$ oscillates is $\omega_{j}$, when the damping rates and Lamb shifts (to be derived below) are much less than all the $\omega_{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

[TABLE]

The regime of weak damping is defined as the regime in which we are close enough to the limit in which $\min_{j}(\omega_{j})/\max_{j}(\gamma_{j})\rightarrow\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_{\mbox{\scriptsize RWA}}$. For the details of these approximations we refer the reader to Wiseman and Milburn (2010); Breuer and Petruccione (2007); Gardiner and Zoller (2010)). The result is the following expression for the evolution of the density matrix of the system in the interaction picture:

[TABLE]

where

[TABLE]

is the interaction between the system and the bath that appears in $H_{\mbox{\scriptsize RWA}}$, in the interaction picture. The operator $\rho^{\mbox{\scriptsize I}}(t)$ is the density matrix of the system in the interaction picture, $\rho_{\mbox{\scriptsize B}}(0)$ is the initial density matrix of the bath, and $\mbox{Tr}_{\mbox{\scriptsize 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_{\mbox{\scriptsize R}}^{\mbox{\scriptsize I}}(t)$ and $H_{\mbox{\scriptsize R}}^{\mbox{\scriptsize I}}(s)$ into the expression above, one of the terms we obtain is

[TABLE]

where

[TABLE]

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$.

Now substituting $G(\omega)$ into $K$ and rearranging we obtain

[TABLE]

It is useful to define

[TABLE]

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

[TABLE]

Here, for compactness, we have defined

[TABLE]

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

We note that the decay rates $\gamma_{j}$ will be

[TABLE]

and the Lamb shifts will be

[TABLE]

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 Appendix B. So long as the spectral density does not decrease with $\omega$, and $\Omega>\omega_{j}$, the Lamb shift $\Delta_{j}$ can be expected to be greater than the damping rate $\gamma_{j}$ (this is true for the Ohmic bath, see below).

The master equation we have derived in Eq.(24), the Bloch-Redfield equation, 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 $\Gamma_{j}=\Gamma_{k}$. In this case the last two lines of Eq.(24) 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_{\mbox{\scriptsize L}}=\hbar D^{\dagger}D$ with $D=\sqrt{I_{j}}\Sigma_{j}+\sqrt{I_{k}}\Sigma_{k}$ and the Lindblad damping term is $\mathcal{D}[\sqrt{\gamma_{j}}\Sigma_{j}+\sqrt{\gamma_{k}}\Sigma_{k}]$.

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 $\omega$. (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 Collett and Gardiner (1984); Gardiner and Collett (1985).) If the spectrum is flat then $\Gamma_{j}=\Gamma_{k}$ for all values of $\omega_{j}$ and $\omega_{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 $\omega$ 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, $\omega_{j}$ and $\omega_{k}$, are different, every term in the last two lines of Eq.(24) will oscillate at the difference frequency $\Delta\omega_{jk}=|\omega_{j}-\omega_{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.(24), which is the non-degenerate master equation. How large does $\Delta\omega_{jk}$ need to be to eliminate the last two lines of Eq.(24)? 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 $\Delta_{j}$ and $\gamma_{j}$, respectively, and those on the last two lines are

[TABLE]

So the terms on the last two lines are eliminated when $\min_{jk}\Delta\omega_{jk}\gg\max_{lm}M_{lm}$ and $\min_{jk}\Delta\omega_{jk}\gg\max_{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.(24), we only need concern ourselves with the regime

[TABLE]

(since we have assumed $O_{jk}\geq M_{jk}$, the regime $\Delta\omega_{jk}\lesssim\min_{jk}O_{jk}$ automatically includes the regime $\Delta\omega_{jk}\lesssim\min_{jk}M_{jk}$). Outside of this regime, the last two lines will be eliminated by the oscillations at the detuning frequency $\Delta\omega_{jk}$.

We now recall that $\gamma_{j}$ and $\Delta_{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

[TABLE]

Now if $\Gamma_{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 $\Gamma_{j}\approx\Gamma_{k}$. More specifically, we consider systems for which the spectral density satisfies

[TABLE]

for which the more precise statement is

[TABLE]

since this implies that $\Gamma(\omega)$ also satisfies the same “slow variation” conditions. Under the condition in Eq.(33) we have

[TABLE]

while at the same time allowing

[TABLE]

The relation (36) 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.(24) because i) when $|\omega_{j}-\omega_{k}|\lesssim O_{jk}$ these replacements are well justified, and ii) when $|\omega_{j}-\omega_{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

[TABLE]

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

[TABLE]

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

[TABLE]

in which $\phi_{j}=\arg[g_{j}]$ as defined below Eq.(10). The term $H_{\mbox{\scriptsize L}}\equiv-\hbar D^{\dagger}D$ is the Lamb-shift Hamiltonian. If we expand it out we see that, so long as $\sigma_{j}^{\dagger}\sigma_{k}\not=0$, the upper levels of the transitions $j$ and $k$ are coupled together via the bath:

[TABLE]

where the phases are given by

[TABLE]

The decay rates $\gamma_{j}$ are given in Eq.(26) and are determined by the value of the spectral density $J(\omega)$ only at transition frequency $\omega_{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 $\Omega$, in which $J(\omega)\propto\omega$ (we choose to define $J(\omega)=\omega/\Omega^{2}$) the Lamb shifts are

[TABLE]

We see that so long as $\Omega$ is larger than $2\omega_{j}$ the Lamb shift is larger than the damping rate by at least a factor of $\Omega/\omega_{j}$. Recall that the Lamb shifts are required to be much less that the transition frequencies. Using the expression for $\Delta_{j}$ above, we have

[TABLE]

Thus to satisfy the condition $\Delta_{j}\ll\omega_{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 $\Delta_{k}/\Delta_{j}\approx 1$, when $\Delta\omega=|\omega_{j}-\omega_{k}|\lesssim\Delta_{j}$. Denoting $\Delta_{j}$ by $\Delta(\omega_{j})$, and writing $\omega_{k}=\omega_{j}+\Delta_{j}$, we have

[TABLE]

Since the master equation is already derived under the conditions that $\Omega\gg\omega_{j}$, and $\omega_{j}\gg\Delta_{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

[TABLE]

in which

[TABLE]

Here $T$ is the temperature of the bath and $k_{\mbox{\scriptsize 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(\omega)$ is multiplied by $n_{T}(\omega)$, so the new terms give a new integral for which we have to calculate the principle value. This integral is

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

and

[TABLE]

with

[TABLE]

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, the secular approximation, 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.(38), 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 $10^{3}\Delta_{\mbox{\scriptsize max}}$ in which $\Delta_{\mbox{\scriptsize max}}$ is the maximum Lamb shift among all the transitions. If we denote the transition operators in the first set by $\sigma^{(1)}_{j}$, with $j=1,\ldots,N_{1}$, and those in the second by $\sigma^{(2)}_{j}$, with $j=1,\ldots,N_{2}$, then the result of making the secular approximation on the master equation in Eq.(39) is

[TABLE]

with

[TABLE]

IV Accuracy of the Lindblad Master Equation: additional confirmation for the Ohmic bath

It is clear from the derivation in the previous section that the master equation we have obtained, given in Eqs.(39) and (52), will be valid so long as the variation of the Lamb shifts $\Delta(\omega_{j})$ and damping rates $\gamma(\omega_{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(\omega)$ 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}$, $\omega_{j}=\omega_{1}+\gamma_{2}/(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 $|\psi_{0}\rangle=(7i|1\rangle+3|2\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\times 10^{-3}$.

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}$, and $\omega_{2}=\omega_{1}+2\gamma$, $\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}\equiv\omega_{j}/\gamma_{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_{\mbox{\scriptsize 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 $|\psi_{0}\rangle=(i|1\rangle_{1}|0\rangle_{2}+|0\rangle_{1}|1\rangle_{2})/\sqrt{2}$, we show the evolution of the populations for both the master equation and the exact simulations of $H_{\mbox{\scriptsize RWA}}$ in Fig. 3b. The error in the evolution of the master equation over the duration shown in Fig. 3a is less than $5\times 10^{-3}$.

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 supplemental material 333The supplemental material can be found at xxxxxx, 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.

V Regime of Validity

We have shown that when the spectral density is sufficient flat the B-R equation 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, does the B-R equation deviate 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 (deviates from the exact evolution) at large damping. To explore how the slope of the spectral density affects the master equations we thus thus require a new spectral density, and we use the one depicted in Fig.4. This density, which we will denote by $J_{r}(\omega)$, 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}(\omega)/d\omega=1/\Omega^{2}$). In the middle segment, which includes the transition frequency(ies) of the system, this slope is increased by a factor of $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 Santra et al. (2017); Barnett and Radmore (2003). We can use system identification, already discussed in Section II, 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 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 $\omega_{1}=\omega_{0}-\gamma/2$, $\omega_{2}=\omega_{0}+\gamma/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}(\omega)$ 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^{n}$, with $n=0$ (the Ohmic spectrum) and $n=3,4,\ldots,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 $\log_{2}r$). Note that the evolution of the master equation is determined solely by the damping rate $\gamma(r)=2\pi g^{2}J_{r}(\omega_{0})$. 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 $\gamma(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 $\gamma$ specified by the master equation, increases just as rapidly as the latter.

We use system identification to determine the dimension of the dynamics that generates the exact evolution for each value of $r$. We define this dimension as the number of dynamical variables required to account for $0.999$ of the combined magnitude of the dynamical eigenvalues (see Section II). 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 master equations break down.

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 $|2\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 behaviour 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\times 10^{-12}$.

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}$.

VI 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 $\gamma_{j}(\omega_{j})$, $\Delta_{j}(\omega_{j})$, and $\sigma_{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 McCauley et al. (2020). 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 $|\psi_{-}\rangle\equiv(|1\rangle-|2\rangle)/\sqrt{2}$, which for $\Delta\omega=0$ will not decay since it is the (sub-radiant) dark state Gross and Haroche (1982); Stannigel et al. (2012). We then change the detuning with time as determined by following two functions:

[TABLE]

The function $F_{1}$ is chosen so that the detuning increases gradually, while $F_{2}$ involves rapid changes. In Fig.6 we compare the adiabatic version of the master equation with the exact evolution for the two cases. For $\Delta\omega(t)=F_{1}(t)$ the maximum error of the adiabatic master equation is $3.4\times 10^{-3}$, and for $\Delta\omega(t)=F_{2}(t)$ the maximum error is $1.5\times 10^{-2}$.

As our final example we consider a generalized version of the Landau-Zener transition Wittig (2005); Wubs et al. (2006), 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}$:

[TABLE]

The detuning, $\Delta\omega\equiv\omega_{1}-\omega_{2}$, starts at $-\Delta_{0}$ and increases as a sinusoid through zero to end at $\Delta_{0}$. 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\times 10^{-3}$.

Discussion

We have shown that a slow variation condition on the spectral density is necessary for Markovianity, 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 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. Eastham et al. (2016).

The new master equation unifies the existing Linblad 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.

Acknowledgements

The authors thank Even Jensen for graphic design of Fig.1. This research was supported in part by appointments to the Student and Postgraduate Research Participation Programs at the U.S Army Research Laboratory 444These programs are administered by the Oak Ridge Institute for Science and Education and the Oak Ridge Association of Universities through an interagency agreement between the U.S Department of Energy and USARL.. D.I.B. was supported by Air Force Office of Scientific Research Young Investigator Research Program (FA9550-16-1-0254).

Appendix A 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 Katayama (2005); van Overschee and de Moor (1996); Keesman (2011); Di Ruscio (1995).

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\leq M<J$. Let us denote a state of the full $J$ dimensional system by the vector $\mathbf{v}$, and the subset of $N$ variables in which we are interested by the $N$-dimensional vector $\mathbf{x}$. The map that gives the state of the total $J$-dimensional system at time $\tau$ given an initial state $\mathbf{v}(0)$ we will call $Z(\tau)$ so that $\mathbf{v}(\tau)=Z(\tau)\mathbf{v}(0)$. Defining $Z\equiv Z(\tau)$ we note that $Z(n\tau)=Z^{n}$. We also define the non-square projector $P$ that projects onto the $N$ variables so that $\mathbf{x}=P\mathbf{v}$.

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\tau)P^{\mbox{\scriptsize T}}$ that maps the $N$ variables at time [math] to their values at time $\tau$, for any time $\tau$.

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

[TABLE]

We now note that $H_{0}$ can be written as an outer product $H_{0}=CD^{\mbox{\scriptsize T}}$ of (non-square) matrices given by

[TABLE]

We can now determine $C$ and $D$ by doing a singular value decomposition of $H_{0}$ to give $H_{0}=USV^{\mbox{\scriptsize T}}$. Note that since the smaller dimension of the matrices $C$ and $D$ is smaller than that of $H_{0}$ we expect many of the columns of U and the rows of $V^{\mbox{\scriptsize T}}$ will be zero, as will many of the eigenvalues of $H_{0}$ 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

[TABLE]

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

[TABLE]

we can obtain the evolutionary map for the model, $M$, from $H_{1}$ using

[TABLE]

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.

Appendix B Principle value of an integral

In deriving the master equation in Section III we used the fact that

[TABLE]

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\in(b,c)$), is defined by

[TABLE]

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 $\int_{-b}^{c}(1/x)dx$. Assuming that $b$ and $c$ are positive and $c>b$ we have

[TABLE]

**Competing interests: ** The authors declare no competing financial interests.

Data availability statement: All numerical results (datasets) generated and analyzed in this work can be obtained from the corresponding author upon reasonable request.

Code availability statement: All computer code used to generate and/or analyze the numerical results presented in this work can be obtained from the corresponding author upon reasonable request.