Non-Markovian incoherent quantum dynamics of a two-state system
aa r X i v : . [ qu a n t - ph ] D ec Non-Markovian incoherent quantum dynamics of a two-state system
M. H. S. Amin and Frederico Brito
D-Wave Systems Inc., 100-4401 Still Creek Drive, Burnaby, B.C., V5C 6G9, Canada
We present a detailed study of the non-Markovian two-state system dynamics for the regime ofincoherent quantum tunneling. Using perturbation theory in the system tunneling amplitude ∆, andin the limit of strong system-bath coupling, we determine the short time evolution of the reduceddensity matrix and thereby find a general equation of motion for the non-Markovian evolutionat longer times. We relate the nonlocality in time due to the non-Markovian effects with theenvironment characteristic response time. In addition, we study the incoherent evolution of asystem with a double-well potential, where each well consists several quantized energy levels. Wedetermine the crossover temperature to a regime where many energy levels in the wells participatein the tunneling process, and observe that the required temperature can be much smaller than theone associated with the system plasma frequency. We also discuss experimental implications of ourtheoretical analysis.
I. INTRODUCTION
It is difficult to overemphasize the importance of thedissipative dynamics of a two-state system (TSS). In gen-eral, standing as a first hand approximation of a muchrather complex level structure, the model of a TSS cou-pled to a dissipative environment has been successfullyapplied to several physical systems. Indeed, the dissipa-tive TSS dynamics is the paradigm for the study of super-conducting devices containing Josephson junctions, two-level atoms in optical cavities, electron transfer in bio-logical and chemical systems and semiconductor quan-tum dots, to name just a few.Despite its simplicity, the description of the TSS dis-sipative dynamics imposes great theoretical challenges,especially when considering non-Markovian processes.This is the case for the analysis of the environment low-frequency noise spectrum, since the long-lived featureof its fluctuations does not allow for a “memoryless”bath (Markov) approximation. In the context of a weakTSS-bath coupling, theoretical efforts have been madeto quantify the low-frequency effect for both spin-boson and 1 / f noise models.Furthermore, it has been largely demonstrated thatlow-frequency noise plays important role in the deco-herence process of superconducting devices containingJospehson junctions. Since those devices areseen as promising candidates to the physical implemen-tation of a quantum bit, this subject has rapidly grownin interest and several studies on describing the micro-scopic origin and characterizing the low-frequency noisein such devices have already been put forward.
Understanding the evolution of a TSS also plays animportant role in understanding the performance ofan adiabatic quantum computer in the presence ofnoise. This is because for many hard problemsthe bottleneck of the computation is passing through apoint where the gap between the ground state and firstexcited state is very small. Near such an energy anticross-ing, the Hamiltonian of the system can be truncated intoa two-state Hamiltonian and in the regime of strong coupling to the environment the two-state results dis-cussed in this paper may be directly applied.Here, following a previous work, we put forward adetailed study of the TSS dissipative dynamics in thepresence of low-frequency noise, for the regime of strongTSS-bath coupling. We show that, for the regime of smalltunneling amplitude ∆, dephasing takes place much ear-lier in the evolution, leading the system to incoherentquantum dynamics. We employ such a property to de-rive equations that describe the non-Markovian evolutionof the system’s density matrix.The paper is organized as follows. In section II, wepresent the system Hamiltonian and a formal solutionfor the time evolution operator. Assuming second or-der perturbation theory in ∆, we calculate in section IIIthe short-time dynamics of the system reduced densitymatrix elements. Section IV, presents a discussion re-garding the non-Markovian behavior of the system whenthe environment is in equilibrium. We determine condi-tions under which the system reaches the detailed balanceregime. Section V provides a systematic derivation of anequation of motion for the system evolution, which ingeneral is non-local in time. We also discuss regimes inwhich the equations governing the diagonal part of thedensity matrix become t -local. Considering a double-wellpotential, section VI puts forward the analysis of intra-and interwell transitions and situations where a classicalmixture of states participate in the quantum tunnelingprocess. Finally, section VII presents our concluding re-marks. II. SYSTEM HAMILTONIAN
We start by considering an open two-state system withHamiltonian H = −
12 [∆( t ) σ x + ǫ ( t ) σ z ] − σ z Q + H B , (1)where Q is an operator acting on the environment de-scribed by the Hamiltonian H B .In order to determine the system evolution operator U ( t , t ), we proceed through two simple steps. First, wewrite the state vector of the system Hamiltonian (1) as | ψ ( t ) i = e iH B t | ϕ ( t ) i . (¯ h = k B = 1, through this paper.)Thus, one finds that the state vector | ϕ ( t ) i evolves intime according to i ∂∂t | ϕ ( t ) i = [ H ( t ) + V ( t )] | ϕ ( t ) i , where H ( t ) = − ǫ ( t ) σ z − σ z Q ( t ) , V ( t ) = −
12 ∆( t ) σ x , (2)and Q ( t ) = e iH B t Qe − iH B t . The environment is assumedto feature fluctuations following Gaussian distribution,therefore all averages can be expressed in terms of thecorrelation function or its Fourier transform, the spectraldensity: S ( ω ) = Z ∞−∞ dt e iωt h Q ( t ) Q (0) i , (3)hence we do not need to specify H B . The next step is to make use of the interaction pic-ture, considering V ( t ) as the perturbation. The statevector in the interaction picture is defined by | ϕ I ( t ) i ≡ U † ( t ) | ϕ ( t ) i , and any operator ˆ O is transformed byˆ O I ( t ) = U † ( t ) ˆ OU ( t ), with U ( t ) = T e − i R t H ( t ′ ) dt ′ = T exp (cid:26) i σ z Z t [ ǫ ( t ′ ) + Q ( t ′ )] dt ′ (cid:27) , (4)where T denotes the time ordering operator. Now, thestate evolution is determined by the interaction potential H I ( t ) = −
12 ∆( t )˜ σ x ( t ) , (5)where ˜ σ x ( t ) = U † ( t ) σ x U ( t ). The time evolution opera-tor in the interaction representation reads U I ( t , t ) = T e − i R t t H I ( t ) dt . (6)Finally, we can write a formal solution for the completetime evolution operator as U ( t , t ) = T e − i R t t H ( t ) dt = e − iH B t U ( t ) U I ( t , t ) U † ( t ) e iH B t . (7)In this paper, we are interested in the strong couplingregime in which the r.m.s. value of the noise W ≡ p h Q i = (cid:18)Z ∞−∞ dω π S ( ω ) (cid:19) / , (8)is much larger than the tunneling amplitude: W ≫ ∆.Physically, W is basically the uncertainty in the energybias ǫ ( t ) and therefore represents the broadening of theenergy levels. In the above regime, consequently, the broadening of the energy levels is larger than the mini-mum gap and therefore the gap will not be well-defined.On the other hand, as we shall see, for the case of lowfrequency noise, W represents the dephasing rate of thesystem. Thus, the above regime is a limit in which thequbit loses quantum coherence before it can tunnel, i.e.,the dynamics is incoherent. III. DENSITY MATRIX CALCULATION
We would like to study the evolution of the reduceddensity matrix. Let ρ SB ( t ) denote the total density ma-trix of the system plus bath. We therefore have ρ SB ( t ) = U ( t, ρ SB (0) U † ( t,
0) (9)= e − iH B t U ( t ) U I ( t, ρ SB (0) U † I ( t, U † ( t ) e iH B t . The system reduced density matrix is defined by ρ ( t ) =Tr B [ ρ SB ( t )], where Tr B [ ... ] means averaging over all en-vironmental modes. We assume that the density matrixat t = 0 is separable, i.e., ρ SB (0) = ρ (0) ⊗ ρ B , where ρ B = e − H B /T is the density matrix of the environment,which we assume to be in equilibrium at temperature T . Under the assumption of separability of the initialdensity matrix, we consider that the system evolutionimmediately follows an initialization in a definite state,implemented, e.g., through a state measurement.If ∆ is the smallest energy scale in the problem, wecan approximate U I ( t,
0) by performing a perturbationexpansion in ∆, which up to second order reads U I ( t, ≈ i Z t dt ′ ∆( t ′ )˜ σ x ( t ′ ) − Z t Z t ′ dt ′ dt ′′ ∆( t ′ )∆( t ′′ )˜ σ x ( t ′ )˜ σ x ( t ′′ ) . (10)If the time interval t is not small enough to make theabove integrals small, i.e., t > ∼ / ∆, the higher orderterms in ∆ must be considered in the expansion. A. Off-diagonal elements of ρ To zeroth order in ∆, we have U I ( t,
0) = 1, therefore ρ SB ( t ) = e − iH B t U ( t ) ρ SB (0) U † ( t ) e iH B t . (11)For this case, since [ U ( t ) , σ z ] = 0, the σ z populationsof the system reduced density matrix ρ are constants ofmotion. Therefore, in the representation of the eigen-states of σ z , σ z | i = − | i ( σ z | i = | i ), only the off-diagonal elements of ρ present dynamics, which, due tothe coupling to environment, decay in time. This caseconstitutes the one of a pure dephasing dynamics. It hasbeen subject of interest for many areas, where severalapproaches have been used to calculate the off-diagonalelements of ρ . Few examples are the (a) spin-bosonmodel assuming a power-law behavior for the spectraldensity of the bath ; (b) spin-fermion model , and(c) spin-two-state fluctuators system . Here, we considera bosonic bath, but do not have to specify the form ofthe bath spectral density. To quantify this decay, let uswrite the reduced density matrix as ρ ( t ) = X i,j =0 , ρ ij ( t ) | i ih j | . (12)We find for the off-diagonal element ρ ( t ) = Tr B [ h | U ( t ) ρ SB (0) U † ( t ) | i ] = ρ (0) × e − i R t ǫ ( t ′ ) dt ′ (cid:28) ←−T e − i R t Q ( t ′ ) dt ′ T e − i R t Q ( t ′ ) dt ′ (cid:29) , (13)where h ... i ≡ Tr B [ ρ B ... ] and ←−T represents the reversetime ordering operator. We expand the exponentials,group those in the same order, take the average of eachterm, and bring them back to the exponent. Because ofthe Gaussian nature of the environment, it is sufficientto expand up to second order in Q . Assuming the envi-ronment is in equilibrium, one finds (cid:28) ←−T e − i R t Q ( t ′ ) dt ′ T e − i R t Q ( t ′ ) dt ′ (cid:29) == e − R t dt ′ R t dt ′′ h Q ( t ′ ) Q ( t ′′ ) i = e − R dω π R t dt ′ R t dt ′′ e iω ( t ′′− t ′ ) S ( ω ) . (14)Thus, using (14) in (13), we obtain ρ ( t ) = e − i R t ǫ ( t ′ ) dt ′ × exp (cid:26) − Z dωπ S ( ω ) sin ( ωt/ ω (cid:27) ρ (0) . (15)This equation represents a complicated decay rate, whichis in general not a simple exponential function of t . How-ever, in two limits it can be simplified. First, for the caseof white noise, i.e., S ( ω ) = S (0), it gives ρ ( t ) = e − i R t ǫ ( t ′ ) dt ′ − S (0) t ρ (0) . (16)Which leads to dephasing rate 1 /T = S (0), as expectedfor white noise.Another interesting limit is when S ( ω ) is dominatedby low frequencies so that one can use sin x ≈ x to get ρ ( t ) = e − i R t ǫ ( t ′ ) dt ′ − W t ρ (0) , (17)where W is the energy level broadening given by (8).The decay is now a Gaussian, whose width determinesthe dephasing rate, 1 /T φ = W . For the case of 1/f noise,where the cutoff of S ( ω ) is not sharp enough, one gets alogarithmic correction to the above equation . B. Diagonal elements of ρ The evolution of the diagonal part of the density ma-trix happens in a time scale much larger than 1 /W . Thecomplete evolution is given by ρ ( t ) = Tr B [ U ( t ) U I ( t, ρ (0) ρ B U † I ( t, U † ( t )] . (18)Let us assume the initial condition ρ (0) = | ih | andtry to calculate ρ ( t ). To zeroth order in ∆, we have ρ ( t ) = 0 as expected, thus we find that the first nonzerocontribution to ρ ( t ) comes from the first-order term in∆ of (10): ρ ( t ) ≈ Z t dt Z t dt ∆( t )∆( t ) × Tr B [ h | ˜ σ x ( t ) | i ρ B h | ˜ σ x ( t ) | i ]= 14 Z t dt Z t dt ∆( t )∆( t ) × Tr B [ h | U † ( t ) U ∗ ( t ) | i ρ B h | U † ( t ) U ∗ ( t ) | i ]= 14 Z t dt Z t dt ∆( t )∆( t ) e i R t t ǫ ( t ′ ) dt ′ × (cid:28) ←−T e i R t Q ( t ′ ) dt ′ T e i R t Q ( t ′ ) dt ′ ←−T e − i R t Q ( t ′ ) dt ′ T e − i R t Q ( t ′ ) dt ′ (cid:29) , (19)where in the second equality we have used ˜ σ x ( t ) = U † ( t ) σ x U ( t ) = U † ( t ) U ∗ ( t ) σ x . One can calculate the ex-pectation value by expanding the exponentials and keep-ing only the the second order terms. The last two linesof (19) become1 + 12 Z t dt ′ Z t t dt ′′ h Q ( t ′ ) Q ( t ′′ ) i + 12 Z t t dt ′ Z t dt ′′ h Q ( t ′ ) Q ( t ′′ ) i . (20)Substituting the inverse Fourier transformation h Q ( t ′ ) Q ( t ′′ ) i = R dω π e − iω ( t ′ − t ′′ ) S ( ω ), we find1 + Z dω π S ( ω ) ω [ e iω ( t − t ) − i (sin ωt − sin ωt )]= 1 + Z dω π S ( ω ) ω (cos ωτ − − i Z dω π S ( ω ) ω (sin ωτ − ωτ ωτ ′ ) , (21)where τ = t − t and τ ′ = ( t + t ) / S ( ω ) is dominated by lowfrequency noise such that for all relevant modes ωτ ≪ ωτ and cos ωτ in (21) to get ρ ( t ) ≈ Z t dτ ′ Z ˜ t − ˜ t dτ ∆( τ ′ + τ τ ′ − τ × e − W τ / − i (cid:16) ǫ p ( τ ′ ) τ − R τ/ − τ/ ǫ ( τ ′ + t ′ ) dt ′ (cid:17) , (22)where ˜ t = min[2 τ ′ , t − τ ′ )], W is given by (8), and ǫ p ( t ) ≡ Z dω π S ( ω ) ω (1 − cos ωt ) . (23)Equation (22) conveys the non-locality in time, expectedfor a non-Markovian environment. If within time τ ∼ /W , ǫ ( t ) and ∆( t ) do not change much (or even if ∆( t )is a fast but linear exponential function), we can write(22) as ρ ( t ) ≈ Z t dτ ′ ∆ ( τ ′ ) Z ˜ t − ˜ t dτ e i [ ǫ ( τ ′ ) − ǫ p ( τ ′ )] τ − W τ / . (24)Therefore, for t < ∼ / ∆( t ), we find the leading term forsystem population rate change given by˙ ρ ( t ) ≈ ∆ ( t )4 Z t − t dτ e i [ ǫ ( t ) − ǫ p ( t )] τ − W τ / . (25)If t > /W , due to gaussian envelope of the integrand,we can extend the integration limits of (25) to ±∞ , ob-taining ˙ ρ ( t ) ≈ Γ p e − [ ǫ ( t ) − ǫ p ( t )] / W , (26)with the peak value of the functions given byΓ p ≡ r π W . (27)It is worth recalling that for times t > ∼ / ∆, Eq. (10) doesnot represent a fair approximation to U I ( t, ǫ p ( t ) behavior. IV. MACROSCOPIC RESONANT TUNNELING
Should ǫ p be constant in time and the Hamiltonian (1)be time independent, one could directly read (26) as anapproximation for the equation of motion˙ ρ ( t ) ≈ Γ − ρ ( t ) − Γ + ρ ( t ) , (28)since the off-diagonal elements of ρ ( t ) become negligiblefor times t > ∼ /W . The rate Γ − (Γ + ) then represents the | i → | i ( | i → | i ) system transition rate. Thus, for theregime 1 /W < ∼ t < ∼ / ∆( t ), the evolution is described by(26). The same argument holds when ǫ p ( t ) is a functionof time, but in that case the tunneling rates will be timedependent: Γ ± ( t ) = Γ p e − [ ǫ ± ǫ p ( t )] / W . (29) An experimental realization of such a tunneling processin a macroscopic quantum device such as a supercon-ducting flux qubit is called macroscopic resonant tunnel-ing (MRT). The tunneling rates Γ ± are therefore simpleshifted Gaussian functions described by (29). An imme-diate consequence of (23) is that the shift ǫ p vanishes fora classical noise, for which S ( ω ) is symmetric. Therefore,a nonzero value of ǫ p is a signature for quantum natureof the noise source.If the environmental source is in equilibrium at tem-perature T , then the symmetric and antisymmetric (infrequency) parts of the noise intensity are related by thefluctuation-dissipation theorem: S s ( ω ) = S a ( ω ) coth (cid:16) ω T (cid:17) (30)Therefore the fluctuation-dissipation theorem relates W and ǫ p ( t ), which are functions of S s and S a respectively.Let us first define ǫ p = P Z dω π S ( ω ) ω , (31)with P representing principal value integral. In the caseof low-frequency noise, when all the relevant frequenciesare small on the scale of temperature T , i.e., ω ≪ T , onecan write coth( ω/ T ) ≃ T /ω . In that case (8), (30),and (31) yield W = 2 T ǫ p . (32)One therefore finds ǫ p ( t ) = ǫ p − P Z dω π S ( ω ) ω cos ωt. (33)The effect of the last term depends on how small or large t is with respect to the time response, τ R ∼ ω − c , ofthe environment. Here, ω c represents the characteristicenergy of the environment. To understand this let usconsider different regimes. A. Large ω c (short τ R ) limit If ω c is large compared to 1 /t , where t is the typi-cal time scale of interest, then the integral in (33) cov-ers many oscillations of the cosine function and thereforevanishes. In that case ǫ p ≈ ǫ p = W T , (34)consequently being independent of t . For a time indepen-dent Hamiltonian, Eq. (29) then yieldsΓ ± ( ǫ ) = Γ p e − [ ǫ ± ǫ p ] / W , (35)in agreement with Ref. 24. It is easy to see thatΓ − ( ǫ )Γ + ( ǫ ) = e ǫ/T , (36)which (in the limit ∆ →
0) is the detailed balance (Ein-stein) relation. Therefore, the transition rates (35)support thermal equilibrium distribution of the systemstates, which is a natural consequence of the fast envi-ronmental response.
B. Small ω c (long τ R ) limit If the environment’s response is slow compared to timescale of the problem, i.e., ω c ≪ /t , the cosine function in(33) will be close to 1 at all times, making ǫ p ≈
0, again independent of t . For a time independent Hamiltonian,therefore, we getΓ − = Γ + = Γ p e − ǫ / W . (37)Such transitions obviously do not satisfy the detailed bal-ance relation and do not lead to equilibrium distribution.Indeed, an environment in ω c → Q that doesnot vary much during the evolution and has a Gaussiandistribution: P ( Q ) = e − Q / W √ πW . (38)In small ∆ regime, the tunneling rate from state | i i tostate | j i can be calculated using the Fermi Golden ruleΓ i → j = 2 π |h i | V | j i| δ ( E i − E j ) , (39)where V = ∆ σ x / Q , one findsΓ − ( Q ) = Γ + ( Q ) = π ∆ δ ( ǫ + Q ) (40)Averaging over all possibilities of Q, we findΓ − = Γ + = π ∆ Z dQP ( Q ) δ ( ǫ + Q )= Γ p e − ǫ / W , (41)which is the same as (37). C. General ω c ( τ R ) regime In general, away from the above two limits, ǫ p ( t ) isa time dependent function given by (33). The explicitfunctionality depends on the spectral density S ( ω ), espe-cially on its characteristic frequency ω c . To see this, letus assume a simple spectral density S ( ω ) = 2 ηω [1 + ( ω/ω c ) ] (cid:18) − e − ω/T (cid:19) , (42) for which analytical solutions is possible. Substituting(42) in (33), we find ǫ p ( t ) = Z dω π η (1 − cos ωt )[1 + ( ω/ω c ) ] = ηω c − e − ω c t (1 + tω c )](43)We can therefore write ǫ p ( t ) = ǫ p (1 − e − ω c t (1 + tω c )) = (cid:26) , ω c t ≪ ǫ p , ω c t ≫ , (44)which yields the above two limiting results in the ap-propriate regimes with an exponential crossover betweenthe two limits. Indeed, the above behavior of ǫ p ( t ), i.e.,the crossover from 0 to ǫ p within time scale ∼ /ω c , isgeneric regardless of the functional detail of ǫ p ( t ). Thetime scale τ R ∼ /ω c represents the response time of theenvironment to an external perturbation. If t ≫ τ R , thenthe environment has enough time to enforce equilibriumto the system, resulting in ǫ p = ǫ p , which is requiredfor detailed balance (i.e., equilibrium) condition. On theother hand, if t ≪ τ R , the environment cannot respondquickly to the system and the equilibrium relation is notexpected. In that case, we find ǫ p = 0, i.e., the environ-ment behaves as a classical noise. In the next section weshall see how such behavior results in time-nonlocality ofthe equation of motion. V. NON-MARKOVIAN EQUATION OFMOTION
Equation (26) gives the short time (1 /W < ∼ t < ∼ / ∆)evolution of the diagonal part of the density matrix. Assoon as t becomes comparable to ∆, higher order correc-tions become important and the second order perturba-tion used in Eq. (26) becomes insufficient. Instead ofintroducing higher order corrections which is a cumber-some task, in this section we take a different path: Wewrite a general equation of motion expected for the evolu-tion of the density matrix for a system like ours and findits parameters in such a way that it agrees with Eq. (26)for short times.In general, the equation of motion for the evolutionof the density matrix is nonlocal in time, reflecting thenon-Markovian nature of the environment. Since the off-diagonal elements vanish very quickly (within t ∼ /W ),for time scales larger than 1 /W , one can write the dy-namical equations only in terms of the diagonal elementsof ρ . Generally, for a non-Markovian dynamics the equa-tion of motion for ρ ( t ) depends on the history˙ ρ ( t )= Z t −∞ dt ′ [ K − ( t, t ′ ) ρ ( t ′ ) − K + ( t, t ′ ) ρ ( t ′ )] , (45)where K ± ( t, t ′ ) are nonlocal integration kernels. Let usfrom now onwards consider a time-invariant Hamiltonianfor which˙ ρ ( t )= Z t −∞ dt ′ [ K − ( t − t ′ ) ρ ( t ′ ) − K + ( t − t ′ ) ρ ( t ′ )] . (46)If the system starts the evolution from state | i at time t = t , the short time evolution is described by˙ ρ ( t ) ≈ Z tt dt ′ K − ( t − t ′ ) = Z t − t dτ K − ( τ ) . (47)This should agree with (26), therefore Z t − t dτ K ± ( τ ) = Λ ± ( t − t ) θ ( t − t ) . (48)where we have defined functionsΛ ± ( t ) ≡ Γ p e − [ ǫ ± ǫ p ( t )] / W . (49)The presence of the θ -function is necessary to ensurecausality to the system dynamics, since we assume thatthe evolution follows a state initialization at t . Takingthe derivative of both sides of (48), we find K ± ( τ ) = ∂ Λ ± ( τ ) ∂τ θ ( τ ) + Λ ± ( τ ) δ ( τ ) . (50)Notice that for constant transition rates Λ ± ( τ ) = Γ ± ,Eq. (50) leads to˙ ρ ( t ) = Γ − ρ ( t ) − Γ + ρ ( t ) , (51)which, as expected, is t -local.In the limit ω c →
0, where the change in Λ ± happenson a time scale (1 /ω c ) much larger than the time evo-lution considered here, the time derivative in (50) canbe neglected and one obtains (51) with transition ratesΓ ± = Λ ± (0) = Γ p e − ǫ / W , with ǫ p ( t ) = 0, as expectedfor a static noise.On the other hand, in the ω c → ∞ limit, variations ofΛ ± ( t ) happen in a very short time, hence ∂ Λ ± ( τ ) /∂τ → t > ∼ τ R ∼ /ω c . Therefore the t ′ -integration in (46)is basically between t − τ R and t . If within this shortrange ρ ( t ′ ) does not change much, one can bring it outsidethe integral. In that case, (46) leads to (51) with Γ ± =Λ ± ( t → ∞ ) = Γ p e − ( ǫ ± ǫ p ) / W , with ǫ p ( t ) = ǫ p , whichis expected in the detailed balance regime.Both of the above regimes led to t -local equations forthe diagonal part of the density matrix. However, for fi-nite ω c , in general, one gets a nonlocal equation in time.If the system evolution is slow compared to the timescale τ R ∼ /ω c , one can substitute the Taylor expan-sion ρ ij ( t ′ ) = ρ ij ( t ) + ( t ′ − t ) ˙ ρ ij ( t ) into (46) obtaining˙ ρ ( t ) = Λ − ( t ) ρ ( t ) − Λ + ( t ) ρ ( t )+ ˙ ρ ( t ) Z tt dt ′ ∂ Λ( t − t ′ ) ∂t ( t − t ′ ) , (52)where Λ( t ) = Λ − ( t ) + Λ + ( t ). Solving for ˙ ρ ( t ), one finds(51) with transition ratesΓ ± = Λ ± ( ∞ )1 − R ∞ dτ τ ∂ Λ( τ ) /∂τ = Λ ± ( ∞ )1 − R ∞ dτ [Λ( ∞ ) − Λ( τ )] , (53) where in the last step we have used integration by parts.The integral limit was taken to infinity, since the inte-grand very quickly vanishes for τ > ∼ /ω c . All the nonlo-cal behavior is captured in the denominator of (53). Theintegrand (53) is maximum at τ = 0, but very quicklyvanishes within τ ∼ /ω c , hence R ∞ dτ [Λ( ∞ ) − Λ( τ )] ∼ [Λ( ∞ ) − Λ(0)] /ω c , leading toΓ ± ≈ Λ ± ( ∞ )1 − [Λ( ∞ ) − Λ(0)] /ω c . (54)Using (49), we obtainΛ( ∞ ) − Λ(0) = Γ p ( e − ( ǫ − ǫ p ) / W + e − ( ǫ + ǫ p ) / W − e − ǫ / W )= 2Γ p e − ǫ / W (cid:16) e − ǫ p / W cosh ǫ T − (cid:17) . (55)Therefore, to the lowest order in Γ p /ω c , we getΓ ± ( ǫ ) ≈ Γ p e − ( ǫ ± ǫ p ) / W (cid:26) p ω c e − ǫ / W (cid:16) e − ǫ p / W cosh ǫ T − (cid:17) (cid:27) . (56)The magnitude and the position of the peak of Γ − ( ǫ ) aregiven by (to the lowest order in Γ p /ω c )Γ peak ≈ Γ − ( ǫ p ) ≈ Γ p (1 + Γ p /ω c ) ,ǫ peak ≈ ǫ p (cid:18) p ω c e − ǫ p / W (cid:19) . (57)The peak value is enhanced by the nonlocal effects. Thepeak position is also shifted, but by a very small amountdue to the exponential suppression. Notice that the peakbecomes asymmetric around its center due to the nonlo-cality.The nonlocal corrections to the transition rates becomenegligible when Γ p ≪ ω c . Also, observe that Γ p is ap-proximately the peak value of the transition rate (27).Therefore, nonlocality becomes important only when themaximum transition rate Γ p is of the order of or largerthan ω c , or equivalently, when the response time ( τ R ) ofthe environment is comparable or longer than the systemtransition time ( ∼ / Γ p ). VI. MRT IN A DOUBLE-WELL POTENTIAL
So far we have studied incoherent tunneling in an ide-alized two state system. However, for most realistic sys-tems, the two state model is only an approximation ofa more complicated multi-level problem. An example ofsuch cases is a system in which the classical potentialenergy has a double-well structure and the kinetic partof the Hamiltonian provides quantum tunneling betweenthe two wells. Experimental implementation of sucha system is possible using superconducting flux qubits,which have been studied considerably both theoreticallyand experimentally . Especially,MRT measurements have been performed both betweenground states as well as excited states of the wells.
In such a double-well system, the energy within eachwell is quantized, with energy level distributions depen-dent on the bias energy between the wells. In general, inthe presence of the environment, a system initialized inone of those levels can experience two types of evolution:intra- and interwell dynamics.The intrawell dynamics are transitions within a singlewell, e.g., when a system excited within a well relaxes to alower energy level in the same well by exchanging energywith the environment. Thus, in this case, the systemdynamics is confined in just one well of the potential,with no tunneling to the opposite well.It is also possible for the system, depending on thetunneling amplitude between the two states, to tunnelto an energy level in the opposite well, leading thus toan interwell dynamics. If the evolution of the system isconfined to the ground states of the two wells and it liesin the incoherent tunneling regime, then the formalismdeveloped herein can describe such an evolution in fulldetail. This, however, is not the only type of evolutionpossible for a double-well system. Here, we also considerpossibilities that the evolution involves the excited states.
A. Tunneling between ground states
At low enough temperatures, the system can only oc-cupy the lowest energy states within the wells. In sucha case, tunneling can occur between those energy levelsif the levels become in resonance. The probability of thesystem being found in state | i at time t is given by ρ ( t ).For a time independent system initialized in state | i , inthe limit Γ p ≪ ω c , ρ ( t ) is the solution of (51). Sucha t -dependence can be measured experimentally and isusually an exponential function with initial value 0 andfinal value given by the equilibrium distribution. Accord-ing to (51), the initial slope of ρ ( t ) gives the transitionrate: Γ − = ˙ ρ (0). Likewise, if the system is initial-ized in state | i , one can extract Γ + in a similar way.Plotting the resulting transition rates versus bias ǫ , oneobtains the tunneling resonant peaks. By fitting the ex-perimental data to the shifted Gaussian line-shapes (35)the parameters ǫ p , W , and Γ p can be extracted and from(27), ∆ can be obtained. Such a procedure, performed inRef. 38, successfully confirmed our theory especially therelation (32) between W and ǫ p .If the transition rate Γ p becomes comparable to theenvironment’s characteristic energy ω c , the t -local equa-tion (51) will not be adequate to describe the evolutionof the system. However, if the nonlocality effect is small,one can still use the same equation but with Γ ± definedby (54), hence (56). In such a case, the peak will notbe symmetric around its center, with an asymmetry thatdependens on ∆. Experimental observation of such an asymmetry is an indication of time-delayed response ofthe environment and may provide information about ω c .It should be reminded that a presence of high frequencymodes in S ( ω ) may also lead to deviation from a symmet-ric Gaussian MRT peak but such an effect is independentof ∆ hence could be easily distinguished from the abovenonlocal effects.Another interesting type of experiment is the Landau-Zener transition in which ǫ is a linear function of timeduring the evolution. For that type of evolution, again inthe Γ p ≪ ω c regime, one can still use (51) but with a timedependent ǫ . Such a procedure was proved successful inproviding accurate description of the experimental datafor flux qubits in Ref. 40.It should be mentioned that the tunneling rate ∆ in ourformalism may not be independent of ǫ as assumed here.In practice, as the double-well potential is tilted, it notonly affects the relative positions of the energy levels inthe two wells but also affects the matrix elements betweenthem. Such dependence is weak for a small bias, butas ǫ becomes large the effect of modulation of ∆ mightbecome visible. B. Tunneling to or between excited states
If the energy tilt is large enough so that the groundstate of the initial well becomes in resonance with an ex-cited state of the opposite well, tunneling to the excitedstate can occur. Alternatively, one may initialize thesystem in an excited state in the initial well, via e.g., mi-crowave excitation, and make the system tunnel betweentwo excited states. It is therefore important to under-stand how such a tunneling can be described within thepresent theory. One can generalize the arguments of theprevious section to calculate the tunneling rate. In thiscase, we need to add intrawell relaxations to the picture.In Ref. 24, it was shown that the tunneling rate fromstate | i i in the left well to state | j i in the right well isgiven byΓ ij ( ǫ ) = ∆ ij Z ∞−∞ dt e i ( ǫ − ǫ p ) t − γ ij | t |− W t , (58)where ǫ is the bias energy with respect to the resonancepoint between | i i and | j i , ∆ ij is the tunneling amplitudebetween the two states, and γ ij = ( γ i + γ j ) /
2, with γ i be-ing the intrawell relaxation rates corresponding to state | i i . If one of the states is the ground state in its own well,then its corresponding intrawell relaxation rate is zero.The transition rate becomes a convolution of Lorentzianand Gaussian functions:Γ ij ( ǫ ) = ∆ ij γ ij √ πW Z ∞−∞ dǫ ′ e − [ ǫ ′ − ǫ p ] / W ( ǫ − ǫ ′ ) + γ ij = r π ij W Re (cid:20) w (cid:18) ǫ ± ǫ p + iγ ij √ W (cid:19)(cid:21) , (59)where w ( x ) = e − x [1 − erf( − ix )] = 2 e − x √ π Z ∞ ix e − t dt (60)is the complex error function. In the limit γ ij →
0, theshifted Gaussian line-shape (35) is recovered. In the op-posite limit, γ ij ≫ W , the peak becomes a Lorentzianwith width γ ij . C. Multi-channel tunneling
So far, we have investigated the dynamics of a defi-nite single tunneling event between the wells. However,as the system’s temperature increases, one should expectthe increase of probability of thermal occupation of theexcited states of each well. Under such conditions, itbecomes unknown what single tunneling event will takeplace. Consequently, when predicting the effective tun-neling rate between wells, one has to take into accountthe statistics of occupation of excited states and their re-spective tunneling probabilities to the opposite well, inan ensemble average. The net of this thermally assisteddynamics is a multi-channel tunneling, which leads to anincrease of the measured tunneling rate. As we shall see,due to the fast increase of the tunneling amplitude ∆ ij between excited states | i i and | j i , T does not need to betoo large for this process to become non-negligible. Forsimplicity we consider zero bias ( ǫ = 0) situation in whichthe two potential wells are in resonance.Let ∆ n and Γ n ± denote the tunneling amplitude andtransition rate between the n -th energy levels in the op-posite wells, and Γ ± the total transition rates between thewells. In thermal equilibrium, the occupation probabil-ity of the n -th state is given by Boltzmann distribution: P n = e − E n /T / P i e − E i /T . ThereforeΓ ± ( ǫ ) = X n P n Γ n ± ( ǫ ) . (61)At small enough T , one can assume P n ≈ e − E n /T (for n > E n = E n − E is the relative energy ofstate | n i compared to the ground state ( n = 0). If γ ij ≪ W for the low-lying energy levels, we may neglect γ ij in(59) and all Γ n will have the same Gaussian functionalform, leading toΓ − ( ǫ ) = X n e − E n /T r π n W e − [ ǫ − ǫ p ] / W , = r π eff ( T ) W e − [ ǫ − ǫ p ] / W , (62)where ∆ eff = ∆ X n ≥ ∆ n ∆ e − E n /T / . (63) Therefore, the net contribution from tunneling events in-volving excited states can be seen as a renormalizationof the tunneling amplitude between wells. Since usually∆ n ≫ ∆ , such contribution becomes important even attemperatures much smaller than the plasma frequency ω p ≡ E . The crossover temperature T co can be ob-tained by requiring (∆ / ∆ ) e − ω p /T ∼
1, such that thecontribution from the first excited state becomes impor-tant: T co = ω p / ∆ ) . (64)Typically ∆ is a few orders of magnitude larger than ∆ and therefore T co can be an order of magnitude smallerthan ω p . High frequency modes of environment may alsorenormalize the tunneling amplitude resulting in a T -dependent ∆ eff even at T < T co . Such a T -dependence istypically much weaker and a crossover to the exponentialdependence in (63) should be observable. VII. CONCLUSIONS
We have shown a systematic procedure to determinethe evolution of a two-state system in the regime of inco-herent quantum dynamics. Considering a second orderperturbation theory in the system bare tunneling rate ∆,and a Gaussian distribution for the environment fluctu-ations, we have determined the short time evolution ofthe system reduced density matrix elements.Under the assumption of high integrated noise W , i.e.,a system-bath strong coupling regime, we verify that,indeed, dephasing process takes place early in the systemevolution, which sets 1 /W as the smallest time scale ofthe evolution, justifying the claim of having a systemwith incoherent dynamics.As for the system populations, we have seen that, ingeneral, one should expect complex non-Markovian dy-namics. We were able to clearly demonstrate how thenon-Markovian evolution can be related to the time re-sponse of the environment, τ R . Indeed, we have veri-fied that for time scales t ≫ τ R , the system follows thedetailed balance dynamics. On the other hand, if theenvironment response is very slow, i.e., t ≪ τ R , the sys-tem sees a static (classical) noise source. In addition,by investigating the equation of motion for the reduceddensity matrix, we have demonstrated how one can sim-plify the non-Markovian effects by introducing modifiedtransition rates for the dynamical equations.Finally, we have inspected the intra- and interwelltransition possibilities inside a double-well potential, andquantified how the multi-channel process can lead to anenhancement of the system tunneling. We have deter-mined the condition for this process to take place, andestimated the crossover temperature which can be an or-der of magnitude smaller than the system plasma fre-quency ω p .Some of the predictions of our theory have already beenconfirmed experimentally. More experiments, how-ever, are necessary especially to confirm our descriptionof non-Markovian dynamics. A simple measure of theasymmetry of the MRT peak in large ∆ regime could beindicative of nonlocality in t . As described in Sec. V, suchan asymmetry should be ∆ dependent and should dis-appear at small ∆. A ∆-independent asymmetry couldresult from high frequency components of the environ-mental noise that make small ωτ expansion in (21) fail.Moreover, a T -dependent measure of the tunneling ratescan reveal the renormalization of the effective tunnel- ing amplitude ∆ due to high frequency noise and thecrossover temperature T co to the multichannel tunnelingregime as described in section VI. Acknowledgments
We would like to thank D.V. Averin, A.J. Berkley, R.Harris, J. Johansson and T. Lanting for useful discussionsand a critical reading of this manuscript. A.J. Leggett, S. Chakravarty, A.T. Dorsey, M.P.A. Fisher,A. Garg and W. Zwerger, Rev. Mod. Phys. , 1 (1987). U. Weiss,
Quantum Dissipative Systems , 2nd ed. (Singa-pore: World Scientic, 1999). Y. Makhlin, G. Sch¨on and A. Shnirman, Rev. Mod. Phys. , 357 (2001). C. Tannoudji, C. Dupont-Roc and J. Grynberg,
Atom Pho-ton Interaction: Basic Processes and Applications (NewYork: Wiley, 1992). A. Garg, J. Onuchic and V. Ambegaokar, J. Chem. Phys. , 4491 (1985). R. Hanson, L.P. Kouwenhoven, J.R. Petta, S. Tarucha andL.M.K. Vandersypen, Rev. Mod. Phys. , 1217 (2007). D.P. DiVincenzo and D. Loss, Phys. Rev. B , 035318(2005). J. Schriefl, Yu. Makhlin, A. Shnirman and G. Sch¨on, NewJ. Phys. , 1 (2006). G. Burkard, Phys. Rev. B , 125317 (2009). O. Astafiev, Yu.A. Pashkin, Y. Nakamura, T. Yamamotoand J.S. Tsai, Phys. Rev. Lett. , 267007 (2004). M. M¨uck, M. Korn, C.G.A. Mugford, J.B. Kycia and J.Clarke Appl. Phys. Lett. , 012510 (2005). F. Yoshihara, K. Harrabi, A. O. Niskanen, Y. Nakamuraand J.S. Tsai, Phys. Rev. Lett. , 167001 (2006). K. Kakuyanagi, T. Meno, S. Saito, H. Nakano, K. Semba,H. Takayanagi, F. Deppe, and A. Shnirman, Phys. Rev.Lett. R.C. Bialczak, R. McDermott, M. Ansmann, M. Hofheinz,N. Katz, E. Lucero, M. Neeley, A.D. O’Connell, H. Wang,A.N. Cleland and J.M. Martinis, Phys. Rev. Lett. ,187006 (2007). R.H. Koch, D.P. DiVincenzo, and J. Clarke, Phys. Rev.Lett. , 267003 (2007). L. Faoro and L.B. Ioffe, Phys. Rev. Lett. , 227005(2008). T. Lanting, A.J. Berkley, B. Bumble, P. Bunyk, A. Fung,J. Johansson, A. Kaul, A. Kleinsasser, E. Ladizinsky, F.Maibaum, R. Harris, M.W. Johnson, E. Tolkacheva andM.H.S. Amin, Phys. Rev. B , 060509(R) (2009). S. Sendelbach, D. Hover, A. Kittel, M. M¨uck, J.M. Mar-tinis and R. McDermott, Phys. Rev. Lett. , 227006(2008). E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lund-gren, and D. Preda, Science , 472 (2001). S. Ashhab, J. R. Johansson, and Franco Nori, Phys. Rev.A , 052330 (2006). M.H.S. Amin, P.J. Love, and C.J.S. Truncik, Phys. Rev.Lett. , 060503 (2008). M.H.S. Amin, D.V. Averin, and J.A. Nesteroff, Phys. Rev.A , 022107 (2009). M.H.S. Amin, C.J.S. Truncik, and D.V. Averin, to appearin Phys. Rev. A, eprint arXiv:0803.1196. M.H.S. Amin and D.V. Averin, Phys. Rev. Lett. ,197001 (2008). Recently, it has been argued that the noise produced by asystem of two-state fluctuators can exhibit non-Gaussiandistribution, playing important role in the qubit-fluctuatorstrong coupling limit. See, e.g., Ref. 26. Y. M. Galperin, B. L. Altshuler, J. Bergli and D. V. Shant-sev, Phys. Rev. Lett. , 097009 (2006); Y. M. Galperin,B. L. Altshuler, J. Bergli, D. Shantsev and V. Vinokur,Phys. Rev. B , 064531 (2007); J. Bergli, Y. M. Galperinand B. L. Altshuler, New J. Phys. , 025002 (2009). G. M. Palma, Kalle-Antti Suomine and A. K. Ekert, Proc.R. Soc. Lond. A , 567 (1996). J. H. Reina, L. Quiroga and N. F. Johnson, Phys. Rev. A , 032326 (2002). Alex Grishin, Igor V. Yurkevich and Igor V. Lerner, Phys.Rev. B , 060509 (2005). Roman M. Lutchyn, Lukasz Cywinski, Cody P. Nave andS. Das Sarma, Phys. Rev. B , 024508 (2008) J.R. Friedman, V. Patel, W. Chen, S.K. Tolpygo and J.E.Lukens, Nature , 43 (2000). I. Chiorescu, Y. Nakamura, C.J.P.M. Harmans and J.E.Mooij, Science , 1869 (2003). G. Burkard, R.H. Koch, and D.P. DiVincenzo, Phys. Rev.B , 064503 (2004). T. Hime, P.A. Reichardt, B.L.T. Plourde, T. L. Robertson,C.-E. Wu, A.V. Ustinov and J. Clarke, Science , 1427(2006). R.H. Koch, G.A. Keefe, F.P. Milliken, J.R. Rozen, C.C.Tsuei, J.R. Kirtley and D. P. DiVincenzo, Phys. Rev. Lett. , 127001 (2006). J. H. Plantenberg, P. C. de Groot, C.J.P.M. Harmans andJ.E. Mooij, Nature , 836 (2007). R. Rouse, S. Han, and J. E. Lukens, Phys. Rev. Lett. 75,1614 (1995). R. Harris, M. W. Johnson, S. Han, A. J. Berkley, J. Johans-son, P. Bunyk, E. Ladizinsky, S. Govorkov, M. C. Thom,S. Uchaikin, B. Bumble, A. Fung, A. Kaul, A. Kleinsasser,M. H. S. Amin and D. V. Averin, Phys. Rev. Lett. ,117003 (2008). D.A. Bennett, L. Longobardi, V. Patel, W. Chen, D.V.Averin and J.E. Lukens, Quant. Inf. Process. , 217 (2009). J. Johansson, M.H.S. Amin, A.J. Berkley, P. Bunyk, V.Choi, R. Harris, M.W. Johnson, T.M. Lanting, S. Lloyd and G. Rose, Phys. Rev. B , 012507 (2009). Yu. Makhlin and A. Shnirman, Phys, Rev. Lett.92