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.
Export citation and abstract BibTeX RIS
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.
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. Thunström 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 Schrödinger 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 Schrödinger 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 Schrödinger 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 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 , where we have assumed units such that ℏ = 1.
Let and denote the free system and bath unitary evolution operators (T+ denotes time ordering), and define U0(t,t') = US(t,t')⊗UB(t,t'). Likewise, let denote the joint Schrödinger picture system–bath unitary evolution operator. We transform to the interaction picture with respect to HS(t) and HB, by defining , which, together with the interaction picture density matrix , satisfies
We restrict the use of tilde variables to refer to variables in the interaction picture. The time-dependent interaction picture Hamiltonian is related to its Schrödinger picture counterpart via:
where we have defined the time-dependent system and bath operators:
2.2. Born–Markov approximation
Writing the formal solution of equation (3b) as
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
We make the standard Born approximation assumption, that we can decompose the density matrix as where χ(t), which expresses correlations between the system and the bath, is small in an appropriate sense and can hence be neglected from now on [30, 34].7 Thus the equation of motion reduces to:
where we defined the two-point correlation functions:
In equation (9), we have assumed for simplicity (but without loss of generality) that , 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:
For notational simplicity we will denote by 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. . More precisely, we shall require throughout that
As we show in appendix B, if τB ≪ 1/g, then we can apply the Markov approximation to each of the four summands in equation (9), i.e., replace by :
where 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)2]. An explicit derivation of the upper bound on the error due to this approximation can be found in appendix B. The resulting Markovian equation cannot resolve the dynamics over a time scale shorter than τB.
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 and . Our goal is to express these integrals in terms of the spectral-density matrix
the standard quantity in master equations. It is convenient to replace the one-sided Fourier transform in the spectral-density matrix by a complete Fourier transform. Thus we split it into a sum of Hermitian matrices,
where we show in appendix C that γ and S are given by
If in addition the correlation function is analytic in the strip between τ = −i β and 0, then it follows that the Fourier transform of the bath correlation function satisfies the detailed balance condition
For a proof see appendix D.
It is important to note that the KMS detailed balance condition imposes a restriction on the decay of the correlation function. To see this, note first that , where we used equation (11). Now assume that equation (12) would have to hold for all n. It would follow that
Thus all derivatives of γab(ω) would have to be finite at ω = 0.
On the other hand, it follows from the KMS condition (18) that
so that in the limit as ω approaches zero from below or above,
This shows that in principle γ'aa(ω) may be discontinuous at ω = 0. Indeed, the continuity condition γaa'(0−) = γaa'(0+) implies, from the KMS condition recast as equation (21), that
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(ω) = ω2/(1 − e−βω) for ω ⩾ 0 and γab(ω) = ω2/(e−βω − 1) for ω ⩽ 0. Both examples violate equation (19) for n ⩾ 2. Note, moreover, that when this happens, it follows from equation (19) that is divergent, so that we must conclude that 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)−2, 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
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 Δ, i.e.
where ε0(t) and ε1(t) are the ground and first excited state energies of HS(t). We shall arrive at master equations of the following general form:
Here is the unitary evolution superoperator including the Lamb shift correction, is the dissipative superoperator in the fully adiabatic limit and is the dissipative superoperator with leading order non-adiabatic corrections.
Let
where s ≡ t/tf is the dimensionless time. To ensure that generates adiabatic evolution to leading order we shall require the standard adiabatic condition
In order to derive equation (25), the three superoperators are ordered in perturbation theory, in the sense that
where the norm could be chosen as the supoperator norm, i.e. the largest singular value (see appendix A). We may then assume that
Combining the O(τB) due to the integral on the rhs of equation (13) (as already remarked there) with the g2 prefactor from equation (9), we have
To ensure the first inequality in equation (28) we thus require
The non-adiabatic correction is of order , and when it appears in it is multiplied by the same factor as , i.e. we have
To ensure the second inequality in equation (28) thus amounts to the adiabatic condition, equation (27). All this is added to the condition
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 τB. If we require the change in the basis to be small on the time scale of the bath τB, we must have:
We justify this claim in section 4.1 and appendix F.
Note that the adiabatic condition, equation (27), implies , while equation (34) can be written as . Putting this together thus yields
Our other inequalities (equations (31) and (33)) can be summarized as
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) {εa(t)}, i.e. HS(t)|εa(t)〉 = εa(t)|εa(t)〉, and Bohr frequencies ωba(t) ≡ εb(t) − εa(t). As shown in appendix E.1 we can then write the system time evolution operator as:
where is the Berry connection. The correction term 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βUS(t − τ,0) with an appropriate adiabatic approximation, which will allow us to take this operator outside the integral. To see how, note first that
We now make two approximations: firstly, as per equation (37a) we replace US(t,0) by UadS(t,0); secondly, we replace U†S(t,t − τ) by ei τHS(t), an approximation justified by the appearance of the short-lived bath correlation function inside the integrals we are concerned with. Thus, we write
and find the bound on the error due to dropping Θ(t,τ) in appendix F. Let
Thus, to the same level of approximation
where
and Γαβ(ωba(t)) is the spectral-density matrix defined in equation (14). Similarly, we have for the other term
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:
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 Uad†S(t,0)AαUadS(t,0). If this is done, we obtain the double-sided adiabatic interaction picture master equation
It is convenient to transform back into the Schrödinger picture. Using (equation (6)) implies that . Hence, using equation (46), we find the one-sided adiabatic Schrödinger picture master equation
where
This form of the master equation has not appeared in previous studies of adiabatic master equations. If we again use the adiabatic approximation for US(t,0), i.e. replace US(t,0)Πab(0)U†S(t,0) by Uad†S(t,0)Πab(0)UadS(t,0), we obtain the double-sided adiabatic Schrödinger picture master equation
where
Comparing to equation (25), the first term in equations (48) and (50) is , while the second is .
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.
4.3. Master equation in the adiabatic limit with rotating wave approximation: Lindblad form
In order to arrive at a master equation in Lindblad form, we can perform a 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) ≈ Uad†S(t,0)AαUadS(t,0). This yields
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:
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 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:
The resulting equation becomes:
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.
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
where
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 UadS(t,t') replaced everywhere by the first order term UadS(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
i.e. with the ≪ replaced by a mere ≲. However, we would still like to perform the approximation of equation (40), in the sense that . This still requires
which is what allows the use of ei τHS(t) in this approximation (as shown in appendix F). Taken together, these two conditions are weaker than equation (35), which we can rewrite as .
Recall that to ensure the validity of the Markov approximation and , we also had to demand the inequalities in equations (31) and (33); these can be now summarized as , 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
We couple this spin-system to a bath of harmonic oscillators, with bath and interaction Hamiltonian
where b†k and bk are, respectively, raising and lowering operators for the kth oscillator with natural frequency ωk, 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
Recall that our analysis assumed that the bath is in thermal equilibrium at inverse temperature β = 1/(kBT), and hence is described by a thermal Gibbs state . We show in appendix H that this yields
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 HLS (equation (55)), i.e. equation (65b), can be neglected, since the relative order of S and HS is g2τB/Δ, and indeed we have assumed g2τ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 , the integrand tends to a constant. Secondly, S should be compared to γ, as both are of the same order g2τB, 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 + ωLSba. 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.
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
by inverse Fourier transform of equation (16a). The result is
where ψ(m) is the mth Polygamma function (see appendix I for the derivation). We first assume
and consider an expansion in large βωc:
followed by an expansion in large τ/β:
This expansion reveals the two independent time-scales that are relevant for us. First, there is the time scale τB associated with the exponential decay (corresponding to the true Markovian bath), given by:
then the time scale associated with non-Markovian corrections (the power-law tail):
For sufficiently large ωc, these two time scales capture the two types of behavior found in , as illustrated in figure 2(a).
Download figure:
Standard imageThe transition between the exponential decay and the power-law tail occurs at a time τtr given by , or equivalently by
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 = −θ/n, and rewriting θne−θ = a ≡ 2/(βωc) as yey = −a1/n/n, whose solution is . However, for our purposes the following observations will suffice. We seek a Markovian-like solution, where τtr is large compared to the thermal timescale set by β, i.e. we are interested in the regime where θ ≫ 1. In this case we can neglect θ2 compared to eθ, and approximate the solution to equation (74) by θ ∼ ln(βωc/2). Thus
This agrees with the first term of the asymptotic expansion , 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:
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
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.
5.3. Numerical results
For concreteness, we take giω = 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: .
We focus on the N = 8 site ferromagnetic chain with parameters:
where we pin the first spin in order to break the degeneracy in the ground state of the classical Ising Hamiltonian. The system is initialized in the Gibbs state:
To help the numerics, we truncate our system to the lowest n = 18 levels (out of 256), rotating the density matrix into the instantaneous energy eigenbasis at each time step. The error associated with this is small as long as our evolution does not cause scattering into higher n states, as we have checked. The forward propagation algorithm used is an implicit second order Runge–Kutte method called TR-BDF2 with an adaptive time step [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 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.
Download figure:
Standard imageIn 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].
Download figure:
Standard image5.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 HXS, i.e. the state |0〉 ≡ ⊗Nj=1| − 〉j, where with energy ε0(0) = −NA(0), and where |↓〉j,|↑〉j are the +1,−1 eigenstates of σzj (computational basis states of the jth spin or qubit). The lowest lying energy states are then the N-fold degenerate states with a single flip of one of the spins in the x-direction, i.e. |i〉 ≡ ⊗i−1j=1| − 〉j| + 〉i⊗Nj=i+1| − 〉j, with energy ε1(0) = −(N − 2)A(0). Therefore the gap between the ground state and the first excited states is:
which is at least twice as large as our kBT ≈ 2.6 GHz, almost until A(tc) = B(tc) ≈ 5 GHz at tc ≈ 0.35tf (see figure 1). Noting that σzi| ± 〉i = | ± 〉i, we can write the following relations in terms of the ground and first excited states:
Noting that σzi|j〉 is a doubly excited state , 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
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 × 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 ρ00 = 〈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
Download figure:
Standard image5.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 equation with coefficient ρjj. Therefore, we find that the rate of scattering from the ground state to excited states is given by
whereas the relaxation rate is given by
Therefore, as we emerge from the gapped phase, we have ρ00 ≫ ρii, so scattering into excited states dominates over relaxation into the ground state. This explains the loss of fidelity in phase 2.
Download figure:
Standard image5.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 ρii coupled to ρ00 (i.e. N is an exponentially large fraction of the dimension of the system Hilbert space), then ρ00 decreases exponentially. By 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 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.
Download figure:
Standard image5.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 σz basis, and so the off-diagonal components of the Lab,i operators vanish (or become very small), i.e. Aiab(t) = 〈εa(t)|σzi|εb(t)〉∝δab (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)) between the instantaneous system density matrix and the instantaneous Gibbs state , where 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 σx and σy 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.
Download figure:
Standard image6. 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 . 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|:
Therefore ∥A|ψ〉∥ ⩽ ∥A∥∞ for all normalized states |ψ〉, and ∥A∥∞ ⩽ ∥A∥1.
All unitarily invariant norms satisfy submultiplicativity:
The norms of interest to us are also multiplicative over tensor products and obey an ordering:
In particular, ∥AB∥1 ⩽ ∥A∥∞∥B∥1.
Another useful fact is that the partial trace is contractive, i.e.
for any operator X acting on the Hilbert space . 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 by and the extension of the upper integration limit to infinity, i.e.
where
and where the ellipsis denotes the three other summands appearing in equation (13). The Markov approximation is more accurate the smaller the error terms Δ1 and Δ2.
Consider first the Δ1 term. We shall show that it is of order g2τ3B. For simplicity we consider only the first of its four summands (the bounds for the other three are identical to that for the first):
where we used the triangle inequality, ∥Aα(t)∥∞ = ∥A∥ = 1 and ∥XY ∥1 ⩽ ∥X∥∞∥Y ∥1 (see appendix A). We can now use equation (9) to upper-bound the time derivative:
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
where we used and equation (12).
This is to be compared to the term we use after the Markov approximation:
Comparing equations (B.5) and (B.6), we see that the relative error is O[(gτB)2].
Next consider the Δ2 term. We have
The last integral converges provided 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 > τtr ∼ βln(βωc), the bath correlation function . In this case, then, we can set t = τtr and bound
which tends to zero as the cutoff tends to infinity at a fixed finite temperature (even if the lower integration limit is replaced by a constant), as expected in the weak coupling limit assumed in our work. Note that the infinite temperature limit, where equation (B.8) diverges, is incompatible with weak coupling, and requires the so-called singular coupling limit [22, 30].
Appendix C.: Properties of the spectral-density matrix Γαβ(ω)
Introducing the Fourier transform pair
and using the property that
where denotes the Cauchy principal value9, we have, using equation (14),
in agreement with equation (15), where
in agreement with equation (16).
Note that:
When ρB commutes with HB (which we have assumed), we further have
Thus the spectral-density matrix satisfies:
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:
Note that using the same technique it also follows that
To prove that this implies the frequency domain condition, equation (18), we compute the Fourier transform:
To perform this integral we replace it with a contour integral in the complex plane, , 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(τ)Ba(0)〉 is analytic in the open strip (0,−i β) and is continuous at the boundary of the strip [48]), so that
where is the integrand of equation (D.3), and the integral is the same as in equation (D.3). After making the variable transformation τ = −x − i β, where x is real, we have
Assuming that 〈Ba(±∞ − iβ)Bb(0)〉 = 0 (i.e. the correlation function vanishes at infinite time), we further have , and hence we find the result:
which, together with equation (D.3), proves equation (18).
Download figure:
Standard imageAppendix 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.
E.1. Adiabatic perturbation theory
To consider adiabatic corrections, we recall that US satisfies
where
Define the 'adiabatic intertwiner' W(t,t'):
where the 'intertwiner Hamiltonian' is
and where P0(t) ≡ |ε0(t)〉〈ε0(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
Now define V via the 'adiabatic interaction picture' transformation:
along with
Note that the time dependence of and is entirely in the energy eigenvalue and not in the eigenstates. Then V obeys the Schrödinger equation:
where
When the evolution is nearly adiabatic is a perturbation, so that we consider a solution of equation (E.8) for V of the form:
with the zeroth order solution associated with the purely adiabatic evolution, including the geometric phase:
Differentiating with respect to t, we have
where
Plugging the equation (E.10) expansion into equation (E.8), we obtain, to first order in 1/tf:10
Note that . Therefore, integrating, and using
we can write the solution as:
Therefore, the system evolution operator can be written as:
where the correction term can be made appropriately dimensionless in a more careful second order analysis, and where
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 (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 Schrödinger equation at tf, e.g. [6, 10–12].
Starting from equation (E.18) and taking matrix elements we have
Changing variables to dimensionless time s = t/tf, μba(t,t') becomes , and
where now involves a dimensionless time integration. Using the fact that , we can integrate equation (E.20) by parts as
Now note that for non-degenerate energy eigenstates, differentiation with respect to t of HS(t)|εa(t)〉 = εa(t)|εa(t)〉, and substitution of s = t/tf in equation (26), directly yields the relation:11
Substitution into equation (E.21) yields, with the help of equation (E.22),
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
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.
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) ≡ |ε0(t)〉〈ε0(t)| and . Recalling equation (E.17), we have
Using equation (E.11b) we find , so that up to a phase P0(t)Q(t,t')UadS(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.
Letting τ = t − t' > 0 and expanding around t, we then have, using equations (E.1) and (E.3),
Using the projector properties and [HS(t),P0(t)] = 0 (similarly for Q0), this simplifies to
and, after dropping the corrections, to
Inserting and into equation (E.28b), and subtracting it from equation (E.28a), we find, using once more, the desired result:
which, together with equation (E.13), proves equation (E.4b).
Appendix F.: Short time bound
We wish to bound the error associated with neglecting Θ in equation (40), i.e. we wish to bound
Using the fact that the operator satisfies:
we can write the formal solution for as:
Therefore we can bound:
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 ∥Q(t,0)∥∞ (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 , and, dropping the subdominant O(t−2f), we can write
We can now bound the error term in equation (42b). Let X(t,τ) ≡ ei τHS(t)UadS(t,0), so that US(t − τ,0) = X(t,τ) + Θ(t,τ). We can then write
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 and recalling that ∥Aα∥∞ = 1. First we assume that equation (12) applies. Then:
where in the last inequality we used the fact that if x ⩽ 2 then , and if x ⩾ 2 then again , with , in order to avoid having to extend equation (12) to higher values of n. In all, then, the approximation error in equation (42b) is .
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 τtr, followed by power-law decay. In this case the terms in (F.6b) can instead be bounded as follows:
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 τtr ∼ β ln(βωc) ≫ β (equation (75)), and that . Thus . It follows that we can safely ignore the term in equation (F.8) provided equation (77) is satisfied. The analysis of the term in (F.6c) does not change this conclusion.
Appendix G.: Derivation of the Schrödinger picture adiabatic master equation in Lindblad form
Starting from equation (53) and performing a transformation back to the Schrödinger picture, along with a double-sided adiabatic approximation, yields
Similarly, the other terms yield
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):
where in the last equality we used the freedom to interchange α and β under the summation sign. Likewise, the terms in equation (G.2b) yield
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).
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:
Recalling that the bath operators , we proceed to calculate 〈Bi(t)Bj(0)〉 = 〈B†i(t)Bj(0)〉 assuming that the bath is in a thermal Gibbs state , where is the partition function. We begin by writing:
The time evolution operator acting on the eigenstates simply produces a phase, so we focus on the operator Bj's role.
Plugging this result in, we find:
The only terms that are non-zero are the middle two, giving us:
Using the fact that:
we can write:
We can simplify our expression to
We can replace the sum with an integral:
where f(ω) is a measure of the number of oscillators at frequency ω and J(ω) is the bath spectral function. Our final result is then:
where we have assumed that oscillators with the same ωk value interact with the ith spin with the same interaction strength, i.e.
Plugging our result into equation (16), we find:
where Θ 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,
and compute the bath correlation function by inverse Fourier transform of equation (16a):
where we changed variables to x = βω.
The Polygamma function is defined as , where is the Gamma function. The Polygamma function may be represented for ℜ(z) > 0 and m > 0 as , so that in particular, for m = 1, we have
We can rewrite equation (I.2)
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
The gap at any given time in the evolution is
Since there are only two states, there are only three γ(ω) terms to calculate (assuming gω = g):
where we have taken the limit ωc → ∞ 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 σx basis, with ground state | − 〉 and excited state | + 〉, and do not change such that:
The resulting equations for the density matrix elements for small times are then
where we have used that 〈 + |σz| + 〉 = 〈 − |σz| − 〉 = 0, 〈 + |σz| − 〉 = 1 and . 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−11. 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 σ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:
where we have again assumed that the initial state is the thermal state and the gap is large.
Footnotes
- 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 ρS(t) [34]. The smallness of χ(t) is consistent with any correlations decaying very rapidly, in our case on the timescale τB 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 ρS(t) is pure then χ(t) must vanish.
- 8
- 9
By definition, , where f belongs to the set of smooth functions with compact support on the real line .
- 10
To see that this is an expansion in powers of 1/tf, transform to the dimensionless time variable s = t/tf.
- 11
Differentiating (where ), we have . Taking matrix elements gives , which yields equation (E.22) provided |εb(t)〉 is an eigenstate non-degenerate with |εa(t)〉.