Dissipative quantum mechanics beyond Bloch-Redfield: A consistent weak-coupling expansion of the ohmic spin boson model at arbitrary bias
Carsten J. Lindner Herbert Schoeller





Institut für Theorie der Statistischen Physik, RWTH Aachen, 
52056 Aachen, Germany and JARA - Fundamentals of Future Information Technology

[Email: ]schoeller@physik.rwth-aachen.de
Institut für Theorie der Statistischen Physik, RWTH Aachen, 
52056 Aachen, Germany and JARA - Fundamentals of Future Information Technology


The ohmic spin boson model at finite bias ϵ and tunneling Δ is an 
important model to study relaxation and decoherence properties of qubits coupled weakly 
to a dissipative bosonic environment. 
Fault-tolerant quantum computation requires the investigation of small errors beyond 
Bloch-Redfield theory. Using perturbation theory and the real-time renormalization group (RG) method 
we present a consistent zero-temperature weak-coupling expansion for the
time evolution of the reduced density matrix ρ(t) out of an uncorrelated but otherwise arbitrary initial 
state. We show that standard Born approximation schemes calulating the effective Liouvillian 
of the kinetic equation up to first order in the dimensionless coupling constant α 
are not sufficient to account for various important corrections one order beyond the Bloch-Redfield solution:
(1) The resummation of all secular terms ∼(Γ t)^n is necessary to obtain the correct 
exponential decay of all terms of the time evolution with decay rate Γ or Γ/2, 
together with the correct 
preexponential functions. We show that this is only possible by a correct analytic continuation 
of the Fourier transform L(E) of the effective Liouvillian into the lower half of the
complex plane. (2) The resummation of all logarithmic terms at high and low energies leads to a 
renormalized tunneling Δ̃ and to preexponential functions 
of logarithmic and power-law form. This is achieved by solving closed differential equations for the
derivatives of L(E) set up by the real-time RG method. (3) The fact that two eigenvalues of L(E) 
are close to each other by O(Γ) requires degenerate perturbation theory for times Γ t∼ O(1),
where certain terms of the Liouvillian in O(α^2) are needed to calculate the 
stationary state and the time evolution of the non-oscillating purely decaying terms up to O(α). 
Solving the real-time RG equations perturbatively we present a renormalized perturbation theory
with analytical results for arbitrary bias covering the whole crossover regime from small times 
Ω t≪ 1 to large times Ω t≫ 1, where Ω=√(ϵ^2+Δ̃^2) denotes 
the renormalized level splitting, and compare to other results published in the literature. We
find that branch cuts starting at pole positions lead to rather slowly varying logarithmic time
dependencies for the preexponential functions, whereas pure branch cuts give rise to a strong
time dependence of the preexponential functions in the crossover regime with a leading long time tail 
∼α/(Ω t) at finite bias, besides other terms ∼α/(Ω t)^2 well-known from 
the zero-bias case. For exponentially large times, we use the real-time RG method for a 
non-perturbative resummation of all logarithmic terms ∼ (αln(Ω t))^n  and present
a numerical solution for arbitrary bias. 
We find an interesting power-law behaviour (1Ω t)^2αϵ^2Ω^2 
with a bias-dependent power-law exponent appearing in the preexponential functions of the oscillating terms.
This power-law has to be contrasted to the one at exponentially small times where the power-law exponent
crosses over to 2α. We discuss that the complexity to calculate one order beyond Bloch-Redfield 
is rather generic and applies also to other models of dissipative quantum mechanics. 




05.10.Cc, 05.30.-d, 05.30.Jp, 73.23.-b



Dissipative quantum mechanics beyond Bloch-Redfield: A consistent
  weak-coupling expansion of the ohmic spin boson model at arbitrary bias
    Herbert Schoeller
    
===========================================================================================================================================





§ INTRODUCTION



The study of the dynamics of two-state quantum systems coupled weakly to a dissipative bath is
a fundamental problem of nonequilibrium statistical mechanics that has become of increasing 
importance due to possible future technological applications in quantum information
processing. To realize scalable and fault-tolerant quantum computation very low error
thresholds are needed which requires an understanding beyond Markov approximation schemes and
lowest order perturbation theory in the coupling to the bath. As a generic model for a bosonic bath 
the spin-boson model has been proposed <cit.> and its dynamical properties have been studied
with various methods <cit.>. This model consists of two levels with level
spacing (bias) ϵ, coupled by a direct tunneling term Δ, and each level is linearly coupled  
to an ohmic bosonic bath. In the case of zero tunneling, the spin-boson model can be solved exactly 
and the stability of
surface-code error correction against realistic dissipation has been studied recently 
<cit.>. However, for finite tunneling and the most important case 
of an ohmic coupling to the bath, we will show in this work that a consistent weak-coupling expansion 
beyond the Bloch-Redfield Markov approximation is still lacking at low temperatures. We will discuss 
various subtleties to obtain a consistent perturbative expansion of the time evolution in the
dimensionless coupling constant α, requiring an essentially non-perturbative treatment in a 
certain sense, not yet accounted for completely in various previous publications on the ohmic spin boson
model. Most importantly, our analysis shows that the systematic calculation of errors to the 
Bloch-Redfield result is generically very complex for all models of dissipative quantum mechanics, 
involving many details of the underlying model, and is not specific to the ohmic spin boson model.

The concrete form of the time evolution depends crucially on the form of the density of states 
of the bath and the energy dependence of the coupling constants g_q of the local system to the 
various bath modes ω_q, conveniently taken together in the spectral density (sometimes also 
called hybridization function) J(ω)∼∑_q g_q^2δ(ω-ω_q). However, for a 
flat spectral density (on the scale of the typical energy scales of the local system) or for special 
cases as the ohmic spin boson model, where J(ω)∼αω, a rather generic discussion 
of the typical form of the time evolution is possible and has been provided in 
Ref. pletyukhov_schuricht_schoeller_10, based on the real-time renormalization group method 
(RTRG), see Rfs. schoeller_09,schoeller_14 for reviews. Starting at time t=0 from an arbitrary 
initial state ρ_0=ρ(t=0) of the local system without any initial system-bath correlations 
(i.e. the bath is assumed 
to be an infinitely large system in (grand)canonical equilibrium), the time evolution for the reduced density
matrix ρ(t) of the local system consists of a sum of terms each of them being ∼F_i(t)e^-iz_it,
with z_i=Ω_i-iΓ_i, i.e. exponentially decaying with decay rate Γ_i>0, oscillating 
with frequency Ω_i, together with a non-exponential preexponential function F_i(t), typically 
depending logarithmically or as some power-law on time. The case Γ_i=0 is exceptional and occurs only
for systems with quantum critical points, where the scaling behavior is not cut off by any decay rate. 
For the ohmic spin boson model there are three modes of a purely decaying mode z_0=-iΓ and two
oscillating modes with z_±=±Ω-iΓ/2. 
This form already suggests where the complexity of calculating the time evolution beyond the lowest
order Markovian Bloch-Redfield theory appears. Bloch-Redfield considers only the leading order term
where F_i is basically a constant of O(1), independent of the coupling α to the bath.
However, there are additional terms to each matrix element of the 2× 2-matrix ρ(t), where 
F_i(t)∼α, also containing an exponential function, usually different from the one of the
Bloch-Redfield term. Expanding this exponential in Γ∼α leads to an ill-defined 
perturbation expansion, since terms ∼α(Γ t)^n appear, which all become of O(α) already on 
time scales of the inverse decay rate (and even diverge for time going to infinity). Therefore, for a 
consistent calculation of the O(α) correction to Bloch-Redfield on time scales where the exponential
damping is still moderate, it is necessary to resum these terms in all orders of perturbation theory to get
the correct exponential behavior. We note that these so-called secular terms (sometimes also called 
van Hove singularities <cit.>) are usually only discussed when expanding   
the exponentials of the Bloch-Redfield terms in α, but similarly also appear in higher
order terms, which are more subtle. Technically, they can all be incorporated by expressing the
perturbative expansion for the effective Liouvillian L(E) in Fourier space not in terms of the
bare Liouvillian but in terms of the full Liouvillian again by taking all self-energy insertions into
account. Within the diagrammatic expansion developed in 
Rfs. schoeller_09,schoeller_14,pletyukhov_schoeller_12,kashuba_schoeller_13 it can
be seen that this is possible in all orders of perturbation theory. This allows for a convenient 
analytic continuation of L(E) into the lower half of the complex plane, from which the position
of all non-zero singularities z_i (poles and branching points) of the Fourier transform 
ρ(E)=i/(E-L(E))ρ_0 can be determined self-consistently, leading to the effect that all
z_i acquire a finite imaginary part -iΓ_i.

In connection with the ohmic spin boson model at zero bias, the 
occurrence of exponentials in the O(α)-correction to Bloch-Redfield has recently been noted and 
corrected in Rfs. slutskin_etal_11,kashuba_schoeller_13. Similar considerations have been performed  
close to α∼1 2, see Rfs. egger_etal_97,kennes_etal_13. For finite bias, 
a Born approximation has been used in Ref. divincenzo_loss_05 to calculate perturbatively 
one order beyond Bloch Redfield, missing the exponentials in those corrections. In this paper 
we will present a perturbative calculation at arbitrary bias including all exponentials and, moreover, 
show that the resummation of secular terms is also important to obtain the correct energy scales in
logarithmic terms of preexponential functions. Furthermore, we will calculate all terms of the
time evolution for an arbitrary initial state of the local system, whereas in 
Ref. divincenzo_loss_05 only the time evolution of the Pauli matrix in z-direction has
been calculated for an initial state without any spin in x- and y- direction.

Besides secular terms proportional to powers of time, there are further subtleties in the calculation
of the time evolution, even in the case where potential logarithmic terms can be treated perturbatively. 
A generic feature of the reduced density matrix ρ(E)=i/(E-L(E))ρ_0 in Fourier space is that there
occurs one singularity at E=0 (determining the stationary state from L(0^+)ρ_st=0) and
a pure decay pole at E=z_0=-iΓ∼ O(α). These two singularities are close to each other 
within the expansion parameter α, and leads to the generic feature that two eigenvalues of 
L(E) are close to each other by O(α). Therefore, degenerate perturbation theory is necessary
for the zero and purely decaying modes, and the calculation of the corresponding projectors on the 
eigenstates of L(E) up to O(α) requires the knowledge of the Liouvillian at least up to 
O(α^2). This fact has already been mentioned at the end of Ref. divincenzo_loss_05, where the 
stationary state was calculated up to O(α) and the influence on the time evolution for the
purely decaying mode was indicated. This again is a generic problem for all models of dissipative
quantum mechanics and shows that lowest order Born approximation is not sufficient to account for all
first-order corrections to Bloch-Redfield. 
In this paper we will show that the special algebra of the ohmic spin model allows for
a simplification of this problem such that the results of 
Ref. divincenzo_loss_05 for the stationary case up to O(α) can be used to calculate also
all terms in O(α) for the time evolution of the purely decaying mode. 

The ohmic spin boson model (and similar many other models with a rather structureless spectral density of
states) has further problems in perturbation theory arising from logarithmically divergent integrals at
high and low energies, which have to be treated by renormalization group. At high energies logarithmic 
divergencies ∼αlnD/Λ_c occur, where D denotes the finite band width and
Λ_c∼max{1/t,Ω} is some high-energy cutoff determined by the largest energy scale 
of the system. For large D a non-perturbative resummation
of all powers of such terms is required. In the short-time regime 1/t≫Ω, this leads to 
well-known terms ∼ 1/(D t)^2α which can also be obtained from
the noninteracting blip approximation (NIBA) <cit.>. For the most important
regime of times which are not exponentially small or large, where |αln(Ω t)|≪ 1, 
we will show in this paper that the logarithmic terms at high energies can be 
incorporated into a renormalized tunneling Δ̃=Δ(Ω/D)^α, where 
Ω=√(ϵ^2+Δ̃^2) is the renormalized Rabi frequency of the local system,
leading also to a renormalized decay rate Γ=παΔ̃^2/Ω. We note that the
correct cutoff scale is Λ_c=Ω and not Δ̃ as first pointed out 
in Ref. divincenzo_loss_05,
where the logarithmic correction was calculated perturbatively in α. Furthermore, we
will show in this paper how the unrenormalized tunneling occurring in various terms of perturbation 
theory has to be replaced by the renormalized one. This is quite nontrivial since both Δ
and Δ̃ appear in the final solution. We will achieve this goal by solving the RTRG equations
perturbatively with the result of a renormalized propagator containing Z-factors with 
Z=Δ̃^2/Δ^2. Subsequently, we will apply renormalized perturbation theory to calculate 
the time evolution analytically in the whole crossover regime from small times Ω t ≪ 1
to large times Ω t≫ 1 with |αln(Ω t)|≪ 1 such that logarithmic terms in time can 
be treated perturbatively. We find that the leading order terms in the preexponential functions
stem from branch cuts starting at a pole position of ρ(E) giving rise to constant terms together with terms 
∼αln(Ω t) showing a rather weak logarithmic time dependence. In contrast, branch cuts 
starting at branching points unequal to the poles of ρ(E) lead to crossover functions with a
strong time dependence of the preexponential functions which all can be expressed by the exponential integral. 
Interestingly, for large times Ω t≫ 1, we find that the leading order terms fall off 
∼α/(Ω t) for finite bias, in contrast to the unbiased case, where all terms fall off 
∼α/(Ω t)^2.  

After having got rid of logarithmic terms at high
energies, one is still left with logarithmic terms at low energies ∼αlnΩ t. If the latter 
can be treated perturbatively, the solution for the time evolution one order beyond Bloch-Redfield  
follows from the above mentioned renormalized perturbation theory with the proper replacement of 
Δ by Δ̃. 
However, for intermediate couplings α∼ 0.1-0.2 or for the case of high bias 
ϵ≫Δ̃ (where the decay rate Γ≪αΩ is very small), it turns out 
that higher order terms ∼(αlnΩ t)^n with n>1 become important already for times
scales t∼ 1/Γ. In these cases a non-perturbative resummation is also necessary for the logarithmic 
terms at low energies to determine the first correction to Bloch-Redfield consistently. The only 
available method up to date to achieve such a resummation is the RTRG method 
<cit.>,
which can account simultaneously for logarithmic terms at high and low energies in all orders 
to determine the time evolution of models of dissipative quantum mechanics in the weak coupling regime. 
The idea is not to consider the perturbative expansion of the effective Liouvillian L(E) but of 
the second derivative d^2dE^2L(E), 
together with a proper resummation of self-energy insertions
and vertex corrections. This leads to a set of closed differential equations for the effective 
Liouvillian and the effective vertices, which are well-defined in the limit D→∞ and contain 
no secular terms and logarithmic divergencies at low and high energies. 
Therefore, the r.h.s. of these differential equations are a well-defined series in α and can be 
truncated systematically. We will consider the RG equations in leading order and solve them 
numerically for the ohmic spin boson model at arbitrary bias. Most importantly, we find for the
leading order terms of the preexponential functions of the oscillating modes a power-law behaviour 
∼ 1/(Ω t)^2αϵ^2Ω^2 for exponentially large times, where the power-law
exponent interpolates between 1 for ϵ=0 and 2α for 
ϵ≫Δ̃. The bias-dependent power-law exponent 2αϵ^2Ω^2 has 
also been proposed in Ref. slutskin_etal_11 but we stress that it is only correct for
very large times and we will show that, for small times, other logarithmic contributions appear
which lead to a complicated crossover to a power-law
∼ 1/(Ω t)^2α for exponentially small times. As already mentioned in 
Ref. kashuba_schoeller_13, the determination of the correct long-time behavior of 
preexponential functions depends crucially on the vertex renormalization not taken into account in any
previous work. At zero bias this has lead to a correction of the NIBA-result <cit.> and
we stress that all our results for nonzero bias presented in this paper can as well only be derived correctly by
including the vertex renormalization. 

The paper is organized as follows. In Section <ref> we
introduce the ohmic spin boson model and the kinetic equation to calculate the time dynamics. 
We provide the perturbative expansion of the effective Liouvillian in Fourier space and explain 
its analytic structure together with the one of the reduced density matrix. We also provide the
propagator in renormalized perturbation theory which will be derived later in 
Section <ref>. In Section <ref> we will explicitly 
calculate the time dynamics in various time regimes. We review the exact solution at zero tunneling
and the Bloch-Redfield solution in Sections <ref> and <ref>
as a reference. In Section <ref> we present the results from renormalized perturbation 
theory and determine the time evolution in the regime of small times in Section <ref>
and in the whole regime where time is not exponentially small or large in 
Section <ref>.
In Section <ref> we use the RTRG method to account for all logarithmic terms nonperturbatively. 
In Section <ref> we set up the RG equations for the ohmic spin boson model and show 
in Sections <ref> and <ref> how the propagator has to be changed 
to account for all logarithmic renormalizations from high energies. The numerical solution of the
RG equations containing also all logarithmic renormalizations at low energies
will be presented in Section <ref>. We close with a summary of our results 
in Section <ref> and discuss their relevance for other models of dissipative quantum
mechanics. We use the unit ħ=1 throughout this paper.




§ MODEL, KINETIC EQUATION AND LIOUVILLIAN



In this Section we introduce the model under consideration and set up the kinetic 
equation to determine the time dynamics of the local reduced density
matrix. In addition we provide the perturbative solution for the effective Liouvillian
in Fourier space. This form is very helpful to understand the proper analytical continuation 
into the lower half of the complex plane and the correct procedure to avoid the occurrence of 
secular terms. Furthermore, we will present the perturbative determination of the decay poles.





 §.§ Model



The Hamiltonian for the spin boson model consists of a local 2-level system
(described by Pauli matrices σ_i) coupled linearly to a bosonic bath with 
energy modes ω_q>0

    H_tot    = H + H_bath + V  ,
    
    H    = ϵ 2σ_z - Δ 2σ_x  ,
    
    H_bath    = ∑_q ω_q a_q^† a_q  ,
    
    V    = 1 2σ_z ∑_q g_q(a_q+a_q^†)  ,

where ϵ denotes the bias, Δ the tunneling, and the coupling to the bath
is described by the coupling parameters g_q. We note that by a convenient spin rotation the 
coupling to the bath can always be chosen in the z-direction and the y-axis can be chosen 
perpendicular to the local spin in the Hamiltonian (the expectation value of the local spin
will of course get all components as function of time). The parameters ϵ, Δ and
g_q are real to guarantee hermiticity of H_tot (please note that the sign convention for Δ 
is sometimes chosen differently in the literature). For convenience we choose Δ,ϵ>0 
which again can always be achieved by an appropriate spin rotation.

The microscopic details of the modes ω_q and the coupling constants g_q enter 
the time dynamics of the local system only via the energy dependence of the spectral density

    J(ω) = π∑_q g_q^2δ(ω-ω_q)  ,

which, for the ohmic spin boson model, is parametrized as

    J(ω) = 2παωθ(ω)J_c(ω)  ,

where α is a dimensionless coupling constant and J_c(ω) is a high-energy 
cutoff function needed since frequency integrals diverge logarithmically at high energies
for all terms in the perturbative series in α. In this paper we choose a
Lorentzian cutoff function (in contrast to exponential cutoffs ∼ e^-ω/D often
used in the literature)

    J_c(ω) = D^2 D^2 + ω^2  ,

where D denotes the band width. This choice is taken to simplify frequency integrals and 
influences only some prefactors of non-logarithmic terms but not the scaling behavior. The 
ohmic spin boson model in weak coupling is defined by the condition α≪ 1 such that 
a perturbative expansion in α makes sense. 

Since we will also work in a basis where the local Hamiltonian H is diagonal we introduce
the unitary transformation

    U = U^† = U^-1 = 1√(2Ω_0)([ -v_-  v_+;  v_+  v_- ])  ,

where v_±=√(Ω_0∓ϵ) and 

    Ω_0 = √(ϵ^2 + Δ^2)

denotes the bare level splitting (Rabi frequency) of the local system. With this unitary 
transformation we get U H U^† = 1 2Ω_0σ_z, i.e. the eigenvalues 
±Ω_0/2 with corresponding eigenvectors given by the two columns of U^†=U. 



 §.§ Kinetic equation



We aim at calculating the time dynamics of the reduced density matrix of the local system

    ρ(t) = Tr_bathρ_tot(t)

with an initial state for the total density matrix 

    ρ_tot(t=0) = ρ_0ρ_bath^eq

factorizing into an arbitrary
initial state ρ_0=ρ(t=0) for the local system and an equilibrium canonical distribution
ρ_bath for the bath. For simplicity we set temperature T=0 in the following. 
Using standard projection operator <cit.>, path integral <cit.>,
or diagrammatic <cit.> techniques one can show that ρ(t) can
be determined from a formally exact kinetic equation

    iρ̇(t) = L_0ρ(t) + ∫_0^t dt' Σ(t-t')ρ(t')  ,

where L_0=[H,·] and Σ(t-t') are superoperators acting on operators. The first term on
the r.h.s. describes the time evolution from the von Neumann equation of the isolated local system,
whereas the second term contains the dissipative kernel Σ(t-t') leading to irreversible time
dynamics into a stationary state ρ_st=lim_t→∞ρ(t). The various methods
described in Rfs. projection_operator,weiss_12,schoeller_09,schoeller_14 just differ in the
technique how to calculate this kernel in perturbation theory in α. Since all quantities are
only defined for positive times, we define the Fourier transform as for retarded correlation functions
(for convenience we use the same symbol for the Fourier transform)

    ρ(E) = ∫_0^∞ dt  e^iEtρ(t)  ,  Σ(E) = ∫_0^∞ dt  e^iEtΣ(t)  ,

which are well-defined analytic functions in the complex plane for all E with positive 
imaginary part (a proper analytic continuation into the lower half of the complex plane will 
be discussed later). From (<ref>) we obtain the formal solution in Fourier space as

    ρ(E) = i E-L(E)ρ_0  ,

where L(E)=L_0+Σ(E) denotes the effective Liouvillian in Fourier space with matrix elements
L_s_1 s_2,s_1^' s_2^' (s denote the states of the local system). The Liouvillian has
the two important properties <cit.>

    TrL(E)· = 0  ,  
    L(E)^c = -L(-E^*)  ,

where Tr denotes the trace over the local system and the c-transform is defined by 
L(E)^c_s_1 s_2,s_1^' s_2^'=L(E)_s_2 s_1,s_2^' s_1^'^*.
From these properties one can show the conservation of probability Trρ̇(t)=0 and the
hermiticity of the density matrix ρ(t)^†=ρ(t) <cit.>.

Once L(E) is known, the
time dynamics can be calculated from inverse Fourier transform as 

    ρ(t) = i 2π∫_C dE e^-iEt E-L(E)ρ_0  ,

where C is a straight line in the complex plane lying slightly above the real axis, i.e.
E=x+iη, with η=0^+ and x running from x=-∞ to x=+∞ (the precise form of 
C in the upper half is not important since ρ(E) is an analytic function there). We 
note that we have used the Fourier and not the Laplace transform (defined by e^-Et in 
(<ref>)) since it makes the analogy to the analytic properties of retarded correlation
functions more transparent.    

As pointed out in detail in 
Rfs. pletyukhov_schuricht_schoeller_10,schoeller_09,schoeller_14,kashuba_schoeller_13 the most 
elegant way to determine the integral over C is to close the integration contour in the
lower half of the complex plane and to use a convenient analytic continuation of L(E) and ρ(E)
into the lower half of the complex plane, such that all branch cuts point into the direction of the
negative imaginary axis and start at the branching points z_i=Ω_i-iΓ_i. For the ohmic
spin boson model we note that ρ(E) has one isolated pole at E=0 determining the stationary state from 

    L(0^+)ρ_st = 0  ,

together with three branch cuts starting at the branching points (or poles)

    z_0 = -iΓ  ,   z_±=±Ω-iΓ/2  ,

with Γ>0, whereas L(E) has only branch cuts starting at z_0 and z_± without any poles. 
If we denote the eigenvalues of the 4× 4-matrix L(E) by γ_i(E) with i=st,0,±, the
pole positions of the propagator 1/(E-L(E)) follow from γ_i(z_i)=z_i and it follows from 
(<ref>) that one eigenvalue must be zero and -γ_i(E)^* are the eigenvalues of
L(-E^*). Thus, -γ_i(-E^*)^* must be also an eigenvalue of L(E), leading to

    γ_st    = 0  ,
    γ_0(E)    = -γ_0(-E^*)^*  ,
    γ_+(E)    = -γ_-(-E^*)^*  .

As a consequence, the pole z_0 is purely imaginary and z_+=-z_-^*, in accordance with (<ref>).

Using the diagrammatic technique
of Rfs. schoeller_09,schoeller_14,pletyukhov_schoeller_12,kashuba_schoeller_13 
one can derive the analytic features in all orders of perturbation theory 
but it is illustrative to study them already from the perturbative solution for L(E) up to O(α), which
will be presented in the next Section.



 §.§ Liouvillian in perturbation theory



With the help of the diagrammatic technique used in Ref. kashuba_schoeller_13 for the
ohmic spin model at zero bias, we calculate the Liouvillian up to O(α) in 
Appendix <ref>. Denoting the two states of the local system by i=1,2
(corresponding to the original Hamiltonian H in (<ref>)) and using the sequence (11,22,12,21) to 
numerate the matrix elements of superoperators, we find:

    L(E)    = L_0 + Σ_a(E) + Σ_s = L_a(E) + Σ_s  ,
    
    L_0    = ([    0 Δτ_-; Δτ_- ϵσ_z ])  ,
    Σ_s    =  iπαΔ([   0   0; τ_+   0 ])  ,
    Σ_a(E)    = α∑_i=0,±F_i(E)M_i  ,  
    M_i=([    0    0;    0 M̂_i ])  ,
    M̂_0    = 2Δ^2Ω_0^2τ_-  ,
    M̂_±    = τ_+ ± ϵΩ_0σ_z  + ϵ^2Ω_0^2τ_-  ,

where τ_±=1 2(1±σ_x) and

    F_i(E)    = (E-λ_i(E))L_i(E)  ,
    L_i(E)    = ln-i(E-λ_i(E)) D  .

Here λ_i(E) are the important functions

    λ_0(E)    = -αΔ^2Ω_0∑_σ=±σL_σ(E)  ,
    λ_±(E)    = ±(Ω_0+αΔ^2Ω_0L_0(E))  ,

which determine the position of the poles (<ref>) of the resolvent 1/(E-L(E)) 
(and also of ρ(E) due to (<ref>)) by solving the self-consistent equations

    z_i = λ_i(z_i)  .

This can be seen from the derivation in Appendix <ref>, where the
λ_i(E) are defined as the eigenvalues of the Liouvillian L̃_Δ(E), defined by

    L̃_Δ(E) = Z'(E)L_Δ(E)  ,   Z'(E) = 1 1-L'(E)  ,

where L_Δ(E) and L'(E) follow from the decomposition

    L_a(E)    = L_Δ(E) + E L'(E)
        = L_0 + Σ_Δ(E) + E L'(E)  ,

with

    Σ_Δ(E)    = -α∑_i=0,±λ_i(E)L_i(E)M_i  ,
    
    L'(E)    = α∑_i=0,±L_i(E)M_i  .

This decomposition is very
helpful since it exhibits the purely logarithmic superoperators L_Δ(E) and L'(E), together with the 
terms linear in E. The eigenvalues of L(E) and L̃_Δ(E) are different but the relation
(note that Σ_s L_a=0)

    1 E-L(E)    = 1 E-L_a(E)(1+Σ_s1 E)
       = 1 E-L̃_Δ(E)Z'(E)(1+Σ_s1 E)  ,

shows that the poles of the two resolvents 1/(E-L(E)) and 1/(E-L̃_Δ(E)) are the same, i.e.
the solutions z_i of the self-consistent equations (<ref>) provide indeed the 
nonzero poles of the resolvent 1/(E-L(E)).

Most importantly, we see from the perturbative result (<ref>)-(<ref>) that z_i 
are not only the poles of the local density matrix in Fourier space but at the same time determine the
branching points of the logarithmic functions L_i(E), i.e. determine the starting points for 
the branch cuts of L(E) in the lower half of the complex plane. The logarithm in Eq. (<ref>) is
the natural logarithm with a branch cut on the negative real axis, i.e. the branch cut w.r.t. the 
Fourier variable E points into the direction of the negative imaginary axis, a choice which will be most
convenient for an analytical determination of the branch cut integral in the long time limit, see
Section <ref>. The fact that the branching points of all logarithmic terms are
the same as the pole positions of the local density matrix is a very important observation and can be 
shown to hold in 
all orders of perturbation theory by using the diagrammatic method developed in
Rfs. pletyukhov_schuricht_schoeller_10,schoeller_09,schoeller_14,kashuba_schoeller_13, see
also some remarks in Appendix <ref>. 
Obviously, for this property it is very important to keep the functions λ_i(E) in the argument of
the logarithm and not to expand L_i(E) in α. As already mentioned in 
Ref. schoeller_14 in all detail, such an expansion leads to secular terms (1/E)^n for 
the Liouvillian, e.g. for the expansion of αF_0(E) one obtains

    αF_0(E)    = α (E-λ_0(E))ln-iE D - αλ_0(E) 
        + 1 2αλ_0(E)^21 E + O(α^4)  .

We note that secular terms start at O(α^3) due to the factor E-λ_0(E) in front of the
logarithm. Therefore, even in a calculation up to O(α^2) one can not see the occurrence of secular
terms in L(E). 
The power of these secular terms increases with increasing order in α and, therefore, have to be
resummed nonperturbatively. They appear directly in the effective Liouvillian L(E) and have to be 
distinguished from secular terms appearing by expanding the resolvent 1/(E-L_0-Σ(E)) in Σ(E).
The resummation of the latter are responsible to obtain the correct exponential behavior 
of the leading order Bloch-Redfield terms for the time evolution, whereas the ones in L(E) have to be 
resummed to obtain the exponential part of all correction 
terms to the Bloch-Redfield solution. Essentially, the fact that logarithmic functions in all
orders of perturbation theory appear always in the form of L_i(E) is due to the property that
all bare propagators of the local system can be replaced by full propagators without any double counting,
see Appendix <ref>. As a consequence the exact eigenvalues of L̃_Δ(E)
appear in the perturbative series and not the bare ones. This fact is very important to notice in
order to find the correct nonanalytic features in the lower half of the complex plane. E.g. by calculating 
F_0(E) only by the first term on the r.h.s. of (<ref>) one obtains a logarithm which has a 
branch cut starting at the origin leading to a term of the time evolution which is not exponentially 
decaying. The expansion (<ref>) is only well-defined for E∼Ω_0, i.e. on time scales
t∼ 1/E ∼ 1/Ω_0, where the solution is just oscillating and the decay has not yet set in.   
In this regime the perturbative solution of Ref. divincenzo_loss_05 can be used but 
not for larger time scales describing the crossover to the regime of exponential decay. 

We note that the perturbative solution (<ref>)-(<ref>) for L(E) can only be used
when the logarithmic terms are small enough, i.e. the condition 

    α|ln-i(E-λ_i(E)) D| ≪ 1

should hold. This is obviously not fulfilled when E approaches the branching point z_i or is too far
away from it. Only the RG method presented in Section <ref> is capable of resumming the logarithmic 
terms in all orders to find the correct scaling behavior for large E or E close to z_i. The condition 
(<ref>) can be reformulated in terms of time by replacing E-λ_i(E)→ 1/t 
leading to

    α|ln(Dt)| ≪ 1  ,

showing that the perturbative theory can not be used to calculate the time evolution for exponentially small or 
large times. However, as we will see in Section <ref> these regimes can be studied as well by using the
RTRG method. 

As a consequence one should also not be concerned by the fact that the solution of the self-consistent
equations (<ref>) with (<ref>) and (<ref>) is 
ill-defined due to the singularity of the logarithm. For times in the regime
(<ref>) we need the functions λ_i(E) only in the typical regime 
(<ref>). Using z_0∼ O(α) and z_±=±Ω_0+O(α), 
this means that for |E-z_0|∼α^nΩ_0 (with some integer n>0) we can replace λ_0(E) by

    λ_0(E) ≈ -αΔ^2Ω_0∑_σ=±σln-i(-σΩ_0) D = -iΓ_1  ,

with

    Γ_1 = παΔ^2Ω_0  ,

up to an error of O(nα^2lnα). Up to the same error,
for |E-z_±|∼α^nΩ_0, we can replace λ_±(E) by

    λ_±(E) ≈ ±(Ω_0+αΔ^2Ω_0ln-i(±Ω_0) D)
     = ±Ω_1-iΓ_1/2  ,

with 

    Ω_1 = Ω_0 - αΔ^2Ω_0lnDΩ_0  .

Therefore, we conclude from the perturbative expansion that the solution of 
(<ref>) is given by

    z^(1)_0    = -iΓ_1 + O(α^2lnα)  ,
    
    z^(1)_±    = ±Ω_1-iΓ_1/2 + O(α^2lnα)  .

In Section <ref> we will resum all logarithmic renormalizations
∼(αln(Ω/D))^n from high energies and show that Ω_1 has to be 
replaced by the renormalized Rabi frequency Ω, which has the same form as Ω_0
but the bare tunneling Δ has to be replaced by the renormalized tunneling 
Δ̃ 

    Ω    = √(ϵ^2 + Δ̃^2)  ,
    Δ̃    = Δ(Ω D)^α = Δ(√(ϵ^2+Δ̃^2) D)^α  .

We note that the low-energy scale cutting off the logarithmic terms in this expression is set by
Ω but not  by the renormalized tunneling as has been stated e.g. in Ref. weiss_12. 
This was already mentioned in Ref. divincenzo_loss_05, where the oscillation 
frequency has been calculated perturbatively up to the first logarithmic term, as given by 
Eq. (<ref>). Furthermore, we note, that besides the logarithmic terms there can be 
other regular terms ∼α^n which depend on the specific high-energy cutoff function
under consideration. The logarithmic terms however are universal, i.e. do not depend on the
specific form of the high-energy cutoff function. This will be explained in Section <ref>. 

Inserting the propagator (<ref>) in (<ref>) and using the
perturbative result (<ref>-<ref>) for L(E), one can systematically 
determine the time dynamics one order beyond Bloch-Redfield using the scheme presented in 
Section <ref>. However, this calculation can be easily improved by using
renormalized perturbation theory, where the renormalized tunneling has to be used at appropriate
places and renormalizations of Z-factors are important. Therefore, although the 
explicit results will be derived later on in Section <ref> using the RTRG method, 
we already state the results for the propagator in the next subsection such that we can use
it in Section <ref> for the perturbative calculation of the time dynamics. 




 §.§ Liouvillian in renormalized perturbation theory



In Section <ref> we will show how the propagator 1/(E-L(E)) has to be slightly modified to
account for all logarithmic renormalizations from high energies. There are two different kinds of 
logarithmic terms, one involving powers of αln(D/Ω) (which can be resummed in the
renormalized tunneling (<ref>)), the other containing powers of logarithmic terms 
αln(Ω t) in time. The latter can be treated perturbatively provided that time is not
exponentially small or large. This defines the regime which we call the regime of times in 
the non-exponential regime 

    |αln(Ω t)| ≪ 1  ,

which corresponds in Fourier space to the regime

    |αln-i(E-z_i)Ω| ≪ 1  .

This is the regime where renormalized perturbation theory can be applied. 
In Section <ref> we will show that in this regime the propagator 
can be written as 

    1 E-L(E)    ≈ 1 E-L̃_a(E)Z^'(1+Σ_s1 E)  ,

with

    Z^'    = ([ 1 0; 0 Z ])  ,  
    Z = Δ̃^2Δ^2
    L̃_a(E)    = L̃_0 + Σ̃_a(E)   ,
    L̃_0    = ([     0  Δτ_-; ZΔτ_-  ϵσ_z ])  ,
    Σ̃_a(E)    = α∑_i=0,±F_i(E)M_i  ,  
    M_i=([    0    0;    0 M̂_i ])  ,
    M̂_0    = 2Δ̃^2Ω^2τ_-  ,
    M̂_±    = τ_+ ± ϵΩσ_z  + ϵ^2Ω^2τ_-  ,

where F_i(E) is defined by (<ref>) with

    λ_0(E)    = -αΔ̃^2Ω∑_σ=±σL_σ(E)  ,
    λ_±(E)    = ±(Ω+αΔ̃^2ΩL_0(E))  ,

and 

    L_i(E)    = ln-i(E-λ_i(E))Ω  .

In comparison to the unrenormalized perturbation theory (<ref>-<ref>) we see that the
renormalized Rabi frequency and the renormalized tunneling appear in Σ̃_a(E) and λ_i(E)
instead of the bare ones and the band width D is replaced by Ω in the logarithmic function 
L_i(E). In addition, L_0 and the propagator
get a renormalization from the Z'-matrix containing the Z-factor Z=Δ̃^2/Δ^2. 
In Section <ref> we will see that Z can be obtained from a poor man scaling equation 
for Z(E)=(-iE/D)^2α cut off at E=iΩ. Our result 
shows that renormalized perturbation theory is not obtained by just replacing 
Δ→Δ̃ defining a local system with a renormalized tunneling. Instead, the
Liouvillian L̃_0 is no longer hermitian, i.e. can essentially be not expressed as a 
commutator with a renormalized local Hamiltonian. 

Since the solutions of z_i=λ_i(z) again define the positions of the poles of the propagator, 
the logarithmic renormalizations from high energies lead in analogy to (<ref>) and (<ref>)
to the renormalized pole positions

    z_0    = -iΓ + O(α^2)  ,
    
    z_±    = ±Ω-iΓ/2 + O(α^2)  ,

with 

    Γ = παΔ̃^2Ω  .


In Section <ref> we will also discuss the regimes of exponentially small or large
times where the condition (<ref>) fails and higher powers of logarithmic 
terms have to be resummed by a proper RG method for the ultraviolett regime (small times or
large energies) and the infrared regime (large times or energies close to the pole positions). 
Although this regime is certainly of minor interest to quantum information processing it is
of high interest from a theoretical point of view since various power-laws appear which are
qualitatively very different in the ultraviolett and infrared regime. Furthermore, these
power laws are not only of academic interest in unrealistic time regimes since they become 
clearly visible for moderate α∼ 0.05 and, moreover, second
order terms ∼(αln(Ω t))^2 can become of order α already for time
scales t∼ 1/Γ where the decay is still moderate depending on the ratio of 
Ω/Δ̃. Using (<ref>) we find for t∼ 1/Γ

    (αln(Ω t))^2∼α  ⇔  ΩΔ̃∼√(πα)  e^1/(2√(α))  ,

leading e.g. to Ω/Δ̃∼ 4 for α=0.05. This are quite realistic values
showing that higher powers of logarithmic terms contribute significantly on the same level as
corrections ∼α to the Bloch-Redfield solution. Although the terms ∼αln(Ω t)
are the leading order terms in this regime, the second order terms ∼(αln(Ω t))^2 are
clearly visible in the time dynamics of the preexponential functions showing a significant 
deviation from a straight line plotted logarithmically as function of ln(Ω t), see 
Section <ref>. Thus, for the spin boson model at finite
bias, the systematic calculation of corrections to Bloch-Redfield is quite subtle and requires
an analysis of higher-order terms beyond O(α) for the Liouvillian for various reasons.

In Section <ref> we will see that the resummation of logarithmic terms in time is very
complicated in the infrared regime and requires a careful solution of the full RG equations, which
we will perform numerically. In contrast, the resummation of time-dependent logarithmic terms in the 
ultraviolett regime is quite straightforward since, for large energies, the energy scales of the
local system do not play an important role and can be treated perturbatively. Therefore, we state
here also the result for the propagator in the regime of small times defined by

    1 D ≪ t ≪ 1Ω  ,

corresponding to the regime of large energies 

    Ω ≪  |E|  ≪  D   .

We note that resumming all logarithmic terms ∼(αln(E/D))^n or ∼(αln(Dt))^n leads to a 
universal result for the time evolution in the regime |αln(Ω t)|∼ 1 and t≫ 1/D 
(where all corrections of O(α) and O(1/(Dt) can be neglected), in contrast to the 
non-universal regime t≲ 1/D, where
bare perturbation theory in α can be used to determine ρ(t) and the result depends crucially
on the shape of the high-energy cutoff function J_c(ω).   

For large energies E∼ 1/t ≫Ω, we neglect all terms 
of relative order αΩ/E∼αΩ t in L̃_Δ(E)
and Z'(E) and find in Section <ref> that

    L̃_Δ(E)    ≈ L̃_0(E)(1 + O(αΩ/E))  ,
    
    Z^'(E)    ≈ ([    1    0;    0 Z(E) ])
    (1 + O(αΩ/E))  ,

with

    L̃_0(E)    = ([         0      Δτ_-; Δ Z(E)τ_-      ϵσ_z ])  ,
    
    Z(E)    = (-iE D)^2α  .

Since Σ_s/E∼αΩ/E can also be neglected in
(<ref>) we find for the propagator the approximation

    1 E-L(E) ≈ 1 E - L̃_0(E)([    1    0;    0 Z(E) ])  .

In Section <ref> we will see that the form for Z'(E) results from a poor man scaling
equation cut off at the largest energy scale E which corresponds to 1/t in time space. If E 
becomes of the order Ω the Z-factor is cut off at E=iΩ, leading to the Z-factor
(<ref>) used in the regime where time is not exponentially small or large.   

The form (<ref>) can be used in the whole regime Ω≪ E≪ D, irrespective
of whether E is exponentially large or not. Thus, we can also use it in the regime where
|αln(-iE/Ω)|≪ 1, where we can expand Z(E) as

    Z(E) = Δ̃^2Δ^2 (1+2αln-iEΩ)  ,

and, after a straightforward calculation, one finds that the propagator (<ref>)
at high energies obtains the same form in leading order in α and Ω/E as the propagator 
(<ref>) in the regime of non-exponentially large energies. 





§ TIME DYNAMICS



In this section we will present the time dynamics of the local density matrix analytically in the
regimes of small times (including the case of exponentially small times) and for the regime of
times which are not exponentially small or large, where renormalized perturbation theory can be applied using
the propagator presented in Section <ref>. The exact solution
for zero tunneling and the lowest order Bloch-Redfield solution will be rederived 
in Sections <ref> and <ref> for reference. 
In Section <ref> we will present renormalized perturbation theory 
to show how the Bloch-Redfield solution has to be modified, together with the systematic calculation of the 
next correction in O(α). For the most interesting regime of times which are not exponentially 
small or large we note that our analytic solution has never been obtained correctly 
in the literature before.




 §.§ Exact solution at zero tunneling



For zero tunneling the time dynamics can be calculated exactly even for an arbitrary 
spectral density and finite temperatures <cit.>. In this case the
local Hamiltonian H=σ_zϵ/2 decouples from the rest and the coupling to the bath 
can be eliminated by a unitary transformation shifting the field operators of the bath

    H_tot    = H + e^σ_z χH_rese^-σ_zχ + c  ,
    χ    = ∑_q g_q (a_q + a_q^†)  ,

with an unimportant constant c=∑_qω_q g_q^2 dropping out for the time dynamics. After a 
straightforward calculation, the time dynamics for the diagonal and non-diagonal matrix elements of
ρ(t) follows as

    ρ(t)_σσ    = ρ(0)_σσ  ,
    ρ(t)_σ,-σ    = e^-iσϵ t⟨ e^2(χ(t)-χ)⟩_resρ(0)_σ,-σ  ,

where σ=±≡ 1,2 denotes the two local states, 
χ(t) is the Heisenberg picture w.r.t. H_res, and ⟨⋯⟩_res
denotes the expectation value w.r.t. to the canonical equilibrium distribution of the reservoir. 
Calculating this average by standard means gives the following result for the expectation values of
the Pauli matrices of the local system

    ⟨σ_x⟩(t)    = e^-h(t) {cos(ϵ t)⟨σ_x⟩(0)
    -sin(ϵ t)⟨σ_y⟩(0)}  ,
    ⟨σ_y⟩(t)    = e^-h(t) {sin(ϵ t)⟨σ_x⟩(0)
    +cos(ϵ t)⟨σ_y⟩(0)
    }  ,
    ⟨σ_z⟩(t)    = ⟨σ_z⟩(0)  ,

with

    h(t) = 1π∫_0^∞ dω J(ω)(1+2 n(ω))1-cos(ω t)ω^2  ,

where n(ω) is the Bose distribution function which vanishes at zero temperature. Thus, at zero
temperature, we get for the ohmic case 

    h(t) = 2α∫_0^∞ dω J_c(ω) 1-cos(ω t)ω  ,

which contains a logarithmic divergence at large ω. Therefore, in the limit Dt≫ 1, we get
the result

    h(t) ≈ 2α(γ + ln(Dt))  ,

where γ is Euler's constant. This leads to the universal power-law 

    e^-h(t)    ≈ (1-2αγ)(1 Dt)^2α
       = (1-2αγ)Δ̃^2Δ^2(1Ω t)^2α

for the time dynamics, where we have written the factor in front up to O(α) in order to compare
it later on to our perturbative solution for arbitrary tunneling. For the second form we have used
(<ref>) to write the result independent of D parametrizing it by the ratio of the renormalized 
tunneling to the unrenormalized one (which is finite even in the limit of zero tunneling). 

As one can see the result (<ref>) contains a resummation of all powers of 
logarithmic terms ∼(αln(Dt))^n and, thus, can only be obtained from the RG procedure presented in
Section <ref>. It will turn out that it holds even at finite tunneling Δ̃≪ϵ, 
provided that the condition Ω t≫ 1 ≫Γ t holds. 




 §.§ Bloch-Redfield solution


The easiest way to derive the Bloch-Redfield solution is to insert
(<ref>) in (<ref>) and use the spectral 
decomposition of the Liouvillian L̃_Δ(E). This gives the formally exact expression

    ρ(t)    = i 2π∑_i=st,0,±∫_C dE e^-iEt E-λ_i(E)·
       · P_i(E)Z'(E)(1+Σ_s1 E)ρ_0  .

Here, λ_i(E) are the eigenvalues of L̃_Δ(E) and P_i(E) are the 
corresponding projectors. These quantities can be calculated by solving for the right and left
eigenstates of L̃_Δ(E)

    L̃_Δ(E)|x_i(E)⟩    = λ_i(E)|x_i(E)⟩  ,
    ⟨x̅_i(E)|L̃_Δ(E)    = ⟨x̅_i(E)|λ_i(E)  ,
    
    P_i(E)    = |x_i(E)⟩⟨x̅_i(E)|  .

The projectors fulfill the property

    P_i(E) P_j(E) = δ_ijP_i(E)  ,  ∑_i P_i(E) = 1  .

We note that the eigenvalues are complex since the superoperator L̃_Δ(E) is a
non-hermitian matrix. One of the eigenvalues is zero (denoted by i=st) 
and the corresponding right/left eigenstates are exactly known in all orders of perturbation 
theory 

    λ_st    = 0   ,
    
    |x_st⟩    = 1√(2)([ 1; 1; 0; 0 ])  ,  ⟨x̅_st| = 1√(2)([ 1 1 0 0 ])
    
    P_st    = ([ τ_+   0;   0   0 ])  .

This can be seen from the matrix structure (<ref>-<ref>), which holds in all
orders of perturbation theory, see Appendix <ref> for the proof. 
We note that the right eigenstate |x_st(E)⟩ for E=0^+ does not give the stationary state 
ρ_st, following from (<ref>), since the eigenstates of 
L̃_Δ(E) and L(E) are different.

The eigenvalues λ_i(E) for i=0,± have already been provided in
perturbation theory up to O(α) in (<ref>) and (<ref>). 
Since P_stZ'(E)=P_st and ⟨x̅_st|Σ_s=0 we note that
the second term involving Σ_s contributes only for ist. 

The Bloch-Redfield solution is obtained by taking P_i(E)Z'(E)≈ P_i^(0) in lowest order in 
α (which is independent of E) and taking the Markovian approximation
λ_i(E)≈λ_i(z_i)=z_i, which again neglects O(α) contributions from the
residua and further corrections arising from possible branch cuts starting at z_i. The pole positions
are taken from (<ref>), (<ref>) and z_st=0. This gives the result

    ρ^(0)(t)    = i 2π∑_i=st,0,±∫_C dE e^-iEt E-z_i^(1)
    P_i^(0)(1+Σ_s1 E)ρ_0
        = (e^-iz^(1)t-1)1 z^(1)_0 P_0^(0)Σ_S ρ_0 + P_stρ_0
        + ∑_i=0,± e^-iz^(1)_it P_i^(0)ρ_0  .

The first term on the r.h.s. arises from the pole at E=0 
from the term Σ_s/E. It is of 
O(1) since 1/z^(1)_0=i/Γ_1∼ 1/α, in contrast
to the contributions from 1/z^(1)_±=1/(±Ω_1-iΓ_1/2)∼ O(1) 
which lead to an O(α) correction to ρ(t).

The projectors in lowest order are the ones for L_0. We note that there is no problem with
degenerate perturbation theory for the two eigenvalues λ_st=0 and λ_0∼α
(requiring in general a knowledge of L̃_Δ up to O(α) to calculate P_st 
and P_0 in lowest order) since the projector P_st is exactly known from (<ref>) 
in all orders of perturbation theory such that P_0^(0)=1-P_st-P_+^(0)-P_-^(0) can be used. 
The projectors for L_0 can be most easily obtained by transforming the matrix L_0 to the basis
of the exact eigenstates of H, which, by using the unitary matrix (<ref>), is described 
by the unitary transformation (A_0)_ij,kl=U_ikU_jl^* leading to

    A_0 = A_0^† = A_0^-1 = 1Ω_0([  Ω_0τ_++ϵτ_-     -Δσ_zτ_+;     -Δσ_zτ_- -ϵτ_+-Ω_0τ_- ]).

In the new basis L_0 is given by

    A_0 L_0 A_0^† = ([      0      0;      0 Ω_0σ_z ])  ,

and the projectors obviously follow from

    A_0 P_stA_0^†    = ([ τ_+   0;   0   0 ])  ,
    
    A_0 P_0^(0) A_0^†    = ([ τ_-   0;   0   0 ])  ,
    
    A_0 P_σ^(0) A_0^†    = ([           0           0;           0 1 2(1+σσ_z) ])  .

Transforming back with the matrix A_0 to the original basis one obtains straightforwardly
the result 

    P_0^(0)    = 1Ω_0^2([    ϵ^2τ_- -Δϵσ_zτ_+; -Δϵσ_zτ_-    Δ^2τ_+ ])  ,
    
    P_σ^(0)    = 1 2Ω_0^2([            Δ^2τ_-          Δϵσ_zτ_+;          Δϵσ_zτ_- ϵ^2τ_+ + Ω_0^2τ_- ]) +
       + σ 2Ω_0([    0 Δτ_-; Δτ_- ϵσ_z ])  .


Inserting (<ref>), (<ref>), (<ref>) and
(<ref>) in (<ref>) we obtain the Bloch-Redfield solution. Using the
formulas (<ref>) and (<ref>) for the pole positions, we can 
decompose the time evolution of the Pauli matrices generically as

    ⟨σ_α⟩(t)    = ⟨σ_α⟩_st + 
    F^0_α(t)e^-Γ_1 t + 
       + F^c_α(t)e^-Γ_1 2 tcos(Ω_1 t) + 
    F^s_α(t)e^-Γ_1 2 tsin(Ω_1 t)  ,

with α=x,y,z. F_α^0,c,s(t) denote the preexponential functions, which 
become time independent in Bloch-Redfield approximation 

    ⟨σ_x⟩_st    = ΔΩ_0  ,  ⟨σ_y⟩_st = 0   ,  ⟨σ_z⟩_st = -ϵΩ_0  ,
    
    F^0_x    = -⟨σ_x⟩_st-ΔΩ_0⟨σ_z^'⟩_0  ,  
    
    F^0_y    = 0   ,  
    
    F^0_z    = -⟨σ_z⟩_st+ϵΩ_0⟨σ_z^'⟩_0  ,
    
    F^c_x    = -ϵΩ_0⟨σ_x^'⟩_0  ,
    
    F^c_y    = ⟨σ_y⟩_0   ,
    
    F^c_z    = -ΔΩ_0⟨σ_x^'⟩_0  ,
    
    F^s_x    = -ϵΩ_0⟨σ_y⟩_0  ,  
    
    F^s_y    = -⟨σ_x^'⟩_0  ,
    
    F^s_z    = -ΔΩ_0⟨σ_y⟩_0  ,

where

    σ_x^'    = 
    -1Ω_0(ϵσ_x + Δσ_z)  ,
    σ_y^'    = -σ_y  ,
    σ_z^'    = 1Ω_0(ϵσ_z - Δσ_x)

are the Pauli spin operators in the basis where the local Hamiltonian is diagonal.

For later reference, we also state the form of the Bloch-Redfield solution
in the regime of small times where Ω_0 t≪ 1. Expanding the exponentials up to linear
order in Ω_1 t and neglecting Γ_1t,(Ω_1-Ω_0)t∼αΔ t  
we obtain

    ⟨σ_x⟩(t)    = ⟨σ_x⟩_0 - ϵ t⟨σ_y⟩_0  ,
    ⟨σ_y⟩(t)    = ⟨σ_y⟩_0 + ϵ t⟨σ_x⟩_0
    + Δ t⟨σ_z⟩_0  ,
    ⟨σ_z⟩(t)    = ⟨σ_z⟩_0 
    - Δ t⟨σ_y⟩_0  .






 §.§ Renormalized perturbation theory



Using the propagators provided in Section <ref> we will now apply
renormalized perturbation theory to calculate the modification of the Bloch-Redfield solution in
lowest order in α (but including all logarithmic corrections ∼(αln(Dt))^n and
∼(αln(D/Ω))^n from high energies in all orders) together with the first systematic
correction in O(α) to the Bloch-Redfield solution. Since renormalized perturbation theory
can only be applied analytically in the regimes of small times or times which are not exponentially
small or large, we will restrict our analysis to these two regimes and find that the two solutions
coincide for small but not exponentially small times, so that also the crossover between these two
regimes can be described with our analytic results, providing a systematic analytic solution beyond
Bloch-Redfield in the most interesting regime for quantum information where the exponential decay 
has not yet destroyed the time dynamics completely. Only the regime of very large times where higher 
powers of logarithmic terms like α^2ln^2(Ω t) become important is not treated analytically 
and will be presented in Section <ref> via a numerical solution of the RG equations. 




  §.§.§ Small times



For small times Ω t ≪ 1 but still in the universal regime t≫ 1/D we 
take the form (<ref>) for the propagator and, since
E∼ 1/t≫ϵ,Δ, can expand the resolvent up to first order in L̃_0(E)

    1 E - L̃_0(E) ≈ 1 E+1 EL̃_0(E)1 E  .

In this way we keep all terms ∼L̃_0/E∼ϵ t,Δ t, which, for Ω t∼α,
can be of the same order as the first correction ∼α to the Bloch-Redfield result.

Inserting (<ref>) and (<ref>) in
(<ref>), using the integrals

    I_1(t)    = i 2π∫_C dE e^-iEtZ(E) E
       = sin(2πα) 2παΓ(1+2α)(1 Dt)^2α
    
    I_2(t)    = i 2π∫_C dE e^-iEtZ(E) E^2
       = -it 1-2αI_1(t)  ,

where Γ(x) denotes the Gamma function with 
Γ(1+x)=1-γ x+O(x^2), and neglecting all terms of O(α^2), 
O(αϵ t) and O(αΔ t), we find for the local density matrix the 
following result in the short time limit

    ρ(t)    = ([ 1 0; 0 0 ])ρ_0
     + 
       + (1 Dt)^2α([               0        -iΔ tτ_-;        -iΔ tτ_- 1-2αγ - iϵ tσ_z ])ρ_0  .

For the time dynamics of the Pauli matrices this gives

    ⟨σ_x⟩(t)    = (1 Dt)^2α{(1-2αγ)⟨σ_x⟩_0 - ϵ t⟨σ_y⟩_0}  ,
    ⟨σ_y⟩(t)    = (1 Dt)^2α{(1-2αγ)⟨σ_y⟩_0 + .
       . +ϵ t⟨σ_x⟩_0
    + Δ t⟨σ_z⟩_0}  ,
    ⟨σ_z⟩(t)    = ⟨σ_z⟩_0 
    - (1 Dt)^2αΔ t⟨σ_y⟩_0  .

We see that this solution contains a power law arising from a resummation of all leading logarithmic terms 
∼(αln(Dt))^n, which appears also in the NIBA approximation <cit.>. We note
that it is not allowed to set t=0 since this result is only valid for t≫ 1/D, i.e. terms  
∼ϵ/D,Δ/D≪ϵ t, Δ t are neglected. 

We can study the short-time solution in two different regimes, the one for exponentially small times
|αln(Ω t)|∼ 1 where we can neglect all terms ∼Δ t and ∼ϵ t, and
the one for small but not exponentially small times |αln(Ω t)|≪ 1 where 
only terms of O(α), O(αln(Ω t)), O(Δ t) and O(ϵ t) need to be considered.
Using (<ref>) we obtain for exponentially small times 

    ⟨σ_x⟩(t)    = Δ̃^2Δ^2(1Ω t)^2α
    (1-2αγ)⟨σ_x⟩_0   ,
    ⟨σ_y⟩(t)    = Δ̃^2Δ^2(1Ω t)^2α
    (1-2αγ)⟨σ_y⟩_0   ,
    ⟨σ_z⟩(t)    = ⟨σ_z⟩_0   ,

and for small but not exponentially small times

    ⟨σ_x⟩(t)    = Δ̃^2Δ^2{(1-2α(γ+ln(Ω t))⟨σ_x⟩_0 
    - ϵ t⟨σ_y⟩_0}  ,
    ⟨σ_y⟩(t)    = Δ̃^2Δ^2{(1-2α(γ+ln(Ω t))⟨σ_y⟩_0 + .
       . +ϵ t⟨σ_x⟩_0
    + Δ t⟨σ_z⟩_0}  ,
    ⟨σ_z⟩(t)    = ⟨σ_z⟩_0 
    - Δ̃^2Δ t⟨σ_y⟩_0  .


For zero tunneling Δ=0 the short-time solution is consistent with the exact solution
(<ref>-<ref>) where we set cos(ϵ t)≈ 1
and sin(ϵ t)≈ϵ t. In contrast, the Bloch-Redfield solution
(<ref>-<ref>) at small times
misses all powers of logarithmic terms αln(D/Ω) (resummed in Δ̃) and
αln(Ω t) together with the O(α) corrections for 
⟨σ_x⟩(t) and ⟨σ_y⟩(t). 

In the next section we will show that our analytic solution for times which are not exponentially 
small or large coincides with 
(<ref>-<ref>) in the regime of small 
but not exponentially small times. This shows that by combining the
solution (<ref>-<ref>) for small times with the solution
of the next section we have an analytic and systematic result one order beyond Bloch-Redfield covering the whole 
time regime from Ω/D≪Ω t≪ 1 up to times Ω t≫ 1 which are not exponentially large
(i.e. |αln(Ω t)|≪ 1). 




  §.§.§ Times in the non-exponential regime



We now study the regime of times which are not exponentially small or large defined by the
condition |αln(Ω t)|≪ 1. Here we can use the propagator in the form presented in 
(<ref>-<ref>) and apply renormalized perturbation 
theory to study the modification of the Bloch-Redfield result and to calculate the next correction
in O(α). We want to determine analytically the whole crossover regime from Ω t≪ 1 
up to Ω t≫ 1 provided that time is not exponentially small or large such that
all logarithmic terms ∼|αln(Ω t)|≪ 1 can be treated perturbatively and 
are on the same level as terms ∼α. In particular this
includes the long-time regime where decay sets in such that Γ t∼ 1 or
Ω t∼ 1/α≫ 1. In this long-time regime we have to be very careful not to expand the
resolvent 1/(E-L̃_0-Σ̃_a(E)) in Σ̃_a(E)∼αΩ∼ 1/t∼ |E-z_i| since
|E-z_i| sets the scale of the lowest order term in the denominator of the resolvent for the pole contributions.
This would be only allowed in the regime Ω t≲ 1 but can not be used to study
the crossover to the long-time regime. Furthermore, in order to calculate systematically 
the first correction to the Bloch-Redfield result in the long-time regime Ω t∼ 1/α,
it is also necessary to discuss carefully terms 
∼α^2Ω∼α 1/t∼α|E-z_i| in Σ̃_a(E). As we will see this requires
a knowledge of certain terms in O(α^2) of the Liouvillian but it will turn out that the contributions
of these terms to the time dynamics of ρ(t) can all be related to the stationary solution up to 
O(α^2) which can be calculated quite efficiently in equilibrium via the partition function, 
see Ref. divincenzo_loss_05.

To account for all these subtleties systematically we proceed as follows. Since we know in all orders
of perturbation theory that the non-analytic features of the propagator are an isolated pole at 
E=z_st=0 together with branch cuts starting at E=z_i, i=0,±, pointing in the direction of the 
negative imaginary axis, we can decompose the time dynamics of ρ(t) in four contributions

    ρ(t) = ρ_st + ∑_i=0,±ρ_i(t)  ,

with

    ρ_i(t)    = i 2π∫_C_i dE e^-iEt1 E-L̃_0-Σ̃_a(E)·
       · Z'(1+Σ_s1 E)ρ_0  ,

where C_i is a curve in the complex plane encircling clockwise the non-analytic feature 
around E=z_i (i.e. an isolated pole for i=st and a branch cut at E=z_i-ix, x>0, 
for i=0,±). Since the zero eigenvalue of L̃_0+Σ̃_a(E) is projected out by the
projector P_st given (in all orders of perturbation theory) by (<ref>), we get 
for the stationary state 

    ρ_st    = P_stρ_0  - 1L̃_0+Σ̃_a(0)Z'Σ_sρ_0
       = 1 2([ 1; 1; 0; 0 ])
     - iπαΔ ZL̃_0+Σ̃_a(0)1 2([ 0; 0; 1; 1 ])  ,

where we have taken (<ref>) and (<ref>) for Σ_s and Z',
respectively, and have used the normalization Trρ_0=1. For i=0,± we obtain with
E=-ix±η (η=0^+)

    ρ_i(t) = F_i(t) e^-iz_it  ,

with the preexponential operator given by 

    F_i(t)    = 1 2π∫_0^∞ dx  e^-xt·
       ·{1 E-L̃_0-Σ̃_a(E)|_E=z_i-ix+η - (η→ -η)}·
       · Z'(1+Σ_s1 z_i-ix)ρ_0  .

Due to the exponential part e^-xt in (<ref>), we get x∼ 1/t. 
The eigenvalues of L̃_0 are either zero or ±Ω (see below),
and Σ̃_a(E)∼αΩ. Thus, for times Ω t≲ 1, Σ̃_a(E)
is a small correction in the denominator and we can expand the resolvent in Σ̃_a(E).
However, for times Ω t∼ 1/α or |E-z_i|∼αΩ, 
Σ̃_a(E)∼αΩ∼ 1/t becomes of 
the same order as x∼ 1/t and the expansion is no longer valid. To cover the crossover to this regime
as well we leave the important term Σ̃_a(z_i)∼αΩ in the
denominator which is essential for the correct position of the poles, and expand only in

    Σ̃_a(E)-Σ̃_a(z_i)    = 
       
    = αF_i(E)M_i + α∑_j=0,±
     j i(F_j(E)-F_j(z_i))M_j
       ≈ α(E-z_i){ln-i(E-z_i)ΩM_i + 
    ∑_j=0,±
     j idF_j dE(z_i)M_j}
       ∼ α(E-z_i) ∼ α t ≪ 1 t ∼ x  ,

where we have used the form (<ref>) and λ_i(z_i)=z_i 
(see (<ref>)), together with the fact that F_j(E) 
can be expanded around E=z_i for j i. 

Therefore, a systematic expansion of the resolvent up to O(α) valid in the whole 
non-exponential time regime is provided by 

    1 E-L̃_0-Σ̃_a(E)   ≈1 E-L̃_0-Σ̃_a^i +
       
    + 1 E-L̃_a^iδΣ̃_a(E)
    1 E-L̃_a^i  ,

where we have defined

    L̃_a^i = L̃_0 + Σ̃_a^i  ,  Σ̃_a^i = Σ̃_a(z_i)  ,

and

    δΣ̃_a^i(E) = Σ̃_a(E)-Σ̃_a(z_i)  .


To complete the justification of this perturbative expansion (<ref>), we 
finally prove that the order of E-L̃_0-Σ̃_a^i with E=z_i-ix is always 
larger than x∼ 1/t in the regime Ω t≳ 1. To show this
we denote the eigenvalues of L̃_a^i by γ̃_j^i, with i,j=st,0,±. 
The lowest order values are given by the eigenvalues of the real but non-hermitian Liouvillian
L̃_0, which can be diagonalized by the transformation 

    A = A^-1 = 1Ω([  Ωτ_++ϵτ_-   -Δσ_zτ_+; -Δ Zσ_zτ_- -ϵτ_+-Ωτ_- ])  ,

which is the analog of (<ref>) but with Ω_0→Ω and the Z-factor 
in the lower non-diagonal resulting in a non-unitary matrix. In this basis L̃_0 is given by

    A L̃_0 A  = ([    0    0;    0 Ωσ_z ])  ,

i.e. two eigenvalues are zero and two are identical to ±Ω in lowest order in α. 
Σ̃_a^i will shift these eigenvalues by O(αΩ) such that, together with 
the symmetry relations (<ref>-<ref>), we get 

    γ̃_st^i    = 0  ,  γ̃_0^0 = z_0   ,  γ̃_σ^σ = z_σ  ,
    γ̃_0^σ   = -(γ̃_0^-σ)^* = O(αΩ)  ,
    γ̃_-σ^σ   = -(γ̃_σ^-σ)^* 
    = -σΩ + O(αΩ)  ,
    γ̃_σ^0    = -(γ̃_-σ^0)^*
    = σΩ + O(αΩ)  ,

We note that (<ref>) holds exactly in all orders of perturbation theory since
1/(E-L_a(E))=(1/(E-L̃_a(E)))Z' with Z' given by (<ref>) can be viewed as 
the definition of L̃_a(E) and, therefore, the pole positions of the two resolvents 
1/(E-L_a(E)) and 1/(E-L̃_a(E)) must be exactly the same. 

For E=z_i-ix, (<ref>-<ref>) leads to 

    |E-γ̃_j^i| = 
    x
       for j=i
    
    |-ix+z_i|
       for j=st
    
    |-ix±Ω + O(αΩ)|
       for j i,st

i.e. for x∼ 1/t≲Ω to the desired result 

    |E-γ̃_j^i| ≳ x  .


Using the expansion (<ref>) it is now straightforward to write down
the various terms for the time dynamics of ρ_i(t). Denoting the projectors on the eigenstates
of L̃_0+Σ̃_a^i by P̃_j^i, with
i,j=st,0,±, we get
for i=0,±

    ρ_i(t)    = i 2π∫_C_i dE e^-iEt{1 E-γ̃_i^iP̃_i^i  + 
        
    + ∑_j,j'=0,±1 E-γ̃_j^i 1 E-γ̃_j'^iP̃_j^iδΣ̃_a^i(E)P̃_j'^i}·
       · Z'(1+Σ_s1 E)ρ_0  .

Here, we have used for the first term (in the first bracket) on the r.h.s. 
that the other projectors P̃_j^i with j i lead to an
analytic function on the curve C_i with zero integral. Furthermore, due to the matrix 
structure (<ref>) of δΣ̃_a^i(E) and the form (<ref>) of 
P̃_st^i, we get 
P̃_st^iδΣ̃_a^i(E)=δΣ̃_a^i(E)P̃_st^i=0 
and only the terms with j,j'st contribute to the second term (in the first bracket) on the
r.h.s. Furthermore, we note that we can omit all analytic terms
∼ (E-z_i)^2 for δΣ̃_a^i(E) since they lead to analytic contributions on the 
curve C_i in (<ref>) with zero integral. Thus, we  
can use the form (<ref>) for δΣ̃_a^i(E). 
Inserting this form and leaving out all analytic
functions on C_i, we can split ρ_i(t) obviously in pole and pure branch cut contributions

    ρ_i(t) = ρ_i^p(t) + ρ_i^bc(t)  ,

with 

    ρ_i^p(t)    = ρ_i^p1(t) + ρ_i^p2(t) + ρ_i^p3(t)  ,
    ρ_i^p1(t)    = P̃_i^i Z'(1+Σ_s1 z_i)ρ_0  e^-iz_it  ,
    ρ_i^p2(t)    = α∑_j=0,±
    j idF_j dE(z_i)
    P̃_i^i M_j P̃_i^i·
       · Z'(1+Σ_s1 z_i)ρ_0 e^-iz_it  ,
    ρ_i^p3(t)    = αi 2π∫_C_i dE e^-iEt1 E-z_iln-i(E-z_i)Ω·
       ·P̃_i^i M_i P̃_i^i Z'(1+Σ_s1 E)ρ_0

for the pole contributions and

    ρ_i^bc(t)    = αi 2π∫_C_i dE e^-iEt
    (E-z_i)ln-i(E-z_i)Ω·
       ·∑_j,j'=0,±
    (j,j')(i,i)1 (E-γ̃_j^i)(E-γ̃_j'^i)P̃_j^i M_i P̃_j'^i Z'(1+Σ_s1 E)ρ_0

for the pure branch cut contributions. We note that the terms involving Σ_s/z_i are very
important for (<ref>) to calculate the terms in O(1) and O(α) 
consistently since

    Σ_s1 z_0    = -ΔΩΔ̃^2(1-Γ^(2)Γ^(1))
    ([   0   0; τ_+   0 ]) + O(α^2)  ,
    Σ_s1 z_σ    =  iσπαΔΩ([   0   0; τ_+   0 ]) + O(α^2)  ,

where we have used (<ref>) for Σ_s and expanded the pole position z_0=-iΓ
in α by using

    Γ = Γ^(1)+Γ^(2)+O(α^3)  ,  Γ^(1) = παΔ̃^2Ω  ,

where we have taken (<ref>) for Γ^(1). This shows that also second order terms
∼α^2 are needed for the Liouvillian to calculate the pole position z_0 up to second order
needed to get all terms in O(α) for the purely decaying mode of the time evolution. A similiar 
term will also occur for the stationary state (see below). 

In contrast, for the two other pole contributions (<ref>) and (<ref>)
the terms involving Σ_s are only needed for the purely decaying mode i=0 and it is sufficient to
take z_0 up to O(α). For the branch cut contribution (<ref>) the term with Σ_s 
can be left out since it leads to a contribution in O(α^2). Furthermore, in 
(<ref>-<ref>) the projectors P̃_j^i and
the eigenvalues γ̃_j^i can be
replaced by their values in lowest order, all other terms contribute in O(α^2). Only for 
the first pole contribution (<ref>) the projector P̃_i^i is needed up to O(α).
Denoting the projectors in lowest and in first order in α by
P̃^(0)i_j and P̃^(1)i_j, we show in Appendix <ref>
by a straightforward calculation that the projectors transformed with the matrix A 
(see (<ref>)) are given by

    A P̃^(0)i_0 A    = ([ τ_-   0;   0   0 ])  ,
    
    A P̃^(0)i_σ A    = 1 2([      0      0;      0 1+σσ_z ])  ,
    
    A P̃^(1)0_0 A    = iπαΔ̃^2ϵΔΩ^2([      0 1 Zτ_-;    τ_-      0 ])  ,
    
    A P̃^(1)σ_σ A    = 1 4iπσαΔ̃^2Ω^2([       0       0;       0 τ_+-τ_- ])  ,

Taking the projectors in lowest order, the number of terms contributing to 
(<ref>-<ref>) is considerably reduced due to 

    P̃^(0)0_j M_0 P̃^(0)0_j'    0 ⇔ j,j' 0   , 
    P̃^(0)σ_j M_σP̃^(0)σ_j'    0 ⇔ j,j' -σ  ,
    P̃^(0)σ_σ M_-σP̃^(0)σ_σ   = 0   .

As a consequence we get

    ρ_0^p1(t)    = {P̃^(0)0_0 Z'(1+iΣ_s1Γ^(1)
    (1-Γ^(2)Γ^(1)))
       
    + P̃^(1)0_0 Z'}ρ_0  e^-iz_0t  ,
    ρ_0^p2(t)    = α∑_σ=±dF_σ dE(z_0)
    P̃^(0)0_0 M_σP̃^(0)0_0·
       · Z'(1+iΣ_s1Γ^(1))ρ_0 e^-iz_0t  ,
    ρ_0^p3(t)    = 0  ,
    ρ_0^bc(t)    = α∑_σ,σ'=±i 2π∫_C_0 dE e^-iEtln-i(E-z_0)Ω·
       ·(E-z_0) (E-σΩ)(E-σ'Ω)P̃^(0)0_σ M_0 P̃^(0)0_σ' Z'ρ_0  ,

and

    ρ_σ^p1(t)    = {P̃^(0)σ_σ
    Z'(1+σΩΣ_s) + 
       
    + P̃^(1)σ_σ Z'}ρ_0  e^-iz_σ t  ,
    ρ_σ^p2(t)    = αdF_0  dE(z_σ)
    P̃^(0)σ_σ M_0 P̃^(0)σ_σ·
       · Z'(1+iΣ_s1Γ^(1))ρ_0 e^-iz_σ t  ,
    ρ_σ^p3(t)    = αi 2π∫_C_σ dE e^-iEt1 E-z_σln-i(E-z_σ)Ω·
       ·P̃^(0)σ_σ M_σP̃^(0)σ_σ Z'ρ_0  ,
    ρ_σ^bc(t)    = αi 2π∫_C_σ dE e^-iEtln-i(E-z_σ)Ω·
       ·{E-z_σ E^2P̃^(0)σ_0 M_σP̃^(0)σ_0 + 
       
    +1 E(P̃^(0)σ_0 M_σP̃^(0)σ_σ + 
    P̃^(0)σ_σ M_σP̃^(0)σ_0) }
    Z'ρ_0  .


The first term on the r.h.s. of (<ref>) and (<ref>) leads 
to the Bloch-Redfield result modified by the Z-factor. All other contributions to the time
evolution are corrections in O(α). All energy integrals can be calculated from

    i 2π∫ dE e^-iEt1 E-z_iln-i(E-z_i)Ω =
       
    = -(γ + ln(Ω t))e^-iz_i t
       i 2π∫_C_i dE e^-iEt1 E-aln-i(E-z_i)Ω =
       
    = e^-iz_i t∫_0^∞ dy e^-y1 y-i(a-z_i)t =
       
    = e^-iz_i tH((a-z_i)t)  ,
       i 2π∫_C_idE e^-iEtE-z_i (E-a)^2ln-i(E-z_i)Ω =
       
    = e^-iz_i t∫_0^∞ dy e^-yy (y-i(a-z_i)t)^2 =
       
    = e^-iz_i tH̃((a-z_i)t)  ,

where γ is Euler's constant, a z_i, and H(x) and H̃(x) 
can be expressed via the exponential integral E_1(-ix) 

    H(x)    = e^-ixE_1(-ix)  ,
    H̃(x)    = (1-ix)H(x)-1  .

It is important to note that, for the energy integrals occurring in (<ref>) and 
(<ref>), the imaginary part of (a-z_i)t is ∼ -iΓ t and can be neglected
in H((a-z_i)t) and H̃((a-z_i)t) (i.e. leading to higher orders in α) compared to 
the real part of (a-z_i)t which is given by ±Ω t. This
holds even in the case Γ t∼ 1, as can be seen from the integrals (<ref>) and
(<ref>). In contrast, for the exponential function e^-iz_i t it is not possible
to expand in the imaginary part of z_i for Γ t∼ 1. As a consequence, only the crossover
functions H(±Ω t) and H̃(±Ω t) will appear for the branch cut integrals.

Finally, the derivatives of F_i(E) can be obtained from (<ref>)

    d F_i dE(E) = 1+ln-i(E-λ_i(E))Ω+O(α)  ,

which gives

    d F_0 dE(z_σ)    = 1-iσπ 2 + O(α)  ,
    d F_σ dE(z_0)    = 1+iσπ 2 + O(α)  .


Using all these relationships together with the form of the various matrices, one can
straightforwardly evaluate (<ref>-<ref>) and calculate the expectation
values of the Pauli matrices. Decomposing the time dynamics according to 
(<ref>) in the various modes, we get for the preexponential functions the
following final result for the time dynamics in the non-exponential time regime 


    F_x^0(t)     = 
    -⟨σ_x⟩_st - (1+2αΔ̃^2Ω^2)Δ̃^2ΔΩ⟨σ̃_z⟩_0
     + αΔ̃^4ϵΔ^2Ω^3{(H_t^' - H̃_t^')⟨σ̃_x⟩_0 + 
    (π + H̃_t^'')⟨σ_y⟩_0}  ,
    
    
    F_y^0(t)     = παΔ̃^2ϵΔΩ^2(1+⟨σ̃_z⟩_0)
     + αΔ̃^4Δ^2Ω^2{H̃_t^''⟨σ̃_x⟩_0+
    (H_t^' + H̃_t^')⟨σ_y⟩_0}  ,
    
    
    F_z^0(t)     = 
    -⟨σ_z⟩_st + (1+2αΔ̃^2Ω^2)ϵΩ⟨σ̃_z⟩_0
     - παΔ̃^2ϵ^2ΔΩ^3⟨σ_y⟩_0
     + αΔ̃^4ΔΩ^3{(H_t^' - H̃_t^')⟨σ̃_x⟩_0 + 
    H̃_t^''⟨σ_y⟩_0}  ,
    
    
    F_x^c(t)     = 
    -f_tΔ̃^2ϵΔ^2Ω⟨σ̃_x⟩_0
     - αΔ̃^4ϵΔ^2Ω^3{2 H_t^'⟨σ̃_x⟩_0 - 
    (π 2 + 2 H_t^'')⟨σ_y⟩_0} -  2αΔ̃^2ΔΩ^3
    (Δ̃^2 H̃_t^' + ϵ^2 H_t^')⟨σ̃_z⟩_0 
      ,
    
    
    F_y^c(t)     = 
    f_tΔ̃^2Δ^2⟨σ_y⟩_0
     - αΔ̃^2ϵΔΩ^2(π + 2 H_t^'⟨σ̃_z⟩_0)
     - π 2αΔ̃^4Δ^2Ω^2⟨σ̃_x⟩_0
      ,
    
    F_z^c(t)     = 
    -f_tΔ̃^2ΔΩ⟨σ̃_x⟩_0
     - π 2αΔ̃^4ΔΩ^3⟨σ_y⟩_0
     +  2αΔ̃^2ϵ^2ΔΩ^3
    ( H_t^'⟨σ̃_x⟩_0 + H_t^''⟨σ_y⟩_0)
     -  2αΔ̃^2ϵΩ^3
    (H_t^'-H̃_t^')⟨σ̃_z⟩_0
      ,
    
    F_x^s(t)     = 
    -f_tΔ̃^2ϵΔ^2Ω⟨σ_y⟩_0
     + παΔ̃^2ϵ^2ΔΩ^3 +  2αΔ̃^4ϵΔ^2Ω^3
    ( H_t^''⟨σ̃_x⟩_0- H_t^'⟨σ_y⟩_0)
     +  2αΔ̃^2ΔΩ^3
    (ϵ^2 H_t^''+Δ̃^2H̃_t^'')⟨σ̃_z⟩_0
      ,
    
    F_y^s(t)     = 
    -f_tΔ̃^2Δ^2⟨σ̃_x⟩_0
     - παΔ̃^4Δ^2Ω^2⟨σ_y⟩_0
     -  2αΔ̃^2ϵΔΩ^2 H_t^'⟨σ̃_z⟩_0
      ,
    
    F_z^s(t)     = 
    -f_tΔ̃^2ΔΩ⟨σ_y⟩_0
     + αΔ̃^2ϵΩ^3{π + 2(H_t^''-H̃_t^'')⟨σ̃_z⟩_0} -  2αΔ̃^2ϵ^2ΔΩ^3
    (H_t^''⟨σ̃_x⟩_0 - H_t^'⟨σ_y⟩_0)
      ,


where we have defined the quantities

    H_t^'    = ReH(Ω t) = 1 2∑_σ=±H(σΩ t)  ,
    
    H_t^''    = ImH(Ω t) = -i 2∑_σ=±σ H(σΩ t)  ,
    H̃_t^'    = ReH̃(Ω t) = 1 2∑_σ=±H̃(σΩ t)  ,
    H̃_t^''    = ImH̃(Ω t)
     = -i 2∑_σ=±σH̃(σΩ t)  ,
    σ̃_x    = 
    -1Ω(ϵσ_x + Δσ_z)  ,
    σ̃_z    = 1Ω(ϵσ_z - 
    Δ̃^2Δσ_x)  ,
    
    f_t    = 1 + αΔ̃^2Ω^2 
    - 2α(γ + ln(Ω t))ϵ^2Ω^2  .

We note that the operators σ̃_x and σ̃_z can not be interpreted as
the Pauli spin operators in the basis where the local Hamiltonian with
Δ→Δ̃ is diagonal since both Δ and Δ̃ appear
in the definition in a subtle way. Only if the renormalization of the tunneling is neglected, 
these operators are identical to the Pauli spin operators defined in 
(<ref>-<ref>).

The stationary values ⟨σ_α⟩_st of the Pauli matrices follow from

    ⟨σ_x⟩_st    = Δ̃^2ΔΩ(1-Γ^(2)Γ^(1)) + 2αΔ̃^4ΔΩ^3  ,
    ⟨σ_y⟩_st    = 0  ,
    ⟨σ_z⟩_st    = -ϵΩ(1-Γ^(2)Γ^(1)) - 2αΔ̃^2ϵΩ^3  .

This can be obtained from (<ref>) via the spectral decomposition of 
L̃_a^st=L̃_0+Σ̃_a(0). Denoting the eigenvalues and projectors of this
Liouvillian by γ̃_j^st and P̃_j^st, with j=st,0,±, we show in 
Appendix <ref> that we get in analogy to 
(<ref>-<ref>) and (<ref>-<ref>)

    γ̃_st^st   = 0  ,
    γ̃_0^st   = 
    -i(Γ^(1)(1-2αΔ̃^2Ω^2)+Γ^(2))+O(α^3)  ,
    γ̃_σ^st   = σΩ + O(αΩ)  ,

and 

    A P̃^st_st A    = ([ τ_+   0;   0   0 ])  ,
    
    A P̃^(0)st_0 A    = ([ τ_-   0;   0   0 ])  ,
    
    A P̃^(0)st_σ A    = 1 2([      0      0;      0 1+σσ_z ])  ,
    
    A P̃^(1)st_0 A    = iπαΔ̃^2ϵΔΩ^2([      0 1 Zτ_-;    τ_-      0 ])  .

Inserting the spectral decomposition in (<ref>) we get up to O(α)

    ρ_st    = 1 2([ 1; 1; 0; 0 ])
    -iΓ^(1) 2Δ{Ωγ̃_0^stP̃_0^(0)st+
       
    +iΩΓ^(1)P̃_0^(1)st
    +∑_σ=±σP̃_σ^(0)st}([ 0; 0; 1; 1 ])  .

Inserting (<ref>) and (<ref>-<ref>) we find that the sum of the
last two terms on the r.h.s. is zero and we get for the stationary density matrix up to O(α)
the final result

    ρ_st    = 1 2([ 1; 1; 0; 0 ])-
    (1+2αΔ̃^2Ω^2-Γ^(2)Γ^(1))1 2Ω([       ϵ;      -ϵ; -Δ̃^2/Δ; -Δ̃^2/Δ ])  ,

which leads to the result (<ref>-<ref>) for the stationary values of
the Pauli matrices.

To calculate the ratio Γ^(2)/Γ^(1) we need an analysis of all second order terms 
∼α^2 of the Liouvillian L(0^+) to get the stationary state up to O(α). This
goes beyond the scope of this paper. However, in Ref. divincenzo_loss_05, such an
analysis has been performed in bare perturbation theory (i.e. using the unrenormalized tunneling) 
with the result (note that we slightly changed the result such that it is valid for a Lorentzian
cutoff function in the bath)

    ⟨σ_x⟩_st    = ΔΩ_0+αΔ^3Ω^3
    +αΔΩ_0^3(Δ^2+2ϵ^2)lnΩ_0 D  .

and it was shown that this agrees with the result from the partition function proving the Ergoden
hypothesis up to O(α). This result is consistent with (<ref>) if we take

    Γ^(2)Γ^(1) = αΔ̃^2Ω^2  ,

such that our final result for the stationary values reads

    ⟨σ_x⟩_st    = Δ̃^2ΔΩ
    +αΔ̃^4ΔΩ^3  ,
    ⟨σ_y⟩_st    = 0  ,
    ⟨σ_z⟩_st    = -ϵΩ
    -αΔ̃^2ϵΩ^3  .


We note that the terms involving Γ^(2) cancel out for the full time dynamics of ρ(t)
in the limit Γ t≪ 1, where the exponential e^-Γ t≈ 1. This is a generic feature
since, in this time regime, |E-L_0|≫Γ, and bare perturbation theory can be used to expand
the resolvent 1/(E-L_0-Σ(E)) in Σ(E), without any need of the Liouvillian up to second
order in α to calculate all terms of the time dynamics up to O(α). Therefore, it is of
no surprise that the time-dependent terms involving Γ^(2) can be related to corresponding terms
of the stationary state.

We now discuss our central result (<ref>-<ref>) and compare it with 
the literature. The leading order term is consistent with the Bloch-Redfield solution 
(<ref>-<ref>), provided one neglects the renormalization of the tunneling. 
Our result shows that the renormalized tunneling appears in a subtle way which can not be
obtained by just replacing Δ→Δ̃. There is a Z-factor renormalization
Z=Δ̃^2/Δ^2 for F^c_x,y and F^s_x,y, and terms ∼Δ or ∼Δ^3 
in the Bloch-Redfield solution are replaced by √(Z)Δ=Δ̃^2/Δ and 
Z^2Δ^3=Δ̃^4/Δ, respectively.

The most interesting correction in O(α) is the slowly varying logarithmic term 
αln(Ω t) appearing in the function f_t multiplying the leading order terms of 
F_α^c/s. We note that the correct energy scale in this logarithmic term is the renormalized
Rabi frequency Ω and not the Lamb-shift Ω-Ω_0 as it was 
obtained in Ref. divincenzo_loss_05. As was already mentioned in 
Section <ref> via Eq. (<ref>) the crucial point is not to neglect the 
O(α) contributions in the logarithmic functions. E.g. if one considers the integral 
(<ref>) for z_i=z_+=Ω-iΓ/2 and neglects all O(α) 
contributions in the argument of the logarithm by setting 
ln(-i(E-z_+)/Ω)≈ln(-i(E-Ω_0)/Ω_0) one obtains 

    i 2π∫ dE e^-iEt1 E-z_+ln-i(E-Ω_0)Ω_0 =
        = ln-i(z_+-Ω_0)Ω_0e^-iz_+ t
    + H((z_+-Ω_0)t)e^-iΩ_0 t  ,

which is obviously quite different from the exact result -(γ + ln(Ω t))e^-iz_+ t 
not only because of the incorrect exponential appearing in the second term on the r.h.s. (which
is just oscillating with the unrenormalized Rabi frequency) but also due to the incorrect
preexponential functions of both terms involving the energy scale of the Lamb shift 
δΩ=Ω-Ω_0. This shows that the resummation of secular terms
contained in logarithmic contributions of the Liouvillian is not only 
important to get the correct exponential part of the time dynamics but
also to obtain the correct preexponential functions. Only in the limit Γ t,δΩ t≪ 1,
where |E-Ω_0|≫δΩ,Γ, it is allowed to neglect secular terms by disregarding
the O(α) terms in the argument of the logarithm. In this case one can use the approximation 
H((z_+-Ω_0)t)≈ -γ-ln(-i(z_+-Ω_0)t) and 
e^-iz_+t≈ e^-iΩ_0t in (<ref>) leading to

    i 2π∫ dE e^-iEt1 E-z_+ln-i(E-Ω_0)Ω_0 ≈
        = -(γ + ln(Ω_0 t))e^-iΩ_0 t  ,

with the correct logarithmic time dependence involving the Rabi frequency and not the
Lamb shift.

For large times Ω t∼ 1/α≫ 1, where the damping is still moderate due to 
Γ t∼ O(1), the logarithmic term ∼αln(Ω t) is
the most important correction to the leading order terms. In this regime, the functions 
H_t and H̃_t lead only to very small contributions and fall off according to

    H^'_t    = 1 (Ω t)^2 + O(1(Ω t)^4)  ,
    
    H^''_t    = 1Ω t + O(1(Ω t)^3)  ,
    H̃^'_t    = -1 (Ω t)^2 + O(1(Ω t)^4)  ,
    H̃^''_t    = O(1(Ω t)^3)  .

The pure branch cut contributions arising from H_t and H̃_t are the only terms 
showing a significant time dependence whereas the logarithmic terms are slowly varying in time. The most 
important term is the one arising from H^''_t which falls off only ∼ 1/(Ω t). It arises
only in the finite bias case for the modes F^c_x/z and F^s_x/z and has never been reported before.
The standard case treated in the literature <cit.> 
is the calculation for the time dynamics of the Pauli matrix in z-direction at zero bias 
for the initial condition ⟨σ_z⟩_0=1 and ⟨σ_x/y⟩_0=0. 
In this case and for Ω t≫ 1 our solution reduces to

    ⟨σ_z⟩(t) ≈ 
    (1+α)cos(Δ̃t)e^-Γ 2t 
    - 2α1 (Δ̃ t)^2e^-Γ t  ,

which, up to the missing exponential for the second term on the r.h.s., agrees with the 
NIBA result <cit.> and the result obtained from the Born approximation
<cit.> (where also the residuum has been calculated for the first term on
the r.h.s.). In Refs. slutskin_etal_11,kashuba_schoeller_13 the correct exponential 
has been obtained for the second term. The important new result for finite bias is that, besides the 
appearance of many other terms falling off ∼α/(Ω t)^2, there are new terms falling
off ∼α/(Ω t). For Ω t∼ 1/α, this are terms in O(α^2) and 
thus of the same order as other constant terms ∼α^2 or slowly varying logarithmic terms
∼α^2ln^2(Ω t) not covered by our analytic solution in the non-exponential regime. 
However, the terms ∼α/(Ω t) 
are consistent in the sense that they determine the leading behavior of those contributions which
show a significant time dependence. In contrast, terms ∼α/(Ω t)^2 are inconsistent 
in this sense, since for finite bias there will be other strongly varying terms ∼α^2/(Ω t)
of the same order which we have not calculated. Keeping only the consistent terms falling off 
∼α/(Ω t) we obtain for large times Ω t≫ 1:


    F_x^0(t)     = 
    -⟨σ_x⟩_st - (1+2αΔ̃^2Ω^2)Δ̃^2ΔΩ⟨σ̃_z⟩_0
     + παΔ̃^4ϵΔ^2Ω^3⟨σ_y⟩_0
      ,
    
    
    F_y^0(t)     = παΔ̃^2ϵΔΩ^2(1+⟨σ̃_z⟩_0)
      ,
    
    
    F_z^0(t)     = 
    -⟨σ_z⟩_st + (1+2αΔ̃^2Ω^2)ϵΩ⟨σ̃_z⟩_0
     - παΔ̃^2ϵ^2ΔΩ^3⟨σ_y⟩_0
      ,
    
    
    F_x^c(t)     = 
    -f_tΔ̃^2ϵΔ^2Ω⟨σ̃_x⟩_0
     + π 2αΔ̃^4ϵΔ^2Ω^3⟨σ_y⟩_0
     + 2αΔ̃^4ϵΔ^2Ω^31Ω t⟨σ_y⟩_0
      ,
    
    
    F_y^c(t)     = 
    f_tΔ̃^2Δ^2⟨σ_y⟩_0
     - παΔ̃^2ϵΔΩ^2 - π 2αΔ̃^4Δ^2Ω^2⟨σ̃_x⟩_0
      ,
    
    F_z^c(t)     = 
    -f_tΔ̃^2ΔΩ⟨σ̃_x⟩_0
     - π 2αΔ̃^4ΔΩ^3⟨σ_y⟩_0
     +  2αΔ̃^2ϵ^2ΔΩ^31Ω t⟨σ_y⟩_0
      ,
    
    F_x^s(t)     = 
    -f_tΔ̃^2ϵΔ^2Ω⟨σ_y⟩_0
     + παΔ̃^2ϵ^2ΔΩ^3 +  2αΔ̃^4ϵΔ^2Ω^31Ω t⟨σ̃_x⟩_0
     +  2αΔ̃^2ϵ^2ΔΩ^31Ω t⟨σ̃_z⟩_0
      ,
    
    F_y^s(t)     = 
    -f_tΔ̃^2Δ^2⟨σ̃_x⟩_0
     - παΔ̃^4Δ^2Ω^2⟨σ_y⟩_0
      ,
    
    F_z^s(t)     = 
    -f_tΔ̃^2ΔΩ⟨σ_y⟩_0
     + παΔ̃^2ϵΩ^3 +  2αΔ̃^2ϵΩ^31Ω t⟨σ̃_z⟩_0
     -  2αΔ̃^2ϵ^2ΔΩ^31Ω t⟨σ̃_x⟩_0
      .



For zero bias ϵ=0 and large times Ω t≫ 1, we keep the
leading terms falling off ∼α/(Ω t)^2 and obtain with the help of f_t=1+α, 
⟨σ_x⟩_st=(1+α)Δ̃Δ,
⟨σ_y⟩_st=⟨σ_z⟩_st=0, 
⟨σ̃_x⟩_0=-ΔΔ̃⟨σ_z⟩_0 and 
⟨σ̃_z⟩_0=-Δ̃Δ⟨σ_x⟩_0 the result


    F_x^0(t)     = -(1+α)Δ̃Δ +  (1+2α)Δ̃^2Δ^2⟨σ_x⟩_0   ,  
    
    F_y^0(t)  = 0   ,  
    
    F_z^0(t)  = -2α1 (Δ̃ t)^2⟨σ_z⟩_0
      ,
    
    
    F_x^c(t)     = -2α1 (Δ t)^2⟨σ_x⟩_0
      ,  
    
    F_y^c(t)  = 
    (1+α)Δ̃^2Δ^2⟨σ_y⟩_0
     + π 2αΔ̃Δ⟨σ_z⟩_0
      ,  
    F_z^c(t)  = 
    (1+α)⟨σ_z⟩_0
     - π 2αΔ̃Δ⟨σ_y⟩_0
      ,
    
    F_x^s(t)     = 0   ,  
    F_y^s(t)  = 
    (1+α)ΔΔ̃⟨σ_z⟩_0
     - παΔ̃^2Δ^2⟨σ_y⟩_0
      ,  
    F_z^s(t)  = 
    -(1+α)Δ̃Δ⟨σ_y⟩_0
      ,


which agrees with the result obtained in Ref. kashuba_schoeller_13, except that we
have also calculated all time-independent corrections for the preexponential functions in O(α) here. 

For very large times in the 
exponential region where αln(Ω t)∼ O(1), our result is no longer valid and the
RG treatment presented in Section <ref> is needed to sum
up all powers of such logarithmic terms (determining the power law exponent in 
leading order in α). As we discuss later on, the main result is that the 
function f_t has to be replaced by the power-law

    f_t → (1Ω t)^2αϵ^2Ω^2(1-2αγϵ^2Ω^2+αΔ̃^2Ω^2)  .

with a power-law exponent depending on the bias. This power-law exponent is consistent with the one predicted 
in Ref. slutskin_etal_11 where Ω was replaced by the unrenormalized 
Rabi frequency Ω_0. However, in this reference, many terms in higher order in
Δ/Ω_0 and α have been neglected and a consistent
RG analysis was lacking whether additional logarithmic terms appear which can
change the power-law exponent e.g. from 2αϵ^2/Ω^2 to 2α. It turns out
that this analysis depends crucially on the time regime under consideration. Whereas, for exponentially
large times, it turns out that the power-law exponent is indeed 2αϵ^2/Ω^2 for the
oscillating modes, a completely different result appears for exponentially small times with
a power-law exponent given by 2α, see (<ref>) and 
(<ref>). There is a complicated crossover between these two power laws
since the real part of the functions H_t^' and H̃_t^' contain additional 
logarithmic terms for small times Ω t≪ 1, see Eqs. (<ref>) and 
(<ref>) below. Only via our consistent RG treatment presented in 
Section <ref> one can be sure to include all terms of the leading logarithmic series
providing the correct power-law exponents in O(α) for exponentially small and large times,
together with the correct crossover behavior in the non-exponential regime. 

One can check that in the limit of small but not exponentially small times our solution 
(<ref>-<ref>) is consistent with 
(<ref>-<ref>). The logarithmic terms
are a result of a combination of logarithmic terms arising from the terms 
∼αln(Ω t) appearing explicitly in (<ref>-<ref>) 
and those arising from the functions H_t and H̃_t, which, for small argument, 
can be expanded as

    H^'_t    = -γ - ln(Ω t) + O(Ω t)  ,
    
    H^''_t    = π 2 + O(Ω t)  ,
    H̃^'_t    = -γ - ln(Ω t) - 1 + O(Ω t)  ,
    H̃^''_t    = π 2 + O(Ω t)  .

Inserting this expansion in (<ref>-<ref>) and neglecting all terms
∼αΩ t (with or without a logarithm) we obtain

    ⟨σ_x/y⟩_st + F^0_x/y + F^c_x/y ≈ 
       ≈ Δ̃^2Δ^2{1-2α(γ + ln(Ω t))}⟨σ_x/y⟩_0   ,
       ⟨σ_z⟩_st + F^0_z + F^c_z  ≈ ⟨σ_z⟩_0   ,
       F^s_x Ω t  ≈  -Δ̃^2Δ^2ϵ t
    ⟨σ_y⟩_0   ,
       F^s_y Ω t  ≈ Δ̃^2Δ^2{ϵ t ⟨σ_x⟩_0 + Δ t ⟨σ_z⟩_0 }  ,
       F^s_z Ω t  ≈  -Δ̃^2Δ^2Δ t ⟨σ_y⟩_0  .

Inserting this result in (<ref>), expanding the 
exponential functions up to linear order in Ω t and again neglecting all terms
∼αΩ t, we obtain precisely the 
expansion (<ref>-<ref>) for
small but not exponentially small times, showing that we cover the correct crossover behavior by combining the 
solutions (<ref>-<ref>) for small or exponentially
small times with (<ref>-<ref>) in the non-exponential regime. 

For moderate times Ω t∼ O(1) the logarithmic terms are of the same order as all other
corrections in O(α). In this regime our full solution (<ref>-<ref>)
is needed to calculate all terms one order beyond Bloch-Redfield. In this case the 
time dependence of the preexponential functions is governed by a complicated combination of slowly 
varying logarithmic terms and terms arising from the functions H_t and H̃_t containing
the exponential integral via (<ref>) and (<ref>). 

Finally, we note that our solution in the non-exponential regime at zero tunneling Δ=0 is 
fully consistent with the exact solution at zero tunneling presented in 
(<ref>-<ref>) 
and (<ref>). It is even fully reproduced if we use the 
replacement (<ref>) for the function f_t. 
 




§ REAL-TIME RENORMALIZATION GROUP



In this section we will present the real-time renormalization group approach to calculate
the Liouvillian L(E) beyond perturbation theory by including the leading logarithmic series at
low and high energies. This provides the basis for the renormalized perturbation theory in the
non-exponential regime together with the calculation of the time dynamics for exponentially small or 
large times, see Sections <ref> and <ref>.



 §.§ RG equations



The leading order RG equations to determine the Liouvillian L(E) for the ohmic spin boson model 
have been derived in Ref. kashuba_schoeller_13. Using the definitions (<ref>), 
(<ref>) and (<ref>), they read


    d dEL̃_Δ(E)    = 2α∑_i Z'(E) G(E) P_i(E) Z'(E) G(E) 
    L̃_Δ(E) - λ_i(E)  E - λ_i(E)  ,
    d dEZ'(E)    = 2α∑_i Z'(E) G(E) P_i(E) Z'(E) G(E) Z'(E)
    1  E - λ_i(E)  ,
    d dEG(E)    = 2α∑_ij G(E) P_i(E) Z'(E) G(E) P_j(E) Z'(E) G(E) 
    L_i(E) - L_j(E)  λ_i(E) - λ_j(E)  ,


with L_i(E)=ln((-i)(E-λ_i(E))/Ω), see (<ref>). For i=j, we get
(L_i(E) - L_j(E))/(λ_i(E) - λ_j(E)) = -1/(E-λ_i(E)) for 
the last factor on the r.h.s. of (<ref>).
Here, λ_i(E) and P_i(E) are the eigenvalues and projectors of L̃_Δ(E),
respectively, defined in (<ref>-<ref>). 
The RG equations for L̃_Δ(E) and Z'(E) are coupled to the RG equation for the 
vertex G(E), which is also a 4× 4-matrix. They can be solved along a certain path 
in the complex plane with the following initial condition at E=iD

    L̃_Δ(iD)    = L_0   ,
    
    Z'(iD)    = 1  ,
    
    G(iD)    = ([   0   0;   0 σ_z ])  .
 
In the limit of large D, this initial condition can also be taken for E=ω-iD, with 
some real ω≪ D. Fixing the real parameter ω, the RG equations are solved numerically 
along the path E=ω-iΛ, starting at Λ=D. For ω=0 and Λ→ 0, 
we obtain the Liouvillian L(0^+) from which we can obtain the stationary state ρ_st via 
(<ref>). For ω=±η (η=0^+) and Λ < -Γ, we obtain a jump
between L(-iΛ+η) and L(-iΛ-η) indicating the branch cut starting at
z_0=-iΓ. Similarly, for ω=σΩ±η and Λ < -Γ/2, there will
be jump between the two solutions L(σΩ-iΛ+η) and L(σΩ-iΛ-η),
corresponding to the two branch cuts starting at z_σ=σΩ-iΓ/2. In this way one
can determine numerically the positions of the branching points z_i (i=0,±) of the Liouvillian, 
together with fixing the branch cuts along the direction of the negative imaginary axis. This is an 
important advantage of the real-time RG method since it uses the complex Fourier variable E as flow
parameter and, via a numerical solution of the RG equations along a certain path in the complex plane, 
allows for an elegant analytical continuation of retarded quantities into the lower half of the
complex plane. Once the RG equations have been solved in this way the time dynamics 
can be calculated from (<ref>) and (<ref>) 
with the preexponential operator given in analogy to (<ref>) by 

    F_i(t)    = 1 2π∫_0^∞ dx  e^-xt·
       ·{1 E-L̃_Δ(E)Z'(E)|_E=z_i-ix+η - (η→ -η)}·
       · (1+Σ_s1 z_i-ix)ρ_0  .

Due to the exponential part e^-xt in the integrand, the calculation of the integrals is numerically 
very stable, which is the basic reason why the direction of the branch cuts is chosen along 
the negative imaginary axis.

As explained in detail in Refs. kashuba_schoeller_13,schoeller_14, the RG equations
(<ref>-<ref>) contain the leading logarithmic series, i.e. power-law exponents
can be calculated reliably up to O(α). All other terms on the r.h.s. of the RG equations are
of higher order in α but scale as function of E in the same way as the leading term, i.e. for
large energies as 1/E and for small energies at most as 1/(E-λ_i(E)). Furthermore, the r.h.s. of
the RG equations is universal, i.e. well-defined in the limit D→∞, the band width D of 
the bath enters only via the initial value E=iD. Therefore, the RG equations provide a well-defined set
of differential equations which can be systematically truncated at order α, including secular terms 
and the leading logarithmic series to all orders.  

We note that the vertex renormalization is very essential. At large energies we show in 
Section <ref> that the vertex obtains an important Z-factor renormalization, which 
is important to determine the correct time dynamics at exponentially small and intermediate times in the
non-exponential regime. For energies close to one of the branching points z_i, i=st,±, the
vertex renormalization is very different, leading to completely different power-laws for the 
preexponential functions of the time dynamics compared to the one at exponentially small times. This
has also been discussed in Ref. kashuba_schoeller_13 for the zero-bias case, where it was
shown that G(E) does no longer renormalize for E close to z_i. At finite bias, this is quite 
different and leads to the power-law (<ref>), as will be demonstrated in 
Section <ref> by a numerical solution of the RG equations. 

Furthermore, we note that the vertex renormalization has to be taken with care when solving the
RG equations along a certain branch cut. For i j, the difference L_i(E)-L_j(E) 
can contain a constant imaginary term due to the jump of the various logarithm across the branch cuts. Although
this constant term does not lead to a logarithmic contribution and contributes inconsistent terms in
O(α^2) for the Liouvillian, it leads to a term ∼α |E-z_i|/Ω for the vertex, 
which becomes of O(1) for |E-z_i|∼Ω/α. This effect is inconsistent and is canceled 
by other contributions from higher orders. However, this effect can easily be avoided by just omitting
this constant imaginary term for the vertex renormalization when E=z_i-ix±η with x>Ω. 
This is consistent with the analytic solution of the RG equations up to O(α) for large energies and 
energies in the non-exponential regime, as shown in Sections <ref> and <ref>. 
  




 §.§ Large energies



For large energies, where Ω≪ |E| ≪ D, we will show in this section that L̃_Δ(E) and
Z'(E) are given by (<ref>-<ref>). For |E|≫|λ_i(E)|,
the RG equations (<ref>-<ref>) can be approximated by


    d dEL̃_Δ(E)    = 2α∑_i Z'(E) G(E) P_i(E) Z'(E) G(E) 
    L̃_Δ(E) - λ_i(E)  E
        = 2α{Z'(E) G(E) Z'(E) G(E) L̃_Δ(E) -  Z'(E) G(E) L̃_Δ(E) Z'(E) G(E)}1 E  ,
    d dEZ'(E)    = 2α Z'(E) G(E) Z'(E) G(E) Z'(E) 1  E   ,
    d dEG(E)    = -2α G(E) Z'(E) G(E) Z'(E) G(E) 1  E  .


With the initial conditions (<ref>-<ref>) these RG equations can be solved 
by the ansatz 

    L̃_Δ(E)    = ([          0       Δτ_-; Δ Z(E) τ_-       ϵσ_z ])
      ,
    
    Z'(E)    = ([    1    0;    0 Z(E) ])
      ,
    
    G(E)    = g(E)
    ([   0   0;   0 σ_z ])
      .

Inserting this ansatz into the r.h.s. of the RG equations we find

    d dEL̃_Δ(E)    = 2α Z(E)^2 g(E)^2 1 E([          0          0; Δ Z(E) τ_-          0 ])
      ,
    d dEZ(E)    = 2α Z(E)^3 g(E)^2 1 E  ,
    d dEg(E)    = -2α Z(E)^2 g(E)^3 1 E  ,

showing that (<ref>) is consistent with the ansatz 
(<ref>). Furthermore we find

    d dE Z(E)^2 g(E)^2    = 0  ,
    d dE Z(E)    = 2α Z(E) 1 E  ,

with the solution

    Z(E)^2 g(E)^2    =  1  ,
    
    Z(E)    = (-iE D)^2α  .

In conclusion, (<ref>), (<ref>) and 
(<ref>) prove the form (<ref>-<ref>) of 
the Liouvillian at large energies which was used in Section <ref> to calculate the
time dynamics for small times. 





 §.§ The non-exponential regime



In the non-exponential regime (<ref>), where αln(-i(E-λ_i(E))/Ω)≪ 1,
the RG equations can be solved perturbatively around the solution 
(<ref>-<ref>) at high energies evaluated at 
E=iΩ, see Ref. schoeller_09 for details. Denoting the latter by

    L̃_0    = ([       0    Δτ_-; Δ Z τ_-    ϵσ_z ])
      ,
    
    Z^'    = ([ 1 0; 0 Z ])
      ,   Z = (Ω D)^2α = Δ̃^2Δ^2  ,
    
    G    = g
    ([   0   0;   0 σ_z ])  ,   g = 1/Z
      ,

the solution of the RG equations in the non-exponential regime can be written up to O(α) as

    L̃_Δ(E)    = L̃_0  +
       
     + 2α∑_i Z' G P_i^(0) Z' G 
    (L̃_0-λ_i(E))L_i(E)  ,
    
    Z'(E)    = Z'  + 2α∑_i Z' G P_i^(0) Z' G Z'L_i(E)  ,

with L_i(E)=ln(-i(E-λ_i(E))/Ω) defined in
(<ref>). Our convention for the notation of Z' and G is chosen
such that when no argument E is written, we implicitly take the high-energy solution 
evaluated at E=iΩ, given by (<ref>-<ref>). 
P_i^(0) denote the projectors of L̃_0 in lowest order in α, 
which are given by (<ref>-<ref>)

    A P^(0)_0 A    = ([ τ_-   0;   0   0 ])  ,
    
    A P^(0)_σ A    = 1 2([      0      0;      0 1+σσ_z ])  ,

where the matrix A is defined in (<ref>). For large |E|≫|λ_i(E)| but not 
exponentially large the solution (<ref>-<ref>)
is consistent with the result at large energies, given by (<ref>), 
(<ref>) and (<ref>), when expanded in αln(-iE/Ω).
As a consequence it is straightforward to see that (<ref>-<ref>) is
indeed the solution of the RG equations in the non-exponential regime up to O(α) since the
differential equation is fulfilled and the boundary condition at large energies is reproduced.

Using (<ref>) and (<ref>), we find from the solution 
(<ref>-<ref>) the following result for the Liouvillian L_a(E) up to
O(α)

    L_a(E)    = E  - 1 Z'(E-L̃_0+Σ̃_a(E))
    Σ̃_a(E)    = 2α∑_i Z' G P_i^(0) Z' G F_i(E)  ,

with F_i(E)=(E-λ_i(E))L_i(E) defined in (<ref>). Inserting the algebra
of the various matrices by using (<ref>-<ref>) and 
(<ref>-<ref>), we obtain (<ref>) for the 
propagator together with (<ref>) for Σ̃_a(E). 

To check that (<ref>-<ref>) are indeed the
eigenvalues of L̃_Δ(E), we transform (<ref>) with the matrix A 
and replace λ_i(E) by the eigenvalues of L̃_0. After a straightforward calculation 
we obtain

    AL̃_Δ(E)A   = 
    Ω([   0   0;   0 σ_z ])+
    2αΔ̃^2ΩL_0(E) 
    ([      0      0;      0 τ_-σ_z ])
       
    -αΔ̃^2Ω∑_σσL_σ(E)
    ([           τ_-             0; ϵΔ (σ_z+σ)τ_-             0 ]) .

As a consequence the four eigenvalues are given up to O(α) by

    λ_st(E)    = 0  ,
    λ_0(E)    = -αΔ̃^2Ω∑_σσL_σ(E)  ,
    λ_σ(E)    = σ(Ω+αΔ̃^2ΩL_0(E))  ,

in agreement with (<ref>-<ref>).



 §.§ Exponentially large times



For exponentially large times, where higher powers in αln(Ω t) become significant and can
no longer be treated in lowest order to analyze the corrections to Bloch-Redfield, we need a solution 
of the RG equations exponentially close to the branching points z_i. Analytically, such an analysis 
is very complicated for arbitrary bias but can be done at zero bias, see Ref. kashuba_schoeller_13.
For arbitrary bias, we have studied the numerical solution of the RG equations and will present a fit to an
analytical ansatz in this section. 

The case of zero bias ϵ=0 has been studied in Ref. kashuba_schoeller_13 
by using the real-time RG method. The main result was that the result 
(<ref>-<ref>) for large times still holds for 
exponentially large times, except for F_x^c(t), which obtains an additional function s_0(t)

    F_x^c(t) = -2αs_0(t) (Δ t)^2⟨σ_x⟩_0  ,

with

    s_0(t) = (1 (1+αln(Ω t))(1-ln(1+αln(Ω t))))
      ,

such that the complete solution for ⟨σ_x⟩(t) reads

    ⟨σ_x⟩(t)    = ⟨σ_x⟩_st(1-e^-Γ t)
     + (1+2α)Δ̃^2Δ^2e^-Γ t⟨σ_x⟩_0
       
    - 2αs_0(t) (Δ t)^2cos(Ω t) e^-Γ t/2⟨σ_x⟩_0  ,

with ⟨σ_x⟩_st=(1+α)Δ̃^2Δ^2. However,
this result is not very important since, at zero bias, the importance of higher orders in
αln(Ω t) for the preexponential function shows only up for Γ t≫ 1, where the
exponential damping leads to a negligible result for the time dynamics. Only for α∼ 1,
the estimation in (<ref>) shows that higher powers of logarithmic terms are
important for times where the damping is moderate. Only from an academic point of view, where the
preexponential function can be studied separately, exponentially large times are also interesting
at zero bias and the function s_0(t) can be identified. As already discussed in 
Ref. kashuba_schoeller_13, we note that there is no change of the power-law exponent of
the 1/t^2 parts, in particular for the time dynamics of ⟨σ_z⟩(t), see
(<ref>), in contrast to the NIBA solution which predicts an incorrect 
power-law exponent 2-2α <cit.>.

As discussed in detail via the estimation (<ref>), the importance of higher orders 
in αln(Ω t) changes significantly for large bias.
At arbitrary bias, we have checked numerically that a power-law appears for the leading-order 
term of the oscillating modes in the regime of very large times, with a bias-dependent exponent 
2αϵ^2/Ω^2.  
E.g. Fig. <ref> shows the numerical solution for the pole contribution of F_z^p,c(t) 
(i.e. the first term on the r.h.s. of (<ref>)) for various values of the bias. 
For large times, the logarithm of this contribution shows indeed a straight line as function of 
ln(Ω t) with a slope given by -2αϵ^2/Ω^2

    ln(F_z^p,c(t)) = -2αϵ^2Ω^2ln(Ω t) + const  ,

where the constant term on the r.h.s. is independent of time but depends on the bias.  







§ SUMMARY



In this work we have presented the solution for the time dynamics of the ohmic spin boson 
model at finite bias by systematically expanding one order beyond Boch-Redfield. 
Using real-time RG and perturbation theory we have set up a renormalized 
perturbation theory to study analytically the whole time regime from exponentially small 
(Ω t∼ e^-1/α) up to large times (Ω t≫ 1). For very large times we used the
real-time RG method to sum up the leading logarithmic series in αln(Ω t). As a result
we obtained several interesting features for the time dynamics: (1) We showed how both the 
unrenormalized (Δ) and renormalized tunneling (Δ̃) enter the time dynamics 
and that it is not possible to account for the renormalization by using a local Hamiltonian 
with a renormalized tunneling. As in Ref. divincenzo_loss_05 we found that the 
renormalized Rabi frequency enters as high-energy cutoff scale to determine Δ̃.
(2) We found that all terms of the time evolution are exponentially damped by summing up all
secular terms ∼ (Γ t)^n. This results from a self-consistent perturbation theory in analogy
to the one presented in Ref. slutskin_etal_11. (3) For the preexponential functions of
the oscillating modes and in the non-exponential time regime we found logarithmic terms
∼αln(Ω t) containing the renormalized Rabi frequency as energy scale together with
terms falling off as α/(Ω t). (4) We showed that some correction terms in O(α) to 
Bloch-Redfield require an analysis of the Liouvillian up to second order in α. We were able
to calculate these terms by relating them to the stationary density matrix. (5) By resumming the
leading logarithmic series in αln(Ω t) in all orders of perturbation theory we found
for the preexponential functions of the oscillating modes an interesting crossover from
a power-law ∼ 1/(Ω t)^2α at exponentially small times to a power-law   
∼ 1/(Ω t)^2αϵ^2Ω^2 at exponentially large times. The latter has also
been proposed in Ref. slutskin_etal_11 but the logarithms determining the crossover to
the power-law at small times have not been discussed there. 

We have identified three important reasons why it is not sufficient to calculate the kernel of
the kinetic equation up to first order in the coupling to the bath to obtain all terms of the first correction 
to the Bloch-Redfield result. We now discuss why these issues are quite generic and are expected to occur
also for other models of dissipative quantum mechanics. 

First, for times of the order of the inverse decay rate t∼Γ^-1, where damping is still  
moderate, the distance of the Fourier variable E to some of the poles z_i of the
propagator is proportional to the decay rate |E-z_i|∼Γ. In this case perturbation theory is 
quite subtle since the denominator E-L(E) of the propagator is of O(Γ). The kernel Σ(E) 
can no longer be considered as a small
correction compared to E-L_0 and can not be expanded up to the numerator. We solved this problem
by expanding all analytic parts of Σ(E) around E=z_i and keeping Σ(z_i) in the denominator
whereas all other higher terms of the Taylor expansion are at least of O(α^2) and can be taken as
a small correction. The non-analytic terms of Σ(E) are more subtle and are some function 
f_i(E-z_i) when E is close to z_i, where f_i(E) ∼α is a non-analytic function 
with branch cut on the negative imaginary 
axis. For the ohmic spin boson model we get f_i(E-z_i)∼α (E-z_i)ln(-i(E-z_i))∼α^2 
such that it can be considered as a small correction. For dissipative quantum models with logarithmic
divergencies at high and low energies it is typical that Σ(E) has a logarithmic form,
see e.g. the Kondo model <cit.> or the interacting resonant level model
<cit.>, see Ref. schoeller_14 for a review. For weak coupling problems and E close
to z_i, Σ(E) contains either logarithmic terms ∼ln(-i(E-z_j)) with branching points 
z_j z_i (i.e. are analytic and can be expanded around E=z_i) 
or are proportional to (E-z_i)ln(-i(E-z_i)) (such that they vanish at
E=z_i). Terms ∼Γln(-i(E-z_i)) with a constant energy scale in front diverge at E=z_i and
are typical for strong coupling problems like e.g. the Kondo model. Most importantly, even for weak coupling
problems, it is never allowed to expand any part of Σ(E) in α by setting 
z_i=z_i^(0)+δ z_i, where z_i^(0) are the pole positions without the bath and 
δ z_i∼ O(α) denotes the correction from the bath, since (E-z_i^(0))/δ z_i is a parameter 
of O(1). Thus, for any model of dissipative quantum mechanics, it is very dangerous to use a naive 
perturbative expansion of the kernel in the coupling to the bath. The positions z_i of the branching points 
of L(E) (or poles of the propagator) should be kept non-perturbatively in a self-consistent way by
using the full propagator and not the bare one between the vertices, as also emphasized in
Ref. slutskin_etal_11. For a non-interacting bath described by a quadratic form 
H_bath=∑_qω_q a_q^† a_q in the field operators, the diagrammatic technique 
developed in Ref. schoeller_09 shows that all bare propagators can be replaced by
full ones without any double counting such that a systematic self-consistent perturbation theory can be
set up. Whether this is also possible for more complicated baths like e.g. spin baths is an open 
question. 

Secondly, we have seen that degenerate perturbation theory is generically needed since the decay poles 
z_i=-iΓ_i and the stationary pole z_st=0 of the propagator are close to each other within
the decay rate Γ_i∼α. Therefore, second order terms are needed for the Liouvillian 
to calculate the stationary state and all terms of the time evolution of the purely decaying modes
up to first order in α. Again this problem occurs only
for times of the order of the inverse decay rate, since for small times |E-L_0|∼ 1/t is 
much larger than Γ and can be considered as the largest term in the denominator of the
propagator such that the full kernel Σ(E) can be expanded up to the numerator. Thus, for two state
models with one purely decaying and two oscillating modes, the complicated terms in the stationary state and the 
purely decaying mode arising from the second order terms of the Liouvillian, must generically cancel for
small times. This simplifies the calculation of those terms for the purely decaying mode since they
can be expressed via the stationary state which, for the equilibrium case, can be easily calculated 
up to first order in α via the partition function. 
This strategy has been taken over in this work by using the stationary state 
calculated in Ref. divincenzo_loss_05 up to O(α). However, for generic models with more than
two local states, several purely decaying modes can occur and the problem of degenerate perturbation theory
can no longer be solved by just calculating the first correction to Bloch-Redfield of the stationary state. 

Whereas the two aforementioned issues are important to be considered for the calculation of the first correction 
to Bloch-Redfield on all time scales, there are further problems with weak coupling expansions in the
regimes of exponentially small or large times. They arise for problems of dissipative quantum mechanics 
with logarithmic divergencies at high and low energies like the ohmic spin boson model, the interacting
resonant level model, quantum dot models and the Kondo model. They have to be 
treated by an appropriate renormalization group method like the RTRG method <cit.>. 
For weak coupling problems, where the renormalized vertices stay small in the whole complex plane, the
RG equations can be truncated systematically such that logarithmic terms are summed up
non-perturbatively in leading or subleading order.  
Whereas logarithmic divergencies at high energies can be incorporated in renormalized parameters from poor man
scaling equations, logarithmic divergencies at low energies close to the branching points z_i are quite
subtle and require a full solution of the RG equations. For models of dissipative quantum mechanics without
logarithmic divergencies this issue is not important. 






§ ACKNOWLEDGMENTS

This work was supported by the Deutsche Forschungsgemeinschaft via RTG 1995. 
We thank D. Loss and D. diVincenzo for valuable discussions. 






§ LIOUVILLIAN IN PERTURBATION THEORY



Here we calculate the Liouvillian up to first order in the coupling α to the bath by using
the diagrammatic technique developped in Ref. kashuba_schoeller_13, where an expansion
in the coupling to the bath is used together with the application of Wick's theorem to integrate out
the phonon bath. In this reference it is shown for the ohmic spin boson model that the 
kernel Σ(E)=Σ_s+Σ_a(E) can be split into two parts,
one stemming from the symmetric and one from the antisymmetric part of the Bose distribution function of the
bath. At zero temperature this leads to Eq. (<ref>) with Σ_s given by (<ref>). The 
antisymmetric part Σ_a(E) involves only the antisymmetric part of the Bose distribution n(ω)
of the bath 

    n_a(ω) = 1 2(n(ω)-n(-ω))  = 1 2sign(ω)

since n(ω)=-θ(-ω) at zero temperature. The lowest order diagram for




Σ_a(E) is shown in Fig. <ref>, where the green line indicates the contraction 
between the bath field operators which involves the spectral density (<ref>)
of the bath and the antisymmetric part of the Bose distribution function via

    γ_a(ω) = 2αω J_c(ω) n_a(ω)
     = α|ω|D^2 D^2+ω^2  .

Using the diagrammatic rules the diagram is translated as

    Σ_a(E) = ∫ dωγ_a(ω) G R_a(E+ω) G   ,

where G is the bare vertex given by (<ref>) and R_a(E)=1/(E-L_a(E)) is the local
propagator of the antisymmetric part only. 
To approximate the ω-dependence of R_a(E+ω), we exhibit the logarithmic
parts by using the decomposition (<ref>) and use the spectral decomposition
(<ref>-<ref>) of L̃_Δ(E)

    R(E+ω)    = 1 E+ω-L̃_Δ(E+ω)Z'(E+ω)
       
    = ∑_i 1 E+ω-λ_i(E+ω)P_i(E+ω) Z'(E+ω)  .

Neglecting the ω-dependence of the logarithmic functions λ_i(E+ω), 
P_i(E+ω) and Z'(E+ω) (leading to higher orders in α), and using the
integral (defined for Im(E)>0 and analytically continued into the lower half of
the complex plane by choosing the branch cut along the direction of the
negative imaginary axis)

    ∫ dω |ω| D^2 D^2 + ω^21 E+ω    = D^2 D^2 + E^2 2E ln-iE D
        2E ln-iE D  ,

where ln(z) is the natural logarithm with branch cut on the negative real axis,
we find from (<ref>)

    Σ_a(E) = 2α∑_i F_i(E) G P_i(E)Z'(E) G   ,

with F_i(E) defined in (<ref>-<ref>). Taking the
projectors P_i(E) and Z'(E) in lowest order, 
given by (<ref>-<ref>) and Z'(E)^(0)=1,
and inserting (<ref>) for G, we find the result (<ref>-<ref>).

We note that the non-analytic features of Σ_a(E) in the lower half of the complex plane are
located at E=z_i-ix, 0<x<∞, where z_i are the positions of the poles of R_a(E). This holds 
exactly and can be shown in all orders of perturbation theory <cit.>. E.g. for the 
lowest order diagram (<ref>) we can see that this holds even when we do not use
any approximation for the ω-dependence of R_a(E+ω). Closing the integration contour
in the upper half and noting that R_a(E+ω) is an analytic function there and γ_a(E) has
non-analytic features only on the imaginary axis, we find the result

    Σ(E)    = i∫_0^∞ dx {γ_a(ix+0^+)-γ_a(ix-0^+)}·
       · G R_a(E+ix) G   .

Since R_a(E+ix) has a pole at E+ix=z_i, we find that Σ(E) is non-analytic for
E=z_i-ix with 0<x<∞. A similar proof can be used to show this in all orders of perturbation theory,
see Ref. schoeller_09,schoeller_14. 

Furthermore, we note that the matrix structure

    Σ_a(E) = ([       0       0;       0 Σ̂_a(E) ])

holds in all orders of perturbation theory. This is due to the fact that the bare vertices G have
the same structure, see (<ref>), and for each diagram all intermediate propagators are
sandwiched between two vertices. A consequence of this matrix structure is that
the projector P_st on the zero eigenvalue of L̃_Δ(E) is exactly known and given
by (<ref>). 





§ PROJECTORS FOR L̃_0+Σ̃_A^I



To calculate the projectors of the matrix L̃_0+Σ̃_a^i up to O(α), 
we first set up the matrix Σ̃_a^i=Σ̃_a(z_i) by setting E=z_i in  
(<ref>) and use

    F_0(0) ∼ O(α)  ,  F_0(z_0) = F_σ(z_σ) = 0   ,
       F_σ(0) , F_σ(z_0) , F_0(z_σ) 
     = -iπ 2Ω + O(α)   ,
       F_-σ(z_σ) = 2σΩln2-iπΩ+O(α)   ,

This gives for Σ̃_a^i transformed with the matrix A
(see (<ref>)) up to O(α) the result

    AΣ̃_a^st A    = AΣ̃_a^0 A = -iπα1Ω([     Δ̃^2 τ_-     Δϵτ_-σ_z; Δ̃^2Δϵσ_zτ_-          ϵ^2 ])   ,
    
    AΣ̃_a^σ A    = -iπαΔ̃^2Ω([   0   0;   0 τ_- ])  + 
        + α1Ωa_σ([          Δ̃^2τ_-     Δϵτ_-(σ_z-σ); Δ̃^2Δϵ(σ_z-σ)τ_-      ϵ^2(1-σσ_z) ]) 
      ,

where a_σ=2σln2-iπ. The transformed Liouvillian AL̃_0A is given by
(<ref>). Due to the matrix structure of Σ̃_a(E), one 
projector is exactly known (in all orders of perturbation theory, see (<ref>) and
(<ref>))

    A P̃^i_st A    = ([ τ_+   0;   0   0 ])  .

Using usual perturbation theory it is straightforward to calculate the projectors 
P̃_σ^st,0,σ in zero and first order in  α as

    A P̃^(0)i_σ A    = 1 2([      0      0;      0 1+σσ_z ])  ,
    
    A P̃^(1)st,0_σ A    =
       
    = -1 2iπαΔ̃^2ϵΔΩ^2([              0 1 Zτ_-(1+σσ_z);    (1+σσ_z)τ_-              0 ])  ,
    
    A P̃^(1)σ_σ A    = 1 4iπσαΔ̃^2Ω^2([       0       0;       0 τ_+-τ_- ])  .

Using ∑_j=st,0,±AP̃_j^iA=1, we find

    A P̃^(0)i_0 A    = ([ τ_-   0;   0   0 ])  ,
    
    A P̃^(1)st,0_0 A    = iπαΔ̃^2ϵΔΩ^2([      0 1 Zτ_-;    τ_-      0 ])  .

This proves (<ref>-<ref>) and (<ref>-<ref>). 
We note that although degenerate perturbation theory is needed to calculate P̃_0^i up
to first order in α, we do not need any second order terms in α for the Liouvillian
since the projector P̃_st^i is exactly known. This is a particular advantage for the
spin boson model.

To derive the formula (<ref>) for the eigenvalue γ̃_0^st of
L̃_a(0) up to second order in α, we relate it to the eigenvalue 
γ̃_0^0=z_0=-i(Γ^(1)+Γ^(2)+O(α^3)) of L̃_a(z_0). We first note 
that, due to the matrix structure of AL̃_a^st,0A 
(see (<ref>) and (<ref>)), 
the second order contribution to the eigenvalues γ̃_0^st,0 is not influenced by
the non-diagonal blocks of AL̃_a^st,0A and arises only from the upper left block.
Denoting this block by (AL̃_a^st,0A)_11 we expand

    L̃_a(0)_11 = L̃_a(z_0)_11 - dΣ̃_a dE(0)_11z_0 + O(α^3)
      .

and use (<ref>) together with 
d dEF_σ(0)=1+iσπ 2+O(α) to get

    L̃_a(0)_11 = L̃_a(z_0)_11 + 
    2iΓ^(1)αΔ̃^2Ω^2τ_-  + O(α^3)  ,

such that 

    γ̃_0^st = z_0  + 2iΓ^(1)αΔ̃^2Ω^2 + O(α^3)
      ,

which proves (<ref>).





99

leggett_87
  A.J. Leggett, S. Chakravarty, T.A. Dorsey, M.P.A. Fisher, A. Garg, and W. Zwerger,
Rev. Mod. Phys. 59, 1 (1987).
  
weiss_12
  U. Weiss, Quantum Dissipative Systems (World Scientific Publishing Company, Singapore, 2012).

grifoni_haenggi_98
M. Grifoni and P. Haenggi, Phys. Rep. 304,229 (1998).

hutter_loss_13
A. Hutter and D. Loss, Phys. Rev. A 89, 042334 (2014).

novais_mucciolo_13
E. Novais and E.R. Mucciolo, Phys. Rev. Lett. 110, 010502 (2013).

pletyukhov_schuricht_schoeller_10
M. Pletyukhov, D. Schuricht, and H. Schoeller, Phys. Rev. Lett. 104, 106801 (2010).

schoeller_09 
H. Schoeller, Eur. Phys. J. Spec. Top. 168, 179 (2009).

schoeller_14
H. Schoeller, “Dynamics of open quantum systems”, in 
Lecture Notes of the 45th IFF Spring School 
“Computing Solids - Models, ab-initio methods and supercomputing” (Forschungszentrum Jülich, 2014).

pletyukhov_schoeller_12
M. Pletyukhov and H. Schoeller, Phys. Rev. Lett. 108, 260601 (2012).

van_hove_loss
L. van Hove, Physica 21, 517 (1955); D. Loss, Physica A 139, 505 (1986).

slutskin_etal_11
A. A. Slutskin , Europhys. Lett. 96, 40003 (2011).

kashuba_schoeller_13
O. Kashuba and H. Schoeller, Phys. Rev. B 87, 201402(R) (2013).

divincenzo_loss_05
D.P. DiVincenzo and D. Loss, Phys. Rev. B 71, 035318 (2005).

egger_etal_97
R. Egger, H. Grabert, and U. Weiss, Phys. Rev. E 55, R3809 (1997).

kennes_etal_13
D. M. Kennes, O. Kashuba, M. Pletyukhov, H. Schoeller, and V. Meden, 
Phys. Rev. Lett. 110, 100405 (2013).

projection_operator
H. Mori, Progress of Theoretical Physics 33, 423 (1965);
H. Grabert, Projection Operator Techniques in Nonequilibrium Statistical Mechanics, 1st ed., 
Springer Tracts in Modern Physics, Vol. 95 (Springer, Berlin, 1982) p. 166;
E. Fick and G. Sauermann, The Quantum Statistics of Dynamic Processes
(Springer Series in Solid-State Sciences 86, Springer-Verlag, Berlin, 1990);
C. W. Gardiner, Quantum Noise
(Springer Series in Synergetics 56, Springer-Verlag, Berlin, 1991);
R. Zwanzig, Nonequilibrium Statistical Mechanics, 3rd ed. (Oxford University Press, New York, 2001) p. 240.

schoeller_reininghaus_09
H. Schoeller and F. Reininghaus, Phys. Rev. B 80, 045117 (2009);
ibid. Phys. Rev. B 80, 209901(E) (2009).

RTRG_irlm
C. Karrasch, S. Andergassen, M. Pletyukhov, D. Schuricht, L. Borda, V. Meden, and H. Schoeller, 
Europhys. Lett. 90, 30003 (2010);
S. Andergassen, M. Pletyukhov, D. Schuricht, H. Schoeller, and L. Borda, Phys. Rev. B 83, 205103 (2011).










