Quantum adiabatic Markovian master equations

We develop from first principles Markovian master equations suited for studying the time evolution of a system evolving adiabatically while coupled weakly to a thermal bath. We derive two sets of equations in the adiabatic limit, one using the rotating wave (secular) approximation that results in a master equation in Lindblad form, the other without the rotating wave approximation but not in Lindblad form. The two equations make markedly different predictions depending on whether or not the Lamb shift is included. Our analysis keeps track of the various time and energy scales associated with the various approximations we make, and thus allows for a systematic inclusion of higher order corrections, in particular beyond the adiabatic limit. We use our formalism to study the evolution of an Ising spin chain in a transverse field and coupled to a thermal bosonic bath, for which we identify four distinct evolution phases. While we do not expect this to be a generic feature, in one of these phases dissipation acts to increase the fidelity of the system state relative to the adiabatic ground state.


I. INTRODUCTION
Recent developments in quantum information processing, in particular theoretical [1] and experimental [2] proposals for adiabatic quantum computation, have generated considerable renewed interest in the old topic of quantum adiabatic dynamics [3,4].While much work has been done on rigorous formulations of adiabatic approximations for closed quantum systems, e.g., Refs.[5][6][7][8][9][10][11][12], adiabatic evolution in open quantum systems is still a relatively unexplored topic.In this regard master equations governing the evolution of a quantum system with a time-dependent Hamiltonian coupled to an external environment or bath are an important tool.
The study of master equations with time-dependent Hamiltonians is certainly not a new topic, going back at least as far as the pioneering work of Davies & Spohn [13], who derived an exact master equation for an adiabatic but infinitesimally weak system-bath interaction.Other, more recent approaches have been attempted, but in each case certain limitations apply.For example, Childs et al. used the Lindblad equation with a time-independent Hamiltonian at each time step as an approximation to the adiabatic evolution of a system with a time-dependent Hamiltonian [14].Sarandy & Lidar derived a phenomenological adiabatic master equation, based on the idea that in the adiabatic limit the dynamical superoperator can be decomposed in terms of independently evolving Jordan blocks [15].This approach is phenomenological in the sense that it does not allow one to derive the various parameters appearing in the final master equation from underlying system and bath Hamiltonians.Approaches based on non-Hermitian effective Hamiltonians, e.g., Refs.[16][17][18], are necessarily also phenomenological.A rigorous phenomenological master equation derivation was given by Joye [19].Oreshkov & Calsamiglia connected open system adiabaticity to the theory of noiseless subsystems [20].Thunström et al. derived a master equation from first principles in the physically reasonable joint limit of slow change and weak open system disturbances, but did not elucidate the relative time and energy scales involved in their approximations [21].Various authors provided derivations for slow periodic Hamiltonians [22][23][24].Such derivations are valuable and can be made rigorous, but the assumption of periodicity can be excessive, especially in the context of adiabatic quantum computation.Bounds on the validity of the adiabatic approximation for open systems, but without master equations, were presented in Ref. [25] (see also Ref. [10]).Various authors derived or studied adiabatic master equations limited to the case of a single qubit, where detailed physical considerations are possible [26][27][28].
Our goal in this work is to derive master equations for adiabatic open system dynamics from first principles, while keeping track of all physical approximations, time-and energy-scales.In this manner we hope to fill a gap in the earlier literature on this topic, and to provide tools allowing for detailed comparisons with experiments satisfying the explicit assumptions behind our approximations.In particular, we derive several Markovian master equations, distinguished by different levels of adiabatic perturbation theory.When we add the rotating wave approximation (sometimes called the secular approximation) we arrive at master equations in Lindblad form [29], for which positivity of the state is guaranteed at all times [30].Our formalism allows for the calculation of non-adiabatic corrections, which we also discuss.We apply our master equations to numerically study the evolution of a transverse field Ising chain coupled to a thermal bath, and discuss generic features of the evolution.
Our basic starting point is the observation that the Markovian and adiabatic limits are fundamentally compatible, in the sense that a Markovian bath is "fast", while an adiabatically evolving system is "slow".As long as the corresponding timescales are appropriately matched, it is possible to derive an adiabatic master equation which is internally consistent.Somewhat more technically, we observe that in the interaction picture with respect to the unperturbed system and bath evolutions, where the bath is evolving "fast" while the system is evolving "slowly", and for sufficiently weak system-bath coupling, it is possible to consistently apply adiabatic perturbation theory to the time-dependent system operators.This insight allows us to derive our adiabatic Markovian master equations.Our work is conceptually similar in its starting point to that of Amin et al. [31] (see also the Supplementary Information of Ref. [2]), but is significantly more general.Some of our results are also conceptually similar to those found by Kamleitner & Shnirman [24] by de Vega et al. [32], and by Yung [33], but are again more general.The master equations we derive are a natural generalization of standard time-independent Hamiltonian results [30,34], and our master equations reduce to the standard results after freezing the time dependence of the system Hamiltonian.
The structure of this paper is the following.We set the stage in Section II, where we define the system-bath model, perform the Born-Markov approximation, and introduce the bath correlation functions and spectral-density.We pay special attention to constraints imposed (via the KMS condition) by baths in thermal equilibrium, and single out the case of the Ohmic oscillator bath.We interrupt the formal development in Section III, where we provide a summary of all the time-and energy-scales appearing in our various approximations, and the inequalities they must mutually satisfy.We then proceed to derive our adiabatic master equations in Section IV.We proceed in several steps.First, in subsection IV A we find the adiabatic limit of the time-dependent system operators.Next, in subsection IV B we find two pairs of master equations in the adiabatic limit, one pair in the interaction picture, the other in the Schrödinger picture.The master equations within a given picture are distinguished by whether we apply the adiabatic approximation once or twice.Next, in subsection IV C, we introduce the rotating wave approximation, which allows us to reduce the Schrödinger picture master equations into Lindblad form.We conclude the formal development by discussing non-adiabatic corrections to our master equations in subsection IV D. We move on to a detailed numerical study of an example in Section V, of a ferromagnetic chain in a transverse field, coupled to an Ohmic oscillator bath at finite temperature.We apply two of our Schrödinger picture master equations, with and without the rotating wave approximation, and show that without inclusion of the Lamb shift term they yield similar predictions.When the Lamb shift is included, however, substantial differences emerge.We discern four distinct phases in the evolution from the transverse field to the ferromagnetic Ising model, which we discuss and analyze.Concluding remarks are presented in Section VI.The paper is supplemented with detailed appendixes where many of the technical details of the derivations and required background are provided, both for ease of flow of the general presentation and for completeness.

A. Model
We consider a general system-bath Hamiltonian where H S (t) is the time-dependent, adiabatic system Hamiltonian, H B is the bath Hamiltonian, and H I is the interaction Hamiltonian.Without loss of generality the interaction Hamiltonian can be written in the form: where the operators A α and B α are Hermitian and dimensionless with unit operator norm, and g is the (weak) system-bath coupling.The joint system-bath density matrix ρ(t) satisfies the Schrödinger equation ρ = −i[H, ρ(t)], where we have assumed units such that = 1.
) denote the free system and bath unitary evolution operators (T + denotes time ordering), and define denote the joint Schrödinger picture system-bath unitary evolution operator.We transform to the interaction picture with respect to H S (t) and H B , by defining Ũ (t, 0) = U † 0 (t, 0)U (t, 0), which, together with the interaction picture density matrix ρ We restrict the use of tilde variables to refer to variables in the interaction picture.The time-dependent interaction picture Hamiltonian HI (t) is related to its Schrödinger picture counterpart via: where we have defined the time-dependent system and bath operators: We are interested in deriving a master equation for the system-only state, where in the second line we used the fact that U B acts only on the bath.

B. Born-Markov approximation
Writing the formal solution of Eq. (3b) as and substituting this solution back into Eq.(3b), we obtain, after tracing over the bath degrees of freedom, the equation of motion for the system density matrix ρS (t) = Tr B ρ(t) We make the standard Born approximation assumption, that we can decompose the density matrix as where χ(t), which expresses correlations between the system and the bath, is small in an appropriate sense and can hence be neglected from now on [30,35,36].Thus the equation of motion reduces to: where we defined the two-point correlation functions: In Eq. ( 9), we have assumed for simplicity that B α 0 = Tr [B α (t)ρ B (0)] = 0, so that the inhomogeneous term in Eq. ( 8) vanishes.Since we assumed that the bath state ρ B is stationary, the correlation function is homogenous in time: For notational simplicity we will denote B αβ (τ, 0) by B αβ (τ ) when this does not lead to confusion.Let us denote the time scale over which the two-point correlations of the bath decay by τ B , e.g., |B αβ (τ )| ∼ exp(−τ /τ B ).More precisely, we shall require throughout that As we show in Appendix B, if τ B 1/g, then we can apply the Markov approximation to each of the four summands in Eq. ( 9), i.e., replace ρS (t − τ ) by ρS (t): where (. . . ) refers to the products of A α and A β operators in Eq. ( 9), and where we have also assumed that t τ B and the integrand vanishes sufficiently fast for τ τ B , so that the upper limit can be taken to infinity.Note that by Eq. ( 12) the integral on the RHS of Eq. ( 13) is of O(τ B ), so that the relative magnitude of the two terms is O[(gτ B ) 2 ].An explicit derivation of the upper bound on the error due to this approximation can be found in Appendix B. The resulting Markovian equation cannot resolve the dynamics over a time scale shorter than τ B .

C. Correlation functions, the KMS condition, and the spectral-density matrix
In computing the terms appearing in Eq. ( 13) we are faced with integrals of the form Our goal is to express these integrals in terms of the spectral-density matrix the standard quantity in master equations.It is convenient to replace the one-sided Fourier transform in the spectral-density matrix by a complete Fourier transform.Thus we split it into a sum of Hermitian matrices, where we show in appendix C that γ and S are given by If we assume not only that the bath state is stationary, but that it is also in thermal equilibrium at inverse temperature β, i.e., ρ B = e −βH B /Z, then it follows that the correlation function satisfies the KMS (Kubo-Martin-Schwinger) condition [30] B If in addition the correlation function is analytic in the strip between τ = −iβ and τ = 0, then it follows that the Fourier transform of the bath correlation function satisfies the detailed balance condition For a proof see Appendix D.
It is important to note that the KMS detailed balance condition imposes a restriction on the decay of the correlation function.To see this, note first that , where we used Eq.(11).Now assume that Eq. ( 12) would have to hold for all n.It would follow that Thus all derivatives of γ ab (ω) would have to be finite at ω = 0. On the other hand, it follows from the KMS condition (18) that so that in the limit as ω approaches zero from below or above, This shows that in principle γ aa (ω) may be discontinuous at ω = 0. Indeed, the continuity condition γ aa (0 − ) = γ aa (0 + ) implies, from the KMS condition recast as Eq. ( 21), that 2γ aa (0) = βγ aa (0) .
On the other hand, Eq. ( 22) tells us that continuity and a lack of divergence at ω = 0 require γ aa (0) to satisfy a condition which prohibits it from being arbitrary.This is indeed the case for the Ohmic oscillator bath discussed in subsection V A, which satisfies Eq. ( 22) with finite γ aa (0).For this case the bath correlation function is exponentially decaying assuming the oscillator bath has an infinite cutoff.However, as we show in subsection V B the Ohmic bath transitions from exponential decay to a power-law tail for any finite value of the cutoff, at some finite transition time τ tr .In this case we find |B αβ (τ )| ∼ (τ /τ M ) −2 , where τ M is a time-scale associated with the onset of non-Markovian effects, and hence we have to relax Eq. (12), and, replace it with

III. TIMESCALES
In this subsection we summarize, for convenience, the relations between the different timescales which shall arise in our derivation.The total evolution time is denoted t f and the minimum ground state energy gap of H S is ∆, i.e., where ε 0 (t) and ε 1 (t) are the ground and first excited state energies of H S (t).We shall arrive at master equations of the following general form: Here L uni = −i[H S (t)+H LS , •] is the unitary evolution superoperator including the Lamb shift correction, L ad diss is the dissipative superoperator in the fully adiabatic limit, and L non-ad diss is the dissipative superoperator with leading order non-adiabatic corrections.Let where s ≡ t/t f is the dimensionless time.To ensure that L uni generates adiabatic evolution to leading order we shall require the standard adiabatic condition In order to derive Eq. ( 25), the three superoperators are ordered in perturbation theory, in the sense that where the norm could be chosen as the supoperator norm, i.e., the largest singular value (see Appendix A).We may then assume that Combining the O(τ B ) due to the integral on the RHS of Eq. ( 13) (as already remarked there) with the g 2 prefactor from Eq. ( 9), we have To ensure the first inequality in Eq. ( 28) we thus require The non-adiabatic correction is of order h ∆ 2 t f , and when it appears in L non-ad diss it is multiplied by the same factor as L ad diss , i.e., we have To ensure the second inequality in Eq. ( 28) thus amounts to the adiabatic condition, Eq. (27).All this is added to the condition for the validity of the Markovian approximation, mentioned in the context of Eq. ( 13), and justified rigorously in Appendix B.
There is one additional, independent time scale we have to concern ourselves with.This is the time scale associated with changes in the instantaneous energy eigenbasis relative to τ B .If we require the change in the basis to be small on the time scale of the bath τ B , we must have: We justify this claim in subsection IV A and Appendix F. Note that the adiabatic condition, Eq. ( 27), implies hτ B ∆t f ∆τ B , while Eq. ( 34) can be written as Putting this together thus yields Our other inequalities [Eqs.(31), (33)] can be summarized as IV. DERIVATION OF ADIABATIC MASTER EQUATIONS

A. Adiabatic limit of the time-dependent system operators
Let us first work in the strict adiabatic limit.We will discuss non-adiabatic corrections in subsection IV D. We denote the instantaneous eigenbasis of H S (t) by {|ε a (t) }, with corresponding real eigenvalues (energies) {ε a (t)}, i.e., H S (t) |ε a (t) = ε a (t) |ε a (t) , and Bohr frequencies ω ba (t) ≡ ε b (t) − ε a (t).As shown in Appendix E 1 we can then write the system time evolution operator as: where U ad S (t, t ) is the "ideal" adiabatic evolution operator.It represents a transformation of instantaneous eigenstate |ε a (t ) into the later eigenstate |ε a (t) , along with a phase where To achieve our goal of arriving at a master equation expressed in terms of the spectral-density matrix, our basic strategy is to replace the system operator ) with an appropriate adiabatic approximation, which will allow us to take this operator outside of the integral.To see how, note first that We now make two approximations: first, as per Eq.(37a) we replace U S (t, 0) by U ad S (t, 0); second, we replace U † S (t, t − τ ) by e iτ H S (t) , an approximation justified by the appearance of the short-lived bath correlation function B αβ (τ ) inside the integrals we are concerned with.Thus, we write and find the bound on the error due to dropping Θ(t, τ ) in Appendix F. Let Neglecting the operator-valued correction term Θ(t, τ ) entirely, we have, upon substituting Eq. (37b) and using where the approximation in is the smallness parameter of the adiabatic approximation, which we have already assumed to be small.The second term, 34)], is new, and its smallness is associated with changes in the instantaneous energy eigenbasis relative to τ B .We are interested in the regime where both terms are small.Thus, to the same level of approximation where and Γ αβ (ω ba (t)) is the spectral-density matrix defined in Eq. ( 14).Similarly, we have for the other term B. Master equations in the adiabatic limit We are now ready to put everything together.Starting from the Born-Markov master equation constructed from Eqs. ( 9) and ( 13), and using the approximations (43) and (45), we arrive at the following one-sided adiabatic interaction picture master equation: Since we used an adiabatic approximation for A β (t − τ ), it makes sense to do the same for A α (t), i.e., to replace the latter with U ad † S (t, 0)A α U ad S (t, 0).If this is done, we obtain the double-sided adiabatic interaction picture master equation It is convenient to transform back into the Schrödinger picture.Using ρS (t) = U † S (t, 0)ρ S (t)U S (t, 0) [Eq.( 6)] implies that ).Hence, using Eq. ( 46), we find the one-sided adiabatic Schrödinger picture master equation where This form of the master equation has not appeared in previous studies of adiabatic master equations.If we again use the adiabatic approximation for U S (t, 0), i.e., replace U S (t, 0)Π ab (0)U † S (t, 0) by U ad † S (t, 0)Π ab (0)U ad S (t, 0), we obtain the doublesided adiabatic Schrödinger picture master equation where Comparing to Eq. ( 25), the first term in Eqs. ( 48) and ( 50) is L uni (t), while the second is L ad diss (t).The master equations we have found so far are not in Lindblad form, and hence do not guarantee the preservation of positivity of ρ S .We thus introduce an additional approximation, which will transform the master equations into completely positive form.

C. Master equation in the adiabatic limit with rotating wave approximation: Lindblad form
In order to arrive at a master equation in Lindblad form, we can perform a secular, or rotating wave approximation (RWA).To do so, let us revisit the t → ∞ limit taken in the Markov approximation in Eq. ( 13).Supposing we do not take this limit quite yet, we can follow the same arguments leading to Eq. (42c), which we rewrite, along with the adiabatic approximation A α (t) ≈ U ad † S (t, 0)A α U ad S (t, 0).This yields We note that µ dc (t, 0) + µ ba (t, 0) = One can now make the argument that when the t → ∞ limit is taken, terms for which the integrand vanishes will dominate, thus enforcing the "energy conservation" condition ω ba = −ω dc .This is a similar rotating wave approximation as made in the standard time-independent treatment, although here, the approximation of phase cancellation is made along the entire time evolution of the instantaneous energy eigenstates.Clearly, in light of the appearance of other terms involving t in Eq. (52b), this argument is far from rigorous.Nevertheless, we proceed from Eq. (52b) to write, in the t → ∞ limit, where we have defined a new index ω such that: Note that the set of {ω}'s involved in the sum ω is evolving in time since it corresponds to differences of the instantaneous energy eigenvalues, but we suppress the time dependence for notational brevity.We show in Appendix G how, by performing a transformation back to the Schrödinger picture, along with a double-sided adiabatic approximation, we arrive from Eq. ( 53) at the Schrödinger picture adiabatic master equation in Lindblad form: where the Hermitian Lamb shift term is and we have defined Since the bath correlations are of positive type, it follows from Bochner's theorem that the matrix γ-the Fourier transform of the bath correlation functions-is also positive [30].Therefore, this Lindblad form for our master equation guarantees the positivity of the density matrix.We emphasize that Eqs. ( 48), (50), and (55) all generalize both the standard Redfield and Lindblad time-independent Hamiltonian results to the case of a time-dependent Hamiltonian in the adiabatic limit. 1 The time-independent result can easily be recovered by simply freezing the time dependence of the Hamiltonian, energy eigenvalues and eigenvectors.To see this explicitly, let us restrict ourselves to the Lindblad case.
The energy eigenvalues and eigenvectors are time independent in this case, so we can replace our time dependent Lindblad operators with time independent ones, and the ω index no longer varies with time: The resulting equation becomes: which is the standard form for the time-independent Lindblad master equation.This should make evident the physical meaning of our derivation.We have systematically generalized the time-independent result such that the Lindblad operators now rotate with the (adiabatically) changing energy eigenstates, which makes them time-dependent.This is a non-trivial difference, as it encodes an important physical effect: the dissipation/decoherence of the system occurs in the instantaneous energy eigenbasis.Equations ( 50) and ( 55) are the two master equations we use for numerical simulations presented later in this paper.We note that Eq. ( 55) appears similar to the Markovian adiabatic master equation found in Ref. [24], but is more general and did not require the assumption of periodic driving.

D. Non-Adiabatic Corrections to the Master Equations
So far, we assumed the adiabatic limit of evolution of the system [see Eq. (37a)].The adiabatic perturbation theory we review in Appendix E 1 allows us to compute systematic non-adiabatic corrections to the master equations we have derived.This perturbation theory is essentially an expansion in powers of 1/t f , and we rederive in Appendix E 1 the well known result [5] that to first order where Thus, to derive the lowest order non-adiabatic corrections to our master equations is a matter of repeating our calculations of subsections IV A and IV C with U ad S (t, t ) replaced everywhere by the first order term U ad S (t, t )V 1 (t, t ), and adding the result to the zeroth order master equations we have already derived.Rather than actually computing these corrections, let us estimate when they are important.
The condition under which the zeroth order adiabatic approximation is accurate is Eq. ( 27), which is now replaced by i.e., with the replaced by a mere .However, we would still like to perform the approximation of Eq. ( 40), in the sense that which is what allows the use of e iτ H S (t) in this approximation (as shown in Appendix F).Taken together, these two conditions are weaker than Eq. ( 35), which we can rewrite as ).
Recall that to ensure the validity of the Markov approximation and L uni L ad diss , we also had to demand the inequalities in Eqs. ( 31), (33); these can be now summarized as max g, g 2 ∆ 1 τ B , to be compared with Eq. (63).

A. The Model
We proceed to use our master equations to study the evolution of the Ising Hamiltonian with transverse field where the functions A(t) and B(t) are shown in Fig. 1, and were chosen for concreteness to describe the interpolation between the transverse field and Ising term in the D-Wave Rainier chip [2].This is a system which begins with the transverse magnetic field H X S turned on while the Ising term H Z S is turned off, and then slowly switches between the two.We couple this spin-system to a bath of harmonic oscillators, with bath and interaction Hamiltonian where b † k and b k are, respectively, raising and lowering operators for the kth oscillator with natural frequency ω k , and g j k is the corresponding coupling strength to spin j.This is the standard pure dephasing spin-boson model [37], except that our system is time-dependent.The resulting form for our operator L [Eq. ( 51)] is Recall that our analysis assumed that the bath is in thermal equilibrium at inverse temperature β = 1/(k B T ), and hence is described by a thermal Gibbs state ρ B = exp (−βH B ) /Z.We show in Appendix H that this yields where only one of the Heaviside functions is non-zero at ω = 0. To complete the model specification, we assume an Ohmic bath spectral function where ω c is a high-frequency cut-off and η is a positive constant with dimensions of time squared.It is often stated that the terms associated with the Lamb shift H LS [Eq.( 56)], i.e., Eq. (67b), can be neglected, since the relative order of S and H S is g 2 τ B /∆, and indeed we have assumed g 2 τ B /∆ 1 [Eq.(31)].However, this argument is misleading for two reasons.First, S can be divergent, as is easy to see in the limit ω c = ∞ for the Ohmic model (68), where for ω max ba ω ba (t), the integrand tends to a constant.Second, S should be compared to γ, as both are of the same order g 2 τ B , and both result in changes to the system relative to its unperturbed state.Indeed, in the interaction picture with respect to H S + H LS [recall Eq. ( 55)], the overall transition rates between states with quantum numbers a and b will depend on the dressed (i.e., shifted) energy gaps ω ba + ω LS ba .The importance of this Lamb shift effect was also stressed by de Vega et al. [32].We analyze the Lamb shift effect in subsection V C.
Finally, we note that although the harmonic oscillators bath with linear coupling to the system provides a γ matrix that satisfies the KMS condition, it is important to note that this model has infrared singularities that destroy the ground state of the total system [38].The KMS condition assumes a stable ground state and stable thermal states, which our underlying spin-boson model violates.However, for the purposes of our work, a γ matrix that satisfies the KMS condition will suffice without too much concern about how it is derived.

B. Correlation function analysis
In light of the subtleties alluded to in subsection II C associated with satisfying the KMS condition, we analyze the different timescales determining the behavior of the Ohmic correlation function in this subsection.Removing the ω dependence from the g's in Eq. (67a) and substituting J(ω) from Eq. (68), it is possible to compute the bath correlation function analytically for the resulting by inverse Fourier transform of Eq. (16a).The result is where ψ (m) is the mth Polygamma function (see Appendix I for the derivation).We first assume and consider an expansion in large βω c : followed by an expansion in large τ /β: This expansion reveals the two independent time-scales that are relevant for us.First, there is the time scale τ B associated with the exponential decay (corresponding to the true Markovian bath), given by: then the time scale associated with non-Markovian corrections (the power law tail): For sufficiently large ω c , these two time scales capture the two behaviors found in B ab (τ ), as illustrated in Fig. 2(a).
The transition between the exponential decay and the power law tail occurs at a time τ tr given by 4π 2 e −τtr/τ B = ( τ M τtr ) 2 , or equivalently by exp(θ) This transcendental equation has a formal solution in terms of the Lambert-W function [39], i.e., the inverse function of f (W ) = W e W , as can be seen by changing variables to y = −θ/n, and rewriting θ n e −θ = a ≡ 2/(βω c ) as ye y = −a 1/n /n, whose solution is θ = −ny = −nW (− 1 n a 1/n ).However, for our purposes the following observations will suffice.We seek a Markovian-like solution, where τ tr is large compared to the thermal timescale set by β, i.e., we are interested in the regime where θ 1.In this case we can neglect θ 2 compared to e θ , and approximate the solution to Eq. ( 76) by θ ∼ ln(βω c /2).Thus This agrees with the first term of the asymptotic expansion W (x) = ln(x) − ln ln(x) + ln ln(x) ln(x) + • • • , which is accurate for x 3 [39].
When βω c 1 is not strictly satisfied, the exponential regime is less pronounced, and τ B is corrected by powers of ω c .By dimensional analysis, the corrections must be of the form: where c n are constants of order one that must be fitted [see Fig. 2(b)].
The implications of this cutoff-induced transition for our perturbation theory inequalities are explored in Appendix F, where we show that a sufficient condition for the theory to hold is which can be interpreted as saying that the cutoff should be the largest energy scale.Equation (79) joins the list of inequalities given in Section III as an additional special condition that applies in the Ohmic case, along with βω c 1.

C. Numerical Results
For concreteness, we take g i ω = g, ω c = 8πGHz, and T = 20mK ≈ 2.6GHz (in units such that = 1; this is the operating temperature of the D-Wave Rainier chip [2]), corresponding to τ B = 0.06ns for the Ohmic model with infinite cutoff, Eq. ( 74), and τ B ≈ 0.07ns for ω c = 8π using Eq.(78).For this value of ω c the transition between the exponential and power law regimes is still sharp (see Fig 2(b)) and occurs at approximately τ tr = 0.33ns.For these values, we satisfy at least one of the cases from Eq. (79), including numerical prefactors: We focus on the N = 8 site ferromagnetic chain with parameters: where we pin the first spin in order to break the degeneracy in the ground state of the classical Ising Hamiltonian.The system is initialized in the Gibbs state: To help the numerics, we truncate our system to the lowest n = 18 levels (out of 256), rotating the density matrix into the instantaneous energy eigenbasis at each time step.The error associated with this is small as long as our evolution does not cause scattering into higher n states, as we have checked.The forward propagation algorithm used is an implicit second order Runge-Kutte method called TR-BDF2 with an adaptive time step [40,41].Figure 3 presents our results for the evolution of the system described in the previous subsection.We computed the overlap between the instantaneous ground state of H S (t) and the instantaneous density matrix predicted by our two master equations (50) (non-RWA) and (55) (RWA).Although our two master equations predict different numerical values for this overlap, the qualitative features of the evolution are the same.We observe a generic feature of four distinct regions of the evolution: a gapped phase (labeled "1" in figure 3), an excitation phase (labeled "2"), a relaxation phase (labeled "3"), and finally a frozen phase (labeled "4").We elaborate on these regions in the following subsection.Furthermore, we observe that the larger t f (therefore more adiabatic) evolution shows a smaller difference between the two master equations for the final fidelity.The results in Fig. 3 illustrate the importance of the Lamb shift.We find that while its effect is small in the RWA, its effect is relatively large in the non-RWA.
In order to study the importance of the relaxation phase, we can study the behavior of the fidelity at t = t f as we change the coupling strength g.We observe (see Fig. 4) that there is a rapid drop in fidelity from the closed system result as soon as the coupling to the thermal bath is turned on, but there is a subsequent steady increase in the fidelity as the coupling strength is further increased.This increase in fidelity is a direct consequence of the higher importance of the relaxation phase (made possible by the increasing coupling strength) in restoring the probability of being in the ground state.However, we also observe that there is a very pronounced difference between the behavior of the results from the two master equations as the coupling is increased.In the case of the Lindblad equation, the fidelity saturates, whereas for the NRWA equation, we see an increase in fidelity and a subsequent violation of positivity.These results bring to light the relative advantages and disadvantages of FIG. 3. Fidelity, or overlap of the system density matrix with the instantaneous ground state along the time evolution for t f = 10µs and ηg 2 /( 2 ) = 1.2 × 10 −4 /(2π).The solid blue curve was calculated using Eq. ( 50) (no RWA), while the dashed red curve was calculated using Eq. ( 55) (Lindblad form, after the RWA).Four phases are indicated: thermal (1), excitation (2), relaxation (3), and frozen (4).Panel (a) includes the Lamb shift terms, while (b) excludes them.
both master equations.For the Lindblad equation, positivity of the density matrix is guaranteed, but it clearly is not capturing the physics associated with the increasing importance of the thermal relaxation that the non-RWA equation captures.However, the non-RWA equation fails to preserve positivity of the density matrix so it is unable to reliably describe the system at higher coupling strength.Others have also noted that the RWA and NRWA can lead to physically different conclusions, e.g., in the context of Berry phases in cavity QED [42].

Phase 1 -the gapped phase
For times sufficiently close to the initial time, the ground state of H S (t) is the ground state of H X S , i.e., the state |0 ≡ ⊗ N j=1 |− j , where |± j = (|↓ j ± |↑ j )/ √ 2 with energy ε 0 (0) = −N A(0), and where |↓ j , |↑ j are the +1, −1 eigenstates of σ z j (computational basis states of the jth spin or qubit).The lowest lying energy states are then the N -fold degenerate states with a single flip of one of the spins in the x direction, i.e., |i ≡ ⊗ i−1 j=1 |− j |+ i ⊗ N j=i+1 |− j , with energy ε 1 (0) = −(N − 2)A(0).Therefore the gap between the ground state and the first excited states is: which is at least twice as large as our k B T ≈ 2.6GHz, almost until A(t c ) = B(t c ) ≈ 5GHz at t c ≈ 0.35t f (see Fig. 1).Noting that σ z i |± i = |∓ i , we can write the following relations in terms of the ground and first excited states: Noting that σ z i |j is a doubly excited state ∀j = i ∈ {1, . . ., N }, we truncate the problem to the ground and first excited states only, so that there are just three types of values of γ(ω) we need to consider: γ(0), γ(ω i0 ), γ(ω 0i ).Recalling the KMS condition, Eq. ( 18), we have which shows that upward transitions are exponentially suppressed relative to downward transitions, by a factor ranging between e −2A(0)/(k B T ) ≈ e −67.4/2.6 ≈ 7 × 10 −12 and e −2A(0.3tf )/2.6 ≈ e −3.8 ≈ 0.02.This explains why for early times (Phase 1) the system hardly deviates from the ground state, which in turn is very close to the thermal state (81).To make this argument more precise, denoting ρ 00 = 0|ρ|0 and ρ ii = i|ρ|i , we can write the effective (truncated to the ground and first excited states) Lindblad equation ( 55) as the simplified rate equations ρii ≈ γ ii (ω i0 ) −ρ ii + e −βωi0 ρ 00 , (85a) where we have assumed that the system is initially in the thermal state (81) and the gap is large (relative to k B T ).A derivation of Eq. ( 85) can be found in Appendix J.We compare our simulation results with the results from the above equations in Fig. 5 and find very good agreement for early times.

Phase 2 -the excitation phase
When A(t) becomes small enough such that β∆(t) ∼ O(1), then the KMS condition no longer suppresses excitations from the ground state to higher excited states (see Fig. 6).If we interpret the master equation as a set of rate equations for the matrix elements of ρ, we can identify the rate of scattering into the ith state from the jth state as being the term in the ρii equation with coefficient ρ jj .Therefore, we find that the rate of scattering from the ground state to excited states is given by whereas the relaxation rate is given by Therefore, as we emerge from the gapped phase, we have ρ 00 ρ ii , so scattering into excited states dominates over relaxation into the ground state.This explains the loss of fidelity in Phase 2.

Phase 3 -the relaxation phase
As the gap begins to grow again and the suppression factor shrinks (see Fig. 6), the KMS condition begins to suppress scattering into higher excited states while allowing relaxation to occur.In our model this causes a resurgence in the overlap with the ground state.Therefore, in this phase, the presence of the thermal bath can actually help to increase the fidelity, as was also observed in Ref. [26,32,43] The excitation and relaxation phases reveal the two competing processes for a successful adiabatic computation.If we spend too long in the excitation phase (or if the gap shrinks too fast relative to t f and the evolution is not adiabatic), the system loses almost all fidelity with the ground state, and the system would have to spend a very long time in the relaxation phase to recover some of that fidelity.
However, we stress that fidelity recovery will not be observed if the population becomes distributed over a large number of excited states in the excitation phase.This would happen, e.g., if when the gap closes there is an exponential number of states close to the ground state, such as in the quantum Ising chain with alternating sector interaction defects [44].To see this more explicitly in the context of our analysis, note from Eq. (85b) that if there is an exponential number N of ρ ii coupled to ρ 00 (i.e., N is an exponentially large fraction of the dimension of the system Hilbert space), then ρ 00 decreases exponentially.By Eq. ( 86) this means that all ρ ii become exponentially small, but not zero.In phase 3, by Eq. ( 87), the relaxation is suppressed as long as the gap is not very large, because the relaxation is proportional to ρ ii , which is exponentially small.This analysis presumes that the system-bath Hamiltonian has non-negligible coupling between the ground and excited state, i.e., that | 0|σ z β |i | > 0 in our model [see Eq. ( J6a)].This suggests another mechanism that can suppress relaxation: the ground state in phase 1 might have a large Hamming distance from the ground state in phase 3. Relaxation is then suppressed simply because the coupling is small.
We might expect that there exists an optimum t f for which the fidelity is maximized by the end of the relaxation phase.However, for the simple case of the spin-chain we considered, we did not observe such an optimum t f .The fidelity continues to grow (albeit slowly) for sufficiently high t f .This is illustrated in Fig. 7.

Phase 4 -the frozen phase
As the gap continues to grow, the relaxation phase ends [notice that the tail in Fig. 6(b) is longer than the actual relaxation phase] and the system's dynamics are frozen in the ground state.This is because H S becomes almost entirely diagonal in the σ z basis, and so the off-diagonal components of the L ab,i operators vanish (or become very small), i.e., A iab (t) = ε a (t)|σ z i |ε b (t) ∝ δ ab [Eq.( 66)] leaving only the diagonal ones.At this point, while off-diagonal elements may continue to decay, the system ground state is no longer affected by the bath; this is the onset of the frozen phase. .Time evolution of the system density matrix using the RWA equation with ηg 2 /( 2 ) = 0.4/(2π) for different t f 's.Solid blue curve is for t f = 10µs, dashed red curve is for t f = 100µs, and dot-dashed green curve is for t f = 1ms.

E. Thermal equilibration
An interesting question is whether the system reaches thermal equilibrium throughout the evolution.To answer this we computed the trace-norm distance [defined in Eq. (A1)] between the instantaneous system density matrix and the instantaneous Gibbs state ρ Gibbs (t) = exp[−βH S (t)]/Z, where Z = Tr [exp (−βH S (t))] is the partition function, for the Ising chain discussed above.The result is shown in Fig. 8.The distance is zero at t = 0 since the system is initialized in the Gibbs state, and then begins to grow slowly as the system transitions from the gapped phase to the excitation phase.Though not generic, the distance decreases as the gap shrinks while the excitation phase becomes the relaxation phase, and the system returns to a near Gibbs state where the gap is minimum (at t/t f ∼ 0.4).As the gap opens up again in the transition from the relaxation phase to the frozen phase, the distance begins to grow and continues to grow throughout the frozen phase.As such, the system is quite far from the Gibbs state at the final time.This feature is significant for adiabatic quantum computation, since it shows the potential for preparation of states which are biased away from thermal equilibrium towards preferential occupation of the ground state.However, we note that on sufficiently long timescales one should expect (from general thermodynamic arguments) terms proportional to σ x and σ y in the system-bath interaction, which we neglected in writing the interaction Hamiltonian H I in Eq. (65), to become important, and to disrupt the frozen phase, allowing the system to fully equilibrate into the Gibbs state.

VI. CONCLUSIONS
Using a bottom-up, first principles approach, we have developed a number of Markovian master equations to describe the adiabatic evolution of a system with a time-dependent Hamiltonian, coupled to a bath in thermal equilibrium.Our master equations systematically incorporate both time-dependent perturbation theory in the (weak) system-bath coupling g, and adiabatic perturbation theory in the inverse of the total evolution time t f .Since we have kept track of the various time-and energy-scales involved in our approximations, higher order corrections (starting at third order in g and second order in 1/t f ) can be incorpo-rated if desired, a problem we leave for a future publication.Using two of our master equations, we studied generic features of the adiabatic evolution of a spin chain in the presence of a transverse magnetic field, and coupled to a bosonic heat bath.We identified four phases in this evolution, including a phase where thermal relaxation aids the fidelity of the adiabatic evolution.We hope that this work will prove useful in guiding ongoing experiments on adiabatic quantum information processing, and will serve to inspire the development of increasingly more accurate adiabatic master equations, going beyond the Markovian limit. where and where the ellipsis denotes the three other summands appearing in Eq. ( 13).The Markov approximation is more accurate the smaller the error terms ∆ 1 and ∆ 2 .Consider first the ∆ 1 term.We shall show that it is of order g 2 τ 3 B .For simplicity we consider only the first of its four summands (the bounds for the other three are identical to that for the first): where we used the triangle inequality, We can now use Eq. ( 9) to upper-bound the time derivative: where the factor of 4 is due to the same number of summands appearing in Eq. ( 9).Substituting this bound back into Eq.(B3) we have where we used | and Eq. ( 12).This is to be compared to the term we use after the Markov approximation: Comparing Eqs.(B5) and (B6), we see that the relative error is Next consider the ∆ 2 term.We have The last integral converges provided |B αβ (τ )| ∼ 1/τ x with x > 1.This is typically the case.E.g., for the Ohmic spin-boson model discussed in subsection V B we show that for t > τ tr ∼ β ln(βω c ), the bath correlation function |B αβ (τ )| ∼ ηg 2 βωc 1 τ 2 .In this case, then, we can set t = τ tr and bound which tends to zero as the cutoff tends to infinity at a fixed finite temperature (even if the lower integration limit is replaced by a constant), as expected in the weak coupling limit assumed in our work.Note that the infinite temperature limit, where Eq. (B8) diverges, is incompatible with weak coupling, and requires the so-called singular coupling limit [22,30].
To prove that this implies the frequency domain condition, Eq. ( 18), we compute the Fourier transform: To perform this integral we replace it with a contour integral in the complex plane, C dτ e iωτ B b (−τ − iβ)B a (0) , with the which, together with Eq. (D3), proves Eq. (18).

Appendix E: Non-Adiabatic Corrections
For completeness we provide a brief review of pertinent aspects of adiabatic perturbation theory, and the derivation of the adiabatic condition relating the total evolution time to the ground state gap.See, e.g., Ref. [5] for additional details.

Adiabatic perturbation theory
To consider adiabatic corrections, we recall that U S satisfies where Define the "adiabatic intertwiner" W (t, t ): where the "intertwiner Hamiltonian" is and where P 0 (t) ≡ |ε 0 (t) ε 0 (t)| is the projection onto the ground state of H S (t).The result (E4b) is by no means obvious and is proven in subsection E 3.
To extract the geometric phase we define Now define V via the "adiabatic interaction picture" transformation: along with H S (t, t ) ≡ HS (t, t ) − HG (t, t ) , (E7c) K(t, t ) = W † (t, t )K(t)W (t, t ) = iW † (t, t )∂ t W (t, t ) . (E7d) Note that the time dependence of HS and H S (t) is entirely in the energy eigenvalue and not in the eigenstates.Then V obeys the Schrödinger equation: where When the evolution is nearly adiabatic Had S (t, t ) is a perturbation, so that we consider a solution of Eq. (E8) for V of the form: with the zeroth order solution associated with the purely adiabatic evolution, including the geometric phase: Differentiating with respect to t, we have where Plugging the Eq.(E10) expansion into Eq.(E8), we obtain, to first order in 1/t f :3 we can write the solution as: Therefore, the system evolution operator can be written as: where the correction term can be made appropriately dimensionless in a more careful second order analysis, and where

Adiabatic timescale analysis
Equation (E17) shows that the first order correction to the purely adiabatic evolution is given by Q.Thus the adiabatic timescale is found by ensuring that the matrix elements of Q are all small.Let us show that a sufficient condition for this is (27)].More rigorous analyses replace the condition with an inequality relating the same parameters to the fidelity between the ground state of H S (t f ) and the solution of the Schrödinger at t f , e.g., [6,[10][11][12].
Starting from Eq. (E18) and taking matrix elements we have where μ now involves a dimensionless time integration.Using the fact that e it f μba (s ,t /t f ) = i t f ω ab (s ) d ds e it f μba (s ,t /t f ) , we can integrate Eq. (E20) by parts as Continued integration by parts will yield higher powers of the dimensionless quantity εa(s)|∂τ H S (s)|ε b (s) . Thus a sufficient condition for the smallness of |Q ab (t, t )| for all a, b is that namely the condition in Eq. ( 27), where we assumed that the minimal Bohr frequency is the ground state gap, ω 1,0 .Of course, this argument is by no means rigorous, and in fact it has been the subject of considerable discussion, e.g., [50][51][52][53][54][55].For rigorous analyses see, e.g., Refs.[5][6][7][10][11][12].For our purposes Eq. ( 27) suffices.
Appendix F: Short Time Bound We wish to bound the error associated with neglecting Θ in Eq. ( 40), i.e., we wish to bound Using that the operator Û (τ ) = e −iτ H S (t) U † S (t, t − τ ) satisfies: we can write a formal solution for Û as: Therefore we can bound: where we used Eq.(E17) and the fact that supoperator norm between two unitaries is always upper bounded by 2 in the second line, and the standard adiabatic estimate to bound Q(t, 0) ∞ (recall subsection E 2).While h of Eq. ( 26) is expressed in terms of a matrix element, a more careful analysis (e.g., Ref. [10]) would replace this with an operator norm.Thus we shall make the plausible assumption that h ∼ t f max t ∈[t−τ,t] ∂ t H S (t ) ∞ , and, dropping the subdominant O(t −2 f ), we can write We can now bound the error term in Eq. (42b).Let X(t, τ ) ≡ e iτ H S (t) U ad S (t, 0), so that U S (t − τ, 0) = X(t, τ ) + Θ(t, τ ).We can then write The first term on the RHS of (F6a) is the approximation we have used in Eq. (42b).The terms in (F6b) and (F6c) can be bounded as follows, using Eq.(F5), the unitarity of X, the fact that ρS ∞ ≤ 1, and recalling that A α ∞ = 1.First we assume that Eq. ( 12) applies.Then: where in the last inequality we used the fact that if x ≤ 2 then [min(2, x)] 2 = x 2 ≤ 2x, and if x ≥ 2 then again [min(2, x)] 2 = 2 min(2, x) ≤ 2x, with x = h ∆ 2 t f + τ 2 h t f , in order to avoid having to extend Eq. ( 12) to higher values of n.In all, then, the approximation error in Eq. (42b Next we recall from the discussion in subsection II C that Eq. ( 12) must, in the case of a Markovian bath with a finite cutoff, be replaced by the weaker condition (23), reflecting fast decay up to τ tr , followed by power-law decay.In this case the terms in (F6b) can instead be bounded as follows: in place of (F7b).A similar modification can be computed for the term in (F6c).To compute the order of Starting from Eq. ( 53) and performing a transformation back to the Schrödinger picture, along with a double-sided adiabatic approximation, yields U S (t, 0) Using that: we can write: We can simplify our expression to e −iω k t g i k g j k + e iω k t−βω k g i k g j k . (H8) We can replace the sum with a integral: where f (ω) is a measure of the number of oscillators at frequency ω and J(ω) is the bath spectral function.Our final result is then: B i (t)B j (0) = ∞ 0 dω J(ω) 1 − e −βω e −iωt g i ω g j ω + e iωt−βω g i ω g j ω . (H10) where we have assumed that oscillators with the same ω k value interact with the i th spin with the same interaction strength, i.e.
Plugging our result into Eq.( 16), we find: where we have used that +|σ z |+ = −|σ z |− = 0, +|σ z |− = 1, and [H S , ρ S ] ≈ 0. Interpreting Eq. (J5) as a rate equation, we see that γ(+∆) is associated with relaxation from the higher energy state to the lower energy state, and γ(−∆(t)) is associated with the excitation from the lower energy state to the higher energy state.Note that at t = 0, we assume that the system is initialized in a thermal state such that ρ ++ /ρ −− = e −β∆(0) ≈ 10 −11 .Since the gap remains relatively large during the gapped phase, the change in the initial thermal state is minimal.

N qubits
For the N -qubit case, we can simply generalize our arguments for the single qubit case.For sufficiently small times the system is diagonal in the σ x basis, and the lowest lying energy states can be considered to be single spin flips in the x direction.Furthermore, the gap between the ground state and the first excited states is ∆(t) = 2A(t) [Eq.( 82)], which is large in our case.We can restrict ourselves to only these low lying energy states since higher excited states will have twice the gap to the ground state, and so their contribution will be further suppressed at short times.Now, using Eq. ( 83) and the KMS condition [Eq.( 84)] we find a similar set of equations as in the single qubit case: ρii ≈ αβ γ αβ (ω i0 ) 0|σ z β |i i|σ z α |0 −ρ ii + e −βωi0 ρ 00 = γ ii (ω i0 ) e −βωi0 ρ 00 − ρ ii (J6a) where we have again assumed that the initial state is the thermal state and the gap is large.

FIG. 4 .
FIG.4.Dependence of the fidelity at t f = 10µs on the coupling strength to the thermal bath is varied using our two master equations.The insets are closeups of the behavior near zero coupling.

FIG. 5 .
FIG.5.Early evolution of diagonal elements of the system density matrix in the thermal regime.Red dots are from using the Lindblad equation (55), while the solid blue curve is from the solution of Eq. (85).

FIG. 6 .
FIG.6.The behavior of the gap and the exponential suppression factor along the time evolution.The dashed line is kBT ≈ 2.6GHz.

FIG. 7
FIG.7.Time evolution of the system density matrix using the RWA equation with ηg 2 /( 2 ) = 0.4/(2π) for different t f 's.Solid blue curve is for t f = 10µs, dashed red curve is for t f = 100µs, and dot-dashed green curve is for t f = 1ms.