Quantum adiabatic Markovian master equations

We consider a general system–bath Hamiltonian

$$\begin{equation} H(t) = H_{\mathrm{S}}(t) + H_{\mathrm{B}} + H_{\mathrm{I}} , \end{equation} \tag{ 1 }$$

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:

$$\begin{equation} H_{\mathrm{I}} = {g} \sum_{\alpha} A_{\alpha} \otimes B_{\alpha} , \end{equation} \tag{ 2 }$$

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 $\dot {\rho } = -\mathrm {i}[H,\rho (t)]$, where we have assumed units such that ℏ = 1.

Let $U_{\mathrm {S}}(t,t')=T_{+}\exp [-\mathrm {i}\int _{t'}^{t}\mathrm {d}\tau \ H_{\mathrm {S}}(\tau )]$ and $U_{\mathrm {B}}(t,t')=\exp (-\mathrm {i}(t-t')H_{\mathrm {B}})$ denote the free system and bath unitary evolution operators (T₊ denotes time ordering), and define U₀(t,t') = U_S(t,t')⊗U_B(t,t'). Likewise, let $U(t,t')=T_{+}\exp [-\mathrm {i}\int _{t'}^{t}\mathrm {d}\tau \ H(\tau )]$ 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 $\tilde {U}(t,0)=U_{0}^{\dagger }(t,0)U(t,0)$, which, together with the interaction picture density matrix $\tilde {\rho }(t)=U_{0}^{\dagger }(t,0)\rho (t)U_{0}(t,0)$, satisfies

$$\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\tilde{U}(t,0) =-\mathrm{i}{\tilde{H}_{\mathrm{I}}(t)\tilde{U}(t,0)},\quad \tilde{U}(0,0)={\mathbbm 1} , \end{equation} \tag{ 3a }$$

$$\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\tilde{\rho}(t) =-\mathrm{i}[\tilde{H}_{\mathrm{I}}(t),\tilde{\rho}(t)],\quad \tilde{\rho}(0)={\rho }(0). \end{equation} \noindent \tag{ 3b }$$

We restrict the use of tilde variables to refer to variables in the interaction picture. The time-dependent interaction picture Hamiltonian $\tilde {H}_{\mathrm {I}}(t)$ is related to its Schrödinger picture counterpart via:

$$\begin{equation} \tilde{H}_{\mathrm{I}}(t)=U_{0}^{\dagger }(t,0){H_{\mathrm{I}}}U_{0}(t,0)= g \sum_{\alpha }A_{\alpha }(t)\otimes B_{\alpha }(t), \end{equation} \noindent \tag{ 4 }$$

where we have defined the time-dependent system and bath operators:

$$\begin{equation} A_{\alpha }(t) =U_{\mathrm{S}}^{\dagger }(t,0)A_{\alpha }U_{\mathrm{S}}(t,0), \end{equation} \tag{ 5a }$$

$$\begin{equation} B_{\alpha }(t) =U_{\mathrm{B}}^{\dagger }(t,0)B_{\alpha }U_{\mathrm{B}}(t,0). \end{equation} \noindent \tag{ 5b }$$

$$\begin{equation} \fl \tilde{\rho}_{\mathrm{S}}(t) \equiv \mathrm{Tr}_{\mathrm{B}}[\tilde{\rho}(t)]= U_{\mathrm{S}}^\dagger(t,0) \mathrm{Tr}_{\mathrm{B}}( U_{\mathrm{B}}^\dagger (t,0)\rho(t) U_{\mathrm{B}}(t,0)) U_{\mathrm{S}}(t,0) \end{equation} \tag{ 6a }$$

$$\begin{equation} \fl\hphantom{\tilde{\rho}_{\mathrm{S}}(t)}= U_{\mathrm{S}}^\dagger(t,0) \mathrm{Tr}_{\mathrm{B}}( U_{\mathrm{B}}(t,0) U_{\mathrm{B}}^\dagger (t,0)\rho(t)) U_{\mathrm{S}}(t,0) =U_{\mathrm{S}}^{\dagger}(t,0)\rho_{\mathrm{S}}(t)U_{\mathrm{S}}(t,0), \end{equation} \noindent \tag{ 6b }$$

Writing the formal solution of equation (3b) as

$$\begin{equation} \tilde{\rho}(t)=\tilde{\rho}(0)-\mathrm{i}\int_{0}^{t}\mathrm{d}\tau \ [\tilde{H}_{\mathrm{I}}(\tau ) , \tilde{\rho}(\tau )], \end{equation} \noindent \tag{ 7 }$$

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 $\tilde {\rho }_{\mathrm {S}}(t)=\mathrm {Tr}_{\mathrm {B}}\tilde {\rho }(t)$

$$\begin{equation} \fl \frac{\mathrm{d}}{\mathrm{d}t}\tilde{\rho}_{\mathrm{S}}(t)=-\mathrm{i}\,\mathrm{Tr}_{\mathrm{B}}[ \tilde{H}_{\mathrm{I}}(t) ,\tilde{\rho}(0)] -\mathrm{Tr}_{\mathrm{B}}\left[ \tilde{H}_{\mathrm{I}}(t),\int_{0}^{t}\mathrm{d}\tau\, [ \tilde{H}_{\mathrm{I}}(t-\tau ),\tilde{\rho}(t-\tau )]\right], \end{equation} \noindent \tag{ 8 }$$

We make the standard Born approximation assumption, that we can decompose the density matrix as $\tilde {\rho }(t)=\tilde {\rho }_{\mathrm {S}}(t)\otimes {\rho }_{\mathrm {B}}+\chi (t)$ 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, 34].⁷ Thus the equation of motion reduces to:

$$\begin{eqnarray} &&\fl \frac{\mathrm{d}}{\mathrm{d}t}\tilde{\rho}_{\mathrm{S}}(t)={g^{2}}\sum_{\alpha ,\beta }\int_{0}^{t}\mathrm{d}\tau \left( A_{\beta }(t\!-\!\tau )\tilde{\rho}_{\mathrm{S}}(t\!-\!\tau )A_{\alpha }(t)\!-\!A_{\alpha }(t)A_{\beta }(t\!-\!\tau )\tilde{\rho}_{\mathrm{S}}(t\!-\!\tau )\right) \mathcal{B}_{\alpha \beta }(t,t\!-\!\tau )+\mathrm{h.c.},\nonumber\\ \end{eqnarray} \tag{ 9 }$$

where we defined the two-point correlation functions:

$$\begin{equation} \mathcal{B}_{\alpha \beta }(t,t-\tau )\equiv \langle B_{\alpha }(t)B_{\beta }(t-\tau )\rangle =\mathrm{Tr}\left[ B_{\alpha }(t)B_{\beta }(t-\tau ){\rho}_{\mathrm{B}}\right] . \end{equation} \tag{ 10 }$$

In equation (9), we have assumed for simplicity (but without loss of generality) that $\langle B_{\alpha }\rangle _{0}=\mathrm {Tr} [ B_{\alpha }(t)\tilde {\rho }_{\mathrm {B}}(0) ] =0$, so that the inhomogeneous term in equation (8) vanishes. Since we assumed that the bath state ρ_B is stationary, the correlation function is homogeneous in time:

$$\begin{equation} \fl \mathcal{B}_{\alpha \beta }(t,t-\tau )=\langle B_{\alpha }(t)B_{\beta }(t-\tau )\rangle =\langle B_{\alpha }(\tau )B_{\beta }(0)\rangle =\langle B_{\beta }(0)B_{\alpha }(\tau )\rangle^*=\mathcal{B}_{\alpha \beta }(\tau ,0). \end{equation} \tag{ 11 }$$

For notational simplicity we will denote $\mathcal {B}_{\alpha \beta }(\tau ,0)$ by $\mathcal {B}_{\alpha \beta }(\tau )$ 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. $|\mathcal {B}_{\alpha \beta }(\tau )|\sim \exp (-\tau /\tau _{\mathrm {B}})$. More precisely, we shall require throughout that

$$\begin{equation} \int_{0}^{\infty }\mathrm{d}\tau \ \tau ^{n}|B_{\alpha \beta }(\tau )|\sim \tau_{\mathrm{B}}^{n+1} ,\qquad {n\in \{0,1,2\}}. \end{equation} \tag{ 12 }$$

As we show in appendix B, if τ_B ≪ 1/g, then we can apply the Markov approximation to each of the four summands in equation (9), i.e., replace $\tilde {\rho }_{\mathrm {S}}(t-\tau )$ by $\tilde {\rho }_{\mathrm {S}}(t)$:

$$\begin{equation} \fl \int_{0}^{t}\mathrm{d}\tau \left( \cdots \tilde{\rho}_{\mathrm{S}}(t-\tau )\cdots \right) \mathcal{B}_{\alpha \beta }(\tau )\approx \int_{0}^{\infty }\mathrm{d}\tau \left( \cdots \tilde{\rho}_{\mathrm{S}}(t)\cdots \right) \mathcal{B}_{\alpha \beta }(\tau )+{ O({\tau_{\mathrm{B}}^{3}g^2})}, \end{equation} \tag{ 13 }$$

where $\left ( \cdots \right )$ refers to the products of A_α and A_β operators in equation (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 equation (12) the integral on the rhs of equation (13) is of O(τ_B), so that the relative magnitude of the two terms is O[(gτ_B)²]. 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.

In computing the terms appearing in equation (13) we are faced with integrals of the form $\int _{0}^{\infty }\mathrm {d}\tau A_{\beta }(t-\tau )\tilde {\rho }_{\mathrm {S}}(t)A_{\alpha }(t)\mathcal {B}_{\alpha \beta }(\tau )$ and $\int _{0}^{\infty }\mathrm {d}\tau A_{\alpha }(t)A_{\beta }(t-\tau )\tilde {\rho }_{\mathrm {S}}(t) \mathcal {B}_{\alpha \beta }(\tau )$. Our goal is to express these integrals in terms of the spectral-density matrix

$$\begin{equation} \Gamma_{\alpha \beta }(\omega )\equiv \int_{0}^{\infty }\mathrm{d}\tau \mathrm{e}^{\mathrm{i}\omega \tau }\mathcal{B}_{\alpha \beta }(\tau ), \end{equation} \tag{ 14 }$$

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,

$$\begin{equation} \Gamma_{\alpha \beta}(\omega) = \frac{1}{2} \gamma_{\alpha \beta} (\omega) + \mathrm{i}\, S_{\alpha \beta}(\omega), \end{equation} \noindent \tag{ 15 }$$

where we show in appendix C that γ and S are given by

$$\begin{eqnarray} &&\gamma_{\alpha \beta}(\omega) = \int_{-\infty}^\infty \mathrm{d} \tau \mathrm{e}^{\mathrm{i} \omega \tau} \mathcal{B}_{\alpha \beta}(\tau ) = \gamma^*_{\beta \alpha}(\omega) {, } \end{eqnarray} \tag{ 16a }$$

$$\begin{eqnarray} &&S_{\alpha\beta}(\omega) = \int_{-\infty}^{\infty} \frac{\mathrm{d}\omega'}{2\pi} \gamma_{\alpha \beta}(\omega') \mathcal{P}\left( \frac{1}{\omega - \omega'} \right) = S^*_{\beta \alpha}(\omega) . \end{eqnarray} \noindent \tag{ 16b }$$

$\rho _{\mathrm {B}}=\mathrm {e}^{-\beta H_{\mathrm {B}}}/\mathcal {Z}$

$$\begin{equation} \langle B_a(\tau)B_b(0)\rangle = \langle B_b(0)B_a(\tau+\mathrm{i}\beta)\rangle . \end{equation} \noindent \tag{ 17 }$$

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

$$\begin{equation} \gamma_{ab}(-\omega) = \mathrm{e}^{-\beta\omega}\gamma_{ba}(\omega). \end{equation} \noindent \tag{ 18 }$$

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 $\vert \mathcal {B}_{ab}(-\tau ) \vert = \vert \mathrm {Tr}[\rho _{\mathrm {B}}U_{\mathrm {B}}(\tau )B_{a}U_{\mathrm {B}}^{\dagger }(\tau )B_{b}]\kern.5pt \vert\kern.5pt = \kern.5pt\vert\kern.5pt \mathrm {Tr}\kern.5pt[B_{a}\kern.5ptU_{\mathrm {B}}^{\dagger }\kern.5pt(\tau )B_{b}\kern.5ptU_{\mathrm {B}}\kern.5pt(\tau )\kern.5pt\rho _{\mathrm {B}}\kern.5pt]\kern.5pt \vert = \vert \kern.5pt\langle B_{a}(0)B_{b}(\tau )\rangle \vert = \vert \langle B_{b}(\tau )B_{a}(0)\rangle \vert = \vert \mathcal {B}_{ba}(\tau ) \vert$, where we used equation (11). Now assume that equation (12) would have to hold for all n. It would follow that

$$\begin{eqnarray} \fl \left\vert \frac{\mathrm{d}^{n}}{\mathrm{d}\omega ^{n}}\gamma_{ab}(\omega )\right\vert_{\omega =0}&=&\left\vert \int_{-\infty }^{\infty }\left. \tau^{n}\mathrm{e}^{\mathrm{i}\omega \tau }\mathcal{B}_{ab}(\tau )\mathrm{d}\tau \right\vert_{\omega=0}\right\vert=\left\vert \int_{-\infty }^{\infty }\tau ^{n}\mathcal{B}_{ab}(\tau )\mathrm{d}\tau \right\vert\nonumber\\ \fl& \leqslant &\int_{0}^{\infty }\tau ^{n}\left(\left\vert \mathcal{B}_{ba}(\tau )\right\vert +\left\vert \mathcal{B}_{ab}(\tau )\right\vert \right) \mathrm{d}\tau \sim 2\tau_{\mathrm{B}}^{n+1}\quad \forall n\in \{0,1,\ldots \}. \end{eqnarray} \noindent \tag{ 19 }$$

Thus all derivatives of γ_ab(ω) would have to be finite at ω = 0.

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

$$\begin{eqnarray} &&\left[ \frac{\mathrm{d}}{\mathrm{d}(-\omega) }\gamma_{ab}(-\omega ) \right]_{\omega \geqslant 0}= \left. \left[ \beta \mathrm{e}^{-\beta \omega}\gamma_{ba}(\omega) - \mathrm{e}^{-\beta \omega}\frac{\mathrm{d}}{\mathrm{d}\omega }\gamma_{ba}(\omega)\right]\right\vert_{\omega \geqslant 0}, \end{eqnarray} \noindent \tag{ 20 }$$

so that in the limit as ω approaches zero from below or above,

$$\begin{eqnarray} &&\gamma_{ab}'(0_-) = \beta \gamma_{ba}(0) - \gamma_{ba}'(0_+) . \end{eqnarray} \noindent \tag{ 21 }$$

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 equation (21), that

$$\begin{equation} 2 \gamma'_{aa}(0) = \beta \gamma_{aa}(0) . \end{equation} \noindent \tag{ 22 }$$

This conclusion can be extended to the entire γ matrix by diagonalizing it first. When equation (22) is not satisfied γ''_aa(0) diverges, so that equation (19) does not hold except for n∈{0,1}. A simple example of this is γ_ab(ω) = c > 0 (constant) for ω ⩾ 0. Another example is a super-Ohmic spectral density γ_ab(ω) = ω²/(1 − e^−βω) for ω ⩾ 0 and γ_ab(ω) = ω²/(e^−βω − 1) for ω ⩽ 0. Both examples violate equation (19) for n ⩾ 2. Note, moreover, that when this happens, it follows from equation (19) that $\int _{0}^{\infty }\tau ^{2}\left (\left \vert \mathcal {B}_{ab}(-\tau )\right \vert +\left \vert \mathcal {B}_{ab}(\tau )\right \vert \right ) \mathrm {d}\tau$ is divergent, so that we must conclude that $|\mathcal {B}_{ab}(\tau )| \gtrsim 1/\tau ^3$ for sufficiently large |τ|, 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 γ'_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 γ_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 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 τ_tr. In this case we find |B_αβ(τ)| ∼ (τ/τ_M)⁻², where τ_M is a time-scale associated with the onset of non-Markovian effects, and hence we have to relax equation (12), and replace it with

$$\begin{equation} \int_{0}^{\tau_{\mathrm{tr}}}\mathrm{d}\tau \ \tau ^{n}|B_{\alpha \beta }(\tau )| \sim \tau_{\mathrm{B}}^{n+1},\quad {n\in \{0,1,\ldots\}}, \end{equation} \tag{ 23a }$$

$$\begin{equation} \int_{\tau_{\mathrm{tr}}}^{\infty}\mathrm{d}\tau \ |B_{\alpha \beta }(\tau )| \sim \int_{\tau_{\mathrm{tr}}}^{\infty}\mathrm{d}\tau \ (\tau/\tau_{\mathrm{M}})^{-2} = \tau_{\mathrm{M}}^2/\tau_{\mathrm{tr}}. \end{equation} \tag{ 23b }$$

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 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:

$$\begin{eqnarray} &&U_{\mathrm{S}}(t,t') =U_{\mathrm{S}}^{\mathrm{ad}}(t,t')+O\left( \frac{h}{ \Delta ^{2}t_{\mathrm{f}}}\right) {, } \end{eqnarray} \tag{ 37a }$$

$$\begin{eqnarray} &&U_{\mathrm{S}}^{\mathrm{ad}}(t,t') =\sum_{a}|\varepsilon_{a}(t)\rangle \langle \varepsilon_{a}(t')|\mathrm{e}^{-\mathrm{i}\mu_{a}(t,t')} { ,} \end{eqnarray} \tag{ 37b }$$

$$\begin{equation} \mu_{a}(t,t')= \int_{t'}^{t}\mathrm{d}\tau \left[\varepsilon_{a}(\tau )-\phi_a(\tau)\right], \end{equation} \tag{ 38 }$$

where $\phi _a(t) = \mathrm {i}\langle \varepsilon _{a}(t)|\dot {\varepsilon }_{a}(t)\rangle$ is the Berry connection. The correction term $O\left ( \frac {h}{\Delta ^{2}t_{\mathrm {f}}}\right )$ 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 A_β(t − τ) = U^†_S(t − τ,0)A_βU_S(t − τ,0) with an appropriate adiabatic approximation, which will allow us to take this operator outside the integral. To see how, note first that

$$\begin{equation} U_{\mathrm{S}}(t-\tau ,0)=U_{\mathrm{S}}(t-\tau ,t)U_{\mathrm{S}}(t,0)=U_{\mathrm{S}}^{\dagger }(t,t-\tau )U_{\mathrm{S}}(t,0) . \end{equation} \tag{ 39 }$$

We now make two approximations: firstly, as per equation (37a) we replace U_S(t,0) by U^ad_S(t,0); secondly, 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 $\mathcal {B}_{\alpha \beta }(\tau )$ inside the integrals we are concerned with. Thus, we write

$$\begin{equation} U_{\mathrm{S}}(t-\tau ,0)=\mathrm{e}^{\mathrm{i}\,\tau H_{\mathrm{S}}(t)}U_{\mathrm{S}}^{\mathrm{ad}}(t,0)+\Theta (t,\tau ), \end{equation} \tag{ 40 }$$

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

$$\begin{eqnarray} &&\mu_{ba}(t,t') \equiv \mu_{b}(t,t')-\mu_{a}(t,t') , \end{eqnarray} \tag{ 41a }$$

$$\begin{eqnarray} &&\Pi_{ab}(t) \equiv |\varepsilon_{a}(t)\rangle \langle \varepsilon_{b}(t)| . \end{eqnarray} \tag{ 41b }$$

$$\begin{equation} \fl \int_{0}^{\infty }\mathrm{d}\tau A_{\beta }(t-\tau )\tilde{\rho}_{\mathrm{S}}(t)A_{\alpha }(t)\mathcal{B}_{\alpha \beta }(\tau ) \end{equation} \tag{ 42a }$$

$$\begin{equation} \fl \quad \approx\int_{0}^{\infty }\mathrm{d}\tau U_{\mathrm{S}}^{\mathrm{ad}\dagger }(t,0)\mathrm{e}^{-\mathrm{i}\tau H_{\mathrm{S}}(t)}A_{\beta}\,\mathrm{e}^{\mathrm{i}\tau H_{\mathrm{S}}(t)}U_{\mathrm{S}}^{\mathrm{ad}}(t,0)\tilde{\rho}_{\mathrm{S}}(t)A_{\alpha }(t)\mathcal{B}_{\alpha \beta}(\tau) \end{equation} \tag{ 42b }$$

$$\begin{equation} \fl \quad =\int_{0}^{\infty }\mathrm{d}\tau \sum_{ab}\mathrm{e}^{-\mathrm{i}\mu_{ba}(t,0)}|\varepsilon_{a}(0)\rangle \langle \varepsilon_{a}(t)|\mathrm{e}^{-\mathrm{i}\tau H_{\mathrm{S}}(t)}A_{\beta}\mathrm{e}^{\mathrm{i}\tau H_{\mathrm{S}}(t)}|\varepsilon_{b}(t)\rangle \langle \varepsilon_{b}(0)|\tilde{\rho}_{\mathrm{S}}(t)A_{\alpha }(t)\mathcal{B}_{\alpha \beta }(\tau ) \end{equation} \tag{ 42c }$$

$$\begin{equation} \fl \quad =\sum_{ab}\mathrm{e}^{-\mathrm{i}\mu_{ba}(t,0)}|\varepsilon_{a}(0)\rangle \langle\varepsilon_{a}(t)|A_{\beta }|\varepsilon_{b}(t)\rangle \langle \varepsilon_{b}(0)|\tilde{\rho}_{\mathrm{S}}(t)A_{\alpha }(t)\int_{0}^{\infty}\mathrm{d}\tau\mathrm{e}^{\mathrm{i}\tau \lbrack \varepsilon_{b}(t)-\varepsilon_{a}(t)]}\mathcal{B}_{\alpha \beta }(\tau ) , \end{equation} \tag{ 42d }$$

$O[\tau _{\mathrm {B}}\min \{1,\frac {h}{\Delta ^2 t_{\mathrm {f}}} + \frac {\tau _{\mathrm {B}}^2\,h}{t_{\mathrm {f}}}\}]$

$\frac {h}{\Delta ^2 t_{\mathrm {f}}}$

$\frac {\tau _{\mathrm {B}}^2\,h}{t_{\mathrm {f}}}$

Thus, to the same level of approximation

$$\begin{eqnarray} &&\fl \int_{0}^{\infty }\mathrm{d}\tau A_{\beta }(t-\tau )\tilde{\rho}_{\mathrm{S}}(t)A_{\alpha }(t) \mathcal{B}_{\alpha \beta }(\tau )\approx \sum_{ab}\mathrm{e}^{-\mathrm{i}\mu_{ba}(t,0)}A_{\beta ab}(t)\Pi_{ab}(0)\tilde{\rho}_{\mathrm{S}}(t)A_{\alpha }(t)\Gamma_{\alpha \beta}(\omega_{ba}(t)),\nonumber\\ \end{eqnarray} \tag{ 43 }$$

where

$$\begin{equation} A_{\alpha ab}(t)\equiv \langle \varepsilon_{a}(t)|A_{\alpha }|\varepsilon_{b}(t)\rangle = A^*_{\alpha ba}(t) , \end{equation} \tag{ 44 }$$

and Γ_αβ(ω_ba(t)) is the spectral-density matrix defined in equation (14). Similarly, we have for the other term

$$\begin{equation} \fl \int_{0}^{\infty }\mathrm{d}\tau A_{\alpha }(t)A_{\beta }(t-\tau )\tilde{\rho}_{\mathrm{S}}(t) \mathcal{B}_{\alpha \beta }(\tau )\approx\sum_{ab}\mathrm{e}^{-\mathrm{i}\mu_{ba}(t,0)}A_{\beta ab}(t)A_{\alpha }(t)\Pi_{ab}(0)\tilde{\rho}_{\mathrm{S}}(t)\Gamma_{\alpha \beta }(\omega_{ba}(t)). \end{equation} \tag{ 45 }$$

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:

$$\begin{equation} \fl \frac{\mathrm{d}}{\mathrm{d}t}\tilde{\rho}_{\mathrm{S}}(t)=g^{2}\sum_{ab}\mathrm{e}^{-\mathrm{i}\mu_{ba}(t,0)}\sum_{\alpha \beta }\Gamma_{\alpha \beta }(\omega_{ba}(t))A_{\beta ab}(t)\left[ \Pi_{ab}(0)\tilde{\rho}_{\mathrm{S}}(t),A_{\alpha }(t) \right] +\mathrm{h.c}. \end{equation} \tag{ 46 }$$

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

$$\begin{eqnarray} &&\fl \frac{\mathrm{d}}{\mathrm{d}t}\tilde{\rho}_{\mathrm{S}}(t)=g^{2}\sum_{abcd}\mathrm{e}^{-\mathrm{i}\left[ \mu_{dc}(t,0)+\mu_{ba}(t,0)\right] }\sum_{\alpha \beta }\Gamma_{\alpha \beta }(\omega_{ba}(t))A_{\alpha cd}(t)A_{\beta ab}(t)\left[ \Pi_{ab}(0)\tilde{ \rho}_{\mathrm{S}}(t),\Pi_{cd}(0)\right] +\mathrm{h.c}.\nonumber\\ \end{eqnarray} \tag{ 47 }$$

It is convenient to transform back into the Schrödinger picture. Using $\tilde {\rho }_{\mathrm {S}}(t)= U_{\mathrm {S}}^{\dagger}(t,0)\rho _{\mathrm {S}}(t)U_{\mathrm {S}}(t,0)$ (equation (6)) implies that $\frac {\mathrm {d}}{\mathrm {d}t}\tilde {\rho }_{\mathrm {S}}(t)=U_{\mathrm {S}}^{\dagger}(t,0)\left ( \mathrm {i}[H_{\mathrm {S}}(t), \rho _{\mathrm {S}}(t)]+\frac {\mathrm {d}}{\mathrm {d}t}{\rho }_{\mathrm {S}}(t)\right ) U_{\mathrm {S}}(t,0)$. Hence, using equation (46), we find the one-sided adiabatic Schrödinger picture master equation

$$\begin{equation} \fl \frac{\mathrm{d}}{\mathrm{d}t}{\rho }_{\mathrm{S}}(t)=-\mathrm{i}[H_{\mathrm{S}}(t),\rho_{\mathrm{S}}(t)]+g^{2}\sum_{ab}\sum_{\alpha \beta }\Gamma_{\alpha \beta }(\omega_{ba}(t))\left[ L_{ab,\beta }'(t)\rho_{\mathrm{S}}(t),A_{\alpha }\right] + \mathrm{h.c.}, \end{equation} \tag{ 48 }$$

where

$$\begin{equation} L_{ab,\beta }'(t)=\mathrm{e}^{-\mathrm{i}\mu_{ba}(t,0)}A_{\beta ab}(t)U_{\mathrm{S}}(t,0)\Pi_{ab}(0)U_{\mathrm{S}}^{\dagger}(t,0). \end{equation} \tag{ 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 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 double-sided adiabatic Schrödinger picture master equation

$$\begin{equation} \fl \frac{\mathrm{d}}{\mathrm{d}t}{\rho }_{\mathrm{S}}(t)=-\mathrm{i}\left[ H_{\mathrm{S}}(t),\rho_{\mathrm{S}}(t)\right] +g^{2}\sum_{\alpha \beta }\sum_{ab}\Gamma_{\alpha \beta }(\omega_{ba}(t)) \left[ L_{ab,\beta }(t)\rho_{\mathrm{S}}(t),A_{\alpha }\right] +\mathrm{h.c.}, \end{equation} \tag{ 50 }$$

where

$$\begin{equation} L_{ab,\alpha }(t)\equiv A_{\alpha ab}(t)|{\varepsilon_a(t)}\rangle\! \langle{\varepsilon_b(t)}| = L^\dagger_{ba,\alpha }(t). \end{equation} \tag{ 51 }$$

Comparing to equation (25), the first term in equations (48) and (50) is $\mathcal {L}_{\mathrm {uni}}(t)$, while the second is $\mathcal {L}_{\mathrm {diss} }^{\mathrm {ad}}(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.

In order to arrive at a master equation in Lindblad form, we can perform a RWA. To do so, let us revisit the t → ∞ 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 A_α(t) ≈ U^ad†_S(t,0)A_αU^ad_S(t,0). This yields

$$\begin{equation} \fl \int_{0}^{t}\mathrm{d}\tau A_{\beta }(t-\tau )\tilde{\rho}_{\mathrm{S}}(t)A_{\alpha }(t) \mathcal{B}_{\alpha \beta }(\tau ) \end{equation}\noindent \tag{ 52a }$$

$$\begin{eqnarray} &&\approx \int_{0}^{t}\mathrm{d}\tau \sum_{abcd}\mathrm{e}^{-\mathrm{i}\left[\mu_{ba}(t,0)+\mu_{dc}(t,0)\right]}|\varepsilon_{a}(0)\rangle \langle \varepsilon_{a}(t)|A_{\beta}|\varepsilon_{b}(t)\rangle \langle \varepsilon_{b}(0)| \tilde{\rho}_{\mathrm{S}}(t) |{\varepsilon_c(0)}\rangle\nonumber\\ && \times\langle \varepsilon_c(t) |A_{\alpha }|{\varepsilon_d(t)}\rangle \langle \varepsilon_d(0)|\mathrm{e}^{\mathrm{i}\,\tau\omega_{ba}(t)}\mathcal{B}_{\alpha \beta }(\tau ) . \end{eqnarray}\noindent \tag{ 52b }$$

$\mu _{dc}(t,0)+\mu _{ba}(t,0) = \int _0^t \mathrm {d}\tau \left [\omega _{dc}(\tau )+\omega _{ba}(\tau )\!-\!(\phi _d(\tau )\!-\!\phi _c(\tau ))+(\phi _b(\tau )\!-\!\phi _a(\tau ))\right ]$

$$\begin{eqnarray} &&\fl \int_{0}^{t}\mathrm{d}\tau A_{\beta }(t-\tau )\tilde{\rho}_{\mathrm{S}}(t)A_{\alpha }(t) \mathcal{B}_{\alpha \beta }(\tau )\approx \sum_{a {\neq} b}A_{\beta ab}(t)A_{\alpha ba}(t)\Pi_{ab}(0) \tilde{\rho}_{\mathrm{S}}(t)\Pi_{ba}(0) \Gamma_{\alpha\beta}(\omega_{ba}(t))\nonumber\\ &&+ \sum_{ab}A_{\beta aa}(t)A_{\alpha bb}(t)\Pi_{aa}(0)\tilde{\rho}_{\mathrm{S}}(t)\Pi_{bb}(0) \Gamma_{\alpha\beta}(0). \end{eqnarray} \tag{ 53 }$$

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 equation (53) at the Schrödinger picture adiabatic master equation in Lindblad form:

$$\begin{equation} \fl \skew3\dot{\rho}_{\mathrm{S}}(t) = - \mathrm{i} \left[ H_{\mathrm{S}}(t) + H_{\mathrm{LS}}(t), \rho_{\mathrm{S}}(t) \right] \end{equation} \tag{ 54a }$$

$$\begin{equation} \fl \qquad + \sum_{\alpha \beta} \sum_{a {\neq} b} \gamma_{\alpha \beta}(\omega_{ba}(t)) \left[ L_{ab, \beta}(t) \rho_{\mathrm{S}}(t) L^\dagger_{a b, \alpha}(t) - \frac{1}{2} \left\{ L^\dagger_{a b, \alpha}(t) L_{a b,\beta}(t), \rho_{\mathrm{S}}(t) \right\}\right] \end{equation} \tag{ 54b }$$

$$\begin{equation} \fl \qquad+ \sum_{\alpha \beta} \sum_{ab} \gamma_{\alpha \beta}(0) \left[ L_{aa, \beta}(t) \rho_{\mathrm{S}}(t) L^\dagger_{b b, \alpha}(t) - \frac{1}{2} \left\{ L^\dagger_{a a, \alpha}(t) L_{b b, \beta}(t), \rho_{\mathrm{S}}(t) \right\} \right] , \end{equation} \tag{ 54c }$$

$$\begin{eqnarray} &&\fl H_{\mathrm{LS}}(t) =\sum_{\alpha \beta} \left[\sum_{a {\neq} b} L^\dagger_{a b, \alpha}(t) L_{a b, \beta}(t) S_{\alpha \beta}(\omega_{ba}(t)) + \sum_{a b} L^\dagger_{a a, \alpha}(t) L_{b b, \beta}(t) S_{\alpha \beta}(0) \right] . \end{eqnarray} \tag{ 55 }$$

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 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 limit⁸. 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:

$$\begin{equation} \sum_{a,b} \to \sum_{\omega} , \quad L_{ab,\alpha} \to A_{\omega, \alpha} = \sum_{\varepsilon_b - \varepsilon_a = \omega} | \varepsilon_a \rangle \langle \varepsilon_a | A_{\alpha} | \varepsilon_b \rangle \langle \varepsilon_b| , \quad \omega_{ba} \to \omega. \end{equation} \noindent \tag{ 56 }$$

The resulting equation becomes:

$$\begin{equation} \fl \skew3\dot{\rho}_{\mathrm{S}}(t) = -\mathrm{i} \left[H_{\mathrm{S}} + H_{\mathrm{LS}}, \rho(t) \right] + \sum_{\alpha \beta} \sum_{\omega} \gamma_{\alpha \beta} (\omega) \left(A_{\omega, \beta} \rho_{\mathrm{S}}(t) A_{\omega, \alpha}^{\dagger} - \frac{1}{2} \left\{A_{\omega, \alpha}^{\dagger} A_{\omega, \beta}, \rho_{\mathrm{S}}(t) \right\} \right), \end{equation} \noindent \tag{ 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.

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/t_f, and we rederive in appendix E.1 the well known result [5] that to first order

$$\begin{equation} U_{\mathrm{S}}(t,t') = U^{\mathrm{ad}}_{\mathrm{S}}(t,t')\left[\mathbb{I}+ V_1(t,t')\right] {, } \end{equation} \tag{ 58 }$$

where

$$\begin{eqnarray} &&V_1(t,t') =- \sum_{a \neq b} |\varepsilon_a(t') \rangle \langle \varepsilon_b(t') |\int_{t'}^t \mathrm{d} \tau \mathrm{e}^{- \mathrm{i} \mu_{ba}(\tau,t')} \langle \varepsilon_a (\tau)|\dot{\varepsilon}_b(\tau) \rangle . \end{eqnarray} \tag{ 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^ad_S(t,t') replaced everywhere by the first order term U^ad_S(t,t')V₁(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

$$\begin{equation} \frac{h}{t_{\mathrm{f}}}\lesssim \Delta^2, \end{equation} \tag{ 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 $U_{\mathrm {S}}(t-\tau ,0)\approx \mathrm {e}^{\mathrm {i}\,\tau H_{\mathrm {S}}(t)}U^{\mathrm {ad}}_{\mathrm {S}}(t,t')\left [\mathbb {I}+ V_1(t,t')\right ]$. This still requires

$$\begin{equation} \frac{h}{t_{\mathrm{f}}}\ll \frac{1}{\tau_{\mathrm{B}}^2}, \end{equation} \tag{ 61 }$$

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 equation (35), which we can rewrite as ${\frac {h}{t_{\mathrm {f}}}} \ll \min (\Delta ^2,\frac {1}{\tau ^2_{\mathrm {B}}})$.

Recall that to ensure the validity of the Markov approximation and $\|\mathcal {L}_{\mathrm {uni}}\| \gg \|\mathcal {L}_{\mathrm {diss}}^{\mathrm {ad}}\|$, we also had to demand the inequalities in equations (31) and (33); these can be now summarized as $\max \left (g,\frac {g^2}{\Delta }\right ) \ll \frac {1}{\tau _{\mathrm {B}}}$, to be compared with equation (61).

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

$$\begin{eqnarray} &&{H_{\mathrm{S}}(t)} = A(t) H_{\mathrm{S}}^X + B(t) H_{\mathrm{S}}^Z {, } \end{eqnarray} \tag{ 62a }$$

$$\begin{eqnarray} &&H_{\mathrm{S}}^X \equiv \sum_{i=1}^N \sigma_i^x { ,} \end{eqnarray} \tag{ 62b }$$

$$\begin{eqnarray} &&H_{\mathrm{S}}^Z \equiv - \sum_{i=1}^N h_i \sigma_i^z + \sum_{i,j=1}^N J_{i j} \sigma_i^z \sigma_j^z, \end{eqnarray} \tag{ 62c }$$

Zoom In Zoom Out Reset image size

We couple this spin-system to a bath of harmonic oscillators, with bath and interaction Hamiltonian

$$\begin{equation} H_{\mathrm{B}} = \sum_{k=1}^\infty \omega_k b_k^\dagger b_k , \quad H_{\mathrm{I}} = \sum_{i=1}^N \sigma_i^z \otimes B_i , \quad B_i = \sum_k g_{k}^i\left(b_k^\dagger + b_k \right) , \end{equation} \tag{ 63 }$$

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 [36], except that our system is time-dependent. The resulting form for our operator L (equation (51)) is

$$\begin{equation} L_{a b, i} = | \varepsilon_a(t) \rangle \langle \varepsilon_a(t) | \sigma^z_{i} | \varepsilon_b(t) \rangle \langle \varepsilon_b(t) | = A_{iab}(t)|{\varepsilon_a(t)}\rangle \langle \varepsilon_b(t)|. \end{equation} \noindent \tag{ 64 }$$

Recall that our analysis assumed that the bath is in thermal equilibrium at inverse temperature β = 1/(k_BT), and hence is described by a thermal Gibbs state $\rho _{\mathrm {B}} = \exp \left ( - \beta H_{\mathrm {B}} \right ) / \mathcal {Z}$. We show in appendix H that this yields

$$\begin{eqnarray} &&\fl \gamma_{ij} (\omega) = \frac{2 \pi J(|\omega|)}{1 - \mathrm{e}^{-\beta| \omega|}} g_{|\omega|}^{i} g_{|\omega|}^{j} ( \Theta(\omega)+ \mathrm{e}^{-\beta |\omega|} \Theta(-\omega) ) \end{eqnarray} \tag{ 65a }$$

$$\begin{eqnarray} &&\fl S_{ i j}(\omega_{ba}(t)) = \int_0^\infty \mathrm{d} \omega \frac{J(\omega)}{1 - \mathrm{e}^{-\beta \omega}} g_\omega^{i } g_\omega^j \left( \mathcal{P} \left( \frac{1 }{\omega_{ba}(t) - \omega} \right) + \mathrm{e}^{- \beta \omega }\mathcal{P} \left( \frac{1}{\omega_{ba}(t) + \omega} \right) \right) , \end{eqnarray} \noindent \tag{ 65b }$$

$$\begin{equation} J(\omega) = {\eta} \omega \,\mathrm{e}^{-\omega/\omega_{\mathrm{c}}} , \end{equation} \tag{ 66 }$$

where ω_c is a high-frequency cutoff and η is a positive constant with dimensions of time squared.

It is often stated that the terms associated with the Lamb shift H_LS (equation (55)), i.e. equation (65b), can be neglected, since the relative order of S and H_S is g²τ_B/Δ, and indeed we have assumed g²τ_B/Δ ≪ 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 = ∞ for the Ohmic model (66), where for $\omega \gg \max _{ba}\omega _{ba}(t)$, the integrand tends to a constant. Secondly, S should be compared to γ, as both are of the same order g²τ_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 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_ba. 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 γ 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 γ matrix that satisfies the KMS condition will suffice without too much concern about how it is derived.

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

$$\begin{equation} \gamma_{ab}(\omega) = \frac{2\pi \omega \,\mathrm{e}^{-|\omega|/\omega_{\mathrm{c}}}}{1-\mathrm{e}^{-\beta\omega}}g^{a} g^{b}, \end{equation} \tag{ 67 }$$

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

$$\begin{equation} \mathcal{B}_{ab}(\tau )= \frac{\eta}{\beta^2} g^{a} g^{b} \left( \psi^{(1)} \left(\frac{1}{\beta\omega_{\mathrm{c}}} + \frac{\mathrm{i}\, \tau}{\beta} \right) + \psi^{(1)} \left(1 + \frac{1}{\beta\omega_{\mathrm{c}}} - \frac{\mathrm{i}\, \tau}{\beta} \right) \right) , \end{equation} \tag{ 68 }$$

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

$$\begin{equation} \beta\omega_{\mathrm{c}} \gg 1, \end{equation} \tag{ 69 }$$

and consider an expansion in large βω_c:

$$\begin{equation} \fl \mathcal{B}_{ab}(\tau) = \frac{\eta}{\beta^2} g^{a} g^{b} \left( - \pi^2 \mathrm{csch}^2\left(\frac{\pi \tau}{\beta} \right) + \sum_{n=1}^{\infty} \frac{\psi^{(n)} \left( 1 - \frac{\mathrm{i}\, \tau}{\beta} \right) + \psi^{(n)} \left( \frac{\mathrm{i}\, \tau}{\beta} \right)}{(n-1)! \left(\beta \omega_{\mathrm{c}}\right)^n}\right) , \end{equation} \tag{ 70 }$$

followed by an expansion in large τ/β:

$$\begin{equation} \mathcal{B}_{ab}(\tau) = \frac{\eta}{\beta^2} g^{a} g^{b} \left( -4 \pi^2 \mathrm{e}^{-\tau/\tau_{\mathrm{B}}} + \frac{1}{(\tau/\tau_M)^2} + O \left( \mathrm{e}^{-2\tau/\tau_{\mathrm{B}}} , \tau^{-3} \right) \right) . \end{equation} \tag{ 71 }$$

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:

$$\begin{equation} \tau_{\mathrm{B}} \overset{\omega_{\mathrm{c}} \rightarrow\infty}{\longrightarrow } \tau'_{\mathrm{B}} = \frac{\beta}{2 \pi} , \end{equation} \tag{ 72 }$$

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

$$\begin{equation} \tau_{\mathrm{M}} = \sqrt{ \frac{2 \beta}{\omega_{\mathrm{c}}}} . \end{equation} \tag{ 73 }$$

For sufficiently large ω_c, these two time scales capture the two types of behavior found in $\mathcal {B}_{ab}(\tau )$, as illustrated in figure 2(a).

Zoom In Zoom Out Reset image size

The transition between the exponential decay and the power-law tail occurs at a time τ_tr given by ${4 \pi ^2} \mathrm {e}^{-\tau _{\mathrm {tr}} / \tau '_{\mathrm {B}}} = (\frac {\tau _{\mathrm {M}}}{\tau _{\mathrm {tr}}})^{2}$, or equivalently by

$$\begin{equation} \frac{\exp(\theta)}{\theta^2} = \frac{1}{2}{\beta \omega_{\mathrm{c}}}, \qquad \theta\equiv 2\pi\frac{\tau_{\mathrm{tr}}}{\beta}. \end{equation} \tag{ 74 }$$

This transcendental equation has a formal solution in terms of the Lambert-W function [38], i.e. the inverse function of f(W) = We^W , as can be seen by changing variables to y = −θ/n, and rewriting θⁿe^−θ = a ≡ 2/(βω_c) as ye^y = −a^1/n/n, whose solution is $\theta =-ny = -nW(-\frac {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 θ² compared to e^θ, and approximate the solution to equation (74) by θ ∼ ln(βω_c/2). Thus

$$\begin{equation} \tau_{\mathrm{tr}} \sim \beta\ln(\beta\omega_{\mathrm{c}}) . \end{equation} \tag{ 75 }$$

This agrees with the first term of the asymptotic expansion $W(x) = \ln (x)-\ln \ln (x) +\frac { \ln \ln (x)}{\ln (x)}+\cdots$, which is accurate for x ≳ 3 [38].

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:

$$\begin{equation} \tau_{\mathrm{B}} = \tau'_{\mathrm{B}} + \sum_{n=1}^{\infty} \frac{c_n}{\omega_{\mathrm{c}}^n \beta^{n-1}}, \end{equation} \tag{ 76 }$$

where c_n 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

$$\begin{equation} \frac{1}{\omega_{\mathrm{c}}\ln(\beta\omega_{\mathrm{c}})} < \min\left\{2\tau_{\mathrm{B}},\frac{\tau_{\mathrm{B}} h}{\Delta^2 t_{\mathrm{f}}} +\frac{\tau_{\mathrm{B}}^3\,h}{t_{\mathrm{f}}}\right\}, \end{equation} \tag{ 77 }$$

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 βω_c ≫ 1.

For concreteness, we take gⁱ_ω = g, ω_c = 8π GHz and T = 20 mK ≈ 2.6 GHz (in units such that ℏ = 1; this is the operating temperature of the D-Wave Rainier chip [2]), corresponding to τ_B = 0.06 ns for the Ohmic model with infinite cutoff, equation (72), and τ_B ≈ 0.07 ns for ω_c = 8π 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 τ_tr = 0.33 ns. For these values, we satisfy at least one of the cases from equation (77), including numerical prefactors: $\frac {1}{\omega _{\mathrm {c}} \ln \left ( \beta \omega _{\mathrm {c}} \right )} < 2 \tau _{\mathrm {B}}$.

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

$$\begin{equation} h_0 = \frac{1}{4} , \quad h_{i>0}=0 , \quad J_{i, i+1} = -1 , \quad i = 0, \ldots, 7 , \end{equation} \tag{ 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:

$$\begin{equation} \rho_{\mathrm{S}}(t = 0) = \frac{\mathrm{e}^{- \beta H_{\mathrm{S}}(0)}}{\mathcal{Z}} . \end{equation} \tag{ 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 H_S(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 t_f (therefore more adiabatic) evolution shows a smaller difference between the two master 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.

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

5.4.1. 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 $|{\pm }\rangle _j = (|{\downarrow }\rangle _j\pm |{\uparrow }\rangle _j)/\sqrt {2}$ with energy ε₀(0) = −NA(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〉 ≡ ⊗ⁱ⁻¹_j=1| − 〉_j| + 〉_i⊗^N_j=i+1| − 〉_j, with energy ε₁(0) = −(N − 2)A(0). Therefore the gap between the ground state and the first excited states is:

$$\begin{equation} \Delta(t) = \varepsilon_1(t) - \varepsilon_0(t) \overset{t \ll t_{\mathrm{f}}}{\approx} [-(N-2) A(t)]-[-N A(t)] = 2 A(t) , \end{equation} \tag{ 80 }$$

which is at least twice as large as our k_BT ≈ 2.6 GHz, almost until A(t_c) = B(t_c) ≈ 5 GHz at t_c ≈ 0.35t_f (see figure 1). Noting that σ^z_i| ± 〉_i = | ± 〉_i, we can write the following relations in terms of the ground and first excited states:

$$\begin{equation} \omega_{i 0} = \Delta = -\omega_{0i}, \quad \sigma_i^z |i \rangle = | 0 \rangle . \end{equation} \tag{ 81 }$$

Noting that σ^z_i|j〉 is a doubly excited state $\forall j\neq i\in \{1,\ldots ,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, equation (18), we have

$$\begin{equation} \gamma(\omega_{0i}) = \mathrm{e}^{- \beta \omega_{i0}} \gamma(\omega_{i0}) , \end{equation} \tag{ 82 }$$

which shows that upward transitions are exponentially suppressed relative to downward transitions, by a factor ranging between e^{−2A(0)/(k_BT)} ≈ e^−67.4/2.6 ≈ 7 × 10⁻¹² and e^{−2A(0.3t_f)/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 ρ₀₀ = 〈0|ρ|0〉 and ρ_ii = 〈i|ρ|i〉, we can write the effective (truncated to the ground and first excited states) Lindblad equation (54) as the simplified rate equations

$$\begin{equation} \skew3\dot{\rho}_{ii} \approx \gamma_{i i} (\omega_{i 0}) \left(- \rho_{ii} + \mathrm{e}^{- \beta \omega_{i 0}} \rho_{00} \right) , \end{equation} \tag{ 83a }$$

$$\begin{equation} \rho_{00} \approx 1 - N \rho_{ii} , \end{equation} \tag{ 83b }$$

where we have assumed that the system is initially in the thermal state (79) and the gap is large (relative to k_BT). 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.

Zoom In Zoom Out Reset image size

5.4.2. 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 figure 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 $\dot {\rho }_{ii}$ equation with coefficient ρ_jj. Therefore, we find that the rate of scattering from the ground state to excited states is given by

$$\begin{equation} \mathrm{Excitation \ rate} \propto \gamma(\Delta(t)) \mathrm{e}^{-\beta \Delta(t)} \rho_{00} , \end{equation} \noindent \tag{ 84 }$$

whereas the relaxation rate is given by

$$\begin{equation} \mathrm{Relaxation \ rate} \propto \gamma(\Delta(t)) \rho_{ii} . \end{equation} \tag{ 85 }$$

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

Zoom In Zoom Out Reset image size

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 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 [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 ρ_ii coupled to ρ₀₀ (i.e. N is an exponentially large fraction of the dimension of the system Hilbert space), then ρ₀₀ decreases exponentially. By equation (84) this means that all ρ_ii 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 ρ_ii, 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 |〈0|σ^z_β|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 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 figure 7.

Zoom In Zoom Out Reset image size

5.4.4. Phase 4—the frozen phase

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)) between the instantaneous system density matrix and the instantaneous Gibbs state $\rho _{\mathrm {Gibbs}}(t) = \exp [-\beta H_{\mathrm {S}}(t)]/\mathcal {Z}$, where $\mathcal {Z}=\mathrm {Tr}\left [\exp \left ( - \beta H_{\mathrm {S}}(t) \right )\right ]$ 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/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 equation (63), to become important, and to disrupt the frozen phase, allowing the system to fully equilibrate into the Gibbs state.

Zoom In Zoom Out Reset image size

To consider adiabatic corrections, we recall that U_S satisfies

$$\begin{equation} \mathrm{i} \partial_t U_{\mathrm{S}}(t,t') = H_{\mathrm{S}}(t) U_{\mathrm{S}}(t,t') , \end{equation} \tag{ E.1 }$$

where

$$\begin{equation} {H}_{\mathrm{S}}(t) = \sum_a \varepsilon_a(t) |{\varepsilon_a(t)}\rangle\langle \varepsilon_a(t)| . \end{equation} \tag{ E.2 }$$

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

$$\begin{equation} W(t,t') \equiv \sum_a |{\varepsilon_a(t)}\rangle\langle \varepsilon_a(t')| = T_+ \exp\left[-\mathrm{i}\int_{t'}^t \mathrm{d}\tau \ K(\tau)\right], \end{equation} \tag{ E.3 }$$

where the 'intertwiner Hamiltonian' is

$$\begin{equation} K(t) \equiv \mathrm{i} [\partial_t W(t,t')]W^\dagger(t,t') = \mathrm{i}\sum_{a} |{\dot{\varepsilon}_a(t)}\rangle\langle \varepsilon_a(t)| \end{equation} \tag{ E.4a }$$

$$\begin{equation} \hphantom{K(t)}= \mathrm{i}[\skew3\dot{P}_0(t),P_0(t)], \end{equation} \tag{ E.4b }$$

and where P₀(t) ≡ |ε₀(t)〉〈ε₀(t)| is the projection onto the ground state of H_S(t). The result (E.4b) is by no means obvious and is proven in section E.3.

To extract the geometric phase we define

$$\begin{equation} H_{\mathrm{G}}(t) \equiv \sum_a \phi_a(t) |{\varepsilon_a(t)}\rangle \langle \varepsilon_a(t)| , \end{equation} \tag{ E.5a }$$

$$\begin{equation} \phi_a(t) = \mathrm{i}\langle {\varepsilon_a(t)}|{\dot{\varepsilon}_a(t)}\rangle, \end{equation} \tag{ E.5b }$$

$$\begin{equation} H'_{\mathrm{S}}(t) \equiv H_{\mathrm{S}}(t)-H_{\mathrm{G}}(t) . \end{equation} \tag{ E.5c }$$

Now define V via the 'adiabatic interaction picture' transformation:

$$\begin{equation} V(t,t') = W^\dagger(t,t')U_{\mathrm{S}}(t,t') , \end{equation} \tag{ E.6 }$$

along with

$$\begin{equation} \tilde{H}_{\mathrm{S}}(t,t') = W^\dagger(t,t') H_{\mathrm{S}}(t) W(t,t') = \sum_a \varepsilon_a(t) |{\varepsilon_a(t')}\rangle \langle \varepsilon_a(t')| , \end{equation} \tag{ E.7a }$$

$$\begin{equation} \tilde{H}_{\mathrm{G}}(t,t') = W^\dagger(t,t') H_{\mathrm{G}}(t) W(t,t') = \sum_a \phi_a(t) |{\varepsilon_a(t')}\rangle \langle \varepsilon_a(t')| , \end{equation} \tag{ E.7b }$$

$$\begin{equation} \tilde{H}'_{\mathrm{S}}(t,t') \equiv \tilde{H}_{\mathrm{S}}(t,t')-\tilde{H}_{\mathrm{G}}(t,t') , \end{equation} \tag{ E.7c }$$

$$\begin{equation} \tilde{K}(t,t') = W^\dagger(t,t') K(t) W(t,t') = {\rm i}W^\dagger(t,t') \partial_t W(t,t'). \end{equation} \tag{ E.7d }$$

Note that the time dependence of $\tilde {H}_{\mathrm {S}}$ and $\tilde {H}'_{\mathrm {S}}(t)$ is entirely in the energy eigenvalue and not in the eigenstates. Then V obeys the Schrödinger equation:

$$\begin{equation} \mathrm{i}\, \partial_t V(t,t') = \tilde{H}_{\mathrm{S}}^{\mathrm{ad}}(t,t') V(t,t') . \end{equation} \tag{ E.8 }$$

where

$$\begin{equation} \tilde{H}_{\mathrm{S}}^{\mathrm{ad}}(t,t') \equiv \tilde{H}_{\mathrm{S}}(t,t') - \tilde{K}(t,t') = W^\dagger(t,t') [H_{\mathrm{S}}(t)-K(t)] W(t,t'). \end{equation} \tag{ E.9 }$$

When the evolution is nearly adiabatic $\tilde {H}_{\mathrm {S}}^{\mathrm {ad}}(t,t')$ is a perturbation, so that we consider a solution of equation (E.8) for V of the form:

$$\begin{equation} V(t,t') = V_0(t,t') \left( {\mathbbm 1} + V_1(t,t') + \cdots \right) , \end{equation} \tag{ E.10 }$$

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

$$\begin{eqnarray} V_0(t,t') &\equiv& T_+ \exp \left[-\mathrm{i}\, \int_{t'}^t \mathrm{d} \tau \tilde{H}'_{\mathrm{S}}(\tau,t')\right] \end{eqnarray} \tag{ E.11a }$$

$$\begin{eqnarray} U^{\mathrm{ad}}_{\mathrm{S}}(t,t') &\equiv& W(t,t') V_0(t,t') = \sum_a | \varepsilon_a(t) \rangle \langle \varepsilon_a (t')| \mathrm{e}^{-\mathrm{i} \int_{t'}^t \mathrm{d} \tau \left[\varepsilon_a(\tau)-\phi_a(\tau)\right] } . \end{eqnarray} \tag{ E.11b }$$

Differentiating with respect to t, we have

$$\begin{equation} \dot{U}^{\mathrm{ad}}_S(t,t') = \dot{W}(t,t')V_0(t,t')+W(t,t')\dot{V}_0(t,t') \end{equation} \tag{ E.12a }$$

$$\begin{equation} \hphantom{\dot{U}^{\mathrm{ad}}_S(t,t')}= -\mathrm{i}\,K(t)W(t,t')V_0(t,t')-\mathrm{i}\,W(t,t')\tilde{H}'_{\mathrm{S}}(t,t')V_0(t,t') \end{equation} \tag{ E.12b }$$

$$\begin{equation} \hphantom{\dot{U}^{\mathrm{ad}}_S(t,t')}= -\mathrm{i}\left[H_{\mathrm{S}}^{\mathrm{ad}}(t)-H_{\mathrm{G}}(t)\right]U^{\mathrm{ad}}_{\mathrm{S}}(t,t') , \end{equation} \tag{ E.12c }$$

where

$$\begin{eqnarray} H_{\mathrm{S}}^{\mathrm{ad}}(t) &\equiv& K(t)+H_{\mathrm{S}}(t) . \end{eqnarray} \tag{ E.13 }$$

Plugging the equation (E.10) expansion into equation (E.8), we obtain, to first order in 1/t_f:¹⁰

$$\begin{equation} \fl \mathrm{i} \partial_t V_1 (t,t') = - V^\dagger_0(t,t') [\tilde{K}(t,t')-\tilde{H}_{\mathrm{G}}(t,t')] V_0(t,t') \end{equation} \tag{ E.14a }$$

$$\begin{equation} \fl \hphantom{\mathrm{i} \partial_t V_1 (t,t')}= -\mathrm{i} U^{\mathrm{ad}\dagger}_{\mathrm{S}}(t,t') \partial_t W(t,t') V_0(t,t') +U^{\mathrm{ad}\dagger}_{\mathrm{S}}(t,t')H_{\mathrm{G}}(t)U^{\mathrm{ad}}_{\mathrm{S}}(t,t') . \end{equation} \tag{ E.14b }$$

Note that $U^{\mathrm {ad}\dagger }_{\mathrm {S}}(t,t')H_{\mathrm {G}}(t)U^{\mathrm {ad}}_{\mathrm {S}}(t,t') = \sum _a\phi _a(t) |{\varepsilon _a(t')}\rangle \langle \varepsilon _a(t')|$. Therefore, integrating, and using

$$\begin{equation} \fl \mu_{ba}(t,t') = \int_{t'}^t \mathrm{d}\tau \, \omega_{ba}(\tau), \quad \omega_{ba}(\tau) = [\varepsilon_b(\tau)-\phi_b(t)] - [\varepsilon_a(\tau) -\phi_a(t)], \end{equation} \tag{ E.15 }$$

we can write the solution as:

$$\begin{eqnarray} &&\fl V_1(t,t') = - \int_{t'}^t \mathrm{d} \tau\ \left[U^{\mathrm{ad}\dagger}_S(\tau,t') \partial_\tau W(\tau , t') V_0(\tau, t') -\sum_a\phi_a(\tau) |{\varepsilon_a(t')}\rangle \langle \varepsilon_a(t')| \right] \nonumber \\ &&\fl \hphantom{V_1(t,t')} =- \sum_{a \neq b} \int_{t'}^t \mathrm{d} \tau \mathrm{e}^{- \mathrm{i} \mu_{ba}(\tau,t')} |\varepsilon_a(t') \rangle \langle \varepsilon_b(t') | \langle \varepsilon_a (\tau)|\dot{\varepsilon}_b(\tau) \rangle . \end{eqnarray} \tag{ E.16 }$$

Therefore, the system evolution operator can be written as:

$$\begin{equation} U_{\mathrm{S}}(t,t') = U^{\mathrm{ad}}_{\mathrm{S}}(t,t') + Q(t,t') U^{\mathrm{ad}}_{\mathrm{S}}(t,t') + O\left(t_{\mathrm{f}}^{-2}\right) , \end{equation} \tag{ E.17 }$$

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

$$\begin{eqnarray} &&\fl Q(t,t') \equiv U^{\mathrm{ad}}_{\mathrm{S}}(t,t') V_1(t,t') U^{\mathrm{ad}\dagger}_{\mathrm{S}}(t,t')\nonumber\\ &&\fl \hphantom{Q(t,t')} = \sum_{a \neq b} \mathrm{e}^{-\mathrm{i} \mu_{ba}(t,t')} \left( - \int_{t'}^t \mathrm{d} \tau \mathrm{e}^{\mathrm{i} \mu_{ba}(\tau,t')} \langle \varepsilon_a (\tau)| \dot{\varepsilon}_b(\tau) \rangle \right) | \varepsilon_a (t) \rangle \langle \varepsilon_b(t)| . \end{eqnarray} \tag{ E.18 }$$

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 $\frac {h}{\Delta ^2 t_{\mathrm {f}}} \ll 1$ (equation (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 equation at t_f, e.g. [6, 10–12].

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

$$\begin{equation} Q_{ab}(t,t') = \mathrm{e}^{-\mathrm{i} \mu_{ba}(t,t')} \left( - \int_{t'}^t \mathrm{d} \tau \mathrm{e}^{\mathrm{i} \mu_{ba}(\tau,t')} \langle \varepsilon_a (\tau)| \partial_\tau | \varepsilon_b(\tau) \rangle \right) . \end{equation} \noindent \tag{ E.19 }$$

Changing variables to dimensionless time s = t/t_f, μ_ba(t,t') becomes $t_{\mathrm {f}} \int _{t'/t_{\mathrm {f}}}^{s} \mathrm {d}s' \, \omega _{ba}(s') = t_{\mathrm {f}} \tilde {\mu }_{ba}(s,t'/t_{\mathrm {f}})$, and

$$\begin{equation} \fl \int_{t'}^t \mathrm{d} \tau \mathrm{e}^{\mathrm{i} \mu_{ba}(\tau,t')} \langle \varepsilon_a (\tau)| \partial_\tau | \varepsilon_b(\tau) \rangle = \int_{t'/t_{\mathrm{f}}}^{s} \mathrm{d}s' \ \mathrm{e}^{\mathrm{i} t_{\mathrm{f}} \tilde{\mu}_{ba}(s'\!,t'/t_{\mathrm{f}})} \langle \varepsilon_a (s')| \partial_{s'} | \varepsilon_b(s') \rangle , \end{equation} \noindent \tag{ E.20 }$$

where $\tilde {\mu }$ now involves a dimensionless time integration. Using the fact that $\mathrm{e}^{ \mathrm{i} t_\mathrm{f} \tilde{\mu}_{ba} (s',t'/t_\mathrm{f}) } = \frac{\rm i}{t_\mathrm{f} \omega_{ab}(s')} \frac{{\rm d}}{{\rm d} s'} \mathrm{e}^{ \mathrm{i} t_\mathrm{f} \tilde{\mu}_{ba} (s',t'/t_\mathrm{f})}$, we can integrate equation (E.20) by parts as

$$\begin{eqnarray} &&\fl \int_{t'/t_\mathrm{f}}^{s} \mathrm{d}s' \ \mathrm{e}^{\mathrm{i} t_\mathrm{f} \tilde{\mu}_{ba}(s',t'/T)} \langle \varepsilon_a (s')| \partial_{s'} | \varepsilon_b(s') \rangle = \frac{\mathrm{i}}{t_\mathrm{f} \omega_{ab}(s')} \mathrm{e}^{\mathrm{i} t_\mathrm{f} \tilde{\mu}_{a b}(s',t'/t_\mathrm{f})} \langle \varepsilon_a(s')|\partial_{s'} | \varepsilon_b(s') \rangle \Big|_{t'/t_\mathrm{f}}^{s} \nonumber \\ &&- \frac{\mathrm{i}}{t_\mathrm{f}} \int_{t'/t_\mathrm{f}}^{s} \mathrm{d}s' \ \frac{\mathrm{e}^{\mathrm{i} t_\mathrm{f} \tilde{\mu}_{a b}(s',t'/t_\mathrm{f})}}{\omega_{ab}(s')}\frac{\rm d}{\mathrm{d} s'} {\langle \varepsilon_a(s')|\partial_{s'} | \varepsilon_b(s') \rangle} . \end{eqnarray} \noindent \tag{ E.21 }$$

Now note that for non-degenerate energy eigenstates, differentiation with respect to t of H_S(t)|ε_a(t)〉 = ε_a(t)|ε_a(t)〉, and substitution of s = t/t_f in equation (26), directly yields the relation:¹¹

$$\begin{equation} \langle {\varepsilon_b(t)}|{\dot{\varepsilon}_a(t)}\rangle = \frac{ \langle \varepsilon_b (t) | \dot{H}_{\mathrm{S}}(t) | \varepsilon_a (t) \rangle }{\omega_{ab} (t)} \sim \frac{h}{\Delta t_{\mathrm{f}}} . \end{equation} \noindent \tag{ E.22 }$$

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

$$\begin{eqnarray*} &&\fl |Q_{ab}(t,t')| \leqslant \frac{|\langle \varepsilon_a (s) | \partial_s H_\mathrm{S}(s) | \varepsilon_b (s) \rangle |}{\omega_{ab}^2(s) t_\mathrm{f}} \Big|_{t'/t_\mathrm{f}}^{t/t_\mathrm{f}} \nonumber + \left| \int_{t'/t_\mathrm{f}}^{t/t_\mathrm{f}} \mathrm{d}s' \ \frac{\mathrm{e}^{\mathrm{i} t_\mathrm{f} \tilde{\mu}_{a b}(s',t'/t_\mathrm{f})}}{\omega_{ab}^2(s') t_\mathrm{f}}\frac{\rm d}{\mathrm{d} s'} {\langle \varepsilon_a (s') | \partial_{s'} H_\mathrm{S}(s') | \varepsilon_b (s') \rangle} \right| . \end{eqnarray*} \noindent$$

Continued integration by parts will yield higher powers of the dimensionless quantity $\frac {\langle \varepsilon _a (s) | \partial _\tau H_{\mathrm {S}}(s) | \varepsilon _b (s) \rangle }{\omega _{ab}^2(s) t_{\mathrm {f}}}$. Thus a sufficient condition for the smallness of |Q_ab(t,t')| for all a,b is that

$$\begin{equation} \frac{\max_{s;a,b} |\langle \varepsilon_a (s) | \partial_s H_{\mathrm{S}}(s) | \varepsilon_b (s) \rangle |}{\min_{s;a,b}\omega_{ab}^2(s) t_{\mathrm{f}}} \ll 1, \end{equation} \noindent \tag{ E.23 }$$

namely the condition in equation (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. [49–54]. For rigorous analyses see, e.g., [5–7, 10–12]. For our purposes equation (27) suffices.

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 P₀(t) ≡ |ε₀(t)〉〈ε₀(t)| and $Q_0(t) \equiv \mathbb {I}-P_0(t)$. Recalling equation (E.17), we have

$$\begin{equation} \fl P_0(t) U_{\mathrm{S}}(t,t') P_0(t') = P_0(t) U^{\mathrm{ad}}_{\mathrm{S}}(t,t')P_0(t') + P_0(t) Q(t,t') U^{\mathrm{ad}}_{\mathrm{S}}(t,t')P_0(t') + O\left(t_{\mathrm{f}}^{-2}\right) . \end{equation} \noindent \tag{ E.24 }$$

Using equation (E.11b) we find $U^{\mathrm {ad}}_{\mathrm {S}}(t,t')P_0(t') = {P_0(t)} \mathrm {e}^{-\mathrm {i} \int _{t'}^t \mathrm {d} \tau \varepsilon _0(\tau ) }$, so that up to a phase P₀(t)Q(t,t')U^ad_S(t,t')P₀(t') = P₀(t)Q(t,t')P₀(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 P₀(t)Q(t,t')P₀(t') = 0. Thus, to first order in 1/t_f the evolution generated by H_S(t) keeps the ground state decoupled from the excited states, i.e.

$$\begin{equation} P_0(t) U_{\mathrm{S}}(t,t') P_0(t') = U^{\mathrm{ad}}_{\mathrm{S}}(t,t') P_0(t') + O\left(t_{\mathrm{f}}^{-2}\right) {,} \end{equation} \tag{ E.25a }$$

$$\begin{equation} Q_0(t) U_{\mathrm{S}}(t,t') Q_0(t') = U^{\mathrm{ad}}_{\mathrm{S}}(t,t') Q_0(t') + O\left(t_{\mathrm{f}}^{-2}\right) {.} \end{equation} \noindent \tag{ E.25b }$$

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

$$\begin{equation} \fl P_0(t) [\mathbb{I} + \mathrm{i}\,\tau H_{\mathrm{S}}(t)][P_0(t) -\tau \skew3\dot{P}_0(t)] = [\mathbb{I} + \mathrm{i}\,\tau H_{\mathrm{S}}^{\mathrm{ad}}(t)][P_0(t) -\tau \skew3\dot{P}_0(t)] + O\left(t_{\mathrm{f}}^{-2}\right) {,} \end{equation} \tag{ E.26a }$$

$$\begin{equation} \fl Q_0(t) [\mathbb{I} + \mathrm{i}\,\tau H_{\mathrm{S}}(t)][Q_0(t) -\tau \skew3\dot{Q}_0(t)] = [\mathbb{I} + \mathrm{i}\,\tau H_{\mathrm{S}}^{\mathrm{ad}}(t)][Q_0(t) -\tau \skew3\dot{Q}_0(t)] + O\left(t_{\mathrm{f}}^{-2}\right) . \end{equation} \noindent \tag{ E.26b }$$

Using the projector properties $P_0 = P_0^2 \Longrightarrow \skew3\dot{P}_0 P_0 + P_0\skew3\dot{P}_0 = \skew3\dot{P}_0$ and [H_S(t),P₀(t)] = 0 (similarly for Q₀), this simplifies to

$$\begin{equation} \fl -\tau P_0(t)\skew3\dot{P}_0(t) + \mathrm{i}\,\tau H_{\mathrm{S}}(t) P_0(t) = -\tau \skew3\dot{P}_0(t) + \mathrm{i}\,\tau H_{\mathrm{S}}^{\mathrm{ad}}(t) P_0(t)+O\left(t_{\mathrm{f}}^{-2}\right) , \end{equation} \tag{ E.27a }$$

$$\begin{equation} \fl -\tau Q_0(t)\skew3\dot{Q}_0(t) + \mathrm{i}\,\tau H_{\mathrm{S}}(t) Q_0(t) = -\tau \skew3\dot{Q}_0(t) + \mathrm{i}\,\tau H_{\mathrm{S}}^{\mathrm{ad}}(t) Q_0(t)+O\left(t_{\mathrm{f}}^{-2}\right) , \end{equation} \noindent \tag{ E.27b }$$

and, after dropping the $O (t_{\mathrm {f}}^{-2})$ corrections, to

$$\begin{equation} \mathrm{i}[H_{\mathrm{S}}^{\mathrm{ad}}(t)-H_{\mathrm{S}}(t)]P_0(t) + \skew3\dot{P}_0(t) P_0(t) = 0, \end{equation} \tag{ E.28a }$$

$$\begin{equation} \mathrm{i}[H_{\mathrm{S}}^{\mathrm{ad}}(t)-H_{\mathrm{S}}(t)]Q_0(t) + \skew3\dot{Q}_0(t) Q_0(t) = 0 . \end{equation} \noindent \tag{ E.28b }$$

Inserting $Q_0(t) = \mathbb {I} - P_0(t)$ and $\skew3\dot{Q}_0(t) = -\skew3\dot{P}_0(t)$ into equation (E.28b), and subtracting it from equation (E.28a), we find, using $\skew3\dot{P}_0 P_0 + P_0\skew3\dot{P}_0 = \skew3\dot{P}_0$ once more, the desired result:

$$\begin{equation} H_{\mathrm{S}}^{\mathrm{ad}}(t) = H_{\mathrm{S}}(t) + \mathrm{i}[\skew3\dot{P}_0(t), P_0(t)], \end{equation} \tag{ E.29 }$$

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

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

$$\begin{equation} H_{\mathrm{S}}(t) = A(t) \sigma^x - B(t) h_z \sigma^z . \end{equation} \noindent \tag{ J.1 }$$

The gap at any given time in the evolution is

$$\begin{equation} \Delta(t) = 2 \sqrt{ A(t)^2 + B(t)^2 h_z^2 } . \end{equation} \noindent \tag{ J.2 }$$

Since there are only two states, there are only three γ(ω) terms to calculate (assuming g_ω = g):

$$\begin{equation} \fl \gamma(0) = \frac{2 \pi g^2}{\beta} , \quad \gamma(+ \Delta(t) ) = 2 \pi g^2 \frac{\Delta(t)}{1 - \mathrm{e}^{- \beta \Delta(t)}} , \quad \gamma(- \Delta(t) ) = 2 \pi g^2 \frac{\Delta(t) \mathrm{e}^{-\beta \Delta( t)}}{1 - \mathrm{e}^{- \beta \Delta(t)}} , \end{equation} \noindent \tag{ J.3 }$$

where we have taken the limit ω_c → ∞ for simplicity. Let us consider the gapped phase of our evolution where t/t_f is small. In this phase, B(t) ≈ 0, so we can safely assume that the energy eigenstates remain diagonal in the σ^x basis, with ground state | − 〉 and excited state | + 〉, and do not change such that:

$$\begin{equation} \langle \pm | \skew3\dot{\rho} | \pm \rangle = \frac{\mathrm{d}}{\mathrm{d}t} \left( \langle \pm | \rho | \pm \rangle \right) \equiv \skew3\dot{\rho}_{\pm \pm} . \end{equation} \tag{ J.4 }$$

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

$$\begin{equation} \fl \skew3\dot{\rho}_{++}(t) = | \langle + | \sigma^z | - \rangle |^2 \left( \gamma(-\Delta(t)) \rho_{--} - \gamma(\Delta(t)) \rho_{++} \right) = \gamma_0(\Delta(t)) \left( \mathrm{e}^{-\beta \Delta(t)} \rho_{--} - \rho_{++} \right) , \end{equation} \tag{ J.5a }$$

$$\begin{equation} \fl \rho_{--}(t) = 1 - \rho_{++}(t), \end{equation} \tag{ J.5b }$$

where we have used that 〈 + |σ^z| + 〉 = 〈 − |σ^z| − 〉 = 0, 〈 + |σ^z| − 〉 = 1 and $\left [ H_{\mathrm {S}}, \rho _{\mathrm {S}} \right ] \approx 0$. Interpreting equation (J.5) 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⁻¹¹. Since the gap remains relatively large during the gapped phase, the change in the initial thermal state is minimal.

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) (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:

$$\begin{equation} \fl \skew3\dot{\rho}_{ii} \approx \sum_{\alpha \beta} \gamma_{\alpha \beta} (\omega_{i 0})\langle 0 | \sigma^z_{\beta} | i \rangle \langle i | \sigma^z_{\alpha} | 0 \rangle \left(- \rho_{ii} + \mathrm{e}^{- \beta \omega_{i 0}} \rho_{00} \right) = \gamma_{i i} (\omega_{i 0}) \left( \mathrm{e}^{- \beta \omega_{i 0}} \rho_{00} - \rho_{ii} \right) \end{equation} \tag{ J.6a }$$

$$\begin{equation} \fl \rho_{00} \approx 1 - N \rho_{ii} , \end{equation} \tag{ J.6b }$$

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