Non-Equilibrium Quantum Dissipation

Non-Equilibrium Quantum Dissipation

Dvira Segal111Corresponding author. Current Address: Chemical Physics Theory Group, Department of Chemistry, University of Toronto, Ontario M5S 3H6, Canada. Email: and David R. Reichman Department of Chemistry, Columbia University, 3000 Broadway, New York, NY 10027    Andrew J. Millis Department of Physics, Columbia University, 538 W 120th St., New York, NY 10027.

Dissipative processes in non-equilibrium many-body systems are fundamentally different than their equilibrium counterparts. Such processes are of great importance for the understanding of relaxation in single molecule devices. As a detailed case study, we investigate here a generic spin-fermion model, where a two-level system couples to two metallic leads with different chemical potentials. We present results for the spin relaxation rate in the nonadiabatic limit for an arbitrary coupling to the leads, using both analytical and exact numerical methods. The non-equilibrium dynamics is reflected by an exponential relaxation at long times and via complex phase shifts, leading in some cases to an ”anti-orthogonality” effect. In the limit of strong system-lead coupling at zero temperature we demonstrate the onset of a Marcus-like Gaussian decay with voltage difference activation. This is analogous to the equilibrium spin-boson model, where at strong coupling and high temperatures the spin excitation rate manifests temperature activated Gaussian behavior. We find that there is no simple linear relationship between the role of the temperature in the bosonic system and a voltage drop in a non-equilibrium electronic case. The two models also differ by the orthogonality-catastrophe factor existing in a fermionic system, which modifies the resulting lineshapes. Implications for current characteristics are discussed. We demonstrate the violation of pair-wise Coulomb gas behavior for strong coupling to the leads. The results presented in this paper form the basis of an exact, non-perturbative description of steady-state quantum dissipative systems.

03.65.Yz, 05.60.Gg, 72.10.Fk, 73.63.-b

I Introduction

Over the past several decades tremendous effort has been put forth to understand the dynamics of a small quantum entity coupled to a thermal bath WeissBook (). Important problems that can be distilled to this form include the interaction between localized magnetic impurities and itinerant electrons (the Kondo problem) Kondo (); Kondobook (), electron transfer in aqueous environments Ulstrup (), and proton tunneling in biomolecules Bell (); Proton (). The study of such quantum dissipative systems cuts across traditional disciplines and impacts fields from biology to quantum information theory WeissBook ().

Our detailed understanding of quantum dissipative systems is essentially confined to problems which involve a single thermal reservoir Leggett (). In this case, traditional measures of dynamical interest are equilibrium correlation functions or simple measures of the decay of one-time quantities when the initial condition is not one of thermal equilibrium for the global system. While important for the understanding of various experimental situations, this latter form of non-equilibrium behavior is well understood and generically takes the form of an asymptotic exponential decay to the thermal equilibrium state or the ground state (at zero temperature) WeissBook (); Leggett (); Saleur ().

A less well understood type of non-equilibrium behavior may manifest when a small quantum system is coupled to more than one reservoir Aditiphonon (); Aditisemi (); Aditispin (); Aditiphase (); Aditicoul (); Zawadowski (); Paaske03 (); Paaske04 (); Paaske (); Kehrein (); AndreiPRL (); Andrei (); Kondononeq (); KondoV (). Here, the generic situation is one of a non-equilibrium steady state, regardless of the initial preparation. Given the fact that this multi-bath scenario is standard for prospective single-molecule devices MolEl (); NitzanRatner () as well as more general problems, it is imperative to understand the fundamental relaxation motifs that emerge in such nontrivial non-equilibrium cases. Recent work raises the question of whether standard tools borrowed from typical equilibrium quantum dissipative systems are useful in the steady state non-equilibrium case. For example, the simple equivalence between bosonic and fermionic baths (as obtained via bosonization Schotte (); Giamarchi (); Hakim ()) is lost in the multi-bath case, while mean-field approaches are fraught with danger due to the fact that a voltage bias may assist tunneling even at zero-temperature, rendering the meaning and stability of Hartree-Fock minima unclear Aditisemi (); Alexandrov (); Komnik (); Nitzan-nano (); Galperin ().

While several recent papers have taken up the task of describing the steady state, non-equilibrium dynamics in different model problems, our goal here is a first step towards a detailed and systematic understanding of dissipative relaxation in the simplest model problems resulting from coupling a small quantum system to several baths, namely generalized spin-boson models Aditiphonon (); Aditisemi (); Aditispin (). It should be noted that the term “spin-boson” is a misnomer; the interesting and relevant case is that of fermionic reservoirs, which dramatically differ from the case of bosonic reservoirs when the system interacts with more than one bath with different chemical potentials. On the other hand, as in the standard spin-boson model, it is the physics of the x-ray edge singularity Nozieres (); Ohtaka (); LevitovM (); Baranger () that forms the fundamental building block of the description of dynamic observables. Here, it is the recently studied non-equilibrium x-ray edge problem Ng (); Combescot (); Braunecker0 (); Braunecker1 (); Muzy2 (); Braunecker2 (); Levitov () that lies at the core of the relaxation behavior of standard correlation functions. The more complex physics of the non-equilibrium edge behavior allows for a richer range of dynamical behavior than in the well-studied equilibrium case.

In this paper we will confine our discussion to calculations that are perturbative in the bare tunneling matrix element of the system, but allow for arbitrarily strong coupling to the leads. We will employ both analytical and numerical techniques to describe the dynamics. The numerical approach involves a computational solution of the non-equilibrium x-ray edge problem that is numerically exact on all relevant time scales. This will allow us to describe the full cross-over behavior from the regime where equilibrium effects dominate, to that where the full non-equilibrium behavior (such as bias-induced dephasing and complex phase shifts) is manifested. This is crucial, since the full frequency dependence of relaxation rates and generalized fluctuation-dissipation ratios depend on the entire time history of the dynamics Aditispin ().

We will demonstrate that interesting behavior occurs in specific parameter regimes that lead to anti-orthogonality effects and bias-induced tunneling. In particular, the bias-induced tunneling regime at zero-temperature may display a very broad Gaussian decay of the polarization at strong system-leads coupling. In this regime, the relaxation behavior shows interesting similarities to the usual high-temperature Marcus (or semiclassical polaron) behavior Holstein (); Mahan (); Marcus (), with potential bias playing the role of temperature, although crucial differences exist that make these analogies imprecise. Lastly, we investigate the crucial question of the accuracy of the pair-wise Coulomb gas decomposition for non-equilibrium steady state systems. We note that the methods discussed in this work form the basis of a numerically exact path-integral description of quantum dissipation in such non-equilibrium problems Next ().

This paper is organized as follows. In section II we describe our model system (the out of equilibrium spin-fermion model). Section III presents an overview of the analytical results for the non-equilibrium dynamics, along with the relation to the non-equilibrium x-ray edge problem, while section IV presents numerical results. In section V we present the implications for the tunneling rate. Our results imply a breakdown of the Coulomb gas picture at intermediate times, described in section VI. In section VII we conclude.

Ii Model

Our model system consists of a biased two state system (spin) coupled to two electronic reservoirs held at different chemical potentials. In what follows we assume that the temperature is zero, and investigate the possibility of voltage activated excitation between the spin states. The extension to non-zero temperature is straightforward, both analytically and numerically. The total Hamiltonian is the sum of three terms:


The spin system consists of a two level system (TLS) (creation operators ) with a bare tunneling amplitude and a level splitting . The reservoir term includes two non-interacting metallic leads , where a non-equilibrium state occurs when the leads have different chemical potentials . The system-bath interaction couples the spin with scattering processes inside the leads (diagonal coupling), and in between each lead (nondiagonal coupling), and we choose conventions such that only one of the spin levels couples to the leads,


Here is the number operator, with as the identity operator Hamilton (). The operator () creates (annihilates) an electron with momentum in the -th lead. In this paper we focus on the model presented in Ref. Ng (), where the momentum dependence of the scattering potential is neglected. System-bath scattering potentials are then given by , where are the Fermi sea indices. Our main conclusions, however, are valid for more general cases.

We assume that the reservoirs have the same density of states , typically modeled using a Lorentzian function


where is a bandwidth parameter. We typically work in the limit of wide bands, , therefore to a good approximation .

Note that we ignore the spin degree of freedom of the reservoir electrons in our discussion. In what follows we refer to the energy difference as a magnetic field, in order to distinguish it from the voltage bias . We also define two auxiliary Hamiltonians that will be useful below


Explicitly, includes the electronic reservoirs and system-bath interaction, given that the subsystem is in the state. The model (2) contains much of the physics of the Kondo model Kondo (), while lacking direct coupling of the reservoir degrees of freedom to spin-flip processes. It also contains the spin-resonant-level model of Ref. Aditispin () with a particular choice of system-bath couplings. We discuss the spin-resonant-level model in more detail in Appendix A.

Crucial parameters of the model are the and scattering phase shifts. In equilibrium the phase shifts are given by Nozieres (); Ng ()


where and () are dimensionless system-bath coupling strengths


and is the density of states at the Fermi energy of the reservoirs. Out-of-equilibrium, , the phase shifts are complex numbers given by Ng (),


Since the reservoirs density of states weakly varies around the Fermi energy, , the phase shifts are approximately energy independent, and are all calculated at the Fermi energy Ng (). For simplicity, throughout the paper we typically consider the case of , and take to be real. We note however that our main results, in particular the appearance of Marcus-type behavior in the non-equilibrium regime at strong coupling, can be rederived using other variants of this model system with no limitations on the strength of the diagonal interactions, as well as for the spin-resonant-level model of Ref. Aditispin (), see Appendix A.

Under these simplifications, the non-equilibrium phase shifts are given by


For the inverse tangent in Eq. (8) has a branch cut, conventionally placed at and . For this special case,


The weak potential limit therefore corresponds to and , so . However, as , diverges, whereas the equilibrium phases are finite.

An important quantity that will be useful below is the sum of the phase shifts squared, . While for the general model Eq. (7) yields complex numbers, when the system is symmetric (), the phase shifts are complex conjugates and is real, although possibly negative. We discuss the implications of this result in section IV.C .

For the sake of completeness and comparison the equilibrium spin-boson model is discussed in Appendix B. In the nonadiabatic limit, this model yields the classical Marcus rate at high temperatures when the system-bath interaction is strong. We analyze the analogous behavior in the non-equilibrium spin-fermion model (2) in section V.

Iii nonadiabatic dynamics

iii.1 Overview

We are interested in the reduced density matrix in the space of occupancy. This is defined in terms of time evolution from an initial condition at time ,


If the parameter in Eq. (2) vanishes, the problem is just electrons in a time-independent potential, and a closed form analytical solution exists. If , may be expressed as an expansion in . Evaluation of any term in the expansion entails solving a problem of electrons in a time-dependent field. In equilibrium, an essentially exact closed-form solution exists, and the main problem is to re-sum the series in Chen (); Chang-Chak (). For non-equilibrium problems an analytical expression is not known. In this paper we present a detailed numerical evaluation of some low order terms in the expansion for . The essential features are revealed by the Golden Rule decay rate, obtained by assuming that (i) at time , and (ii) is the density matrix corresponding to the ground state of with , and (iii) that an expansion to suffices. This level of description is equivalent to the ”non-interaction blip approximation” (NIBA) in the standard spin-boson model WeissBook (); Leggett (); Aslangul () and yields the (nonadiabatic) Fermi Golden Rule for the forward () and backward () transition rates between the spin levels as Mahan (); Lax (); Kubo (); Golosov ()


where denotes time ordering, , and is the ground state energy of the two uncoupled reservoirs. The Hamiltonians are defined in Eq. (4), refers to the real part of the integral, and the trace is performed over the electronic degrees of freedom. For convenience, the term including the energy bias is taken outside of the trace.

The object of our calculation is therefore the correlation function , which should be evaluated for non-equilibrium conditions covering time scales from up to . , an energy of the order of the Fermi seas bandwidth, and the potential drop , specify two inverse time scales in the problem, where we typically work in the limit of . Note that, unlike the equilibrium case, there is no exact analytical approach to calculate valid for all time scales. Approximate analytical approaches and exact numerics may be performed, as discussed below.

iii.2 Short and long-time asymptotics; Non-equilibrium x-ray edge problem

The correlation function [Eq. (11)] is a crucial element in the theory of the x-ray edge problem, an effect originating from the many body response of a Fermi system to the fast switching of a scattering potential, e.g. the creation of a core hole Nozieres (); Ohtaka (). The x-ray edge Hamiltonian is a simplified version of the spin-fermion model, Eq. (2), with a static subsystem that is either empty or populated,


Here , are creation and destruction operators of the core electron, and () creates (destroys) an electron in the -th lead with momentum . The single band x-ray singularity problem was originally solved exactly in the asymptotic limit by Nozieres and De Dominicis (ND) Nozieres ().

In the last ten years there has been a growing interest in understanding the x-ray edge effect in the mesoscopic regime LevitovM (); Baranger () and for non-equilibrium systems Ng (); Combescot (); Braunecker0 (); Braunecker1 (); Muzy2 (); Braunecker2 (); Levitov (), where the core hole couples to more than one Fermi sea at different chemical potentials. Standard equilibrium techniques, e.g. bosonization Schotte (); Giamarchi (); Hakim (), cannot be simply generalized to handle these non-equilibrium systems (see Appendix C). The first to address the non-equilibrium problem was Ng, who generalized the Nozieres-De Dominicis solution to include more than one Fermi sea with different chemical potentials Ng (). Ng demonstrated that the edge singularity could be described by generalized phase shifts which are real for equilibrium systems and complex when the system is driven out-of-equilibrium. Physically, complex phase shifts reflect the finite lifetime of a non-equilibrium system. More recently, Muzykantskii et al. Braunecker1 (); Muzy2 () formally solved the out-of-equilibrium problem using the Riemann-Hilbert approach. The result, given in terms of the scattering matrix, was later generalized to include finite temperature effects Braunecker2 (). An exact formal determinant solution was presented in Ref. Levitov () for the study of tunneling in a non-equilibrium electron gas.

A formal solution for is obtained from the linked cluster theorem (valid also for non-equilibrium problems) Nozieres (); Ng (),


with the matrix Green function in the space of the leads for Eq. (1), but with . For this model solves the Dyson equation


with and the unperturbed Green’s functions


is the reservoir density of states, taken to be the same for the and leads.

In equilibrium, ND showed Nozieres () that this equation can be solved exactly and the coupling constant integral performed, leading to


with an energy of the order of the Fermi sea bandwidth, , and defined in Eq. (5). For the model studied explicitly


Eq. (16) also holds for the non-equilibrium problem at times .

At long times, , the equation was solved by Ng Ng (); see also Braunecker1 (). The coupling constant integral may similarly be performed, leading to


with , and given by Eq. (8). Here () refers to the imaginary part of the phase shift. For the model studied numerically () we have




iii.3 Intermediate time

Although the long time and short time behavior is known essentially exactly, a transparent non-perturbative analytical expression for that encompasses all time scales and coupling strengths has not been developed. Indeed, in this work we argue that at strong coupling () a different functional form dominates at intermediate times where a prominent Gaussian decay emerges, . We first offer a perturbative calculation which suggests this result, and then present exact numerical simulations which prove this behavior. The dominance of the Gaussian behavior at intermediate times translates into a Marcus-type rate in frequency domain, with bias voltage activation (see discussion in section V), instead of temperature activation, as in the classical Marcus rate (see Appendix B).

The correlation function can be evaluated using the cumulant expansion Mahan (). Note that unlike the bosonic case, all cumulants contribute,

where denotes time ordering, , and denotes a cumulant average. The first cumulant yields an energy shift, while the second term is given explicitly by (, )


For details see Appendix D. The sine and cosine integrals are defined as , , and is the Euler-Mascheroni constant.

This expression reproduces the weak coupling limits of the analytical results Eqs. (16) and (18) at short and long times respectively, and provides an interpolation between the two times. In particular, the second line describes how the long-time dissipation term is ”turned on” as increases from a small value to values much greater than unity. The first and last terms describe how the equilibrium orthogonality is turned off as increases: is a function which interpolates between for , and a constant at . Note that in this model the leading logarithmic term at long times is , consistent with the cancellation of logarithmic terms in the long time limit of Eq. (22), i.e. with the absence of a term which ”turns on” the non-equilibrium power law. This cancellation does not necessarily occur at order in other models, e.g. the spin-resonant-level model, see Appendix A.

We can clearly distinguish between three regimes in Eq. (22):


While the first (equilibrium) limit and the third regime are well established in the literature Ng (); Braunecker1 (); Aditispin (); Muzy2 (), the intermediate domain, leading to an interesting new dynamic has not been discussed. In the strong coupling limit the Gaussian behavior may have a dominant effect on the relaxation, as discussed below. We would like therefore to phenomenologically extend the second cumulant expression, Eq. (22), to larger phase shifts (strong coupling).

Perturbative expressions analogous to Eq. (22) motivated Mitra and Millis Aditisemi () to propose an interpolation function constructed by replacing the factors of in the expression above by the exact phase shifts. For the model considered here their procedure leads to


Note that our approximation for the scattering potentials, , implies that there is no Fumi energy shift. At the short time limit, , the factor dies out, leading to the correct equilibrium behavior (16). In contrast, at long times the cosine integral diminishes, which implies that the dynamics is ruled by an exponential decay with a rate constant , [Eq. (19)], modified by a power law term , Eq. (20).

Our numerical results, to be presented below, show that at weak to moderate coupling, , the correlation function and the resulting transition rates are well described by expression (24). However, Eq. (24) is found to be a poor approximation at strong coupling. Instead, at intermediate times we return to Eq. (22) and replace the weak coupling phase shift by the equilibrium strong coupling phase shift, . The physical picture is that on these time scales the phase shifts are essentially still the equilibrium ones. Only at longer times the non-equilibrium dynamics is reflected in the complex phase shifts (8). This conjecture yields


where the equilibrium function is the same as in the zero bias case Nozieres (),


while the non-equilibrium term provides a quadratic time decay


Notice that the prefactor depends only on the scattering potential . The last element in Eq. (25) is the energy shift . We assume that it is given by the equilibrium limit of the Fumi’s theorem,


For the energy shift is zero.

Similarly to expressions (25)-(28), Eq. (24) gives the first correction to the equilibrium result which is proportional to , but in contrast to these equations, the coefficient involves the non-equilibrium exponent, and in fact does not provide a wide regime of behavior. Our numerical simulations, presented below, support expressions (25)-(27) in the broad window . We have not been successful in constructing a general analytical expression, valid on all time scales and coupling strengths. It is possible that consideration of the fourth cumulant may yield some insight here.

Iv Numerics

iv.1 Methods

The fermionic correlation function can be directly calculated by expressing the zero temperature many body average as a determinant of the single particle correlation functions Mahan (); Temp ()


Here , where are the single particle Hamiltonians for the individual conduction electrons. are the single particle eigenstates of , and the determinant is evaluated over the occupied states. is the Fermi energy of the -th reservoir. In our numerical calculations we have used a Lorentzian density of states, with tails that are long enough to eliminate artificial reflections from the boundaries. The Lorentzian function is centered around the equilibrium Fermi energy with a full width at half maximum , Eq. (3). This quantity sets energy and time scales in our simulations. We have typically used for the two reservoirs, and . We also take the diagonal coupling to be zero () in all of our simulations, unless otherwise stated. For these parameters, we have found that for short-time evolution (), even for strong coupling, it is satisfactory to model the fermionic reservoirs using 400 states per bath, where bias is applied by depopulating one of the reservoirs with respect to the other.

We can also employ the renormalization group (RG) method, originally developed by Wilson for the calculation of the thermodynamic properties of the Kondo problem Wilson (), for the numerical solution of the non-equilibrium x-ray edge problem. In equilibrium, Oliveira et al. Oliveira1 (); Oliveira2 () have used the RG technique to calculate the x-ray absorption spectrum. This study can be generalized to include two reservoirs with different chemical potentials by following a three-step procedure: (i) define the conduction bands on a logarithmic scale. (ii) convert the (isolated) reservoir Hamiltonians into semi-infinite tight binding chains, as is done in Refs. Oliveira1 (); Oliveira2 (). In this representation the impurity couples the chains’ first levels. (iii) build the Hamiltonians in the new basis, first including the occupied levels of the and reservoirs, then adding the empty levels. The determinant (29) is performed over occupied levels only.

In equilibrium, the RG technique is highly advantageous over constant/Lorentzian discretization methods, as it converges rapidly to the continuum limit even for gross discretization. For small voltage differences () this method nicely reveals the crossover of from equilibrium to non-equilibrium behavior with increasing bias. In contrast, for large bias the Lorentzian discretization is more convenient, since energies far from the Fermi energy are not well represented within the RG technique. We present a numerical example in Fig. 1, demonstrating the strength of the RG approach over standard linear discretization for systems in equilibrium. The RG technique provides stable dynamics for long times (full line), where constant discretization fails (dotted line), yielding an artificial rise of the correlation function due to discretization errors. The theoretical value of for , nicely agrees with the numerical slope of . Deviations are due to the sharp energy cutoff used at =1, with the conduction band energies extending from to . We also present the results of an RG calculation with a very small voltage drop (dashed line), where linear discretization would require a very fine grid.

In this work we typically focus on systems far from equilibrium, . Since the RG method samples the Fermi sea states predominantly near the Fermi energy, while high energy states are under-represented, we find the Lorentzian discretization to be more convenient.

iv.2 Results:

Representative results are displayed in figure 2. The main plot presents the logarithm of the correlation function at strong coupling for an applied voltage . Three different regimes are clearly identified: a power law decay at short times , see lower left inset (a), an exponential decay at long times , and remarkably, an intermediate regime of approximately Gaussian behavior [upper right inset (b)]. The short and long time behaviors are consistent with the theoretical results. The intermediate time quasi-Gaussian regime is a new finding with important consequences. We analyze the short time dynamics , enlarged in Fig. 2(a), by fitting the data to the analytic expression . This provides an effective bandwidth and a decay constant consistent, within numerical errors, with the theoretically expected . We can also fit the intermediate time behavior, shown in Fig. 2(b), by a Gaussian function which yields the prefactor =0.03.

Figure 3 presents a more detailed examination of the Gaussian behavior, showing that at both short and intermediate times, the data can be well described by the approximate function


with the theoretically predicted short time (equilibrium) exponent . The inset proves that the data follows the same linear trend when plotted as a function of , with a slope of . This value nicely agrees with the constant predicted by Eq. (27), =0.029 ().

Fig. 4 provides more insight by deconstructing the observed time decay of into the equilibrium power law and non-equilibrium Gaussian components. Another important observation deduced from Figs. 2, 3 and 4 is that the correlation function decays to its initial value by the time the exponential decay begins to dominate. This implies that the Gaussian behavior governs the rate constant at strong enough coupling, leading to a voltage activated regime analogous to the high-temperature semiclassical polaron transport regime. We call this ”fermionic Marcus” behavior.

Fig. 5 presents the evolution of the correlation function as coupling strength is varied from weak to strong. All other parameters are the same as in Fig. 2. For all coupling strengths the short time logarithmic and the approximate long time exponential behavior are observed. However, as the coupling strength is increased, increasingly wide intermediate regime is observed. We have verified, by an analysis similar to that shown in the lower left inset of Fig. 2, that the short time behavior is always a power law with the theoretically predicted exponent . Also note that while at weak coupling () the correlation function weakly decays before the turnover to an exponential decay takes place, for very strong coupling, , the dynamics is critically controlled by the Gaussian form, as the correlation function has decayed to zero before the exponential decay takes place. This implies that the resulting decay rate [Eq. (11)] essentially shows different characteristics in these two regimes.

Figure 1: Decay of as computed by both the RG approach and standard linear discretization for . Constant discretization with 200 states per band, (dotted), showing an artificial rise of the correlation function due to discretization errors; logarithmic discretization (RG) with 100 states per band, , (full); logarithmic discretization (RG) with 100 states per band, , (dashed). is a logarithmic scale parameter, where for each conduction energy there are states with energies , . Inset: The finite bias case exhibits an exponential decay at long times. The slope agrees with the theoretical value, Eq. (19). The Fermi sea bandwidth is in all plots with .

Figure 2: The correlation function for =0.24, =0.95, manifesting a power law decay at short times (a) (notice the log-log scale), a Gaussian decay at intermediate times (b), and an exponential decay at long times .

Figure 3: Evidence for the Gaussian decay at intermediate times 1-10 and strong coupling , (full), (dotted) and (dashed-dotted). The inset, which includes the three lines one on top of the other, reveals that , =0.13, with .

Figure 4: Time dependence of the correlation function for , . The fitting function (dotted) with , and compared to exact numerical solution (full). The Gaussian (dashed-dotted) and the power law (dashed) parts of are also displayed separately.

Figure 5: The correlation function for different coupling strengths , manifesting the increasing dominance of intermediate time Gaussian behavior at strong coupling. in all plots.

We now systematically explore the Gaussian decay at intermediate times, and the exponential decay rate at long times, and compare the numerical coefficients and with the theoretical values, Eqs. (19) and (27), respectively. This is done by calculating the correlation function for coupling strengths (see Fig. 5), then extracting both the quadratic intermediate slope and the long time exponential slope . Fig. 6 presents these coefficients showing excellent agreement with the values predicted from the phenomenological ansatz, Eqs. (25)-(27).

Next, in Fig. 7 we examine the crossover to the analytic long time behavior, Eq. (24). We compare the numerical correlation function with two functions: the approximate fitting function defined above, and the long time perturbation theory result of Ref. Aditispin (), given by exponentiating Eq. (24). We refer to this second function as . We see that describes the data well at long times, but that as the coupling strength is increased, the range over which the Gaussian description applies increases. This feature can be qualitatively described by the approximate crossover function

which captures the crossover from a Gaussian dynamics to an exponential decay. An increase of leads to a strong enhancement of , while reaches saturation, resulting in a counterintuitive lengthening of the range of the intermediate Gaussian dynamics with increased .

In summary, we have shown that the crossover between equilibrium () and non-equilibrium () behavior is described by a regime of Gaussian relaxation negligible for weak coupling, but for strong coupling extending over the wide range , with parameters determined by the equilibrium exponents. In Section V we examine the consequences for spin relaxation.

Figure 6: Testing the validity of Eqs. (19) and (27) for describing the long time and intermediate time behavior, respectively. Top: The relaxation rate . Numerical results, calculated from the slope of vs. at long times (squares); analytical results using Eq. (19) (dashed line). Bottom: The coefficient . The numerical slope of vs. at intermediate times (squares); analytical results using Eq. (27) (dashed line). These data were computed with =0.24, and was extracted from the short time dynamics for each value of .

Figure 7: The turnover from a Gaussian decay to an exponential relaxation. Top: Comparison between the numerical correlation function (full) and the fitting functions (dotted) and (dashed) defined in the text, , =0.24. Inset: Exposing the turnover by taking away the short time and intermediate time terms with =0.002, . Bottom: Same with , leading to , . While explicitly includes the Gaussian decay, correct to , the function captures the correct slope at longer times.

iv.3 Orthogonality and anti-orthogonality

We focus next on the power law contribution to Eq. (18). Unlike the standard equilibrium case, where the system always experiences dephasing, Anderson (), in our model the power law term in Eq. (18) acquires a positive exponent , enhancing the correlation function, see Eq. (20). We refer to this situation as an ”anti-orthogonality” effect. For a general system-bath coupling model, Eq. (24) reveals that


with complex, non-equilibrium, phase shifts given by Eq. (7) Ng (). It is clear that in the special limit of zero diagonal interactions (), the phase shifts are purely imaginary and for all values of , leading to . In contrast, for large diagonal coupling we typically find that , which is the standard orthogonality behavior. The anti-orthogonality effect is therefore a footprint of a non-equilibrium situation.

We next turn to a numerically exact exploration of the anti-orthogonality effect. Fig. 8 shows the correlation function at long times when expression (24) holds. We numerically extract the long time slope , and recover the weak power law dependence by multiplying the correlation function by the inverse of the exponential decay. The standard orthogonality effect is presented in panels (a)-(b) for and , for which . When and the anti-orthogonality effect clearly manifests itself with (c)-(d). Interestingly, the correlation function at long times shows a complicated behavior, more complex than that predicted in Eq. (24), as evidenced by the mild deviations from strict power law behavior displayed in Fig. 8(d).

Fig. 9 presents an ”orthogonality-anti-orthogonality” map as a function of the diagonal () and nondiagonal () couplings using the general expressions of Eq. (7) with . We find that for large , , manifesting the standard orthogonality effect. For large nondiagonal interactions typically anti-orthogonality may be observed.

Figure 8: Orthogonality and anti-orthogonality effects in the non-equilibrium system . (a)-(b) =0.3, =0.5, manifesting the standard orthogonality effect (). (c)-(d) =0.8, =0, revealing the anti-orthogonality effect ().

Figure 9: Map of orthogonality [] and anti-orthogonality [] behavior using expression (7) with .

In addition to the long time exponential decay and anti-orthogonality behavior, non-equilibrium dynamics may be reflected in the appearance of complex power law exponents Ng (). When the spin impurity is symmetrically coupled to the two leads (), the phase shifts are complex conjugates, see Eq. (7), and is always real. This situation was discussed above. In contrast, asymmetric systems may acquire a complex coefficient with , a direct outcome of a non-equilibrium situation.

The imaginary contribution to is resolved in Fig. 10. Since at weak coupling the imaginary term is very small, we investigate a strong coupling system with . Motivated by Eq. (24), we assume the generic form . We numerically extract the phase factor , then plot the function for different diagonal coupling strengths. As expected, in symmetric situations, . In contrast, asymmetric systems () reveal an additional decaying contribution which is expected to oscillate at longer times. We did not succeed in fitting to the stretched-oscillatory function , indicating that at strong coupling the dynamics is more involved. Finally, we note that, consistent with our observations above, equilibrium effects dominate up to . Only at longer times begins to deviate from unity due to the emerging influence of the imaginary term .

Though the imaginary term can strongly affect the correlation function , (Fig. 10), its practical contribution to the Golden Rule rate is small. We find that for weak to intermediate coupling, , leading to . On the other hand, for strong coupling, manifests itself only at long times , when the correlation function has essentially decayed to zero.

Figure 10: Resolving the complex part of the power law exponent . , =0 (full); , =0 (dashed); , =0 (dashed-dotted); , =0.2 (dotted); and in all plots.

V Relaxation

v.1 Qualitative discussion

In this section we present calculations of the nonadiabatic relaxation rates defined in Eq. (11). The given physical model corresponds to two states separated by an energy which depends on the bare level splitting and on renormalizations arising from the coupling to the leads. corresponds to the up-scattering rate describing transitions from the lower level to the upper level, while corresponds to down-scattering. In equilibrium at , =0, i.e. there is no up-scattering. At temperature , the detailed balance relation of equilibrium thermodynamics implies . In this section we examine the rates in the non-equilibrium situation. We show that the Gaussian form of the correlation function which occurs at strong coupling has important consequences for the physics.

Before discussing our results in detail, we establish the relevant energy scales. The general expression, Eq. (11), may be written (neglecting overall factors) as


Here is an effective exponent which changes from the equilibrium power , Eq. (17), to the non-equilibrium value [defined as in Eq. (20)], as changes from less than unity to much great then unity. is the physical energy level difference, given by the sum of [Eq. (2)] and the level shift arising from the system-bath coupling, and is an energy scale of the order of the Fermi sea bandwidth.

The naive assumption Aditisemi () is that the only important energy scale is the relaxation rate given by the current flow across the system, . In fact the numerical and analytical results presented in the previous sections indicate that the situation is more subtle. At short times, whereas at long times . The interplay of , which is proportional to coupling strength at weak coupling but saturates at strong coupling [Eq. (27)], and which is proportional to at weak coupling but diverges at strong coupling [Eq. (19)], gives a richer behavior.

Appendix E gives details of an asymptotic analysis of Eq. (33). This analysis reveals that to discuss the relaxation rate one should distinguish strong and weak coupling. In the weak coupling limit, there are two relevant scales, and (the latter multiplied by various factors which are in practice fairly close to unity). For we get the equilibrium down-scattering rate; for we find a nontrivial approximately Lorentzian behavior, and for we reproduce the behavior found by Mitra et al. Aditisemi (). Specifically,


Here is the complete Gamma function. Note that the formulae match at because in weak coupling leading to . In the strong coupling limit, two frequency scales turn out to be important: and . We find


The Gaussian behavior found at intermediate frequency scales is a consequence of the wide regime of behavior found in the time evolution function, and may be roughly understood as the Fourier transform of although as the results of Appendix E show, this argument must be treated with some care.

We call Eqs. (37)-(37) the ”fermionic Marcus rate,” the analogue of the classical Marcus result for spin-boson systems (B12), which holds in non-equilibrium situations at strong coupling. This expression indicates that the voltage activates the absorption rate, similarly to the role of temperature in the bosonic case. The result differs from the bosonic solution (Appendix B) in some important aspects: (i) In the fermionic case the Gaussian decay is modified by a weak power law term. (ii) For bosonic systems the activation factor depends on the temperature as , while for fermionic systems we get a voltage squared activation, . Therefore, there is no simple linear mapping between temperature and voltage drop in the strong coupling regime. We note however that the classical Marcus rate is applicable in the high temperature limit (see Appendix B), while we typically assume here that . Therefore, in our system the energy window for reorganization processes is the bias voltage, rather than the full bandwidth . Thus, we may interpret the factor in the denominator of the Gaussian decay (37) as a reorganization energy of the non-equilibrium fermionic system, , multiplied by the driving force . In contrast, in the equilibrium spin-boson model, reorganization energies are of order of the cutoff frequency, , and the driving force for absorption processes is temperature . Qualitatively, both models then recast to the familiar Marcus-like form, Marcus-comm (). Further, both the fermionic and the standard Marcus behaviors share similar qualitative features such as the existence of an inverted regime, as discussed in the next section.

Figure 11: Golden Rule relaxation rate calculated for the weak coupling limit (, ), (dashed line); (solid line); (dotted line). The inset presents an expanded view of low frequency regime.

Figure 12: Fermi Golden Rule rate, Eq. (33), as a function of frequency. and (). Inset: Voltage activated excitation rate ().

Figure 13: Spin polarization as a function of potential bias for different energy biases (full), (dashed), (dashed-dotted) for , (). Inset: Mean field calculation of the current for the same frequencies as in the main plot.

v.2 Numerics: rates, population and current

We numerically evaluate the integral (33) using the coefficients , , and as determined by the coupling strengths, Eqs. (19), (20) (26) and (27) respectively. For convenience, we disregard the multiplicative factor .

The main panel of Fig. 11 shows on a semi-logarithmic scale the relaxation rate computed numerically for the relatively weak coupling (, , non-equilibrium exponent ) and two choices of chemical potential, and . Also shown as the dashed line is the equilibrium result. The inset shows an expanded view of the small frequency regime, demonstrating the Lorentzian behavior. We clearly observe the three regimes as discussed in Eq. (34): For small frequencies (large bias voltage) the spin levels are approximately degenerate, and the rates are symmetric around (inset). In the opposite limit, the absorption rate is practically zero, while the emission rate approaches the equilibrium limit. In between, a voltage activated excitation behavior is revealed.

We analyze next the strong coupling limit. Fig. 12 shows that the excitation process is activated by a finite potential difference as prescribed by Eq. (37). More quantitatively, the inset verifies that the relationship holds. Similar to classical bosonic Marcus rate Marcus (), an inverted regime appears for the fermionic system. However, in the present case the rate in the inverted regime decays weakly as a power law rather than as a Gaussian. At large frequencies, , equilibrium behavior is observed where approached zero, and becomes insensitive to voltage. We have also calculated the Golden Rule rate using the numerical correlation function (depicted e.g. in Figs. 2 and 3), instead of the approximate analytical function in (33), and find that the results agree perfectly.

We now turn to a study of the spin polarization. In the incoherent tunneling regime, for small tunneling parameter , the populations of the two levels obey a Markovian balance equation


with the absorption and emission rates given by Eq. (33). The polarization , shown in Fig. 13, manifests a transition from a fully polarized system to an unpolarized system as is increased. Typically, we find that the crossover takes place when the energy bias becomes comparable to the bias voltage. While at high frequencies, , , leading to full polarization, at very large bias the emission and absorption rates are comparable, resulting in equal population of the two levels and zero polarization. The Gaussian activation term in Eq. (37) is therefore reflected in the enhancement of polarization with bias voltage.

It is also interesting to note that the electron current through the system, calculated at the level of mean field theory, , is strongly suppressed for weak bias, , see inset of Fig. 13. In contrast, for very large bias, , and the current increases linearly with . Therefore, it is the intermediate regime of that manifests prominent nonlinear current-voltage characteristics, emerging due to the interplay between the Gaussian relaxation and the power-law dynamics. We can compare our results to the weak coupling Bloch-type rate equations of Gurvitz et al. Gurvitz () which yield at long enough times, independent of voltage drop and energy bias. In contrast, Fig. 13 reveals a rich dynamics in the strong coupling regime with a prominent dependence on system energetics and the non-equilibrium conditions.

Vi Beyond : Coulomb Gas behavior

Figure 14: Schematic representation of spin-flip evens on the Keldysh contour. Plotted are examples for particular two, four and six spin-flip processes, respectively.

In this section we discuss a crucial ingredient of the physics of our model that allows for a description beyond the Golden Rule [] level. A formally exact solution for the impurity spin problem (2) can be written by a power series in the tunneling matrix element Chang (). Here, we restrict ourselves to an exact numerical investigation of the electronic correlation functions that appear in this power series to see if the usual ”Coulomb gas” behavior is observed even when the system is out-of-equilibrium. In particular, the reduced density matrix of the spin impurity is given by


with forward and backward time evolution branches. Here is the total density matrix, and the trace is performed over the reservoir electronic states. We decompose the propagators, including all spin-flip events along the time ordered contour, and obtain, e.g. for the spin up population Leggett (); Chang (),