Quantum adiabatic Markovian master equations Academic research paper on "Physical sciences"

lopscience

¡opscience.iop.org

Home Search Collections Journals About Contact us My IOPscience Quantum adiabatic Markovian master equations

This content has been downloaded from IOPscience. Please scroll down to see the full text. View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 140.211.127.19

This content was downloaded on 20/05/2015 at 17:22

Please note that terms and conditions apply.

New Journal of Physics

The open access journal for physics

Quantum adiabatic Markovian master equations

Tameem Albash1,2,6, Sergio Boixo2,3,4, Daniel A Lidar2,4,5 and Paolo Zanardi12

1 Department of Physics and Astronomy, University of Southern California, Los Angeles, CA 90089, USA

2 Center for Quantum Information Science and Technology, University of Southern California, Los Angeles, CA 90089, USA

3 Information Sciences Institute, University of Southern California, Los Angeles, CA 90089, USA

4 Department of Electrical Engineering, University of Southern California, Los Angeles, CA 90089, USA

5 Department of Chemistry, University of Southern California, Los Angeles, CA 90089, USA

E-mail: albash@usc.edu

New Journal of Physics 14 (2012) 123016 (40pp)

Received 9 July 2012 Published 11 December 2012 Online at http://www.njp.org/

doi:10.1088/1367-2630/14/12/123016

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

6 Author to whom any correspondence should be addressed.

Content from this work may be used under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

New Journal of Physics 14 (2012) 123016 1367-2630/12/123016+40$33.00

© IOP Publishing Ltd and Deutsche Physikalische Gesellschaft

Contents

1. Introduction 2

2. Preliminaries 4

2.1. Model..............................................................................4

2.2. Born-Markov approximation......................................................5

2.3. Correlation functions, the Kubo-Martin-Schwinger condition, and the spectral-density matrix......................................................................6

3. Timescales 8

4. Derivation of adiabatic master equations 10

4.1. Adiabatic limit of the time-dependent system operators..........................10

4.2. Master equations in the adiabatic limit............................................11

4.3. Master equation in the adiabatic limit with rotating wave approximation: Lindbladform......................................................................12

4.4. Non-adiabatic corrections to the master equations................................14

5. An illustrative example: transverse field Ising chain coupled to a boson bath 14

5.1. The model..........................................................................14

5.2. Correlation function analysis ......................................................16

5.3. Numerical results..................................................................18

5.4. The four different phases..........................................................20

5.5. Thermal equilibration ..............................................................22

6. Conclusions 24 Acknowledgments 24 Appendix A. Norms and inequalities 24 Appendix B. Markov approximation bound 25 Appendix C. Properties of the spectral-density matrix Гар (w) 26 Appendix D. Proof of the Kubo-Martin-Schwinger condition 27 Appendix E. Non-adiabatic corrections 28 Appendix F. Short time bound 33 Appendix G. Derivation of the Schrodinger picture adiabatic master equation

in Lindblad form 35

Appendix H. Calculations for the spin-boson model 36

Appendix I. Derivation of the Ohmic bath correlation function 37

Appendix J. Derivation of the effective rate equations in the thermal phase 38

References 39

1. 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. [5-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 and 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 [14] 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. Sarandy and Lidar [15] 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. 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. [16-18], are necessarily also phenomenological. A rigorous phenomenological master equation derivation was given by Joye [19]. Oreshkov and Calsamiglia [20] connected open system adiabaticity to the theory of noiseless subsystems. Thunstrom et al [21] 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. Various authors provided derivations for slow periodic Hamiltonians [22-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 [25] (see also [10]). Various authors derived or studied adiabatic master equations limited to the case of a single qubit, where detailed physical considerations are possible [26-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 (RWA) (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 [2]), but is significantly more

general. Some of our results are also conceptually similar to those found by Kamleitner and 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, 35], 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 2, 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 Kubo-Martin-Schwinger (KMS) condition) by baths in thermal equilibrium, and single out the case of the Ohmic oscillator bath. We interrupt the formal development in section 3, 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 4. We proceed in several steps. First, in section 4.1 we find the adiabatic limit of the time-dependent system operators. Next, in section 4.2 we find two pairs of master equations in the adiabatic limit, one pair in the interaction picture, the other in the Schrodinger picture. The master equations within a given picture are distinguished by whether we apply the adiabatic approximation once or twice. Next, in section 4.3, we introduce the RWA, which allows us to reduce the Schrodinger picture master equations into Lindblad form. We conclude the formal development by discussing non-adiabatic corrections to our master equations in section 4.4. We move on to a detailed numerical study of an example in section 5, of a ferromagnetic chain in a transverse field, coupled to an Ohmic oscillator bath at finite temperature. We apply two of our Schrodinger picture master equations, with and without the RWA, 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 6. 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.

2. Preliminaries

2.1. Model

We consider a general system-bath Hamiltonian

where HS (t) is the time-dependent, adiabatic system Hamiltonian, HB is the bath Hamiltonian and HI is the interaction Hamiltonian. Without loss of generality the interaction Hamiltonian can be written in the form:

where the operators Aa and Ba are Hermitian and dimensionless with unit operator norm, and g is the (weak) system-bath coupling. The joint system-bath density matrix p(t) satisfies the Schrodinger equation p = —i[H, p(t)], where we have assumed units such that h = 1.

H (t) = Hs (t) + HB + HI ,

Let US(t, t') = 7+exp[-i/t' dr HS(r)] and UB(t, t') = exp(-i(t -1')HB) denote the free system and bath unitary evolution operators (7+ denotes time ordering), and define U0(t, t') = US(t, t') ® UB (t, t'). Likewise, let U(t, t') = 7+ exp[-i f{, dr H(r)] denote the joint Schrodinger picture system-bath unitary evolution operator. We transform to the interaction picture with respect to HS (t) and HB, by defining U (t, 0) = U} (t, 0)U (t, 0), which, together with the interaction picture density matrix p(t) = u0(t, 0)p(t)U0(t, 0), satisfies

U(t, 0) = -iHi(t)U(t, 0), U(0, 0) = 1,

-p(t) = -i[Hi(t), p(t)], p(0) = p(0). dt

(3a) (3b)

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 Schrodinger picture counterpart via:

Hi(t) = Uj(t, 0)H1U0(t, 0) = g £ Aa(t) ® Ba(t),

where we have defined the time-dependent system and bath operators: A„(t) = US(t, 0)AaUS(t, 0), Ba(t) = UB(t, 0)BaUb(t, 0). We are interested in deriving a master equation for the system-only state, Ps(t) = TrB[p(t)] = US(t, 0)TrB(UB(t, 0)p(t)Ub(t, 0))Us(t, 0)

= US(t, 0)TrB(Ub(t, 0)UB(t, 0)p(t))Us(t, 0) = Us(t, 0)ps(t)Us(t, 0), where in the second line we used the fact that UB acts only on the bath.

(5a) (5b)

(6a) (6b)

2.2. Born-Markov approximation

Writing the formal solution of equation (3b) as

p(t) = p(0) - i dr [Hi(r),p(r)],

and substituting this solution back into equation (3b), we obtain, after tracing over the bath degrees of freedom, the equation of motion for the system density matrix pS (t) = TrB p(t)

Ps(t) = -iTrB[Hi(t), p(0)] - TrB

Hi(t), dr [Hi(t - r),p(t - r)]

We make the standard Born approximation assumption, that we can decompose the density matrix as p(t) = pS (t) ® pB + x(t) where x(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, 34].7 Thus the equation of motion reduces to: d f t

-tps (t) = g2 V / dT (Ap(t-T)ps (t-r) Aa(t)-Aa(t) A pit-r)ps (t-r)) Bap(t, t-r) +h.c.,

dt ap J0

where we defined the two-point correlation functions:

Bap(t, t - T) = < Ba(t) Bp (t - T)) = Tr [ Ba(t) Bp(t - T)pB ] . (10)

In equation (9), we have assumed for simplicity (but without loss of generality) that <Ba)0 = Tr[Ba(t)pB(0)] = 0, so that the inhomogeneous term in equation (8) vanishes. Since we assumed that the bath state pB is stationary, the correlation function is homogeneous in time:

Bap(t, t - T) = <Ba(t)Bp(t - r)) = <Ba(T)Bp(0)) = <Bp(0)Ba(r)T = Bap(r, 0). (11)

For notational simplicity we will denote Bap(r, 0) by Bap(r) when this does not lead to confusion. Let us denote the time scale over which the two-point correlations of the bath decay by rB, e.g. IBap(r)l ~ exp(-r/rB). More precisely, we shall require throughout that

f drrn | Bap(r)\- TB+1, n e{0, 1, 2}. (12)

As we show in appendix B, if rB ^ 1 /g, then we can apply the Markov approximation to each of the four summands in equation (9), i.e., replace pS (t - r) by pS (t):

f' f ~ 3 2 / dr (■ ■ ■ ps(t - T) ■■■) Bap(T) * dr (■ ■ ■ ps(t) ■ ■ ■) Bap(T) + O(r^g2), (13)

00 where (■ ■ ■) refers to the products of Aa and Ap operators in equation (9), and where we have also assumed that t ^ rB and the integrand vanishes sufficiently fast for x ^ xB, so that the upper limit can be taken to infinity. Note that by equation (12) the integral on the rhs of equation (13) is of O (rB), so that the relative magnitude of the two terms is O [(grB)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 rB.

2.3. Correlation functions, the Kubo-Martin-Schwinger condition, and the spectral-density matrix

In computing the terms appearing in equation (13) we are faced with integrals of the form drAp(t - r)ps(t)Aa(t)Bap(r) and drAa(t)Ap(t - r)ps(t)Bap(r). Our goal is to express these integrals in terms of the spectral-density matrix

/» to

Tap(v) dr eiMT Bap(r), (14)

7 When this term is not neglected it can be shown to give rise to an inhomogeneous contribution to the master equation [30, 34] and also to restrict the set of system states for which the master equation can be used, due to the requirement of positivity of ps (t) [34]. The smallness of x(t) is consistent with any correlations decaying very rapidly, in our case on the timescale rB of the bath correlation functions. It is also consistent with our goal of developing a master equation for an adiabatically evolving system, whose state at all times is close to the ground state and hence very nearly pure; clearly if ps (t) is pure then x (t) must vanish.

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,

(«) = 1 Yos («) + i Saj3 («), where we show in appendix C that y and S are given by

dT eimT B„s(r) = y;o(a>),

d«' / 1 \

=s;<».

rTO d«' , / 1

S„ß(«) = / — Yaß(«')P -:

J-oo 2n V« - «'

(16a )

If we assume not only that the bath state is stationary, but that it is also in thermal equilibrium at inverse temperature S, i.e. pB = e-sHb/Z, then it follows that the correlation function satisfies the KMS condition [30]

< Ba (t) Bh (0)) = ( Bh (0) Ba (t + iß)>.

If in addition the correlation function is analytic in the strip between t = —i ^ and 0, then it follows that the Fourier transform of the bath correlation function satisfies the detailed balance condition

Yah (-«) = e ß«Yha («).

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 IBab(—t)| = |Tr[pb% (t) BaUB (t) Bb ] | = | Tr[ BaUB (t) Bb Ub (t)Pb] | = k Ba (0) Bb (t)>| = K Bb (t) Ba (0)>| = IBba(t)|, where we used equation (11). Now assume that equation (12) would have to hold for all n. It would follow that

Yah («)

/to I/» TO

Tn eia>T Bah (T)dT|«=0 = 1/ Tn Bah (T)dT

-TO «/—TO

< Tn (|Bha (T)| + |Bah (T)|) dT - 2^ Vn G {0, 1,.

Thus all derivatives of yab («) would have to be finite at« = 0.

On the other hand, it follows from the KMS condition (18) that

d(-«)

Yah (-«)

ße-ß«Yha («) - e-ß^—Yha («)

so that in the limit as « approaches zero from below or above, Y'b (0-) = SYba (0) - Yba (0+).

This shows that in principle y'aa (œ) may be discontinuous at œ = 0. Indeed, the continuity condition Yaa (0—) = Yaa (0+) implies, from the KMS condition recast as equation (21), that

2Ya'a (0) = PYaa (0). New Journal of Physics 14 (2012) 123016 (http://www.njp.org/)

This conclusion can be extended to the entire y matrix by diagonalizing it first. When equation (22) is not satisfied y^a (0) diverges, so that equation (19) does not hold except for n e {0, 1}. A simple example of this is yab (&) = c > 0 (constant) for & ^ 0. Another example is a super-Ohmic spectral density yab(&) = &2/(1 - e-p&) for & ^ 0 and yab(&) = &2/(e-pM - 1) for & ^ 0. Both examples violate equation (19) for n ^ 2. Note, moreover, that when this happens, it follows from equation (19) that /0to t2 (\Bab(-r) \ + \Bab(r)\) dr is divergent, so that we must conclude that \Bab(r) \ > 1/r3 for sufficiently large \r\, meaning that the correlation function has a power-law tail and, in particular, cannot be exponentially decaying.

On the other hand, equation (22) tells us that continuity and a lack of divergence at & = 0 require y'aa (0) to satisfy a condition which prohibits it from being arbitrary. This is indeed the case for the Ohmic oscillator bath discussed in section 5.1, which satisfies equation (22) with finite yaa (0). For this case the bath correlation function is exponentially decaying assuming the oscillator bath has an infinite cutoff. However, as we show in section 5.2 the Ohmic bath transitions from exponential decay to a power-law tail for any finite value of the cutoff, at some finite transition time rtr. In this case we find \Bap(r)\ — (r/rM)-2, where rM is a time-scale associated with the onset of non-Markovian effects, and hence we have to relax equation (12), and replace it with

f * dr rn\Bap(r)\- rB+1, n e {0, 1,...}, (23a)

f™ f™ 9 9

/ dr \Bap(r)\- dr (t/tm)-2 = rM/rtr. (23b)

J Ttr J Ttr

3. 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 tf and the minimum ground state energy gap of Hs is A, i.e.

A = min s1 (t) - s0(t), (24)

t e[0,tf]

where s0 (t) and s1 (t) are the ground and first excited state energies of Hs (t). We shall arrive at master equations of the following general form:

ps (t) = [Luni (t) + Lddss (t) + Csns-ad (t )]ps (t)- (25)

Here Luni = -i[Hs(t) + HLs, ■] is the unitary evolution superoperator including the Lamb shift correction, Lddss is the dissipative superoperator in the fully adiabatic limit and L^o*-^ is the dissipative superoperator with leading order non-adiabatic corrections. Let

h = max \<Sa (s )\3sHs (s )\sb (s ))\, (26)

s e[0,1];a ,b

where s = t/tf is the dimensionless time. To ensure that Luni generates adiabatic evolution to leading order we shall require the standard adiabatic condition

<< 1. (27)

In order to derive equation (25), the three superoperators are ordered in perturbation theory, in the sense that

IILunill»l|LSL11 » IILS-*II, (28)

where the norm could be chosen as the supoperator norm, i.e. the largest singular value (see appendix A). We may then assume that

IILnII» A »IlLdsII. (29)

Combining the O (tb) due to the integral on the rhs of equation (13) (as already remarked there) with the g2 prefactor from equation (9), we have

IILdssII~ g2tb. (30)

To ensure the first inequality in equation (28) we thus require

¿r «1. (31)

The non-adiabatic correction is of order a^ , and when it appears in L^sT^ it is multiplied by the same factor as Lddss, i.e. we have

IILT*II ~ g2tb^• (32)

To ensure the second inequality in equation (28) thus amounts to the adiabatic condition, equation (27). All this is added to the condition

gtb « 1, (33)

for the validity of the Markovian approximation, mentioned in the context of equation (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 tb. If we require the change in the basis to be small on the time scale of the bath tb, we must have:

— « 1. (34)

We justify this claim in section 4.1 and appendix F.

Note that the adiabatic condition, equation (27), implies ATr « Atb, while equation (34)

can be written as ATr « a^" . Putting this together thus yields

h tb / 1 \

B « mini Atb , — . (35)

Atf V ATb

Our other inequalities (equations (31) and (33)) can be summarized as

gtb « min(1, A/g). (36)

10 IOP Institute of Physics <J>deutsche physikalische Gesellschaft

4. Derivation of adiabatic master equations

4.1. Adiabatic limit of the time-dependent system operators

Let us first work in the strict adiabatic limit. We will discuss non-adiabatic corrections in section 4.4. We denote the instantaneous eigenbasis of HS(t) by {|£a(t))}, with corresponding real eigenvalues (energies) {sa(t)}, i.e. HS(t)lsa(t)) = sa(t)lsa(t)), and Bohr frequencies o>ba (t) = sb (t) — £a (t). As shown in appendix E.1 we can then write the system time evolution operator as:

Us (t, t') = Usad(t, t') + ^—, (37a)

Usad(t, t7) = ^ |sa(t)){sa(t')|e~ißa(t,t'), (37b)

where USJd(t, t') is the 'ideal' adiabatic evolution operator. It represents a transformation of instantaneous eigenstate |sa (t')) into the later eigenstate |sa(t)), along with a phase

fra (t, t') = dT [sa (T) — fa (t)] , (38)

where (t) = i{sa(t)\sa(t)) is the Berry connection. The correction term O ^is derived in appendix E.2.

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 Ap(t — x) = Uj(t — x, 0)ApUS(t — x, 0) with an appropriate adiabatic approximation, which will allow us to take this operator outside the integral. To see how, note first that

Us(t — t, 0) = Us(t — t, t)Us(t, 0) = USf(t, t — t)Us(t, 0). (39)

We now make two approximations: firstly, as per equation (37a) we replace US(t, 0) by USad(t, 0); secondly, we replace Uj(t, t — x) by eixHs(t), an approximation justified by the appearance of the short-lived bath correlation function 8ap(T) inside the integrals we are concerned with. Thus, we write

Us(t — t, 0) = eitHs(t)USad(t, 0) + &(t, t), (40)

and find the bound on the error due to dropping S(t,t) in appendix F. Let

frba(t, t') = ^b(t, t') — ^a(t, t'), (41a)

nab (t) = \&a (t))(eb (t )\. (41b)

Neglecting the operator-valued correction term S(t, x) entirely, we have, upon substituting equation (37b) and using eitHs(t)\sa(t)) = eiT£a(t)\sa(t)), that

dT Aß(t — T)ps (t) Aa(t )Baß(T) (42a)

/ dTUsadt(t, 0)e—iTHs(t)Aß eiTHs(t)U|d(t, 0)Ps(t)Aa(t)Baß(r) (42b) o

11 ЮР Institute of Physics <J>deutsche physikalische Gesellschaft

/»то

= / dr У]e-i^ba(t,0)|efl(0))(^(t)|e-iTHs(t)AßeiTHs(t)(t)){sb(0)|Ps(t)Aa(t)Baß{x) (42c) Jo

/»то

J2e-i"ba (t ,0) К (0))(ee (t )| Ae l£b (t ))(£b (0)|Ps (t) A„(t) dr eiT [eb (t)-*a (t)] B„e(r), (42d)

where the approximation in (42a) is shown in appendix F to be O [tb min{1, ^ + ^}]. The

h A2 tf

first term, Ah2r, is the smallness parameter of the adiabatic approximation, which we have

already assumed to be small. The second term, (mentioned in equation (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

dT Aß(t - T)Ps (t) Aa(t )B«e(T) « ^ e^ba (t ,0) Ae ab (t )nab (0)p S (t) Aa(t )rae(^ba (t)),

f arAe(t - T)Ps(t)Aa(t)tf«e(T) ^ > e

Aaab (t) = (efl (t )| Aa leb (t )> = A*aba (t), (44)

and ra^(^ba (t)) is the spectral-density matrix defined in equation (14). Similarly, we have for the other term

/» TO

/ dr Aa(t)Ae(t - r)ps(t)Bae(T) e^ba(t,0)Aeab(t)Aa(t)ПаЬ(0)^(t^(^ba(t)). (45)

J° ab

4.2. Master equations in the adiabatic limit

We are now ready to put everything together. Starting from the Born-Markov master equation constructed from equations (9) and (13), and using the approximations (43) and (45), we arrive at the following one-sided adiabatic interaction picture master equation:

ddt PS (t) = g2 £ e-i^(t ,0) £ Tae^ba (t)) Aeab (t) [nab (0)Ps (t), A„(t)] +h.C. (46)

ab ав

Since we used an adiabatic approximation for Ae(t — x), it makes sense to do the same for Aa (t), i.e. to replace the latter with U^ (t, 0) Aa U^(t, 0). If this is done, we obtain the double-sided adiabatic interaction picture master equation

dtPS(t) = g2 £ e—i[^dc(t(t,0)] £ r„e(«ba(t)) Aacd(t) Aeab(t) [nab(0)Ps(t), ncd(0)] + h.c.

abcd ав

It is convenient to transform back into the Schrodinger picture. Using ps (t) = Ust(t, 0)ps(t)Us(t, 0) (equation (6)) implies that £Ps(t) = ^(t, 0) (i[Hs(t), Ps(t)] + dPs(t)) US(t, 0). Hence, using equation (46), we find the one-sided adiabatic Schrodinger picture master equation d dt

—Ps(t) = -i[Hs(t), Ps(t)] + g2 £ £ r«e(®ba(t)) [L'ab,ß(t)Ps(t), Aa] + h.c., (48)

ab ав

Lab,p(t) = e—ifXba(t,0)Apab(t)Us(t, 0)nab(0)U<[(t, 0). (49)

This form of the master equation has not appeared in previous studies of adiabatic master equations. If we again use the adiabatic approximation for US (t, 0), i.e. replace US(t, 0)nab(0)US(t, 0) by USadt(t, 0)nab(0)USad(t, 0), we obtain the double-sided adiabatic Schrodinger picture master equation

dtps(t) = — i[Hs(t), PS(t)] + g2 ^^ raj3(vba(t)) [Lab,p(t)ps(t), Aa] +h.c., (50)

aft ab

Laba(t) = Aaab (t )\^a (t ))(eb (t )\ = L^ a(t). (51)

Comparing to equation (25), the first term in equations (48) and (50) is Luni(t), while the second

is Ldss(t).

The master equations we have found so far are not in Lindblad form, and hence do not guarantee the preservation of positivity of pS. We thus introduce an additional approximation, which will transform the master equations into completely positive form.

4.3. Master equation in the adiabatic limit with rotating wave approximation: Lindbladform

In order to arrive at a master equation in Lindblad form, we can perform a RWA. To do so, let us revisit the t ^ to limit taken in the Markov approximation in equation (13). Supposing we do not take this limit quite yet, we can follow the same arguments leading to equation (42c), which we rewrite, along with the adiabatic approximation Aa(t) & US?dt(t, 0)AaUSd(t, 0). This yields

dT Aft(t — t)Ps (t) Aa(t )Baft(T) (52a)

& i dT ^e—i[frba(t,0)+frdc(t,0)]\ £a(0))(sa(t)\Aft\eb(t))(eb(0)\ps(t)\ec(0))

K(ec (t )\Aa\ed (t ))(ed (0)\ei ™ba (t) Baft(x). (52b)

We note that frdc (t, 0) + frba (t, 0) = /J dT [todc (t) + toba (T) — (fod (t)—&(t)) + (fob (T)—foa (t))]. One can now make the argument that when the t ^ to limit is taken, terms for which the integrand vanishes will dominate, thus enforcing the 'energy conservation' condition b = c and a = d, or c = d and a = b (without over-counting the case a = b = c = d). This is a similar RWA as made in the standard time-independent treatment, although here, the approximation of phase cancelation is made along the entire time evolution of the instantaneous energy eigenstates. Clearly, in light of the appearance of other terms involving t in equation (52b), this argument is far from rigorous. Nevertheless, we proceed from equation (52b) to write, in the t ^ to limit,

i dT Aft(t — t)Ps (t) Aa(t )Baft(T) & ^ Aftab (t) Aa ba (t )^b (0)ps (t )nba (0)Taft(toba (t))

j0 a=b

+ ^ Aft aa (t) Aa bb (t )^a (0) Ps (t )nbb (0)^ft (0). (53)

We show in appendix G how, by performing a transformation back to the Schrodinger picture, along with a double-sided adiabatic approximation, we arrive from equation (53) at the Schrodinger picture adiabatic master equation in Lindblad form:

PS(t) = -i [Hs (t) + Hls (t), PS(t)] (54a)

Y Y Yaß(«ba(t)) Lab.ß(t)PS(t)Lab a(t) - 1 {Llb a(t)Lab,ß(t), Ps(t)

aß a =b

+ Y Y ^«^(0) Laa,fi(t)ps(t)L^.„(t) - ^ {L^,a(t)Lbb,p(t), PS(t)}

where the Hermitian Lamb shift term is

(54b) (54c)

Hls (t) = Y

Y Lib.a (t)Lab,ß(t)Saß(«ba (t)) + Y LJa,a(t)Lbb.ß(t)Saß(ö)

Since the bath correlations are of positive type, it follows from Bochner's theorem that the matrix y —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 equations (48), (50) and (54) all generalize both the standard Redfield and Lindblad time-independent Hamiltonian results to the case of a time-dependent Hamiltonian in the adiabatic limit8. 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 sum over the energy eigenvalues can be replaced by the sum over their differences, such that:

, Lab,a ^ = £ X^a I Aa |^b X^b I, &ba ^ (56)

a ,b & &b —£a =&

The resulting equation becomes:

PS(t) = -i [Hs + Hls, p(t)] + £ £ Ya^(&) ( (t)A&a - 1 {A& aA^, ps(t)M , (57)

which is the standard form for the time-independent Lindblad master equation. This should make the physical meaning of our derivation evident. 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 (54) are the two master equations we use for numerical simulations presented later in this paper. We note that equation (54) appears similar to the Markovian adiabatic master equation found in [24], but is more general and did not require the assumption of periodic driving.

8 See for example equation (3.143) in [30] and equation (8.1.33) in [35].

4.4. Non-adiabatic corrections to the master equations

So far, we assumed the adiabatic limit of evolution of the system (see equation (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/tf, and we rederive in appendix E.1 the well known result [5] that to first order

Us(t, t') = USad(t, t') [I + Vi (t, t')], (58)

Vi (t, t') = - £ \Sa (t ')){Sb (t ')\Jt dT e-i'tba ^ (Sa (z)\Sb (t)) .

t') = — 2J \ea(t'))(eb(t')\ I dTe—"tba(z,t')(ea(z)\^(z)). (59)

Thus, to derive the lowest order non-adiabatic corrections to our master equations is a matter of repeating our calculations of sections 4.1 and 4.3 with U|d(t, t') replaced everywhere by the first order term U|d(t, t') V1 (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 equation (27), which is now replaced by

h < A2, (60)

i.e. with the ^ replaced by a mere <. However, we would still like to perform the approximation of equation (40), in the sense that Us(t — z, 0) & eizHs(t)USad(t, t') [I + V(t, t')]. This still requires

- « -, (61) f TB

which is what allows the use of eizHs(t) in this approximation (as shown in appendix F). Taken together, these two conditions are weaker than equation (35), which we can rewrite as h « min(A2, Z2).

Recall that to ensure the validity of the Markov approximation and ||Lni II > II Lddss 11, we also had to demand the inequalities in equations (31) and (33); these can be now summarized as

max ^g, g^ « -1, to be compared with equation (61).

5. An illustrative example: transverse field Ising chain coupled to a boson bath

5.1. The model

We proceed to use our master equations to study the evolution of the Ising Hamiltonian with transverse field

HS (t) = A(t) HX + B (t) HZ, (62a)

30 25 20 15 10 5

L° Tf

0.2 0.4 0.6 0.8 Figure 1. The A(t) and B(t) functions in equation (62a). A(0) = 33.7 GHz.

HSZ s - £ hi of + £ Jj of of, (62c)

i = 1 i, j = 1

where the functions A(t) and B (t) are shown in figure 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 HSX turned on while the Ising term HSZ 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

то N

HB = £ akb\bk, H = £ of ® B, B = £ gk (b\ + bk) , (63)

k=1 i =1

where b\ and bk are, respectively, raising and lowering operators for the kth oscillator with natural frequency , and gJk is the corresponding coupling strength to spin j. This is the standard pure dephasing spin-boson model [36], except that our system is time-dependent. The resulting form for our operator L (equation (51)) is

Labi = |£a (t )>(£a (t )|of I^b (t )><£b (t )| = Aiab (t )|^a (t )><£b (t )|. (64)

Recall that our analysis assumed that the bath is in thermal equilibrium at inverse temperature ^ = 1/(kB T), and hence is described by a thermal Gibbs state pB = exp (-^HB) /Z. We show in appendix H that this yields

2nJ(M) • j ,

Yij(»> = , _ e_wl)g;,|gj.|(e(a) + e-«e(-a») (65a)

Si,0» = J0 d.^¿rf (p (-JL-.) (65b)

where only one of the Heaviside functions is non-zero at & = 0. To complete the model specification, we assume an Ohmic bath spectral function

J (&) = n& e~&/&c, (66)

where &c is a high-frequency cutoff and n is a positive constant with dimensions of time squared.

It is often stated that the terms associated with the Lamb shift HLS (equation (55)), i.e. equation (65b), can be neglected, since the relative order of S and HS is g2zB/A, and indeed we have assumed g2zB/A « 1 (equation (31)). However, this argument is misleading for two reasons. Firstly, S can be divergent, as is easy to see in the limit &c = to for the Ohmic model (66), where for & » maxba &ba (t), the integrand tends to a constant. Secondly, S should be compared to y, as both are of the same order g2zB, and both result in changes to the system relative to its unperturbed state. Indeed, in the interaction picture with respect to HS + HLS (recall equation (54a)), the overall transition rates between states with quantum numbers a and b will depend on the dressed (i.e. shifted) energy gaps &ba + &LS The importance of this Lamb shift effect was also stressed by de Vega et al [32]. We analyze the Lamb shift effect in section 5.3.

Finally, we note that although the harmonic oscillators bath with linear coupling to the system provides a y 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 [37]. 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 y matrix that satisfies the KMS condition will suffice without too much concern about how it is derived.

5.2. Correlation function analysis

In light of the subtleties alluded to in section 2.3 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 equation (65a) and substituting J(&) from equation (66), it is possible to compute the bath correlation function analytically for the resulting

2n&e—\&\/&c b

Yab (&) = 1—e—,& gagb, (67)

by inverse Fourier transform of equation (16a). The result is

= (¿c^)(1+¿c - £)), (68)

where ty(m) is the mth Polygamma function (see appendix I for the derivation). We first assume

P&c » 1, (69)

and consider an expansion in large /5&c:

W(n)( 1 - f) + f<n>(f)

(/ \ œ w\ 1 — — -n2csch2 ( — ) + V —^-^ ' „vP/ | , (70)

\P) t=i (n - 1)! (M )n y

followed by an expans i on in large x/fi :

Bab(x) = j2 gagb(--4n2e-T/TB + j-j-^ + O (e-2T/TB, t-3)) . (71)

ln(S (0)

IOP Institute of Physics <j)Deutsche physikalische gesellschaft ln(£ (t))

t (ns)

Figure 2. (a) An example of the bath correlation function Bab (t) for R = 2.6-1 ns and cc = 32n GHz. The solid black curve is the function ln |Bab (t)|. The blue dashed curve is the function ln4n2 exp(-2nT/R) and the red dot-dashed curve is the function ln2RT-2/cc. (b) Same as in (a), but cc = 8n GHz, and the blue dashed curve is the function ln d0 exp(-(2n/R + c1 cc + c2CR)t) with c1 — -0.039, C2 - -0.004, do - 33.13. (a) C = 32n GHz, (b) C = 8n GHz.

This expansion reveals the two independent time-scales that are relevant for us. First, there is the time scale tb associated with the exponential decay (corresponding to the true Markovian bath), given by:

«c^œ f R

TB TB = 2n ' (72)

then the time scale associated with non-Markovian corrections (the power-law tail):

Tm = J —. (73)

For sufficiently large oc, these two time scales capture the two types of behavior found in Bab (t), as illustrated in figure 2(a).

The transition between the exponential decay and the power-law tail occurs at a time Ttr given by 4n2e-Ttr/tB = (^)2, or equivalently by

exp(6>) K n O Ttr nA\ --— = - R cc, 6 = 2n—. (74)

62 r R v 7

This transcendental equation has a formal solution in terms of the Lambert- W function [38], i.e. the inverse function of f (W) = WeW, as can be seen by changing variables to y = -6/n, and rewriting 6ne-6 = a = 2/(Rcc) as yey = —a1/n/n, whose solution is 6 = — ny = —nW(-1 a1/n). However, for our purposes the following observations will suffice. We seek a Markovian-like solution, where Ttr is large compared to the thermal timescale set by R, i.e. we are interested in the regime where 6 ^ 1. In this case we can neglect 62 compared to e6, and approximate the solution to equation (74) by 6 — ln(Rcc/2). Thus

Ttr - R ln(Rcc). (75)

This agrees with the first term of the asymptotic expansion W(x) = ln(x) - lnln(x) + '^X^ + • • •, which is accurate for x > 3 [38].

When P&c » 1 is not strictly satisfied, the exponential regime is less pronounced, and tB is corrected by powers of &c. By dimensional analysis, the corrections must be of the form:

Tb = rR +

T, "i ß

where cn are constants of order one that must be fitted (see figure 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

< min < 2tb

tb h tB h

&c ln(P&c)

which can be interpreted as saying that the cutoff should be the largest energy scale. Equation (77) joins the list of inequalities given in section 3 as an additional special condition that applies in the Ohmic case, along with P&c » 1.

5.3. Numerical results

For concreteness, we take glM = g, &c = 8n GHz and T = 20 mK ^ 2.6 GHz (in units such that h = 1; this is the operating temperature of the D-Wave Rainier chip [2]), corresponding to tb = 0.06 ns for the Ohmic model with infinite cutoff, equation (72), and tb ^ 0.07 ns for &c = 8n using equation (76). For this value of &c the transition between the exponential and power-law regimes is still sharp (see figure 2(b)) and occurs at approximately Ttr = 0.33 ns. For these values, we satisfy at least one of the cases from equation (77), including numerical

prefactors: < 2tb.

We focus on the N = 8 site ferromagnetic chain with parameters: 1

h 0 = 4, hi >0 = 0, Ju+1 = —1, i = 0,..., 7, (78)

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:

e—P Hs (0)

PS (t = 0) = . (79)

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 [39, 40].

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 HS (t) and the instantaneous density matrix predicted by our two master equations (50) (non-RWA) and (54) (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 tf (therefore more adiabatic) evolution shows a smaller difference between the two master

(s0(t)\p(t)\s0(t)) 1.0

IOP Institute of Physics <j) deutsche physikalische Gesellschaft <e0(0lp(0le0(0>

1.0 tf

0.6 (b)

1.0 tf

Figure 3. Fidelity, or overlap of the system density matrix with the instantaneous ground state along the time evolution for tf = 10 js and ng2/(h2) = 1.2 x 10-4/(2n). The solid blue curve was calculated using equation (50) (no RWA), while the dashed red curve was calculated using equation (54) (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.

{s0(tf)\p(tf)\s0(tf)) 1.0

<fiö(i/)|p(i/)|e0(i/)>

0.8 0.6 0.4 0.2

0.75 0.5

0.00025 0.0005

0.4 0.5

Ö.Ö000 0.0001 0.0002 0.0003 0.0004 0.000 (b)

Figure 4. Dependence of the fidelity at tf = 10 js 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. (a) RWA, (b) non-RWA.

equations for the final fidelity. The results in figure 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 = tf as we change the coupling strength g. We observe (see figure 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 non-RWA equation, we see an increase in fidelity and a subsequent violation of positivity. These results bring to light the relative advantages and disadvantages of 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 non-RWA can lead to physically different conclusions, e.g. in the context of Berry phases in cavity QED [41].

5.4. The four different phases

5.4.1. Phase 1—the gapped phase. For times sufficiently close to the initial time, the ground state of HS(t) is the ground state of HSX, i.e. the state |0) = ®N=iI—)j, where l±)j = (||)j ±

It)j)/V2 with energy 80(0) = — NA(0), and where | |)j, | t)j are the +1, —1 eigenstates of oz (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. Ii) = I—) j |+)i ®]N=i+1 I—) j, with energy 81 (0) = —(N — 2)A(0). Therefore the gap between the ground state and the first excited states is:

A(t) = 61 (t) — 80(t)t& [—(N — 2)A(t)] — [—NA(t)] = 2A(t), (80)

which is at least twice as large as our kB T & 2.6 GHz, almost until A(tc) = B(tc) & 5 GHz at tc & 0.35tf (see figure 1). Noting that atz|±)i = I±)i, we can write the following relations in terms of the ground and first excited states:

<0 = A = —<0i, of Ii ) = I0). (81)

Noting that of | j) is a doubly excited state V j = i e {1,..., N}, we truncate the problem to the ground and first excited states only, so that there are just three types of values of y(w) we need to consider: y(0), y(^i0), y(<v0i). Recalling the KMS condition, equation (18), we have

Y(^0i) = e—P<0 y(<i 0), (82)

which shows that upward transitions are exponentially suppressed relative to downward transitions, by a factor ranging between e—2A(0)/(kBT) & e—67 4/2 6 & 7 x 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 (79).

To make this argument more precise, denoting p00 = (0|p|0) and pii = {i |p|i), we can write the effective (truncated to the ground and first excited states) Lindblad equation (54) as the simplified rate equations

Pii & Yii (<0) {—Pii + e—P<0P00), (83a)

P00 & 1 — Nph , (83b)

where we have assumed that the system is initially in the thermal state (79) and the gap is large (relative to kBT). A derivation of equation (83) can be found in appendix J. We compare our simulation results with the results from the above equations in figure 5 and find very good agreement for early times.

0.002 0.004 0.006 0.008 0.010 tf

1.4x10" -6

1.2x10" -6

1. x 10" -6

8. x 10" -7

6. x 10" -7

4. x 10" -7

2. x 10" -7

0.002 0.004 0.006 0.008 0.010 tf (b)

Figure 5. Early evolution of diagonal elements of the system density matrix in the thermal regime. Red dots are from using the Lindblad equation (54), while the solid blue curve is from the solution of equation (83). (a) Ground state, (b) excited states.

,-jSA (t)

0.4 0.6

1.0 tf

0.4 0.6 (b)

1.0 tf

Figure 6. The behavior of the gap and the exponential suppression factor along the time evolution. The dashed line is kBT ^ 2.6 GHz. (a) The gap, (b) suppression factor.

5.4.2. Phase 2—the excitation phase. When A(t) becomes small enough such that RA(t) ~ O (1), then the KMS condition no longer suppresses excitations from the ground state to higher excited states (see figure 6). If we interpret the master equation as a set of rate equations for the matrix elements of p, we can identify the rate of scattering into the ith state from the jth state as being the term in the pa equation with coefficient pjj . Therefore, we find that the rate of scattering from the ground state to excited states is given by

Excitation rate a y (A(t))e-M(t)p00, (84)

whereas the relaxation rate is given by

Relaxation rate a y (A(t))pu. (85)

Therefore, as we emerge from the gapped phase, we have p00 » pn, so scattering into excited states dominates over relaxation into the ground state. This explains the loss of fidelity in phase 2.

5.4.3. Phase 3—the relaxation phase. As the gap begins to grow again and the suppression factor shrinks (see figure 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 [26, 32, 42].

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 tf 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 [43]. To see this more explicitly in the context of our analysis, note from equation (83b) that if there is an exponential number N of pn coupled to p00 (i.e. N is an exponentially large fraction of the dimension of the system Hilbert space), then P00 decreases exponentially. By equation (84) this means that all pii become exponentially small, but not zero. In phase 3, by equation (85), the relaxation is suppressed as long as the gap is not very large, because the relaxation is proportional to Pn, which is exponentially small. This analysis presumes that the system-bath Hamiltonian has non-negligible coupling between the ground and excited states, i.e. that I(0IojzIi)I > 0 in our model (see equation (J.6a)). 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 tf 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 tf. The fidelity continues to grow (albeit slowly) for sufficiently high tf. This is illustrated in figure 7.

5.4.4. Phase 4—the frozen phase. As the gap continues to grow, the relaxation phase ends (notice that the tail in figure 6(b) is longer than the actual relaxation phase) and the system's dynamics are frozen in the ground state. This is because HS becomes almost entirely diagonal in the oz basis, and so the off-diagonal components of the Labj operators vanish (or become very small), i.e. Aiab(t) = (8a(t)IozI8b(t)) a 8ab (equation (64)) 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.

5.5. 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 equation (A.1))

<Go(Olp(Ö|eö(Ö> 1.0

1.0 tf

Figure 7. Time evolution of the system density matrix using the RWA equation with ng2/(h2) = 0.4/(2n) for different tf's. Solid blue curve is for tf = 10 p,s, dashed red curve is for tf = 100 ^s and dot-dashed green curve is for tf = 1 ms.

IIP(0-PGibbs(0lll

Figure 8. Trace-norm distance between the evolving system density matrix (using the Lindblad equation) and the Gibbs state for an annealing time of tf = 10 ^s and ng2/(h2) = 0.4/(2n).

between the instantaneous system density matrix and the instantaneous Gibbs state pGibbs (t) = exp[-¡5HS(t)]/Z, where Z = Tr [exp (—¡5HS(t))] is the partition function, for the Ising chain discussed above. The result is shown in figure 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/tf ~ 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 ox and oy in the system-bath interaction, which we neglected in writing the interaction Hamiltonian HI in equation (63), to become important, and to disrupt the frozen phase, allowing the system to fully equilibrate into the Gibbs state.

6. 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 tf. 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 / tf) can be incorporated 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.

Acknowledgments

We are grateful to Robert Alicki for extensive and illuminating discussions regarding the relation between the KMS condition and correlation functions. We are also grateful to Mohammad Amin for useful discussions and for reading an early version of the manuscript. This research was supported by the ARO MURI grant W911NF-11-1-0268 and by NSF grant numbers PHY-969969 and PHY-803304 (to PZ and DAL).

Appendix A. Norms and inequalities

We provide a brief summary of norms and inequalities between them, as pertinent to our work. For more details see, e.g., [44-46]. Let |AI = VA^A. Unitarily invariant norms are norms that satisfy, for all unitary U, V, and for any operator A: || UAV ||ui =\\ A ||ui. Examples of unitarily invariant norms are the trace norm

where si (A) are the singular values (eigenvalues of | A|), and the supoperator norm, which is the largest eigenvalue of |A|:

|| AHn= sup y/(fI Af AIf) = max si (A).

If):(f If) = 1 1

Therefore || AIf) K|| A W^ for all normalized states If), and || A W^ A ||1.

25 IOP Institute of Physics <J>deutsche physikalische Gesellschaft

All unitarily invariant norms satisfy submultiplicativity:

II AB ||ui ^ || A|U II B ||ui. (A.3)

The norms of interest to us are also multiplicative over tensor products and obey an ordering: || A ® B| = || Ah||B| i = 1, (A.4a)

IIAB || ui ^ IIA|| OO 11 ^ 11 ui.

In particular, || AB ||i<|| A |U| B ||i.

Another useful fact is that the partial trace is contractive, i.e.

||TrB(x)||i, ||TrA(x)||i < IIX||i,

(A.4b)

for any operator X acting on the Hilbert space HA ® HB. Actually this result extends to other unitarily invariant norms, with a prefactor depending on the dimension of the traced-out Hilbert space [46].

Appendix B. Markov approximation bound

We wish to derive an upper bound associated with the error from the approximation made in equation (13), which involves the replacement of pS (t — t) by pS (t) and the extension of the upper integration limit to infinity, i.e.

:= f dT {Aß(t - t)ps(t - t)Aa(t) + •••} Eaß(x) Jo

^ dT {Aß(t - t)ps(t)Aa(t) + •••} Baß(T) = I + Ai + A2,

:= / dT {Aß(t - t) [ps(t) - ps(t - t)] Aa(t) + • • •} Baß(T) o

dT {Aß(t - T)ps(t - T)Aa(t) + • • •} Baß(T),

and where the ellipsis denotes the three other summands appearing in equation (13). The Markov approximation is more accurate the smaller the error terms A1 and A2.

Consider first the A1 term. We shall show that it is of order g2For simplicity we consider only the first of its four summands (the bounds for the other three are identical to that for the first):

/ dTAß(t - t) [ps(t) - pS(t - T)] Aa(t)Baß(T) o

/ dT||Aß(t -o

T)\\nnT max

= dT T max

Jo t'e[t t ,t]

t ;e[t—t ,t ]

dp s (t') dt'

dp s (t')

dt' IBaß(t)l,

|| Aa(t )UrolBaß(T)l

where we used the triangle inequality, ||Aa(t)||TO = ||A\\ = 1 and ||XY||j ^ ||X||TO||Y||1 (see appendix A). We can now use equation (9) to upper-bound the time derivative:

-tps (t) dt

Ti dr (||Aß(t - r

a.ß Jo

IIPs (t - T)y 11| Aa(t) y to + •••) IBaß(T)\

= 4g2 Tf dr\BaP(T)\,

a.ft J°

where the factor of 4 is due to the same number of summands appearing in equation (9). Substituting this bound back into equation (B.3) we have

dps (t')

dr t max

t 'e[t-r,t ]

/»TO /»t

\Baß(r)\ ^ 4g2T drr \Baß(r)\ max / dr'\Baß(r')\ 1 ß V J0 ß ^-T't]Jo P

< 4gdTT \Eap(x)\xB ^ 4g2rB^]l,

aft aft

where we used /°t dr r\Bap(x' )\ ^ /°TO dr r\Bap(x 7)| and equation (12).

This is to be compared to the term we use after the Markov approximation:

II /» TO /»TO

II / dr Aft(t - T)ftps (t) Aa (t )Baft(T) ^ / dr\Baft(T)l — Tb . \\J° 1 J°

Comparing equations (B.5) and (B.6), we see that the relative error is O[(gtb)2]. Next consider the A2 term. We have

dr Aß(t - r)ps(t - r)Aa(t)Baß(r)

/TO /» TO

dr || Aß(t - r) | to I Ps (t - t)IIi || Aa (t)IUBaß(r)\ = J dr \Baß(r)\.

The last integral converges provided \Bap(T) \ — 1/tx with x > 1. This is typically the case. For example, for the Ohmic spin-boson model discussed in section 5.2 we show that for t > Ttr — ftln(ftvc), the bath correlation function \Bap(T)\ — ftg^Tj. In this case, then, we can set t = Ttr and bound

/ ,o fTO 1

IIA2Ill dr —

ft^O J Ttr

t 2 ft 2 ^cln(ft^c)

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 equation (B.8) diverges, is incompatible with weak coupling, and requires the so-called singular coupling limit [22, 3°].

Appendix C. Properties of the spectral-density matrix Tap(&)

Introducing the Fourier transform pair

fTO d&> ; Baft(T. °) = 0 J-oo 2n

drela>rBaß(r, 0). (C.l)

and using the property that

f dT ei(o-o/)T = n8(o - o') + iP ( ——J , (C.2)

Jo \o - o'J

where P denotes the Cauchy principal value9, we have, using equation (14),

pro pro do'

Taß(o) = dT eior / — e-io'Tyaß(o') Jo J-ro 2n

/ro do' fro 1

— Yaß(o') dT ei(o-°)T = -Yaß(o) + iSaß(o), (C.3)

-ro 2n Jo 2

in agreement with equation (15), where

fro do' , ( 1 \ Saß(o) = —Yaß(o')P[-; , (C.4)

J-ro 2n \o - o'/

in agreement with equation (16). Note that:

Blß(T, o) = < Ba(T) Bß (o)>* = Tr [ Ba(T) Bß(o)ßB ]* = Tr [ Bß(o) Ba(T)pB ] = Bßa(o, T). (C.5) When pB commutes with HB (which we have assumed), we further have

Bßa(o, t) = Tr [Bß(o)UB(T)Ba(o)ÜB(t)pb] = Tr [peUb(t)Bß(o)UB(t)Ba(o)]

= Bßa(-T, o). (C.6)

Thus the spectral-density matrix satisfies:

o ro o ro o o

T*aß(o) = dre-iOTB*aß(T, o) = dz e-ioTBßa(-T, o) = dr elorBßa(x, o)

Jo Jo J-ro

= 2 Yßa (o) - iSßa (o). (C.7)

Appendix D. Proof of the Kubo-Martin-Schwinger condition

The proof of the time-domain version of the KMs condition, equation (17), is the following calculation:

<Ba(T)Bb(o)> = Tr[peUB(T, o)BaUe(t, o)Bb] = ^Tr[Bbe-(ß-ir)HeBae-irHe] (D.1a)

= —Tr[ Bb ei(r+iß) He Ba e-i(r+iß) Hee-ß He] (D.1b) Z

= Tr[peBbUB(r + iß, o)BaUe(r + iß, o)] (D.1c)

= < Bb (o) Ba (T + iß)>. (D.1d)

9 By definition, P (1) [ f ] = limefmx- €] ^XT dx, where f belongs to the set of smooth functions with compact support on the real line R.

Figure D.1. Contour used in proof of the KMS condition.

Note that using the same technique it also follows that

( Ba (t) Bb (0))-( Bb (-T - i ft) Ba (0)).

To prove that this implies the frequency domain condition, equation (18), we compute the Fourier transform:

/TO /»TO

dreiitfT {Ba(r)Bb(0)) = drela>r {Bb(-r — i p)Ba(0)).

-TO J —TO

To perform this integral we replace it with a contour integral in the complex plane, §C dre1MT {Bb(—r — i ft)Ba(0)), with the contour C as shown in figure D.1. This contour integral vanishes by the Cauchy-Goursat theorem [47] since the closed contour encloses no poles (the correlation function {Bb(r)Ba(0)) is analytic in the open strip (0, —i ft) and is continuous at the boundary of the strip [48]), so that

1 (-}=0-l (-} *L (->(->+i

(•••) +\ (•••) +

where (• • •) is the integrand of equation (D.3), and the integral f^ is the same as in equation (D.3). After making the variable transformation r = —x — i ft, where x is real, we have

(...) -

( Bb ( x ) Ba (0)).

Assuming that {Ba(±to — ip)Bb(0)) = 0 (i.e. the correlation function vanishes at infinite time), we further have ft (• • •) = f^ ( • • • ) = 0, and hence we find the result:

dTe№T (Bb(-T - i ft)Ba(0))- e

(bb(t)ba(0))- eftmyba(-oS),

which, together with equation (D.3), proves equation (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., [5] for additional details.

—io x

29 IOP Institute of Physics <J>deutsche physikalische Gesellschaft

E.1. Adiabatic perturbation theory

To consider adiabatic corrections, we recall that Us satisfies

idtUs(t, t') = Hs(t)Us(t, t'), (E.1)

Hs (t) = ^ £a (t )\Sa (t)><Sa (t)\. (E.2)

a | a >< a

Define the 'adiabatic intertwiner' W(t, t):

W(t, t') = ^ \(t)><Sa(t')\ = r+exp \-i f dr K(r)

a L Jt' J

where the 'intertwiner Hamiltonian' is

K(t) = i[9tW(t, 1f)]Wf(t, t') = i^ \ea(t)><Sa(t)\ (E.4a)

= i[Po(t), Po(t)], (E.4b)

and where P0(t) = |e0(t)}{s0(t)| is the projection onto the ground state of HS(t). The result (E.4b) is by no means obvious and is proven in section E.3. To extract the geometric phase we define

HG (t) = ^ (t) | sa (t )}{sa (t )|, (E.5a)

<Pa (t) = \{ea (t )|£a (t)}, (E.5b)

HS (t) = Hs (t) — Hg (t). (E.5c) Now define V via the 'adiabatic interaction picture' transformation:

V(t, t') = Wf (t, t')Us(t, t'), (E.6)

along with

Hs(t, t') = Wf(t, t')Hs(t)W(t, t') = J2 Sa(t)\Ba(t')><£a(t')\, (E.7a)

Hg(t, t') = Wf(t, t')Hg(t) W(t, t') = ^<pa(t)\Sa(t')><Sa(t% (E.7b)

HS(t, t') = Hs(t, t') — Hg(t, t'), (E.7c)

IK(t, t') = Wf(t, t')K(t)W(t, t') = iWf(t, t')dtW(t, t'). (E.7d)

Note that the time dependence of HS and HS (t) is entirely in the energy eigenvalue and not in the eigenstates. Then V obeys the Schrodinger equation:

i dtV(t, t') = HSad(t, t') V(t, t'). (E.8)

HSad(t, t') = (t, t') — ^(t, t') = Wf (t, t')[HS(t) — K(t)] W(t, t'). (E.9)

When the evolution is nearly adiabatic (t, t') is a perturbation, so that we consider a solution of equation (E.8) for V of the form:

V (t, t') = Vo (t, t') (1 + Vi (t, t') + • •• ), (E.10)

with the zeroth order solution associated with the purely adiabatic evolution, including the geometric phase:

Vo (t, t ') = 7+exp

-if dz HS (x, t ' )

(E.11a)

Usad(t, t') = W(t, t')Vo(t, t') = £ \Sa(t)){Sa(t')|e-itidT[ea(x)-*a(t)]. (E.iib)

Differentiating with respect to t, we have

Uf (t, t') = W (t, t') V0 (t, t') + W (t, t') V0 (t, t') (E.12a)

= —i K(t)W(t, t') V0(t, t') — i W(t, t')HS(t, t')V>(t, t') (E.12b)

= —i [Hsad(t) - Hg(t)] Usad(t, t'), (E.12c)

Hsad (t) = K (t) + Hs (t). (E.13) Plugging the equation (E.10) expansion into equation (E.8), we obtain, to first order in 1/tf:10

i9t V1 (t, t') = — V0f (t, t')[K (t, t') — Hg(t, t')] V0(t, t') (E.14a)

= —iUSadt(t, t')8tW(t, t')V0(t, t') + Usadt(t, t')HG(t)USad(t, t'). (E.14b) Note that Usadt(t, t')HG(t)USad(t, t') = J2a (t)lsa(t')){ea(t')|. Therefore, integrating, and using

^ba (t, t') = j dr Mba (r), Wba (r) = [Sb(T) — fyb(t)] - [^a (r) — (t)], (E.15) we can write the solution as:

Vi (t, t') = -J dT

Uf (T, t ')dx W (T, t ') Vo (T, t ' ) - £ (T)\ea (t ')){£a (t )\

= — E f dre-iflbaISa(t')){Sb(t')l{Sa(r)ISb(r)). (E.16)

a=bJt'

Therefore, the system evolution operator can be written as:

Us(t, t') = USd(t, t') + Q(t, tf)USd(t, tf) + O (tf—2), (E.17)

where the correction term can be made appropriately dimensionless in a more careful second order analysis, and where

q (t, t') = usad (t, t') V1 (t, t' )uadt (t, t')

= Ee—iMba(t,t')( — f dreiMba(r,t/){Sa(r)ISb(r))] S(t)){sb(t)l. (E.18)

a=b V Jt' '

10 To see that this is an expansion in powers of 1/ tf, transform to the dimensionless time variable s = t/tf. New Journal of Physics 14 (2012) 123016 (http://www.njp.org/)

E.2. Adiabatic timescale analysis

Equation (E.17) 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 ^ 1 (equation (27)). More rigorous analyses replace the ^ condition with an inequality relating the same parameters to the fidelity between the ground state of HS (tf) and the solution of the Schrodinger equation at tf, e.g. [6, 10-12].

Starting from equation (E.18) and taking matrix elements we have

Qab(t, t') = e-i/ba(t- jf drei/ba<ea(r)\3T^(t)>

(E.19)

Changing variables to dimensionless time s = t/ tf, ¡ba (t, t') becomes tf Jtsf/tf ds' o>ba (s') = tf/ba(s, t'/tf), and

t drei/ba^<sa(r)\3T\Sb(t)>= f ds7 eitf/ba(s'^t'/tf)<sa(s')\3A£b(sf)>, (E.20)

Jt' Jt' / tf

where / now involves a dimensionless time integration. Using the fact that eitf/iba(s ,t /tf) =

tf &ab (s') ds'

i--^eitf//ba(s',tf/tf), we can integrate equation (E.20) by parts as

f ds' eitf/ba(s'/T<sa(s7)\ds(s')> =-^-eitf/ab(s'tf)<ea(s')\dAeb(s')>

Jt'/tf tf^ab (s')

t' / tf

i PS eitf / ab (s ,t / tf) d

— ds7 -—- — <^a (s s )\9s '\sb (s ')>. (E.21)

tfjt' /tf Vab (s') ds7

Now note that for non-degenerate energy eigenstates, differentiation with respect to t of HS(t)\sa(t)> = £a(t)\sa(t)>, and substitution of s = t/tf in equation (26), directly yields the relation:11

, \ | . /f\\ fa (t )\Hs (t )\Sa (t )> h

<sb (t )\e a (t )> =-—--—. (E.22)

^ab (t) - tf

Substitution into equation (E.21) yields, with the help of equation (E.22),

\ Qab (t, t' )\ ^

(s )\3sHs (s )\eb (s )>\

Kb (s )tf

t / tf

t' / tf

t / f eitf /ab (s' ,t' / tf) d

ds7 - — <ea (s' )\9s / Hs (s Ofa (s f)>

/tf vzab(s')tf ds'

Continued integration by parts will yield higher powers of the dimensionless quantity

. Thus a sufficient condition for the smallness of \ Qab (t, t') \ for all a, b is that

<ea (s )\3t Hs (s )\eb (s )> (o2ab (s )tf

maxs;a,b \<ea(s)RHs(s)\eb(s)>\

min^ab toah (s)tf

« 1, (E.23)

a b^ab

11 Differentiating (where X = dtx), we have HS\ea(t)> + HS(t)\ea(t)> = ea(t)\ea(t)> + ea(t)\ea(t)>. Taking matrix elements gives <eb(t)\HS\ea(t)> + eb(t)<eb(t)\ea(t)> = ea(t)8ab + ea(t)<eb(t)\ea(t)>, which yields equation (E.22) provided \eb(t)> is an eigenstate non-degenerate with \ea(t)>.

namely the condition in equation (27), where we assumed that the minimal Bohr frequency is the ground state gap, mi,0. Of course, this argument is by no means rigorous, and in fact it has been the subject of considerable discussion, e.g. [49-54]. For rigorous analyses see, e.g., [5-7, 10-12]. For our purposes equation (27) suffices.

E.3. Intertwiner Hamiltonian

Here we provide an explicit proof of equation (E.4b), a well-known result due to Avron et al [55]. The original proof lacks some detail, so our purpose here is to fill in the gaps, for completeness of our presentation (see also [25] for details of the proof). Define the ground state and orthogonal projectors P0(t) = |s0(t)){s0(t)| and Q0(t) = I — P0(t). Recalling equation (E.17), we have

P0(t)Us(t, tf)P0(tf) = P0(t)Usad(t, tf)P0(tf) + P0(t)Q(t, tf)Usad(t, tf)P0(tf) + O (t—2). (E.24)

Using equation (E.11b) we find Usad(t, tf) P0(tf) = P0(t)e—i ft drS0(t), so that up to a phase P0 (t) Q (t, t')Usad (t, t') P0 (t') = P0 (t) Q (t, t') P0 (t'). However, Q (t, t') (equation (E.18)) contains only off-diagonal terms (since we subtracted the Berry connection in equation (E.5c)), so that P0(t) Q (t, t') P0(t') = 0. Thus, to first order in 1/tf the evolution generated by Hs (t) keeps the ground state decoupled from the excited states, i.e.

P0(t)Us(t, t')P0(t') = Usad(t, t')P0(t') + O (tf—2), (E.25a)

Q0(t)Us(t, tf)Q0(tf) = Usad(t, tf)Q0(tf) + O (tf—2). (E.25b)

Letting r = t — t' > 0 and expanding around t, we then have, using equations (E.1) and (E.3),

P0(t)[I + i rHs(t)][P0(t) — rP0(t)] = [I + i rHsad(t)][P0(t) — rP0(t)] + O (t—2), (E.26a)

Q0(t)[I + i rHs(t)][Q0(t) — r Q0(t)] = [I + i rHad(t)][Q0(t) — r Q0(t)] + O (t—2). (E.26b)

Using the projector properties P0 = P02 P0P0 + P0P0 = P0 and [Hs(t), P0(t)] = 0 (similarly for Q0), this simplifies to

—rP0(t)P0(t) +irHs(t)P0(t) = —rP0(t) +irHsad(t)P0(t) + O (tf—2), (E.27a)

—r Q0(t)Q0(t) + i rHs(t) Q0(t) = —r Q0(t) + i rHad(t) Q0(t) + O (t—2), (E.27b)

and, after dropping the O (tf—2) corrections, to

i[ Hsad (t) — Hs (t)] P0 (t) + P0 (t) P0 (t) = 0, (E.28a)

i[ Hsad (t) — Hs (t)] Q 0 (t) + Q0 (t) Q 0 (t) = 0. (E.28b)

Inserting Q0(t) = I — P0(t) and Q0(t) = — P0(t) into equation (E.28b), and subtracting it from equation (E.28a), we find, using P0P0 + P0P0 = P0 once more, the desired result:

Hsad(t) = Hs(t) +i[P0(t), P0(t)], (E.29)

which, together with equation (E.13), proves equation (E.4b).

33 ¡op Institute of Physics <J>deutsche physikalische Gesellschaft

Appendix F. Short time bound

We wish to bound the error associated with neglecting & in equation (40), i.e. we wish to bound

Ilea ,T)|L = ||^S (t-T, 0)-e T Hs (t Usad (t, 0)|L

= |Usadt (t, 0)e-iTHs(t>Us, (t, t-t)Us(t, 0) - I|U. Using the fact that the operator U(t) = e-iTHs(t)USt(t, t - t) satisfies:

U = -i [Hs(t) - Hs(t - t)] U(t),

we can write the formal solution for U as:

U (t) = I - i Therefore we can bound: Ilea, t)||oo =

i f dt' [Hs(t) - Hs (t -t')] U(t'). Jo

U adt US

(t, 0)Us(t, 0)-i f dt'Usad(t, 0) [Hs(t)-Hs(t -1')] U(t')Us(t, 0) -1 Jo

< min 2, || Q(t, 0)11«, +

f dt' [Hs(t) - Hs (t -1')]

+ O (t-2)

< min {2, —— + -t max ll^Hs(t')IL + O(t-2) , A2tf 2 t'e[t - t ,t ] 1

(F.4a) (F.4b)

(F.4c)

where we used equation (E.17) 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 II Q(t, 0) ||TO (recall appendix E.2). While h of equation (26) is expressed in terms of a matrix element, a more careful analysis (e.g., [10]) would replace this with an operator norm. Thus we shall make the plausible assumption that h ~ tfmaxt/G[t-Tt] ||dt/Hs(t7)||TO, and, dropping the subdominant O (tf-2), we can write

i h t2 h ll&(t »TOI«, ^ min i 2, —— +

A2 tf tf

We can now bound the error term in equation (42b). Let X(t, t) = eiTHs(t)Usad(t, 0), so that Us(t - t, 0) = X(t, t) + e(t, t). We can then write

P to /»TO

/ d tA^(t - T)ps(t)Aa(t)Bafi(t) = / d tXt(t, t)AX(t, t)ps(t)Aa(t)Bafi(t) (F.6a)

p TO p TO

+ / dTXt(t, t)ApS(t, T)ps(t)Aa(t)BaP(t) + / d tXet(t, t)ApX(t, x)~ps(t)Aa(t)Bap(t)

(F.6b) (F.6c)

/ d T0t(t, T)Aß&(t, T)Ps(t)Aa(t)Baß(T). 0

The first term on the rhs of (F.6a) is the approximation we have used in equation (42b). The terms in (F.6b) and (F.6c) can be bounded as follows, using equation (F.5), the unitarity of X, the fact that Hps^ 1 and recalling that \\ Aa ||^= 1. First we assume that equation (12) applies. Then:

dxXf(t, x)Aß©(t, T)PS(t)Aa(t)Baß(x)

dxX©^(t, x)AßX(t, t)Ps(t)Aa(t)Baß(x)

iœ I x h

J dx\\©(t, T)\\œ\Bap(T)\ < min J 2tb+ lb

(F.7a)

pœ pœ

/ dx©f (t ,x) Aß®(t ,x)Ps (t ) Aa(t)Baß(x) ^ dx ||©(t ,x)\\2œ\B^(x)\ oo

< mm\4xB, 2^ + 2^B

xB h ^ xB h

(F.7b)

where in the last inequality we used the fact that if x ^ 2 then [min(2, x)]2 = x2 ^ 2x, and

if x ^ 2 then again [min(2, x)]2 = 2min(2, x) ^ 2x, with x = ^ + ^, in order to avoid having to extend equation (12) to higher values of n .In all, then, the approximation error in

equation (42b) is O[minfe, ^ + ^}].

Next we recall from the discussion in section 2.3 that equation (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 rtr, followed by power-law decay. In this case the terms in (F.6b) can instead be bounded as follows:

/ dxXt(t, x)Aß©(t, t)Ps(t)Aa(t)Baß(x) o

Il f œ

II / dxX(t, x)AßX(t, T)PS(t)Aa(t)Baß(x)

</ dx\©(t,x)\\„,\Baß{x)\ + dx \\©(t,T)\UBaß(T)\

xB h xB h 1 t

< minJ2TB,^ + ^1+2^,

A2 tf tf

in place of (F.7b). A similar modification can be computed for the term in (F.6c). To compute the

order of —, we recall that Ttr ~ / ln(/œc) ^ / (equation (75)), and that tm = ^J22fiJœ~c. Thus

— ~ [Mc ln(^Mc)] It follows that we can safely ignore the — term in equation (F.8) provided

Ttr Ttr

equation (77) is satisfied. The analysis of the term in (F.6c) does not change this conclusion.

Appendix G. Derivation of the Schrodinger picture adiabatic master equation in Lindblad form

starting from equation (53) and performing a transformation back to the schrodinger picture, along with a double-sided adiabatic approximation, yields

Us(t, 0) f d tAß(t - T)ps(t)Aa(t)Baß(T)Ust(t, 0) +h.c. (G.1a)

Aß ab (t) Aa ba (t №b (t )ps (t )Uba (t )Vaß (^ba (t))

+ ^ Ap aa (t) Aa bb (t )naa (t )Ps (t )^b (t )Vap (0) + h.c. (G.1b)

= ^ Labp(t )Ps (t) L I,«« )rap(teba (t)) + ^ Laa p(t )Ps (t) L^^t )Vap(0) +h.c. (G.1c)

a =b ab

similarly, the other terms yield

-Us(t, 0) i d tAa(t)Ap(t - T)Ps(t)Bap(T)Ut(t, 0) (G.2a)

^ L abtU(t) Lab,p (t )Ps (t )rap(toba (t)) - ^ L\ba(t) Laa p(t )Ps (t )^p(0) - h.c. (G.2b)

a=b ab

Using equations (15) and (16) for the spectral-density matrix and its complex conjugation, we

are able to combine the Hermitian conjugate terms, starting from the terms in equation (G.1c):

J2 Taß(toba(t))Labß(t)ps(t)L^t) +h.c. (G.3a)

= ^ \^aß(toba (t)) Lab,ß(t )ps (t) L^t) + Kß^ba (t)) Labta(t )Ps (t) L^t)] (G.3b)

= J2 Yaß(&ba (t)) Labß(t )Ps (t) L ^^t), (G.3c)

where in the last equality we used the freedom to interchange a and p under the summation sign. Likewise, the terms in equation (G.2b) yield

J2 rap(^ba(t))Lla(t)Labp(t)Ps(t) +h.c. (G.4a)

= ^[FapiVba (t)) Llbfa(t) Labp(t )ps (t) + T^ba (t ))Ps (t) L^ hp(t) Lab%a(t)] (G.4b)

= ^ Yap&ba (t)) 1 { Ll^a(t) Lab,p(t), Ps (t)} +i Sap&ba (0)^^(0 Labp(t), Ps (t)],

(G.4c)

where {,} denotes the anticommutator. Combining these results with a similar calculation for the aa and bb terms in equations (G.1c) and (G.2b), finally results in equations (54) and (55).

36 IOP Institute of Physics <J>deutsche physikalische Gesellschaft

Appendix H. Calculations for the spin-boson model

We recall the following basic facts about the bosonic operators appearing in equation (63), where {n} = {nk}k is the set of occupation numbers of all modes:

[bk, 4] = hk>, bk\{n}) = \{nk - 1, n}), b[\{n}) =y/nk + 1 \{nk +1, n}), (H.1a)

HB\{n}) = En}\{n}) = EnkoA \{n}).

(H.lb)

Recalling that the bath operators Bi = kgk (bJ + bk), we proceed to calculate {Bi(t)Bj(0)) = {BJ(t)Bj(0)) assuming that the bath is in a thermal Gibbs state pB = exp (—ftHB) /Z, where Z = Tr (exp (—ftHB)) is the partition function. We begin by writing:

{Bi(t)Bj(0)) = TrB(UB(t, 0)BJUb(t, 0)Bjpb)

= E {{m }|UB (t, 0) BJ l{n }){{n }|Ub (t, 0) Bj |{ p}){{ p}lpb l{m}). (H.2)

{m },{n},{p}

The time evolution operator acting on the eigenstates simply produces a phase, so we focus on the operator Bj's role.

({n}\Bj\{p}) = E (gkVPk + 1^{n},{pk+1,p} + gkypPk^M,{pk-i,. k

Plugging this result in, we find:

(Bi (t) Bj (0)) = — E E ei( E{m}-E{n} )t e-jE{m}

{m},{n},{p} k,k'

(gkSi'V nk + W pk' + 1 ${m },{n ,nk+1}&{n},{p, pk,+1} ${ p}{m}

+ g'kgk'V^^V pk, + 1${m },{n — 1}&{n},{p, pa+1} ${ p}{m} + glkgiV nk +1 Vpk &{m },{n ,nk+1} 8{n},{ p, pk, — !} p}{m} + g'kSkc'V^a*/Pk'&{m},{n,m — 1}${n},{p,pk, — \}${p}{m}) .

The only terms that are non-zero are the middle two, giving us:

(Bt(t)Bj(0)) = EE (mk + 1) (e-ßE{m}ei

( E {m}-E{mk+1} )t + e-ß E{mk+1} e/( E {mk+1}-E {m})

^ . (H.5)

Using the fact that:

TOTO TO 1

E{m} - E{mk+1} = -Ok, Z = n E e-ßmkc = n Y^1

k=1 mk=0 k= 1

_ e-ßOk

we can write:

y (mk + 1) e-ßE{m} = 1 - 1 ^ =-1

^ k ß Z 1 — e

.-ßcOk '

We can simplify our expression to

(B(t)Bj(0)> = £ (e-'«ktg<kgl + em-ßmkgkgj) . (H.8)

We can replace the sum with an integral:

V = V/ dmf (m)S(m - m') = dm J (m), (H.9)

k m j0 j0

where f (m) is a measure of the number of oscillators at frequency m and J (m) is the bath spectral function. Our final result is then:

<B,:(t)Bj (0)> = jH dm (e^'sizl + eimt-ßmg^gi) , (H.10)

where we have assumed that oscillators with the same »k value interact with the ith spin with the same interaction strength, i.e.

g[ = g[r, if »k = »k. (H.11)

Plugging our result into equation (16), we find:

Yj (Vba (t)) = ^ e-ppbfa1 ^a (t) I g La (t) e(Vba (t)) +e-p 1 -a (t) 1 g L (t) j (t) e^ba (t)))

(H.12a)

(»»a(t» = j[TOd»y^ ^ (v^) + g^e-p-p , ».2/»

where e denotes the Heaviside step function.

Appendix I. Derivation of the Ohmic bath correlation function

Here we derive the bath correlation function for the Ohmic oscillator bath with a finite frequency cutoff, equation (68). We start from the simplified expression for the spectral density,

»e-1 » 1

Y(») = 2nng2 »^^-p», (I.1)

and compute the bath correlation function by inverse Fourier transform of equation (16a):

1 f to -B(t) = — d» e-i»TY(»)

2n J-TO

riQ 2 / C0 x ex/(p»c) r to x e-x/(p»c )\

= ™dx e-ixT/p --+ dx e-ixT/p --), (I.2)

p2 U-to 1 - e-x I 1 - e-x / , ' J

where we changed variables to x = p».

The Polygamma function is defined as ty(m\z) = dm^-ln T(z), where T(z) =

/0to e-xxz-1dx is the Gamma function. The Polygamma function may be represented for

^(z) > 0 and m > 0 as \z) = (-1)m+1 /0to ^^ dx, so that in particular, for m = 1, we

f to x e xz

f(1)(z) = J ^Z7dx. (I.3)

We can rewrite equation (I.2)

2 / rœ e-x[-ix/ß+1+1/(ß«c)] pœ ß-x[ix/ß+1/(ß«c)]

B(x)=^{l dx—r-e-x—+1 dx ■ - I (I4a)

= jï tJ3 + 1 + 1J(/œc)) + f(i\i tJ3 + 1J(/a>c))), (I.4b)

which is the desired result.

Appendix J. Derivation of the effective rate equations in the thermal phase

Here we derive the rate equation (83).

J.1. Single qubit

For illustration, let us consider a single qubit such that the Hamiltonian is given by

HS(t) = A(t)ax - B(t)hzaz. (J.1)

The gap at any given time in the evolution is

A(t) = 2yJ A(t )2 + B (t )2 h ^. (J.2)

since there are only two states, there are only three y(<) terms to calculate (assuming g< = g):

2n g2 2 A(t) 2 A(t )e—ftA(t)

Y(0) = -g-. Y(+A(t)) = 2n g 2T-A—At), Y(—A(t)) = 2n g2 ^^' (J.3)

where we have taken the limit <c ^ to for simplicity. Let us consider the gapped phase of our evolution where t/ tf is small. In this phase, B (t) ~ 0, so we can safely assume that the energy eigenstates remain diagonal in the ax basis, with ground state l—) and excited state |+), and do not change such that:

{±|pl±) = d ({±lp|±)) ^ p. (J.4)

The resulting equations for the density matrix elements for small times are then

p++ (t) = l{+laz|—)|2 (Y(—A(t))p—— — Y(A(t))p++) = Y0(A(t)) (e-ftA(t)p—— — p++), (J.5a)

p—— (t) = 1 — p++ (t), (J.5b)

where we have used that {+|az|+) = {—lazl—) = 0, {+|azl—) = 1 and [Hs, ps] & 0. Interpreting equation (J.5) as a rate equation, we see that y(+A) is associated with relaxation from the higher energy state to the lower energy state, and Y(—A(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 p++/p__= e-ftA(0 & 10—n. since the gap

remains relatively large during the gapped phase, the change in the initial thermal state is minimal.

J.2. 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 ax 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 A(t) = 2A(t) (equation (80)), 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 equation (81) and the KMs condition (equation (82)) we find a similar set of equations as in the single qubit case:

Pii & E Yaft(<i0){0lazli){i a|0) (—Pii +e—ft<i0P00) = Yii(<0) (e—ft<i0P00 — Pa) (J.6a)

P00 & 1 — Npu, (J.6b)

where we have again assumed that the initial state is the thermal state and the gap is large.

References

[1] Farhi E, Goldstone J, Gutmann s, Lapan J, Lundgren A and Preda D 2001 Science 292 472

[2] Johnson M W etal 2011 Nature 473 194

[3] Born M and Fock V 1928 Z Phys. 51 165

[4] Kato T 1950 J. Phys. Soc. Japan 5 435

[5] Teufel s 2003 Adiabatic Perturbation Theory in Quantum Dynamics (Berlin: springer)

[6] Jansen s, Ruskai M-B and seiler R 2007 J. Math. Phys. 48 102111

[7] Rigolin G, Ortiz G and Ponce V H 2008 Phys. Rev. A 78 052508

[8] Boixo s, Knill E and somma R D 2009 Quantum Inform. Comput. 9 833

[9] Boixo s, Knill E and somma R D 2010 arXiv:1005.3034

[10] Lidar D A, Rezakhani A T and Hamma A 2009 J. Math. Phys. 50 102106

[11] Boixo s and somma R D 2010 Phys. Rev. A 81 032308

[12] Wiebe N and Babcock N s 2012 New J. Phys. 14 013024

[13] Davies E B and spohn H 1978 J. Stat. Phys. 19 511

[14] Childs A M, Farhi E and Preskill J 2001 Phys. Rev. A 65 012322

[15] sarandy M s and Lidar D A 2005 Phys. Rev. A 71 012331

[16] Garrison J and Wright E 1988 Phys. Lett. A 128 177

[17] Yi X X, Tong D M, Kwek L C and Oh C H 2007 J. Phys. B: At. Mol. Opt. Phys. 40 281

[18] Huang X L, Yi X X, Wu C, Feng X L, Yu s X and Oh C H 2008 Phys. Rev. A 78 062114

[19] Joye A 2007 Commun. Math. Phys. 275 139

[20] Oreshkov O and Calsamiglia J 2010 Phys. Rev. Lett. 105 050503

[21] Thunstrom P, Aberg J and sjoqvist E 2005 Phys. Rev. A 72 022328

[22] Alicki R, Lidar D A and Zanardi P 2006 Phys. Rev. A 73 052311

[23] Pekola J P, Brosco V, Mottonen M, solinas P and shnirman A 2010 Phys. Rev. Lett. 105 030401

[24] Kamleitner I and shnirman A 2011 Phys. Rev. B 84 235140

[25] O'Hara M J and O'Leary D P 2008 Phys. Rev. A 77 042319

[26] Amin M H s, Love P J and Truncik C J s 2008 Phys. Rev. Lett. 100 060503

[27] salmilehto J, solinas P, Ankerhold J and Mottonen M 2010 Phys. Rev. A 82 062112

[28] Dridi G, Guerin s, Jauslin H R, Viennot D and Jolicard G 2010 Phys. Rev. A 82 022109

[29] Lindblad G 1976 Commun. Math. Phys. 48 119

[30] Breuer H-P and Petruccione F 2002 The Theory of Open Quantum Systems (Oxford: Oxford University Press)

[31] Amin M H s, Truncik C J s and Averin D V 2009 Phys. Rev. A 80 022303

[32] de Vega I, Banuls M C and Perez A 2010 New J. Phys. 12 123010

[33] Yung M 2008 Thermal noise on adiabatic quantum computation arXiv:0807.4819

[34] stelmachovic P and Buzek V 2001 Phys. Rev. A 64 062106

[35] Blum K 2010 Physics of Atoms and Molecules: Density Matrix Theory and Applications (Berlin: springer)

[36] Leggett A J, Chakravarty s, Dorsey A T, Fisher M P A, Garg A and Zwerger W 1987 Rev. Mod. Phys. 59 1

[37] Alicki R 2004 Open Syst. Inform. Dyn. 11 53

[38] Weisstein E W Lambert W-Function MathWorld—A Wolfram Web Resource http://mathworld.wolfram.com/

LambertW-Function.html

[39] Coughran Jr W M, Grosse E and Rose D 1983 IEEE Trans. Electron Devices 30 1207

[40] Bank R, Coughran Jr W M, Fichtner W, Grosse E, Rose D and smith R 1985 IEEE Trans. Electron Devices

32 1992

[41] Larson J 2012 Phys. Rev. Lett. 108 033601

[42] Patane D, Amico L, silva A, Fazio R and santoro G E 2009 Phys. Rev. B 80 024302

[43] Reichardt B W 2004 Proc. 36th Annu. ACMSymp. on Theory of Computing (Chicago, IL) (New York: ACM)

[44] Bhatia R 1997 Graduate Texts in Mathematics vol 169 Matrix Analysis (New York: springer)

[45] Watrous J 2004 Lecture Notes on Theory of Quantum Information Online at http://www.cs.uwaterloo.ca/

watrous/lecture-notes/701/all.pdf

[46] Lidar D A, Zanardi P and Khodjasteh K 2008 Phys. Rev. A 78 012308

[47] Mathews J H and Howell R W 2012 Complex Analysis: for Mathematics and Engineering 6th edn (sudbury,

MA: Jones and Bartlett)

[48] Haag R, Hugenholtz N M and Winnink M 1967 Commun. Math. Phys. 5 215

[49] Marzlin K-P and sanders B C 2004 Phys. Rev. Lett. 93 160408

[50] sarandy M s, Wu L-A and Lidar D A 2004 Quantum Inform. Process. 3 331

[51] Tong D M, singh K, Kwek L C and Oh C H 2007 Phys. Rev. Lett. 98 150402

[52] Amin M H s 2009 Phys. Rev. Lett. 102 220401

[53] Comparat D 2011 Phys. Rev. Lett. 106 138902

[54] Yong T 2012 Commun. Theor. Phys. 57 343

[55] Avron J E, seiler R and Yaffe L G 1987 Commun. Math. Phys. 110 33