Dynamics of open quantum systems
Herbert Schoeller





Lecture Notes of the 45^ th IFF Spring
School “Computing Solids - Models, ab initio methods and supercomputing”
(Forschungszentrum Jülich, 2014). All rights reserved. 




















§ INTRODUCTION


An open quantum system open quantum system
consists of a local quantummechanical system of fixed size coupled 
to infinitely large reservoirs in statistical equilibrium via a well-defined interaction, 
see Fig. <ref> for a sketch of the system. The analysis of the time
evolution time evolution of the reduced density matrix of the local system
is of fundamental importance for nonequilibrium statistical mechanics. Of particular interest
is the interplay of quantum coherence in the local quantum system and dissipation generated 
by the reservoirs, which is the reason why this field is called dissipative quantum mechancis
dissipative quantum mechanics. Many topics are here of current interest: 
(1) The development of a microscopic theory for irreversible
time evolution of the local density matrix; (2) The characterization of the reduced dynamics,
in particular by generic features independent of the microscopic details of high-energy
processes; (3) The calculation of typical relaxation and decoherence rates; (4) The analysis of
quantum fluctuations induced by the system-reservoir interaction beyond perturbation theory; 
(5) The analysis of the influence of strong correlations in the local quantum system, induced by 
charging energies; (6) The analysis of the influence of inhomogeneous boundary conditions, induced by different
temperatures and/or chemical potentials of several reservoirs, leading to energy, particle,
and spin currents; (7) The analysis of deviations of the stationary local density matrix from a
grandcanonical one, induced by quantum fluctuations from the system-reservoir coupling or by the
presence of several reservoirs; (8) The development of nonequilibrium renormalization group (RG)
renormalization group
methods capable of resumming logarithmic divergencies occuring in higher-order perturbation
theory in the system-reservoir coupling, either at high energies w.r.t. the band width
of the reservoirs or at low energies w.r.t. the inverse time in the long-time limit; 
(9) The study of non-Markovian dynamics non-Markovian dynamics
leading to additional terms in the time evolution with
unexpected oscillation frequencies and decay rates together with non-exponential time evolution;
(10) The crossover between coherent and incoherent dynamics induced by the sign and size of the
system-reservoir coupling and other tunable parameters. It is the purpose of this
tutorial introduction to present a microscopic theory for the time evolution of open quantum
systems, to discuss some of the above aspects from a generic point of view, and to characterize 
explicitly the reduced dynamics for elementary 2-level quantum systems coupled via energy, 
particle or spin exchange to external reservoirs. 

Although the field of dissipative quantum mechanics has a long history, the field has regained an
enormous interest in the last decades due to its relevance in quantum transport phenomena in 
nanoelectronic systems and quantum information processing, and due to the controlled realization of 
low-dimensional quantum systems in cold atom gases. To describe the time evolution of the reduced
density matrix of the local quantum system microscopically, one starts from the von Neumann equation 
for the total density matrix ρ_tot(t) of the full system (we set ħ=e=k=1)

    i ρ̇_tot(t) = [H_tot(t),ρ_tot(t)] = 
    L_tot(t) ρ_tot(t)  ,

where L_tot(t) is the so-called 
Liouville operator, a superoperator superoperator which acts on
an arbitrary operator via L_tot(t)A=[H_tot(t),A]. The central idea is always 
to integrate out the reservoir degrees of freedom and to set up a formally exact kinetic equation
kinetic equation for the local density matrix 
ρ(t)=Tr_resρ_tot(t), defined by the
trace Tr_res over the reservoir degrees of freedom of the total density matrix.
This kinetic equation has the form

    i ρ̇(t) = ∫_t_0^t dt' L(t,t') ρ(t')  ,

where t_0 is the initial time and L(t,t') is an effective Liouville operator acting only on 
operators of the local quantum 
system. This superoperator contains all the information of the reservoir degrees of freedom and the 
system-reservoir interaction. For a time-translational invariant Hamiltonian, L(t,t')=L(t-t') 
depends only on the relative time difference. The effective Liouvillian effective Liouvillian
L(t,t') is only defined for t>t',
i.e. it acts as a response function relating the density matrix at time t' to the one at the later
time t. This acounts for memory effects
and leads to non-Markovian dynamics. The only assumption needed to derive the kinetic equation 
(<ref>) is the factorization of the total density matrix at the initial time 
t_0 in an arbitrary local part ρ(t_0) and an equilibrium part for the reservoirs

    ρ_tot(t_0) = ρ(t_0) ρ_res^eq  ,  ρ_res^eq = ∏_αρ_α^eq  ,  ρ_α^eq = 1 Z_α e^-(H_α-μ_α N_α)/T_α  ,

where T_α, μ_α, H_α, N_α, and Z_α are the temperature, the chemical 
potential, the Hamiltonian, the particle number, and the partition function of reservoir α, respectively.
However, by changing the Hamiltonian at a certain quench time t_q>t_0 abruptly, other initial
conditions can be realized where system and reservoirs are correlated. 

Various techniques have been developed to calculate the effective Liouvillian L(t,t'). The 
traditional ones are projection operator techniques <cit.> and functional 
integrals <cit.>. Recently, a quantum field theoretical approach has been
developed, which allows for a systematic classification of all processes in all orders of
perturbation theory in the system-reservoir coupling <cit.>. With this method,
it is possible to go beyond bare perturbation theory which is necessary at low temperatures
due to various logarithmic divergencies at high and low energies. The method is capable of  
identifying these logarithmic divergencies very effectively and an RG method in nonequilibrium 
has been set up to resum them. This allows a systematic weak-coupling expansion in
the renormalized coupling constants to be formulated with which the time evolution on all
time scales even when the reservoirs have different chemical potentials or temperatures
can be discussed. This technique has been applied successfully to the Kondo model Kondo model 
<cit.>, the interacting
resonant level model (IRLM) interacting resonant level model
<cit.>, 
and the ohmic spin boson model spin boson model <cit.>. 
In particular, it has turned out that the RG formulation is most effective for the calculation of the time
evolution if the Fourier variable E conjugate to the time t is used as flow parameter, i.e. as
the paramater w.r.t. which derivatives of the various quantities of interest are taken to obtain
differential equations (the so-called RG equations). This technique is called
the E-RTRG E-RTRG method <cit.>. The models treated so far fall
into the special class where the density of states in the reservoirs and the frequency dependence
of the system-reservoir couplings is weak and varies on the scale of the high-energy cutoff D. 
Physically, the high-energy cutoff can either be the band width of the reservoirs or it is some 
internal high-energy scale of the local quantum system, like e.g. charging energies, arising when 
effective models are used by integrating out high-energy processes (e.g. quantum dots in the 
Coulomb blockade regime where charge degrees of freedom can be eliminated, see the lecture B3 by 
T. Costi). For such models it is often possible to find universal physics where the special form of the 
high-energy cutoff function is not important and influences only the value of certain low-energy
scales (e.g. the Kondo temperature for the Kondo model). In such a case the high-energy cutoff D
does no longer occur explicitly. Furthermore, for a wide class of time-translational invariant models 
it turns out that the effective Liouvillian has the form

    L(E) = L_Δ(E) + E L'(E)  ,

where L(E)=∫_0^∞ dt e^iEtL(t) is the Fourier-transform of the response function
L(t-t')=L(t,t')θ(t-t'). In this decomposition L_Δ(E) and L'(E) are slowly varying
logarithmic functions, where L_Δ(E) is proportional to some energy scale Δ of the model 
which can be anything except for the Fourier variable E. This form will be shown by the RG analysis in 
Section <ref> for the concrete models under consideration but it remains an interesting 
question for the future how generic this form is. A large part of this tutorial deals with
the technical details of calculating the appearing functions L_Δ and L'. Before we
do that, we will first investigate the physical consequences for the time evolution in
Section <ref>. We will see that when L(E) has the form (<ref>) 
the time evolution can generically be decomposed as

    ρ(t) = ∑_n F_n(t) e^-iz_n t ρ_t=0  ,

where z_n=±Ω_n-iΓ_n, with Ω_n,Γ_n≥ 0, determine the oscillation
frequencies and decay rates of exponential decay, and F_n(t) are pre-exponential functions,
which typically consist of power-laws ∼ 1/t^k (k=1,2,…) and logarithmic corrections
in the long-time limit t≫ 1/|z_n|. At least one of the exponential scales is zero z_st=0, which
determines the stationary state. 

It is the purpose of the present article to first discuss the generic physics of the time evolution
on the basis of the form (<ref>) of the effective Liouvillian, and with this 
motivation discuss the E-RTRG method for the calculation of L(E) and its decomposition into
(<ref>). Then we will summarize the results for the time evolution of the 
Kondo model, the ohmic spin boson model, and the IRLM. We note that other RG methods 
have been developed recently to discuss the 
time evolution of open quantum systems. The most important ones are the flow-equation method 
<cit.> and the functional RG <cit.>. The latter will be introduced in
the lecture B7 by V. Meden and is a method where one expands systematically in the short-ranged
renormalized interaction parameter present in the local system but not in the system-reservoir coupling, 
i.e. it is complementary to the RTRG technique where arbitrary local interactions can be treated but an
expansion in the renormalized system-reservoir coupling is needed. Besides the analytical RG methods,
there is also an extensive research going on to develop numerical methods to describe the time evolution,
like e.g. time-dependent numerical renormalization group <cit.>, time-dependent density matrix
renormalization group <cit.>, iterative stochastic path integrals <cit.>, and quantum 
Monte Carlo <cit.>. Furthermore, for special models, field-theoretical methods have been used to
find exact results <cit.>.  





§ BASIC MODELS



We start with the description of the basic models under consideration, where the quantum system 
consists of 2 states coupled via spin (Kondo model), charge and potential (IRLM), or energy (spin boson)
fluctuations to a noninteracting environment. The total Hamiltonian is assumed to be time-translational
invariant and consists of three parts

    H_tot = H + H_res + V  ,   H_res = ∑_α H_α  ,  
    H_α = ∑_kσ ϵ_ασ k a^†_ασ ka_ασ k

where H is the Hamiltonian of the local quantum system, V is the system-reservoir interaction, and
H_res describes the noninteracting (fermionic or bosonic) reservoirs with field operators
a_ασ k. α is the reservoir index, σ the channel index (e.g. spin), and
the quantum number k characterizes the energy. For convenience, for given α and σ, we 
will denote by ω=ϵ_ασ k-μ_α the energy of the reservoir states 
measured relative to the chemical potential, and we assume that the relation between ω and k is
unique. As a consequence, the field operators of the reservoirs can be characterized by the multi-index
1≡ηασω, where η=± distinguishes between creation (η=+) and 
annihilation operators (η=-). Depending on the model under consideration, we will define below
convenient forms of the field operators a_1≡ a_ηασ(ω) in continuum notation,
such that the commutation relations read (the upper/lower case refers always to bosons/fermions)

    [a_ασ(ω),a^†_α'σ'(ω')]_∓ = δ_αα' δ_σσ' δ(ω-ω') ρ_ασ(ω)  ,

where [·,·]_∓ denotes the commutator/anticommutator for bosons/fermions.
As defined below the spectral function ρ_ασ(ω) contains the d.o.s. of the 
reservoirs and possibly frequency-dependencies of the system-reservoir couplings. Together with
the commutation relations 

    [a_1,H_α] = -η (ω + μ_α) a_1   ,   
    [a_1,N_α] = -η a_1  ,

it follows that the contraction of two reservoir field operators w.r.t. the equilibrium distribution
is given by

    a_1 a_1'(-20,11) 
        (-22,8)(0,1)3(-22,11)(1,0)12(-10,8)(0,1)3(20,11) 
       = Tr_res a_1 a_1' ρ_res^eq = δ_11̅' ρ_ασ(ω) f^η_α(ω)
     = δ_11̅' {[ η; 1 ]} ρ_ασ(ω) 
    f_α(ηω)  ,

where 1̅≡ -ηασω is defined by reversing the sign of η,
δ_12=δ_η_1η_2δ_α_1α_2δ_σ_1σ_2δ(ω_1-ω_2),
f^+_α(ω)=f_α(ω), f^-_α(ω)=1± f_α(ω), and
f_α(ω)=(e^ω/T_α∓ 1)^-1 is the Bose/Fermi distribution. 

In terms of the continuum field operators, the system-reservoir interaction V is 
generically written as a sum of terms of the form 

    V = 1 n! {[          1; η_1η_2…η_n ]}  :a_n a_n-1… a_1:  g_12… n → 1 n! g_12… n :a_1a_2… a_n:   ,

where n=1,2,… is any integer, imlicit summation/integration is assumed over the multi-indices
i≡η_iα_iσ_iω_i, the operator g_12… n acts only on the 
local system, and :… : denotes normal-ordering w.r.t. to the equilibrium distribution 
(<ref>) of the reservoirs (i.e. in any Wick-decomposition contractions are
not allowed within the normal-ordered expression). We call the operators g_1… n n-point vertex
operators since, together with the corresponding superoperators (<ref>),
they will appear in the diagrammatic technique as vertices with n reservoir lines, see 
Section <ref>. For bosons the two forms for V shown in (<ref>)
are the same. For fermions, the first form is needed for the definition of the vertex operators and, for 
n odd, g_1… n is of fermionic nature and anticommutes with the reservoir field operators. However,
it can be shown <cit.> that, for the calculation of any local observables, the second
form for V can be used and local and reservoir operators can be taken as if they commute.
The vertex operators have the properties 

    g_1… i… j… n = ± g_1… j… i… n  ,  
    g^†_1… n = g_n̅…1̅  .

The first relation can always be achieved by (anti-)symmetrization of the reservoir field operators
within the normal-ordering in (<ref>), whereas the second one is necessary for 
the property V=V^†. In the following we will specify the definition of the continuum 
reservoir field operators a_1, the spectral density ρ_ασ(ω) and the 
vertex operators g_1… n for the various models. 




The Kondo model. In its most basic form the Kondo model describes a local spin-1 2 system coupled via
short-ranged and isotropic exchange couplings to fermionic reservoir spins, see Fig. <ref> for a 
sketch of the system. It is a model system to describe local spin fluctuations. spin fluctuations
For the case of a single channel the Hamiltonian reads

    H = h^(0) S_z   ,   
    V = ∑_αα' J_αα'^(0)√(ρ_α^(0)ρ_α'^(0)) ∑_σσ' 1 2σ_σσ' 
    :∑_k a^†_ασ k∑_k'a_α'σ' k': S  ,

where the isotropic exchange couplings J^(0)_αα'=J^(0)_α'α are real and
dimensionless, σ=(σ_x,σ_y,σ_z) are the Pauli matrices, 
S is the local spin, and h^(0) is the local bare magnetic field. ρ^(0)_α denotes the
d.o.s. of the reservoirs at the Fermi level. The Kondo model can be derived via a Schrieffer-Wolff 
transformation from the single-impurity Anderson model (see the lecture B3 by T. Costi),
in which case the exchange couplings fulfil the relation 

    J^(0)_αα' = 2 √(x_α x_α') J^(0)  ,  ∑_α x_α = 1
      ,   0<x_α<1  ,

where x_α are asymmetry factors weighting the energy broadening of the local level
from reservoir α. Defining the continuum field operator by 
a_1=1√(ρ_α^(0))∑_k δ(ω-ϵ_ασ k+μ_α)a_ηασ k, 
with 1≡ηασω, and
assuming a flat d.o.s. in the reservoirs, we obtain for the spectral function and the vertex operator

    ρ_ασ(ω)    = ρ(ω) = 1ρ_α^(0)∑_k δ(ω-ϵ_ασ k+μ_α) = D^2 D^2+ω^2  ,
    
    g_11'    = 1 2 J_αα'^(0)  σ_σσ'·S  ,  for  η=-η'=+  ,

where D is the band width of the reservoirs and, for convenience, we have chosen a Lorentzian for
the high-energy cutoff function. The case η=-η'=- is obtained from g_11'=-g_1'1. 

The IRLM. The IRLM is a basic model to describe charge and potential fluctuations.
potential fluctuations It consists
of a single fermionic level, which is coupled to fermionic reservoirs via tunneling and
a local Coulomb interaction, see Fig. <ref> for a sketch of the system. Disregarding
the spin, the Hamiltonian is defined by 

    H = ϵ c^† c   ,  
    V = ∑_α t_α√(ρ_α^(0)) ∑_k (a^†_α kc + h.c.)
     + ∑_α U_αρ_α^(0) (c^† c - 1 2) ∑_kk' :a^†_α ka_α k':  ,

where c is the field operator annihilating a particle on the local system, 
t_α are the tunneling matrix elements (in units of 1√(ρ_α^(0))),
U_α denote the dimensionless Coulomb couplings, and ϵ is the bare energy of
the local level. At ϵ=0 the model fulfils particle-hole symmetry. Defining the 
continuum field operators as for the Kondo model, with 1≡ηαω (i.e. omitting the
spin index), we find the same result (<ref>)
for the spectral function, and the vertex operators are given by

    g_1 = t_α {[       c for η=+;     c^† for η=- ].
      ,  
    g_11' = δ_η,-η' δ_αα' η U_α (c^† c -1 2)  .







The spin boson model. The spin boson model describes energy fluctuations, energy fluctuations 
where a 2-level system is
coupled linearly to a phonon bath, see Fig. <ref> for a sketch of the system.
The Hamiltonian is given by

    H_res = ∑_k ω_k a^†_k a_k  ,  
    H = 1 2 ϵ σ_z - 1 2 Δ σ_x   ,  
    V = 1 2 σ_z ∑_k α_k (a_k + a_k^†)  ,

where ϵ and Δ denote the bias and the tunneling of the local
2-level system, respectively. The phonon frequencies ω_k>0 are positive, and the 
equilibrium phonon distribution is characterized by temperature T. The k-dependence of
the real coupling constants α_k and the phonon frequencies ω_k is considered by defining the 
continuum field operators by a_1=∑_k α_k δ(ω-ω_k)a_η k with 1≡ηω. 
This leads to the following spectral function and vertex operator

    ρ(ω) = ∑_k α_k^2 δ(ω-ω_k) = 
    2 α ω (ω D)^s-1θ(ω) D^2 D^2+ω^2  ,  
    g_1 = 1 2 σ_z   ,

where α is a dimensionless coupling constant, and we have again chosen a Lorentzian
high-energy cutoff function with band width D. The special form chosen for ρ(ω) 
describes the ohmic case for s=1 considered in this article, whereas s<1 (s>1) define 
the sub-ohmic (super-ohmic) cases. For the special case Δ=0 the spin boson model can
be solved exactly <cit.> with the result

    ⟨σ_x,y⟩(t) = e^-h(t) ⟨σ_x,y⟩_t=0  ,  ⟨σ_z⟩(t) = ⟨σ_z⟩_t=0  ,

with h(t)=-∫ dω(ρ(ω)/ω)(1-cos(ω t))(1+2f(ω)), where
f(ω)=(exp(ω/T)-1)^-1 is the Bose function.





For the special case of α close to 1 2, the ohmic spin boson model can be mapped
on the IRLM with a single reservoir (with μ=0) <cit.>. The parameters
U and t of the IRLM are related to α and Δ of the ohmic spin boson
model in the following way

    U = 1 - √(2α)  ,  Γ^(0) = 2 π t^2 = Δ^2 D  .

The local occupation ⟨ n⟩(t)=⟨ c^† c⟩(t) of the IRLM is 
related to the expectation value ⟨σ_z⟩(t) of the ohmic spin boson model via

    2 ⟨ n ⟩(t) - 1 = ⟨σ_z⟩(t)  ,

whereas the expectation value ⟨σ_x,y⟩(t) of the spin boson model 
is related to expectation values of highly nonlinear operators involving reservoir degrees of
freedom in the IRLM. The value α=1 2 is of special importance since, at this point,
the time evolution of ⟨σ_z⟩(t) changes from an oscillating one 
(for α<1 2) to a purely decaying one (for α>1 2) 
<cit.>. 
Correspondingly, for the IRLM, this crossover
occurs when the sign of the Coulomb interaction U is changed.




§ KINETIC EQUATION AND TIME EVOLUTION



In this section we aim at discussing the time evolution from a generic point of view
based on the general form (<ref>) of the kinetic equation and the
form (<ref>) of the effective Liouvillian for the case of a time-translational
invariant Hamiltonian. Using L(t,t')=L(t-t') the kinetic equation reads

    iρ̇(t) = ∫_0^t dt' L(t-t') ρ(t')  ,

where, for convenience, we have set the initial time t_0=0.
The reduced density matrix ρ(t) acts only in local space, i.e. has matrix elements
ρ(t)_ss'=⟨ s|ρ(t)|s'⟩, where s and s' are states of the local quantum system.
In contrast, the superoperator L(t) acts on local operators A, i.e. the matrix elements can be written
as L_s_1s_2,s_1^' s_2^'=⟨ s_1 s_2|L(t)|s_1^' s_2^'⟩, where
|ss'⟩=|s⟩⟨ s'| are the basis elements (= operators) in Liouville space 
Liouville space and
⟨ ss'|A=⟨ s|A|s'⟩ are the corresponding dual vectors. The density matrix fulfils
the property of conservation of probability Trρ(t)=1 and is self-adjoint ρ(t)=ρ(t)^†.
It is straightforward to show that the kinetic equation respects these properties if and only if the 
effective Liouville operator fulfils the properties

    Tr L(t) = ∑_s L(t)_ss,·· = 0   ,  
    L(t)^c = - L(t)   ,

where the c-transform is defined by 
L(t)^c_s_1s_2,s_1^' s_2^'=L(t)^*_s_2s_1,s_2^' s_1^' and fulfils the useful property
(L(t)A)^† = L(t)^c A^†. In Fourier space L(E)=∫_0^∞ dt e^iEtL(t) this 
means

    Tr L(E) = 0   ,  
    L(E)^c = - L(-E^*)   ,

or for the quantities L_Δ(E) and L'(E) appearing in the decomposition (<ref>)

    Tr L_Δ(E) = Tr L'(E) = 0  ,  
    L_Δ(E)^c = -L_Δ(-E^*)  ,  
    L'(E)^c = L'(-E^*)  .




With ρ(E)=∫_0^∞ dt e^iEtρ(t), the kinetic equation reads in Fourier space
Eρ(E)-iρ_t=0=L(E)ρ(E) leading to the formal solution 

    ρ(E) = i R(E) ρ_t=0  ,   R(E) = 1 E - L(E)  .

We now investigate the consequences of the generic form (<ref>) of the 
effective Liouvillian L(E). Using inverse Fourier transform, the time evolution can be 
calculated for t>0 from

    ρ(t) = i 2π ∫_-∞+i0^+^∞+i0^+dE 
    e^-iEt R(E) ρ_t=0 = i 2π ∫_-∞+i0^+^∞+i0^+dE 
    e^-iEt R̃(E) Z'(E) ρ_t=0  ,

where we have defined 

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

By convention, Z'(E) is called the Z'-factor operator. 



The last form of (<ref>) is very helpful for the evaluation of the
energy integral because it explicitly exhibits the slowly varying logarithmic functions
L̃_Δ(E) and Z'(E). 
The energy integral ∫ dE is calculated by closing
the integration contour in the lower half of the complex plane and deforming the contour
such that the poles and branch cuts of the integrand are enclosed, see Fig. <ref>.
To identify the singularities of the integrand we use the spectral decomposition of the 
Liouvillian L̃_Δ(E) in terms of its eigenvalues λ_k(E) and 
corresponding projectors P_k(E) 

    L̃_Δ(E) = ∑_k λ_k(E) P_k(E)  .

Since we deal with a non-hermitian superoperator, we have to distinguish the left and
right eigenvectors, which we denote in Dirac notation by |x_k(E)⟩ and ⟨x̅_k(E)|

    L̃_Δ(E) |x_k(E)⟩ = λ_k(E) |x_k(E)⟩  ,  ⟨x̅_k(E)| L̃_Δ(E) = ⟨x̅_k(E)| λ_k(E)  .

The eigenvectors fulfil the orthonormalization condition 
⟨x̅_k(E)|x_k'(E)⟩=δ_kk' and the projectors are given by 
P_k(E)=|x_k(E)⟩⟨x̅_k(E)| with ∑_k P_k(E)=1. 

Due to the condition TrL̃_Δ(E)=0, we obtain either λ_k(E)=0
or Tr |x_k(E)⟩=0. Therefore, the Liouvillian has always an eigenvalue zero,
which we characterize by the index k=st since it corresponds to the stationary state (see below).
The other eigenvalues are numerated by k=0,± 1,± 2, …. Normalizing the eigenvector with 
k=st according to Tr|x_st(E)⟩=1 and using 
⟨x̅_st(E)|=Tr, we get

    Tr |x_st(E)⟩    =    ∑_s ⟨ ss|x_st(E)⟩ = 1  ,  ⟨x̅_st(E)|ss⟩ = 1
    Tr |x_k(E)⟩    =    ∑_s ⟨ ss|x_k(E)⟩ = 0  ,  for   k=0,± 1,± 2,…  .

As a consequence we get 

    P_st(E) = |x_st(E)⟩ Tr  ,

and the property
TrL_Δ(E)=TrL'(E)=0 can also be written as

    P_st(E) Z'(E) = P_st(E)  ,  
    P_st(E) L_Δ(E) = 0  .

Due to the condition L̃_Δ(E)^c=-L̃_Δ(-E^*), the eigenvalues and projectors 
occur always in pairs (except for k=0,st where we define k≡ -k) with 

    λ_-k(E) = -λ_k(-E^*)  ,  
    P_-k(E) = P_k(-E^*)^c  .




Using the spectral representation, the time evolution can be written as

    ρ(t) = i 2π ∑_k ∫_γ dE e^-iEt 1 E-λ_k(E) P_k(E) Z'(E) ρ_t=0  ,

where γ is an integration contour which encloses the lower half of the
complex plane including the real axis. Poles are located at E=z_k^p=λ_k(z_k^p)=±Ω_k-iΓ_k,
with Ω_k,Γ_k≥ 0, where z_st^p=0 is a pole at the origin, see 
Fig. <ref>.
At zero temperature, which we consider from now on, additional nonanalytic features occur from
branch cuts since λ_k(E), P_k(E) and Z'(E) depend
logarithmically via terms ∼ln(D E-z_n) generated by the ultraviolet
divergencies from the high-energy cutoff D (at finite temperature the branch cuts turn
into an infinite number of discrete poles separated by 2π T). From the structure of the
perturbation theory (see below) it can be seen that the singularities z_n are associated
with poles of the resolvents R̃(E_1… n), where 

    E_1… n = E + μ̅_1… n  ,  μ̅_1… n = μ̅_1 +… + μ̅_n  ,  μ̅_1 = η_1 μ_α_1  ,

i.e. are located at z_n=E with E_1… n=E+μ̅_1… n=z_k^p.
Therefore, the singularities z_n=z_k^p-μ̅_1… n are generically given by the 
poles shifted by some linear combination of the chemical potentials of the reservoirs. 

In Section <ref> we will see how L̃_Δ(E) and Z'(E) can be determined from 
differential equations, see Eq. (<ref>),
where we differentiate w.r.t the Fourier variable E. These differential 
equations are defined in the whole complex plane and will be the RG equations of the E-RTRG method. 
E is called the flow paramater and a solution of the RG equations along a certain path is called
the RG flow. The particular advantage
is that these RG equations can be solved along the paths E=z_n+iΛ± O^+, with Λ
real, starting at some high value Λ∼ D down to Λ=-∞. Since no singularities
are present on these paths, it can even be numerically enforced that the branch cuts start at z_n
and point into the direction of the negative imaginary axis. Furthermore, the jump of the Liouvillian 
at the branch cuts can be determined from the difference of the two solutions and the integrals
around the branch cuts can be calculated. The choice that the branch cuts point into the direction
of the negative imaginary axis is very convenient since e^-iEt=e^-iz_n te^-xt 
is exponentially decaying in xt, which allows an analytical discussion of the long-time limit
(see below). Using E=z_n-ix± O^+, the integration around a particular branch cut (including the
case when the branching point is a pole) gives the contribution ρ_n(t)=F_n(t)e^-iz_n tρ_t=0
to the time evolution with

    F_n(t) = 1 2π ∫_0^∞ dx e^-xt {R(z_n-ix+0^+)-R(z_n-ix-0^+)}

such that the total time evolution can be written in the form (<ref>)

    ρ(t) = ∑_n ρ_n(t) = ∑_n F_n(t) e^-iz_n t ρ_t=0  .

For the further evaluation of F_n(t) it is important to distinguish between the cases
when the branching point is a pole or not. We label the contributions from branching poles
z_k^p by F_k^p(t) and ρ_k^p(t) and the others by F_n^b(t) and ρ_n^b(t), such
that (<ref>) reads

    ρ(t)    = ∑_k ρ_k^p(t) + ∑_n ρ_n^b(t)
        = ∑_k F_k^p(t) e^-iz_k^p t ρ_t=0 + ∑_n F_n^b(t) e^-iz_n^b t ρ_t=0  .

Thereby we note that the same singularity z_k^p=z_n^b can appear as a branching pole and as a branching point, 
since a certain term involving λ_k(E) in (<ref>) can have a branch
cut at z_n^b=z_k'^p with k' k.
Generically, for weakly coupled system-reservoir systems, the contributions ρ_n^b(t) are smaller 
since they are proportional to the system-reservoir coupling (see below Eq. (<ref>)). 
However, if
the decay rates occuring in z_n^b are smaller than those ones of z_k^p, the relative order of the
various terms can change as function of time, as discussed e.g. in detail in 
Refs. <cit.>
for the IRLM with positive Coulomb interaction or the ohmic spin boson model for α close but
slightly below the value α=1 2. In the Markovian approximation, only the contributions
ρ_k^p(t) remain and the pre-exponential functions are approximated by constants of O(1). 

Time-evolution regimes.
Using the general expressions (<ref>) and (<ref>),
one can discuss the qualitative form of the time evolution in different time regimes. For
short times t≪ 1/|z_n|, only high frequencies E∼ 1/t≫ |z_n| matter in 
Eq. (<ref>), i.e. the cutoff scales z_n 
in the logarithmic terms are unimportant and can be neglected. Furthermore, to leading order,
we can replace E→ 1/t in the logarithmic parts, and we obtain from 
(<ref>)

    ρ(t) = i 2π ∫_γ dE e^-iEt 1 E-L̃_Δ(1/t) Z'(1/t) ρ_t=0 = e^-iL̃_Δ(1/t)t Z'(1/t) ρ_t=0  .

Expanding the exponential one finds in leading order that the logarithmic dependence of 
Z'(1/t) and L̃_Δ(1/t) at high energies determine the short time behavior. 
This means that the RG equations are cut off at the large 
energy scale E=1/t, which is the poor man scaling regime, where all the cutoff scales z_n
are unimportant. In this regime the time evolution is determined by the scaling of Z'(1/t)
and L̃_Δ(1/t). If, in addition, t≫ 1/D, where D is the high-energy 
cutoff, one obtains universal time evolution in the short-time regime. It means that all 
leading logarithmic divergencies ∼ (αln(Dt))^n have been resummed in the 
functions Z'(1/t) and L̃_Δ(1/t), where α≪ 1 is some 
small dimensionless coupling parameter. Based on this unified picture the universal short-time 
behaviour has been derived in 
Refs. <cit.> for
the Kondo model, the IRLM, and the ohmic spin boson model, in accordance with similiar results
of previous literature. 


For intermediate and long times t≳ 1/|z_n|, we have to study the contributions from
the poles and branch cuts in detail, based on the decomposition (<ref>).
We start with the contributions from the branch cuts starting at a pole z_k^p, which we evaluate by using the form
(<ref>). For the branch cut integral
we set E=z_k^p-ix± 0^+ and replace in leading order λ_k(E)→ z_k^p 
and the logarithmic function P_k(E)Z'(E) by its average
P̅_k(z_k^p-ix)Z̅'(z_k^p-ix) over the branch cut, where 
A̅(E)=1 2(A(E+0^+)+A(E-0^+)). Furthermore, in leading order, we can 
use x→ 1/t in the logarithmic functions. This gives the result

    F_k^p(t) ≈ 1 2π ∫_0^-^∞ dx e^-xt (1 -ix+0^+-1 -ix-0^+) P̅_k(z_k^p-i/t) Z̅'(z_k^p-i/t)
      .

Using 1 -ix+0^+-1 -ix-0^+=2πδ(x), we obtain the following 
contribution to the total time evolution (<ref>)

    ρ_k^p(t) ≈ P̅_k(z_k^p-i/t) Z̅'(z_k^p-i/t) e^-iz_k^p t ρ_t=0  ,

i.e., for z_k^p=±Ω_k-iΓ_k, an exponential one with oscillation Ω_k
and decay rate Γ_k, modulated by a logarithmic function. 
For the special term k=st, where z^p_st=0, 
P_st(E)=|x_st(E)⟩Tr and P_st(E)Z'(E)=P_st(E), 
we get the following contribution to the time evolution

    ρ_st^p(t) ≈ |x_st(-i/t)⟩  ρ_st = |x_st(0)⟩  ,

i.e. we see that for t→∞ one always gets the stationary distribution
ρ_st but, if z_st^p is a branching pole, logarithmic corrections 
can occur for the time evolution which do not decay exponentially. 
We note that for the models discussed here, there is no logarithmic term in the 
diagrammatic series involving the pole z^p_st. In addition, there is no accidental 
pole z_kst^p=0, and therefore the pole at E=0 is isolated and has no attached branch cuts. 

The evaluation of a branch cut starting at a branching point z_n^b which is not a pole is more 
subtle since both λ_k(E) and P_k(E)Z'(E) can be discontinuous and cancellations can occur 
between the two contributions. Therefore, it is more convenient to start from the first expression
of (<ref>) involving the resolvent R(E). Denoting by 
δ A=A_+-A_- the jump across the branch and by A̅=1 2(A_++A_-) the
average value, with A_±=A(E± 0^+)=A̅±1 2δ A, one finds for the 
jump of the resolvent expanding in the small quantity δ L (leading to higher orders in the
renormalized coupling constants)

    δ R(E) = R_+ δ L R_- = 1 E-L̅ δ L 1 E-L̅ + O(δ L^3)  .

Using AB-A̅B̅=1 4δ Aδ B, we get

    1 E-L̅ = 1 E-L+O(δ L^2) = ∑_k1 E-λ_kP_k Z'+O(δ L^2) = ∑_k1 E-λ̅_kP̅_kZ̅'+O(δ L^2)

Inserting this in (<ref>), neglecting O(δ L^3), and approximating
E=z_n^b-ix→ z_n^b-i/t in the
logarithmic functions λ̅_k, P̅_k and Z̅', we get  
the following result for the branch cut integral

    F_n^b(t) ≈ 1 2π ∑_k k'z_k^p,z_k'^p z_n^b ∫_0^∞ dx e^-xt 1 z_n^b-ix-λ̅^n_k P̅^n_k Z̅^' n δ L(z_n^b-ix) 1 z_n^b-ix-λ̅^n_k' P̅^n_k' Z̅^' n ,

where λ̅_k^n=λ̅_k(z_n^b-i/t), P̅_k^n=P̅_k(z_n^b-i/t) and
Z̅^' n=Z̅'(z_n^b-i/t). We have omitted the cases z_k^p=z_n^b or z_k'^p=z_n^b since
we consider a branching point and not a branching pole. Since
λ̅_k^n≈ z_k^p, we can neglect x in the denominators of the resolvents
for times t∼ 1/x ≫ 1/|z_n^b-z_k,k'^p|. In this case, the long-time scaling 
is determined by the scaling of δ L(z_n^b-ix) for small x. Besides additional 
logarithmic corrections (which again can be treated by replacing x→ 1/t), 
we will show in Section <ref> that 

    δ L(z_n^b-ix) ∼ θ(x)

for models with charge fluctuations charge fluctuations (like the IRLM) and 

    δ L(z_n^b-ix) ∼  x θ(x)

for models with spin/orbital or energy fluctuations (like the Kondo and the ohmic spin boson model),
see Eq. (<ref>).
Therefore, if x can be neglected in the resolvents of the integrand of (<ref>),
we obtain (up to logarithmic corrections) ρ_n^b(t)∼ 1/t for charge fluctuations and
ρ_n^b(t)∼ 1/t^2 for spin/orbital and energy fluctuations. For special resonant cases, 
where z_n^b comes close to z_k^p or z_k'^p, one can also define time regimes 
1/|z_n^b|≲ t≪ 1/|z_n^b-z_k,k'^p|, where x dominates in the denominators for certain
values of k or k', leading to different scaling. If x is not neglected in (<ref>),
the integral can also be calculated exactly, leading typically to exponential integrals from which
the whole crossover behaviour from intermediate t∼ 1/|z_n| to long times t≫ 1/|z_n|
can be calculated. 

In the regime of intermediate to long times the cutoff scales z_n are very important. Each term
of the series (<ref>) has to be treated separately, leading to different
scaling of the individual terms (in contrast to the short-time regime, where all exponentials
can be approximated by one and only the sum of all pre-exponential functions matters). As we have
seen above, various functions K(z_n-i/t) with logarithmic scaling occur in the projectors, 
the Z-factors, and the jump of the Liouvillian. In bare perturbation theory, the logarithmic
functions K(E) will contain powers of terms ∼αlnD E-z_m. To get rid of
the high-energy cutoff D, a standard technique is to resum first all leading logarithmic divergencies
∼ (αlnDΛ_c)^n, where Λ_c≳|z_m| is some maximal physical
low energy scale. Technically, this can be achieved by cutting off the RG flow at Λ_c, 
defining renormalized coupling constants α_c at this point, and expanding the full
solution for |E|≲|z_m| in α_c. This is possible if α_c is small, i.e. if
Λ_c is much larger than the strong coupling scale Λ^*, where the coupling constants
become of O(1). As a result, K(E) will contain powers of 
logarithmic terms ∼α_clnΛ_c E-z_m. For E=z_n-i/t the most dangerous case
is n=m, leading to powers in the time-dependent parameter α_t∼α_clnΛ_c t. 
Since α_c≪ 1, this parameter is small α_t≪ 1, unless time is exponentially large.
Therefore, it can be treated perturbatively, leading to logarithmic corrections ∼α_t 
in the pre-exponential functions. This strategy has been 
used in Refs. <cit.> 
to determine the time evolution at
intermediate and long times (but not exponenitally large times) for the Kondo model, the IRLM,
and the spin boson model. 

Finally, the most complicated time regime is the one at 
exponentially large times, where α_t∼ O(1). In this regime, a perturbative treatment
is no longer possible and all powers of α_t are important. These logarithmic divergencies
at low energies are independent of those at large energies and can even arise if there is
no logarithmic divergence at high energies. Their occurence is related to the fact that, concerning
the time evolution, the final cutoff scale at low energies is set by inverse time 1/t and not by 
decay rates. The latter holds only for the calculation of stationary properties, see 
Refs. <cit.>.
The E-RTRG method is unique in the sense that it is also capable of resumming the logarithmic 
divergencies at low energies, provided the renormalized coupling constants remain small when E approaches
one of the singularities z_n. Recently, this has been achieved in a controlled way for the ohmic
spin boson model <cit.>, where deviations from previously predicted scaling behaviour have 
been found. Results for the Kondo model and the IRLM are still under investigation in this regime.
In particular for the Kondo model, the problem is that the renormalized coupling constants become
of O(1) when approaching one of the singularities although they might be small for the calculation
of stationary quantities at E=0. Thus, weak-coupling problems for stationary quantities can turn
into strong-coupling ones for the calculation of the long-time behaviour at exponentially large times. 




§ DIAGRAMMATIC EXPANSION



Effective Liouvillian. In this section we will derive a quantum field theoretical diagrammatic representation of the 
effective Liouvillian by expanding in the system-reservoir interaction V,
following Refs. <cit.>. Although this can
be done for the general case of an explicitly time-dependent Hamiltonian <cit.>, here
we will restrict ourselves to the more simpler case of a time-translational invariant
Hamiltonian. To find a diagrammatic expansion of the effective Liouvillian L(E) in Fourier space, we try to
bring the local density matrix ρ(E) into the form (<ref>). We start from the 
formal solution of the von Neumann equation (<ref>) for the total density matrix, use
the initial condition (<ref>), and obtain by expanding in the system-reservoir interaction

    ρ(E)    = ∫_0^∞ dt e^iEt Tr_res ρ_tot(t)
     = ∫_0^∞ dt e^iEt Tr_res e^-iL_tott ρ_tot(t=0)
        = Tr_res i E - L_tot ρ_t=0 ρ_res^(eq) = Tr_res i E - L^(0) - L_res - L_V ρ_t=0 ρ_res^(eq)
        = i Tr_res R^(0)(E-L_res) ∑_k=0^∞ 
    (L_V R^(0)(E-L_res))^k ρ_t=0 ρ_res^(eq)  ,

where we have defined

    R^(0)(E) = 1 E - L^(0)  ,  
    L^(0) = [H,·]  ,  
    L_res = [H_res,·]  ,  
    L_V = [V,·]  .

Using the form (<ref>) of the system-reservoir interaction, a similiar form can
be derived for the Liouville superoperator L_V

    L_V = 1 n! ∑_p=± G^(0)p… p_1… n 
    :A^p_1… A^p_n:  .

Here, p=± is the so-called Keldysh index, which indicates whether the interaction V arises
from the first or the second part of the commutator L_V b=Vb-bV (b is an arbitrary operator).
A^p_1 are reservoir field superoperators in Liouville space defined by 

    A_1^p b  = σ^p_res {[ a_1 b   p=+; b a_1   p=- ].  ,

and G^(0)p… p_1… n is
a superoperator acting in Liouville space of the local quantum system defined by 

    G^(0)p… p_1… n b = {[   1   n; σ^p   n ]} {[   g_1… n b        p=+; -b  g_1… n        p=- ].  .

σ^p and σ_res^p are convenient sign superoperators which account for
fermionic signs and measure the parity of the fermionic particle number difference N_s-N_s' 
of intermediate states |ss'⟩ in Liouville space via the definition 
(N_s denotes the particle number of state s and ± refers to bosons/fermions)

    σ^+ = 1   ,  σ^-_s_1 s_2,s_1^' s_2^' = δ_s_1 s_1^' δ_s_2 s_2^' 
    (±)^N_s-N_s'  ,

and a corresponding definition for σ^p_res by replacing local states s
by reservoir states. Since the total parity (local system plus reservoirs)
of all intermediate states must be even in Liouville space for fermions (note that it is 
impossible to prepare a nondiagonal matrix element of the total density matrix where the 
total fermionic particle number difference is odd, see Refs. <cit.>
for a detailed discussion and the consequences of this point), we obtain the important property

    σ^p σ^p_res = 1   .

From the definition of the reservoir field superoperators one can straightforwardly derive
how the product :A_1^p… A_n^p: occuring in Eq. (<ref>) acts in Liouville space

    :A_1^p… A_n^p: b = {[       1       n; σ^p_res       n ]} {[  :a_1… a_n: b           p=+; b  :a_1… a_n:           p=- ].  ,

i.e. similiar to G_1… n^(0)p… p but a minus sign is missing for p=-. Taking this
equation together with (<ref>) and using the property (<ref>),
one can easily prove the representation (<ref>) for L_V. 

Most importantly, the reservoir field superoperators are defined such that the usual Wick theorem
can be applied (see Ref. <cit.> for an elegant proof), i.e. the average 
Tr_res A_1^p_1… A_n^p_nρ_res
decomposes into a product of pair contractions and the sum has to be taken over all combinations,
with the usual definition of a fermionic sign to disentangle the various contractions. Using 
(<ref>) a single contraction is given by the expression

    γ_11'^pp'    = A^p_1 A^p'_1'(-20,11) 
        (-26,10)(0,1)3(-26,13)(1,0)16(-10,10)(0,1)3(20,11) 
       = Tr_res A^p_1 A^p'_1' ρ_res^eq
       = δ_11̅' {[  1; p' ]}  ρ_ασ(ω) f^p'η_α(ω)
     = δ_11̅' p' {[ η; 1 ]} ρ_ασ(ω) 
    f_α(p'ηω)  .


Using the form (<ref>) in (<ref>) one can shift all
reservoir field superoperators A_i^p_i to the right by using the
analog of the commutation relation (<ref>) in Liouville space

    A^p_1 L_res = (L_res - η(ω+μ_α)) A^p_1   .

This means that by shifting a certain field superoperator A_1^p through all 
resolvents to the right, we shift all reservoir Liouville operators L_res standing
right to A_1^p by -η(ω+μ_α), where 1≡ηασω. We note that, with
the second form (<ref>) of the interaction, there is no fermionic sign when commuting local and 
reservoir operators. Shifting all reservoir field superoperators to the right and using the notation

    X_1… n = ω̅_1… n + μ̅_1… n  ,  ω̅_1… n    = ω̅_1 + … ω̅_n  ,  ω̅_1 = η_1 ω_1  ,
    μ̅_1… n    = μ̅_1 + … μ̅_n  ,  μ̅_1 = η_1 μ_1  ,

we obtain for (<ref>) the form

    ρ(E)    = i ∑_k=0^∞ Tr_res R^(0)(E-L_res) 
    (1 n!G^(0)) R^(0)(E+X_M_1-L_res) ·
       · (1 n!G^(0)) R^(0)(E+X_M_2-L_res) … 
    (1 n!G^(0)) R^(0)(E+X_M_k-L_res) ρ_t=0 ·
       · (:A… A:) (:A… A:) … (:A… A:) ρ_res^(eq)
       = i ∑_k=0^∞ R^(0)(E) (1 n!G^(0)) R^(0)(E+X_M_1)
     … (1 n!G^(0)) R^(0)(E+X_M_k) ρ_t=0 ·
       · Tr_res {(:A… A:) … (:A… A:)} ρ_res^(eq)  ,

where we have used Tr_resL_res=0 in the last step. Thereby, the set
M_i includes those indices of reservoir field superoperators which were standing
left to the corresponding resolvent in the original expression. As a result
the local and reservoir degrees of freedom have been decoupled and the trace over the reservoir
degrees of freedom can be performed by the application of Wick's theorem in Liouville space.
Since all diagrams give the same contribution when the indices of a particular vertex 
G^(0)p… p_1… n
are permuted, the factor 1 n! is cancelled, except for the case when two vertices are
connected by m contractions, leaving a symmetry factor 1 m!.
This leads to a sum of diagrams which symbolically are translated by the rule

    ρ(E) → i (±)^N_p S (∏γ) 
    R^(0)(E) G^(0) R^(0)(E+X_M_1) … G^(0) R^(0)(E+X_M_k) ρ_t=0  ,

where ∏γ denotes the product over all contractions (<ref>),
N_p is the number of permutations of reservoir field superoperators to disentangle the
fermionic contractions, and S=∏_i m_i! is a symmetry factor arising for the case when
pairs of vertices are connected by m_i contractions. 

The determination of the shift variables X_M_i is simplified by noting that,
according to (<ref>), a single contraction γ_12^p_1 p_2
between A_1^p_1 and A_2^p_2 is only possible for η_1=-η_2, 
α_1=α_2 and ω_1=ω_2. This gives ω̅_12=μ̅_12=0, 
i.e. if the two indices fall both into the same set M_i, there is no contribution 
to the shift X_M_i. As a consequence, the left index 1 of a contraction γ_12^p_1 p_2
will contribute only to those resolvents, which stand between the two field operators A_1^p_1
and A_2^p_2 in the original series. For this reason, the last resolvent in (<ref>)
has no shift X_M_k=0 and is given by R^(0)(E). 

With the diagrammatic rules it is straightforward to translate a particular diagram, 
which we visualize as follows:

    ρ(E)    → i (10,10)
    (5,-20)
        < g r a p h i c s >
    ρ_t=0
    
       = i γ_12^p_1 p_2 γ_36^p_3 p_6 γ_45^p_4 p_5 
    R^(0)(E) G^(0)p_1_1 R^(0)(E_1+ω̅_1) G^(0)p_2_2 
    R^(0)(E)
    
       
    G^(0)p_3 p_4_34 R^(0)(E_34+ω̅_34) G^(0)p_5_5 
    R^(0)(E_3+ω̅_3) G^(0)p_6_6 R^(0)(E) ρ_t=0  ,

where we used the notation E_1… n=E+μ̅_1… n, see (<ref>).
In the diagrams, the green lines are the contractions, the circles denote the vertices, and the 
black lines connecting the vertices represent the resolvents R^(0) describing the dot propagation in 
Fourier space. The indices of the shift variables of a particular resolvent can
be determined by drawing a vertical line at the position of that resolvent and taking the left
indices of all contractions which cut through this line. We note that we do not distinguish
between diagrams which differ only by a permutation of the contractions connected to a certain
vertex, i.e. the permutation of the two green lines connected to the indices 3 and 4 in
the above example does not lead to a new diagram. 

To bring the density matrix ρ(E) into the form (<ref>) and to identify the
effective Liouvillian L(E), we note that each diagram consists of a sequence of connected
blocks, defined by the property that each vertical line will at least hit one contraction,
connected by resolvents R^(0)(E). E.g., the diagram (<ref>) consists
of a sequence of two blocks. Denoting the sum of all connected diagrams by Σ(E), the
diagrammatic series can be written as

    ρ(E)    = i {R^(0)(E) + R^(0)(E) Σ(E) R^(0)(E) + .
       .+ R^(0)(E) Σ(E) R^(0)(E) Σ(E) R^(0)(E) + …} ρ_t=0 = i E - L^(0) - Σ(E) ρ_t=0  .

Comparing to (<ref>), we see that the effective Liouvillian is given by

    L(E) = L^(0) + Σ(E)  ,

and Σ(E) consists of the sum of all connected diagrams with translation rule

    Σ(E) → (±)^N_p S (∏γ)_ con 
    G^(0) R^(0)(E_M_1+ω̅_M_1) … 
    G^(0) R^(0)(E_M_k+ω̅_M_k) G^(0)   ,

where (∏γ)_ con means that only connected diagrams are considered.
E.g. some of the lowest order diagrams of Σ(E) are given by

    Σ(E)    = (10,10)
    (5,-5)
        < g r a p h i c s >
     + (10,10)
    (5,-5)
        < g r a p h i c s >
     + (10,10)
    (5,-5)
        < g r a p h i c s >
     + …

Σ(E) is the dissipative part of the effective Liouvillian, which contains the whole 
information of the coupling to the reservoirs and leads to irreversible time evolution.
In time space we obtain L(t)=L^(0)δ(t-0^+)+Σ(t), such that the kinetic equation
(<ref>) reads

    iρ̇(t) = L^(0) ρ(t) + ∫_0^t dt' Σ(t-t') ρ(t')  .

The first term describes the von Neumann equation in the absence of the reservoirs, whereas
the second one is the dissipative part. The two terms are the analog of the
flow” and the collision” term of quantum Boltzmann equations.

Local observables. From the density matrix ρ(E) the time evolution of all averages of local observables 
can be calculated. The diagrammatic expansion can also be formulated for the calculation of
arbitrary observables containing reservoir degrees of freedom or correlation functions. E.g., if
an observable I of the generic form (<ref>) is taken

    I = 1 n! {[          1; η_1η_2…η_n ]}  :a_n a_n-1… a_1:  i_12… n → 1 n! i_12… n :a_1a_2… a_n:   ,

we define a corresponding superoperator L_I by the anticommutator

    L_I = i 2 [I,·]_+ = 1 n! ∑_p=± I^(0)p… p_1… n 
    :A^p_1… A^p_n:  ,

with 

    I^(0)p… p_1… n b = i 2 {[   1   n; σ^p   n ]} {[  i_1… n b       p=+; b  i_1… n       p=- ].  ,

such that the average can be written as

    ⟨ I ⟩(t) = Tr_tot I ρ_tot(t) = 
    -i Tr Tr_res L_I e^-i L_tott ρ_t=0 ρ_res^eq  .

This expression has a formal similiarity to 

    ρ̇(t)    = -i Tr_res L_tot e^-iL_tott ρ_tot(t=0)
     = -i Tr_res(L^(0)+L_res+L_V) e^-iL_tott ρ_t=0 ρ_res^eq
       = -i L^(0) ρ(t) - i Tr_resL_V e^-iL_tott ρ_t=0 ρ_res^eq  ,

where Tr_resL_res=0 has been used in the last line. Comparing
to the kinetic equation (<ref>), we find

    - i Tr_resL_V e^-iL_tott ρ_t=0 ρ_res^eq = -i ∫_0^t dt' Σ(t-t') ρ(t')  .

Therefore, when applying the same perturbative expansion to (<ref>), we obtain
the result

    ⟨ I ⟩(t) = -i ∫_0^t dt' Tr Σ_I(t-t') ρ(t')
      ,  ⟨ I ⟩(E) = -i Tr Σ_I(E) ρ(E)  ,

with the only difference that the first vertex of the kernel Σ_I(E) has to be the
vertex I^(0) instead of G^(0), i.e. the diagrammatic rule (<ref>) changes to 

    Σ_I(E)    → (±)^N_p S (∏γ)_ con
       
    I^(0) R^(0)(E_M_1+ω̅_M_1) 
    G^(0) R^(0)(E_M_2+ω̅_M_2) … 
    G^(0) R^(0)(E_M_k+ω̅_M_k) G^(0)   .

A prominent example for an observable is the particle current operator flowing from reservoir 
α into the local system defined by

    I_α = -d dt N_α = -i [H_tot,N_α]
     = -i [V,N_α]  .

Inserting the form (<ref>) of V, one finds after some straightforward manipulations 

    (i_α)_1… n = i ∑_k=1^n η_k δ_α_kα g_1… n  ,  
    (I_α)^(0)p… p_1… n = -1 2 ∑_k=1^n η_k δ_α_kα 
    p G^(0)p… p_1… n  .



Using the diagrammatic expansion (<ref>) and (<ref>), 
the kernels Σ(E) and Σ_I(E) can be calculated in perturbation
theory w.r.t. the bare vertices G^(0), and the local density matrix ρ(E) and the 
average ⟨ I ⟩(E) of any observable I follow from (<ref>) 
and (<ref>) in Fourier space. Using inverse Fourier transform the time evolution 
can finally be discussed following the strategy described in Section <ref>.
Stationary quantities are obtained from 

    ρ_st    = lim_t→∞ ρ(t)
     = lim_E→ 0^+ (-iE) ρ(E)
      ⇔   L(E=0^+) ρ_st = 0
    ⟨ I ⟩_st    = lim_t→∞ ⟨ I ⟩(t)
     = lim_E→ 0^+ (-iE) ⟨ I ⟩(E)
     = -i Tr Σ_I(E=0^+) ρ_st  .

Applications of these perturbative schemes for the calculation of transport properties will be discussed  
in the lecture C7 by M. Wegewijs. Similiar schemes have also been developed to
calculate correlation functions <cit.> and to consider explicitly
time-dependent Hamiltonians <cit.>. Concerning the latter first 
applications have considered adiabatic response <cit.> and quantum quenches 
<cit.> for the IRLM.

Analytic properties. From the perturbative expansion one finds that the effective Liouvillian L(E) has
a branch cut on the real axis and is analytic in the upper and lower half of the complex plane. This can
be seen from the resolvents R^(0)(E_M+ω̅_M) since L^(0)=[H^(0),·] is a self-adjoint
superoperator with real eigenvalues and all frequency variables ω̅ are integrated over the real axis. 
The same analytic property holds for the resolvent 
R(E)=1 E-L(E), since, due to (<ref>), we get for any initial density matrix ρ_t=0

    R(E) ρ_t=0    = -i ρ(E) = -i ∫_0^∞ dt e^iEt ρ(t) = 
    -i ∫_0^∞ dt e^iEt Tr_res ρ_tot(t)
        = -i ∫_0^∞ dt e^iEt Tr_res 
    e^-iL_tott ρ_tot(t=0) = Tr_res 1 E - L_tot ρ_tot(t=0)  .

This function can only have a branch cut on the real axis since 
L_tot=[H_tot,·] is a self-adjoint superoperator with real
eigenvalues. To calculate the time evolution we have seen from (<ref>) that the 
integration ∫ dE is slightly above the real axis and has to be closed in the lower half of 
the complex plane (due to t>0). It is very inconvenient to calculate this integral by enclosing 
the branch cut of the integrand on the real
axis due to the rapidly oscillating function e^-iEt in the integrand on the scale 1/t. Therefore,
in analogy to the standard procedure for response functions, one tries to find an appropriate analytic
continuation of the functions L(E) and R(E) into the lower half of the complex plane such that all branch cuts
point into the direction of the negative imaginary axis starting at certain singularities z_n with
Imz_n≤ 0. We achieve this in two steps. First, we will transform the perturbative series for L(E)
into a self-consistent equation by resumming all blocks of connected
diagrams on the propagators connecting the vertices. The diagrammatic representation allows this to
be done in a unique way and, as a result, the bare resolvents R^(0)(E_M+ω̅_M) are replaced 
by the full ones R(E_M+ω̅_M) and no diagrams are allowed with connected sub-blocks without any
free lines on the propagators which we indicate by (∏γ)_ irr

    Σ(E) → (±)^N_p S (∏γ)_ irr 
    G^(0) R(E_M_1+ω̅_M_1) … 
    G^(0) R(E_M_k+ω̅_M_k) G^(0)   .

For ImE>0, all resolvents are analytic functions w.r.t. the integration variables ω̅_i
in the upper half of the complex plane. Therefore, in the second step, we can close all 
integration contours in the upper half and have to enclose only the nonanalytic features arising 
from the spectral function ρ_ασ(ω) and the Bose/Fermi distribution 
f_α(p'ηω) in the contraction γ_11'^pp', defined in 
(<ref>). Thereby we assume that the frequency dependence of the vertices  
G^(0)p… p_1… n can be neglected. Decomposing the Bose/Fermi distribution
in symmetric and antisymmetric parts and using the representation in terms of the Masubara frequencies
ω_n^α=2nπ T_α (ω_n^α=(2n+1)π T_α) for bosons (fermions), 
we can write the contraction in the form 

    γ_11'^pp'    = δ_11̅' p' {[ η; 1 ]} ρ̅_ασ(ω̅)
     {∓ 1 2 + p' (f_α(ω̅) ± 1 2)} = δ_11̅' (p' γ_1^s + γ_1^a)  ,
    γ_1^s    = ∓ 1 2{[ η; 1 ]} ρ̅_ασ(ω̅)  ,
    γ_1^a    = {[ η; 1 ]} ρ̅_ασ(ω̅)
     {f_α(ω̅) ± 1 2}
       = {[ η; 1 ]} ρ̅_ασ(ω̅)  T_α 1 2 ∑_n (1ω̅-iω_n^α + 1ω̅+iω_n^α)  ,

where ρ̅_ασ(ω̅)=ρ_ασ(ω). Thus, after performing
all integrations ∫ dω̅_i, and assuming for the moment that the spectral 
function ρ̅_ασ(ω̅) is an analytic function in the upper half, 
the quantities ω̅_M occuring in the resolvents R(E_M+ω̅_M) will consist of 
a sum of positive Matsubara frequencies i∑_j∈ M|ω_n_j^α_j|. As a consequence, 
the analytic continuation w.r.t. E of this result for L(E) into the lower half of the 
complex plane will lead to nonanalytic features at E_M+i∑_j∈ M|ω_n_j^α_j|=z_k^p, 
where z_k^p are the poles of the
resolvent R(E) after the analytic continuation into the lower half of the complex plane. Since R(E) is
analytic in the upper half, the poles z_k^p have to lie in the lower half, and we find that L(E) has
an infinite series of poles in the lower half located at

    E = z_k^p - μ̅_M - i∑_j∈ M |ω_n_j^α_j|  ,

which, at zero temperature, turn into a series of branch cuts in the direction of the negative imaginary
axis with branching points z_n located at the poles z_k^p shifted by any combinations of the
chemical potentials of the reservoirs 

    z_n = z_k^p- μ̅_M  .

This result has formed the basis for the generic discussion of the time evolution in 
Section <ref>.

Influence of spectral function.
We note that the spectral function 
ρ̅_ασ(ω̅) of the models introduced in Section <ref>
does not change this picture. For quantum dots coupled to Fermi liquid leads, like the Kondo 
model or the IRLM, the spectral function 
ρ̅_ασ(ω̅)=D^2 D^2+ω̅^2 defines just a high-energy
cutoff function with pole at ω̅=iD in the upper half and residuum -iD/2. The contribution
of this pole to the frequency integration leads for D→∞ either to a vanishing or to
a regular contribution in E. For the ohmic spin boson model we get from (<ref>) that 
ρ̅(ω̅)=2α |ω̅|θ(ηω̅)D^2 D^2+ω̅^2,
which has a branch cut on the whole imaginary axis. However, since the vertex g_1=1 2σ_z 
is independent of η, we can sum the contraction γ_11'^pp' over η and η' at
fixed ω̅ and ω̅' (which are the integration variables) and get from
(<ref>) for the case of bosons the effective contraction

    γ_11'^pp' = δ(ω̅+ω̅') p' 2 α ω̅ D^2 D^2+ω̅^2 f(p'ω̅)
     = δ(ω̅+ω̅') (p'γ_1^s + γ_1^a)
      ,

with

    γ_1^s = -α ω̅ D^2 D^2+ω̅^2  ,  γ_1^s = α ω̅ D^2 D^2+ω̅^2 (2f(ω̅) + 1)  ,

where the index 1≡ω̅ contains only the frequency variable. Since ω̅ is
an analytic function, there is no change of the analytic structure of L(E) for the ohmic spin boson model.
For a generic frequency dependence of the spectral function, the analytic structure might change. If
ρ̅(ω̅) has a branch cut in the upper half in the direction of the positive 
imaginary axis starting at Δ_ρ + iγ_ρ, with γ_ρ≥ 0, the position 
(<ref>) of the branching points of L(E) can be shifted by multiples of
-Δ_ρ-iγ_ρ. This can e.g. happen for superconducting leads, where Δ_ρ
corresponds to the superconducting gap and γ_ρ=0. For sub- or super-ohmic
spin boson models there is no change of the analytic properties since the branch cuts of the
spectral function start at the origin. 

Symmetric part of the contraction. 
We note that the part γ_1^s of the contraction (<ref>) involving
the symmetric part of the Bose/Fermi distribution plays a special role. It is the only
part of the contraction γ_11'^pp' which depends on the Keldysh indices via p' and it depends on 
the frequency only via the spectral function. In particular for a spectral function of the form 
ρ̅(ω̅_1)=D^2 D^2+ω̅_1^2, i.e. if it just acts as a high-energy cutoff 
function but has no other special form, the frequency integration ∫ dω̅_1 will involve only
the pole of the spectral function at ω̅_1=iD when closed in the upper half. In the limit
D→∞ this means that this integration gives either zero (if more than one resolvent 
involves ω̅_1) or a constant if this contraction connects two consecutive vertices

    ∫ dω̅_1 D^2 D^2+ω̅_1^2 R(E_1… n+ω̅_1… n) = π D R(E_1… n+ω̅_2… n+iD)  -i π  .

As a result, the symmetric part of the contraction can be integrated out analytically and can be 
incorporated in an effective vertex by taking the two consecutive vertices together to a single one. 
The same can be shown for the ohmic spin boson model <cit.> due to its special algebra,
whereas for more general spectral functions with nonanalytic features in the upper half this is not 
the case. If it holds, one important consequence of this property is that the special pole 
at E_M+ω̅_M=z_st^p=0 of the resolvents R(E_M+ω̅_M) leads to regular contributions 
in the limit D→∞ and does not contribute to the branch cuts of L(E).
The reason is the special form P_st(E)=|x_st(E)⟩ Tr for the projector
of the mode k=st (see Eq. (<ref>)) together with the property

    ∑_p Tr G^(0)p… p_1… n = 0,

which follows straightforwardly from the definition (<ref>). Since the contractions
γ_11'^pp' are independent of the first Keldysh index p (see Eq. (<ref>)),
this means that if the projector P_st(E_M+ω̅_M) is inserted between two 
consecutive vertices, at least one 
of the contractions γ_11'^pp' associated with the right vertex must point into the left direction and 
only its p'-dependent symmetric part p'γ_1^s will contribute

    (10,10)
    (-80,-20)
        < g r a p h i c s >

As shown above this means that this contraction has to connect the two consecutive vertices and the 
frequency integration gives a constant. Thus, the case k=st does not contribute to the
positions (<ref>) of branch cuts for L(E) or R(E). There might be other 
accidental poles z_k^p=0 for k st, like e.g. for multi-channel Kondo models with non-Fermi
liquid behaviour <cit.>, 
but for the models discussed in Section <ref>, this is not the case. Therefore, for these
models, the pole at z_st=0 is isolated, as already stated in Section <ref>
after Eq. (<ref>).


Breakdown of perturbation theory for time-evolution problems.
Finally we note that it is very important to use the perturbative expansion of L(E) in the self-consistent 
form (<ref>) in order to find the right position of the branching
points z_n of the branch cuts of R(E). For the original series (<ref>) involving
the bare resolvent R^(0)(E_M+ω̅_M), the same considerations as above lead to
branch cuts of L(E) starting on the real axis at the value z_n=λ^(0)-μ̅_M,
where λ^(0) is a real eigenvalue of the bare Liouvillian L^(0). This would have a dramatic
effect on the long-time evolution because it leads to non-exponential decay. However, this result
is not correct since perturbation theory is very dangerous in the regime |E-z_n|∼Γ, where
Γ is a typical decay rate. At low frequencies, the resolvents can then become very large of the 
order of the inverse coupling constant, raising serious questions about convergence. In particular 
for the original series (<ref>) involving R^(0) this effect is most dramatic 
since a series of 
connected sub-blocks contains an arbitrary number of resolvents with exactly the same argument, i.e.
the singularity at low frequency appears to an arbitrary power. Such a series is certainly not 
convergent and it is necessary to resum it first to the self-consistent version 
(<ref>) before determining the position of the branching points. E.g.
consider a contribution to the effective Liouvillian of the form (α is some small dimensionless
coupling constant and Δ denotes a typical low-energy scale)

    α (-iE + αΔ) lnD -iE + αΔ = α (-iE + αΔ) lnD -iE - α^2 Δ - 1 2 i α^3 Δ^2 E + O(α^4)  .

The logarithm on the l.h.s. has a branching point at E=-iαΔ but the expanded form gives in 
O(α) and O(α^2) a branching point at E=0. The mistake can only be seen by 
considering higher orders in α, where an infinite series of terms 
with a pole at E=0 is obtained. Due to the factor
-iE+αΔ in front of the logarithm, this artifact is even not visible in O(α^2) but starts
in O(α^3) or higher. The form (<ref>) arises e.g. for the Kondo model and for
the ohmic spin boson model, where we will show in Section <ref> that 
∂∂ EL(E) must be a slowly varying logarithmic function leading to 
typical terms of the form (<ref>), see also Eq. (<ref>). 
For the ohmic spin boson model previous calculations <cit.> have 
predicted terms with non-exponential decay which have been corrected recently 
<cit.>.





§ RENORMALIZATION GROUP



General remarks. At low temperatures the perturbative calculation of the effective Liouvillian can 
break down even at small reservoir-system coupling for two reasons. First, at high energies
(the so-called ultraviolett regime), the frequency integrals are typically logarithmic leading 
to logarithmic contributions 
∼α^k(lnD E-z_n)^l in higher-order perturbation theory, with k≥ l, 
where z_n are the branching points (<ref>) of the resolvent R(E) and α is an
appropriate dimensionless coupling constant. Secondly, even if perturbation theory does
not contain ultraviolett logarithmic divergencies in the limit D→∞, it may contain
logarithmic terms ∼αlnE-z_n E-z_m, which, for E→ z_n, turn into 
the form ∼αlnE-z_n z_n-z_m, which can lead to a breakdown
of perturbation theory at low energies (the so-called infrared regime). Therefore, a method is
needed capable of reorganizing perturbation theory such that all ultravioltett and infrared 
logarithmic divergencies are resummed. Concerning high energies, resumming all logarithmic 
contributions ∼α^k(lnD E-z_n)^l with l=k,k-1,k-2,… is called 
leading order, sub-leading order, sub-sub leading order, etc. approximation (sometimes also referred to 
as 1-loop, 2-loop, 3-loop, etc.). In traditional (so-called poor man
scaling) RG methods <cit.>, one tries to perform this resummation by integrating out 
high-energy scales, i.e. the band width D is successively reduced in infinitesimal steps and the 
physical quantity of interest is kept invariant by renormalizing the coupling constants and other 
energy scales. Provided that the renormalized coupling constants remain small (the so-called 
weak-coupling regime), a well-controlled truncation scheme can be set
up by neglecting higher-order terms in the renormalized couplings. This strategy has also been used 
for calculating stationary quantities of nonequilibrium problems <cit.>, but it 
turns out that for the calculation of the effective Liouvillian L(E) it is very hard to set up a 
systematic truncation scheme. This has been improved by using a high-energy cutoff on the imaginary
axis by cutting off the Matsubara frequencies of the Bose/Fermi distribution function <cit.>
and applied to various models 
<cit.>.
Here we will follow another route by describing the E-RTRG method <cit.>, which is 
unique in the sense that it is capable of dealing with all logarithmic divergencies at high 
and low energies. Technically, 
this is achieved by considering the perturbation theory not for L(E) but for its first or second 
derivative w.r.t. the Fourier variable E together with a proper resummation in terms of effective vertices. 
Whether a first or a second derivative is needed depends on the model under consideration. 
This leads to a series where all frequency integrals converge at high energies and 
the limit D→∞ can be performed in all orders. As a consequence one obtains a universal 
differential equation (called RG equation) for L(E) independent of the specific 
choice of the high-energy cutoff function. Furthermore, the RG equation for L(E) turns out
to be such that the divergence at low energies for E→ z_n is at most ∼1 E-z_n
multiplied with a perturbative series in terms of effective vertices which exists in the limit E→ z_n.
This allows for a systematic solution at low energies as well. Besides the effective Liouvillian 
also effective vertices will appear in the RG equation due to the resummation procedure, for which 
similiar universal RG equations can be derived. Provided that the effective vertices stay small 
(so-called weak coupling problems) the RG equations can be systematically truncated and well-controlled 
universal properties can be determined at high as well as at low energies. According to the discussion
in Section <ref> this allows a well-controlled discussion of the time evolution
at short and long times together with the crossover behaviour. The high-energy cutoff D 
will only appear in the initial condition for the various quantities which are calculated by a 
well-controlled perturbation theory in the bare couplings at E=iD. This procedure has the 
advantage that by construction only the universal properties of the model are obtained, although 
it is also possible to keep D fixed and solve the RG equations for a given high-energy cutoff 
function. Furthermore, the use of a physical scale E as flow parameter of the RG equations has 
the advantage that at each stage of the flow the solution L(E) provides a result for a
physical quantity. Moreover, since E is a complex flow parameter, the flow can be solved on any path in 
the complex plane which is very helpful to find
appropriate analytic continuations of retarded functions into the lower half of the complex plane, even by 
using numerical methods, see also the discussion in Section <ref> after
Eq. (<ref>).

Derivation of the E-RTRG equations.
To illustrate the general strategy for the derivation of the RG equations within the E-RTRG method we consider 
here, for simplicity, a spectral function of the form 
ρ̅(ω̅)=D^2 D^2+ω̅^2, which arises typically for fermionic metallic
reservoirs where the d.o.s. is approximately a constant in the physically relevant energy regime.
Therefore, we consider only the fermionic case in the following.
Furthermore, we assume that only 1- and 2-point vertices occur in the original model and that
the frequency dependence of the bare vertices g_1 and g_12 can be neglected.
This applies to the Kondo model and the IRLM introduced in Section <ref>. For the ohmic spin
boson model, a similiar procedure can be used to derive the RG equations, see Ref. <cit.>.
For such models, one obtains a problem with convergence at high energies if the number of frequency 
integrations is larger or equal to the number of resolvents where the frequencies occur. For models with 
1- and 2-point vertices this means that diagrammatic sub-elements of the form

    (10,10)
    (-50,-0)
        < g r a p h i c s >
    (10,10)
    (-50,-0)
        < g r a p h i c s >
    (10,10)
    (-50,-0)
        < g r a p h i c s >
    
    
       (10,10)
    (-50,-0)
        < g r a p h i c s >
    (10,10)
    (-50,-0)
        < g r a p h i c s >
    (10,10)
    (-50,-0)
        < g r a p h i c s >

lead to problems at high energies and have to be avoided. This can be achieved by taking a single 
or a double derivative w.r.t. E of the resolvents occuring in these diagrams.
Therefore, the idea is to consider a perturbative expansion
for the derivatives ∂∂ EL(E) or ∂^2∂ E^2L(E) and to
resum the series such that no subelements of the form (<ref>) remain. 
The procedure is quite straightforward and we illustrate it for 
the case of a model where only 2-point vertices occur, like e.g. the Kondo model. Here, to 
guarantee convergence we consider two derivatives w.r.t. E of the diagrammatic series 
(<ref>) of L(E) in the self-consistent form. Since the 
E-dependence occurs only in the resovents R(E_M+ω̅_M), we can either take two
derivatives of a single resolvent or two single derivatives of different resolvents. Fixing the
positions of the resolvents we can then resum all remaining diagrams in a unique way such that
the bare 2-point vertices G^(0)pp_12 are replaced by full effective 2-point 
vertices G^p_1 p_2_12(E), which are defined as the sum of all connected diagrams 
with 2 external reservoir lines. With the convention that these two external lines are directed
to the right, it turns out that the energy argument of an effective vertex is identical to
the one of the preceding resolvent, i.e. only the combination 
R(E_M+ω̅_M)G^p_1 p_2_12(E_M+ω̅_M) can occur in the diagrammatic expansion.
Furthermore, for all diagrams
contributing to the effective vertex G_12^p_1 p_2(E), where the two external lines have the 
sequence 21, a fermionic sign has to be added. After this resummation, the diagrammatic series
for ∂^2∂ E^2L(E) up to third order in the effective vertices reads

    1/2∂^2/∂ E^2L(E) = 1/2 1/2 
    -0.5em
        < g r a p h i c s >
     + 
    -0.5em
        < g r a p h i c s >
    
     + 𝒪(G^4)  ,

where the red slash indicates a derivative ∂∂ E of the corresponding resolvent
(two slashes indicate the second derivative ∂^2∂ E^2). This is one of the
central equations in the E-RTRG approach.
Prefactors arsing from the symmetry factor 1 S have explicitly been indicated and all
vertices are full effective 2-point vertices from now on. All frequency integrations are 
convergent even if one neglects the frequency-dependence of the effective vertices (those can
only enhance convergence). Therefore on the r.h.s. of this differential equation we can take
the limit D→∞. This property holds in all orders since, by construction, all
diagrammatic subelements (<ref>) leading to a divergence in the infinite-D
limit have been eliminated by the resummation procedure. To close the equation one can also derive
in the same way a differential equation for the effective 2-point vertex 

    ∂/∂ EG^p_1p_2_12(E) = [
    -0.75em
        < g r a p h i c s >
     - (1↔ 2)]
    
        + 1/2 
    -1em
        < g r a p h i c s >
     + 1/2 
    -1em
        < g r a p h i c s >
    
        + [
    
    -1em
        < g r a p h i c s >
    
     + 
    -1em
        < g r a p h i c s >
    
     - (1↔ 2)]
     + 
    O(G^4)  .

After the limit D→∞ has been taken, the symmetric part (<ref>) of 
the contraction becomes an analytic function and does not contribute to the frequency integration when closing
the integration in the upper half of the complex plane. This means that the Keldysh indices no longer appear
explicitly in the RG equations, i.e. only the effective 2-point vertices averaged over the Keldysh 
indices are needed, which we denote by G_12(E)=∑_p_1 p_2G_12^p_1 p_2(E). This simplifies the
analysis considerably. As a result all contractions can be replaced by the antisymmetric part given by
(<ref>)

    γ_11'^pp' → γ_1^a = f_α(ω̅) - 1 2  ,

where we have already taken the limit D→∞ and integrated out the trivial
part δ_11̅ of all contractions in the RG diagrams. By convention, ω̅ is
always the frequency variable of the left vertex. 

Frequency dependence. 
To calculate the integrals over the internal frequencies in the RG diagrams, it is necessary to
know the frequency dependence of the effective vertices and the Liouvillian. This can be treated
systematically by the formalism. Provided that the bare vertices are frequency-independent, one
finds for the vertices that the diagrammatic series for the difference 
G_12(E)-G_12(E)_ω̅_1=ω̅_2=0 can be resummed by a similiar procedure in terms
of effective 2-point vertices such that the limit D→∞ is well-defined. The reason is 
that at least one resolvent in the original perturbative series must involve the difference 
R(E_M+ω̅_M+ω̅_M_ex)-R(E_M+ω̅_M), where M_ex contains 
some of the external indices {1,2}. Fixing this resolvent and resumming the rest of the diagram in
terms of effective 2-point vertices yields in lowest order the equation 

    -1.em
        < g r a p h i c s >
    
        = 
    -1.em
        < g r a p h i c s >
    
     + 
    -1.em
        < g r a p h i c s >
     - 
    -1.em
        < g r a p h i c s >
     + 𝒪(G^3)  ,

where the filled double dots represent the effective vertices at zero frequency.
This is the second key equation in the E-RTRG approach.
A contraction with an open circle and external frequency ω̅_i indicates that the
resolvent corresponding to the vertical cut at the position of that circle has to be replaced 
by the difference R(E_M+ω̅_M +ω̅_i)-R(E_M+ω̅_M). This difference falls off
∼(ω̅_M)^2 w.r.t. the internal frequency integration variables ω̅_M and, 
therefore, all frequency integrations are convergent in the limit D→∞. For 
the frequency dependence of the Liouvillian L(E_M+ω̅) it turns out that the similiar diagrammatic 
series for the difference L(E+ω̅)-L(E) does not exist in the 
limit D→∞, similiar to the fact that two derivatives are needed for convergence
(see above). Therefore, one defines a discrete version of the second derivative
Δ^2_ω̅L(E) via

    L(E+ω̅) = L(E) + ∂∂ EL(E) ω̅ + Δ^2_ω̅ L(E)
     = L(E) + ∂∂ EL(E) ω̅ + O(G^2)  ,

and finds that Δ^2_ω̅L(E) exists in the limit D→∞ and is 
at least of O(G^2) since it involves second and higher-order derivatives of the Liouvillian.
Neglecting O(G^2) (note that this contributes O(G^4) to the RG equations
(<ref>) and (<ref>)), the resolvents occuring in the RG diagrams 
of (<ref>) and (<ref>) can be written as

    R(E_M+ω̅_M) = 1ω̅_M + χ(E_M) Z(E_M) + O(G^2)  ,

with 

    χ(E) = Z(E) (E - L(E))  ,  
    Z(E) = 1 1 - ∂∂ EL(E)  ,
 
where the RG equation for Z(E) follows from the one for ∂^2∂ E^2L(E) by

    ∂∂ EZ(E) = Z(E) {∂^2∂ E^2L(E)} Z(E)
     ∼ O(G^2)  .

Inserting (<ref>) and (<ref>) into the RG equations (<ref>) and
(<ref>), calculating all frequency integrations and neglecting all terms of O(G^4), 
one obtains a closed set of RG equations for G_12(E)|_ω̅_1=ω̅_2=0 and L(E),
which can be easily solved numerically. These constitute the basic equations of the E-RTRG approach.
The crucial step in the formalism is the parametrization of the frequency dependence,
otherwise a numerical solution would be very
time consuming. Truncating the RG equations at O(G^2) provides the solution up to leading order,
whereas a truncation at O(G^3) includes in addition all sub-leading terms. An important check for the
reliability of the solution is whether these two truncation schemes lead approximately to the same universal
solution. For the nonequilibrium Kondo model at zero magnetic field, the equations have been solved 
in Ref. <cit.> to calculate the stationary conductance with reliable results even
in the strong coupling regime. Similiar RG equations can be 
set up for the IRLM and the spin boson model which have been studied in Refs. <cit.>.


RG equations for the slowly varying parts of the Liouvillian.
We are now ready to show how the decomposition (<ref>) can be derived together with 
RG equations for the slowly varying functions L_Δ(E) and L'(E). First of all, one can see from
the RG equations (<ref>) and (<ref>) that ∂∂ EL(E) and
G_12(E) are slowly varying logarithmic functions. At large E we find from dimensional arguments
that 

    ∂^2∂ E^2L(E) , ∂∂ EG_12(E) ∼ 1 E (1 + O(Δ E))  ,

where Δ is some physical scale except E. For large E, we can neglect the higher orders 
∼ O(Δ E) and we see that, due to the factor 1 E, logarithmic functions are generated 
by integrating over E. For E close to some branching point z_n, we find, that even in the worst case
when all resolvents contain the same branching point, that 
∂^2∂ E^2L(E) and ∂∂ EG_12(E) can at most diverge 
∼1 E-z_n for E→ z_n. As a result, also for E→ z_n, 
∂∂ EL(E) and G_12(E) are slowly varying logarithmic functions of E-z_n. 
This can only be the case if L(E) consists of terms 

    L(E) ∼  (E-z_n) K_n(E-z_n) = -z_n K_n(E-z_n) + E K_n(E-z_n)  ,

where K(E) is a slowly varying function, or, more precisely, z_n=z_k^p-μ̅_M will be replaced 
by λ_k(E_M)-μ̅_M if E is not close to one of the singularities. Therefore, we see
that L(E) can be decomposed in the form (<ref>),

    L(E) = L_Δ(E) + E L'(E)  ,

with slowly varying functions 
L_Δ(E)∼ -z_n K_n(E-z_n) and L'(E)∼ K_n(E-z_n). We note that we used precisely 
this form at the end of Section <ref> in Eq. (<ref>).
It shows that L_Δ(E) and L'(E) have a quite similiar structure. 

We note that the property that ∂∂ EL(E) and G_12(E) are slowly varying 
logarithmic functions can also be seen directly from the original perturbative expansion 
(<ref>) since in all orders of perturbation theory the number of
frequency integrations is identical to the number of resolvents. This leads to logarithmic integrals
at large and low energies even if all resolvents contain the same cutoff scale at low energies. For
the proof it is essential that the perturbation theory is taken in the self-consistent form 
(<ref>) since this leads to the property that all resolvents 
involve a different combination of the frequencies. The same can be shown for the ohmic 
spin boson model where the decomposition (<ref>) holds also in all orders 
of perturbation theory. For models with 1-point vertices and a flat spectral function 
(like e.g. quantum dot models in the charge fluctuation regime), the number of resolvents 
can be arbitrarily larger than the number of frequency integrations. Here, to show the 
logarithmic scaling at low energies in all orders of perturbation theory, it is very important that the 
resolvents do not only have different frequency combinations but many of them have also different 
cutoff scales at low energies. In contrast to models with spin/orbital fluctuations, it turns out
that already the first derivative ∂∂ EL(E) exists in the limit 
D→∞, see e.g. the first diagram of (<ref>). This means that
∂∂ EL(E)∼Γ E-z_n multiplied with 
a well-controlled series with no divergence at high or low energies. This part influences
only the function L_Δ(E) but not E L'(E). The systematic treatment of all orders
in the tunneling for models with charge fluctuations is still an issue of ongoing research.

To find RG equations for L_Δ(E) and L'(E), we try to bring the 
RG equation (<ref>) for ∂^2∂ E^2L(E) into the form

    ∂^2∂ E^2L(E) = ∂∂ EL'(E)
     + ∂∂ E{∂∂ EL_Δ(E)
     + E∂∂ EL'(E)}  ,

such that ∂∂ EL_Δ(E) and ∂∂ EL'(E) can be 
identified and that L_Δ(E) is proportional to some physical scale Δ except E.
For simplicity we show the procedure only up to O(G^2), for O(G^3) see Ref. <cit.>.
Taking only the first term on the r.h.s. of the RG equation (<ref>), replacing the
vertices by the ones at zero frequency via (<ref>), and shifting the two derivatives of
the resolvent 
∂^2∂ E^2R(E_12+ω̅_12)=∂∂ω̅_1∂∂ω̅_2R(E_12+ω̅_12) via two partial integrations
to the contractions, we obtain

    ∂^2∂ E^2L(E)     =   1 2
    (10,10)
    (5,-7)
        < g r a p h i c s >
    
    
    +  ∂∂ E{1 2
    (10,10)
    (5,-7)
        < g r a p h i c s >
    
    }  +   O(G^3)  ,

where a cross at a contraction denotes the derivative 
∂∂ω̅γ_1^a=∂∂ω̅f_α(ω̅),
see (<ref>). The dashed line in the second term indicates that
the resolvent is replaced by the Z'-factor R(E_M+ßω̅_M)→ Z'(E_M+ω̅), 
defined in (<ref>). Therefore, this term is of O(G^3) and can be added without violating
the consistency of the truncation scheme up to O(G^2). The term has been added in such a way that 
when identifying (<ref>) with (<ref>), the derivative 
∂∂ EL_Δ(E) will become proportional to a physical scale Δ. 
Together with the relation

    Z'(E_M+ω̅_M) - E R(E_M+ω̅_M)
        = χ_Δ(E,μ̅_M+ω̅_M) R(E_M+ω̅_M)  ,
    χ_Δ(E,μ̅_M+ω̅_M)    = μ̅_M + ω̅_M - L̃_Δ(E_M+ω̅_M)  ,

which follows from (<ref>) with the definition L̃_Δ(E)=Z'(E)L_Δ(E), we obtain

    ∂∂ E L_Δ(E) = 1 2(10,10)
    (5,-12)
        < g r a p h i c s >
     + O(G^3)
      ,  ∂∂ E L'(E) = 1 2 (10,10)
    (5,-7)
        < g r a p h i c s >
     + O(G^3) ,

where the symbol χ_Δ at the resolvent means that the resolvent multiplied
with χ_Δ(E,μ̅_M+ω̅_M) has to be taken. Obviously, χ_Δ is
proportional to a physical scale, since μ̅_M, ω̅_M and
L̃_Δ(E_M+ω̅) have this property. For ω̅_M this follows
from the fact that the RG equations contain only the derivatives 
∂∂ω̅γ_1^a=∂∂ω̅f_α(ω̅)
of the contractions, such that |ω̅_i|≲ T_α_i. In contrast to the 
RG equation (<ref>) for the full Liouvillian L(E), the RG equations (<ref>)
for L_Δ(E) and L'(E) are first order differential equations. Therefore, the 
differences L_Δ(E_M+ω̅_M)-L_Δ(E_M) and L'(E_M+ω̅_M)-L'(E_M) are
of O(G^2) such that the frequency dependence of the resolvent and χ_Δ entering the
RG equations (<ref>) can be approximated by

    R(E_M+ω̅_M)    = 1ω̅_M + E_M - L̃_Δ(E_M) Z'(E_M)
     + O(G^2)
    χ_Δ(E,μ̅_M+ω̅_M)    = μ̅_M + ω̅_M - L̃_Δ(E_M) + O(G^2)  .

As a consequence, all frequency integrations can be straightforwardly performed such that the
differential-integro equations (<ref>) are converted into differential equations.
E.g., at zero temperature, the two frequency integrations in (<ref>)
are trivial leading to the explicit expression

    ∂∂ EL_Δ(E)    = 1 2 G_12(E) μ̅_12 - L̃_Δ(E_12)  E_12 - L̃_Δ(E_12) 
    Z'(E_12) G_2̅1̅(E_12)  ,
    ∂∂ EL'(E)    = 1 2 G_12(E) 1  E_12 - L̃_Δ(E_12) 
    Z'(E_12) G_2̅1̅(E_12)  ,

together with the RG equation for the vertex which follows from the lowest order term of
(<ref>) as

    ∂∂ EG_12(E) = 
    G_13(E) 1  E_13 - L̃_Δ(E_13) Z'(E_13) G_3̅2(E_13)
     - (1↔ 2)  .



Solution of approximate E-RTRG equations.
The first-order RG equations (<ref>) and (<ref>) for
L_Δ(E) and L'(E) provide the most convenient starting point 
for an analytical solution of the RG equations at least in that regime of the complex plane where 
the effective vertices stay small, see Refs. <cit.> for details. 
The strategy is to solve the RG equations approximately in three different energy regimes by expanding in the
effective vertices but keeping large logarithmic terms (either at large or low energies) to all orders, and 
matching the different solutions to fix the integration constants. Denoting the small dimensionless 
coupling constant by α, we distinguish the following regimes: (1) The regime of high energies
|E|≫|z_n|, where the RG resums all ultraviolett logarithmic terms ∼ (αlnD -iE)^k; 
(2) The regime of intermediate and small energies |E-z_n|≲ O(|z_n|) but E not too close to
the branching points such that one can expand in the small parameter 
α|ln|z_n||E-z_n||≪ 1; (3) The regime of small energies exponentially close to some 
of the branching points, i.e. |E-z_n|≪ O(|z_n|) and α|ln|z_n||E-z_n||∼ 1, where
the RG resums all infrared logarithmic terms ∼(α|ln|z_n||E-z_n||)^k. In particular for
the ohmic spin boson model and the IRLM, we will see in 
Section <ref> that the coupling constant α stays
small in the whole complex plane such that a well-controlled analytical solution is possible for all E,
showing that the resummation of logarithmic terms for high and low energies gives very different results.
For the Kondo problem a weak-coupling solution is only possible for high, intermediate and small energies,
but not for exponentially small energies where the coupling constant α∼ O(1).

Initial conditions.
The initial conditions for the RG flow
at large energies are set up at the value E=iD, where D≫|z_n| is the high-energy cutoff.
The motivation for the choice E=iD lies in the fact that, for D≫|E|≫|z_n|, the bare perturbation 
series for L_Δ(E) and L'(E) contain logarithmic terms ∼ (lnD -iE)^k of all
powers k (we have chosen -iE in the argument, such that the branch cut is directed towards 
the negative imaginary axis). All other terms ∼ (|E| D)^n are neglected since they vanish in the
limit D→∞ and thus do not contribute to the universal solution
which is independent of the cutoff details. Extrapolating
this result up to E=iD has the effect that all logarithmic terms vanish, which sets the initial point
for the universal RG flow. The calculation of the initial values can be done by bare perturbation theory
for D≫|E|≫|z_n| and omitting all logarithmic contributions. For small bare coupling constants it is
sufficient to take the lowest order term if it is universal, otherwise one takes zero for the initial 
condition. We note that this procedure works well to determine the universal initial condition for  
L_Δ(E) and L'(E) at E=iD but fails for the initial condition of L(E) since L(E) contains
terms linear in E which are very large for E=iD. Therefore, for the RG equation (<ref>) one
either has to keep the high-energy cutoff function in the RG equations and start the RG flow at 
|E|≫ D (where one can take the bare values as initial condition), or one has to find a reference
point at low energies where L(E) is known from exact results, see e.g. the solution of the Kondo 
model in strong coupling in Ref. <cit.>.

RG for the Liouvillian discontinuity jumps.
Finally we show how RG equations can be derived for the jump δ L of the Liouvillian at a particular 
branch cut, as this jump is needed to evaluate the branch cut contributions to the
time evolution, see Eq. (<ref>). As described in Section <ref>
the branch cuts of L(E) occur only at zero temperature and can be identified in the perturbative expansion 
(<ref>) by closing all frequency integrations in the upper half of the
complex plane and considering the branch cuts of the Fermi distribution functions on the positive imaginary
axis. This means that the frequencies are shifted to the positive imaginary axis
ω̅_M→ i|ω̅_M|. In leading order, a given branch cut at 
E=z_n-ix± 0^+ is generated by some resolvent which is resonant, i.e. the jump of this 
resolvent across the branch cut becomes a δ-function. If the resolvent contains the 
eigenvalue λ_k(E_M+i|ω̅_M|) the resonance occurs if z_n=z_k^p-μ̅_M.
With E_M=z_n+μ̅_M-ix± 0^+=z_k^p-ix± 0^+, we replace approximately 
λ_k→ z_k^p, P_k→P̅_k(z_k^p-ix) and 
Z'→Z̅'(z_k^p-ix), which gives for the jump of the resolvent the 
following δ-function

    δ R    =    (1 -ix+iω̅_1… n+0^+-
    1 -ix+i|ω̅_M|-0^+) P̅_k(z_k^p-ix) Z̅'(z_k^p-ix)
    
       =    2πδ(|ω̅_M|-x) P̅_k(z_k^p-ix) Z̅'(z_k^p-ix)  .

As expected, the frequency integrals give only a contribution for x>0, since this is
the region where the branch cut starts.
The RG equation for ∂δ L∂ E(z_n-ix) is obtained by fixing
the resonant resolvent together with the resolvent where the E-derivative is taken 
and resumming the rest of the perturbative series in terms of effective 2-point vertices.
If both resolvent are the same, the E-derivative is replaced by a frequency derivative
∂∂ω̅_i and is shifted via partial integration to the 
derivative of some contraction crossing over the resolvent. Thus, we obtain in leading order

    ∂δ L∂ E(z_n-ix) = 
    i∂∂ xδ L(z_n-ix)  =  - 1 2
    (10,10)
    (5,-7)
        < g r a p h i c s >

where the symbol δ at the resolvent means that we replace the resolvent by its jump
δ R given by (<ref>). The initial condition for the RG equation is 
δ L(z_n)=0. At zero temperature both frequency integrations are trivial. The contraction
with the cross gives 
∫ dω̅f_α^'(ω̅){…}=-{…}_ω̅=0,
whereas the contraction without the cross leads to

    ∫ dω̅ (f_α(ω̅)-1 2) {…} = 
    -1 2 ∫ dω̅ sign(ω̅) {…} = 
    -i ∫_0^∞ dω̅ {…}_ω̅→ iω̅  ,

where we have used that the jump of the sign-function at the branch cut is given by 2.
Therefore, (<ref>) gives explicitly

    ∂∂ xδ L(z_n-ix) = 
    -π θ(x) G̅_12(z_k^p-μ̅_12-ix) P̅_k(z_k^p-ix) Z̅'(z_k^p-ix) G̅_2̅1̅(z_k^p-ix)  ,

where we have replaced all vertices that are discontinous across the branch cut by their
average value G̅. Up to the corrections from the weak x-dependence of the 
logarithmic functions on the r.h.s. of this equation, we obtain 

    δ L(z_n-ix) ∼  x θ(x)  ,

giving Eq. (<ref>) used in Section <ref>.
Therefore, if x can be neglected in the 
resolvents of the integrand of (<ref>),
we obtain (up to logarithmic corrections) ρ_t^n,b∼ 1/t^2 in the long-time limit
for all models with spin or orbital fluctuations, like e.g. the Kondo model. Similiar
considerations show that the same holds for the ohmic spin boson model, whereas for models
with charge fluctuations (like the IRLM), one obtains δ L(z_n-ix)∼θ(x).





§ RESULTS



In this section we will discuss the application of the formalism to the models introduced
in Section <ref>. Since the Kondo model has the simplest algebra, we will take
this model as a tutorial example to discuss the solution of the RG equations and the consequences
for the time evolution in all detail. Since the general strategy is always the same, we will
then show briefly the results for the ohmic spin boson model and the IRLM, and will
concentrate on interesting features which are different from the ones for the Kondo model.




 §.§ Kondo model


We consider the nonequilibrium Kondo model at zero magnetic field h^(0)=0 and zero
temperature T=0 for the antiferromagnetic case J^(0)>0. We assume that the 
local spin-1 2 is coupled to several 
reservoirs with chemical potentials μ_α. For the special case of two reservoirs
α≡ L/R≡±, we take μ_α=αV 2, where V denotes
the bias voltage across the system. Following Refs. <cit.> 
our aim is to calculate the time evolution of the local spin ⟨S⟩(t). 


The model is spin rotational invariant and therefore the effective Liouvillian L(E) should 
be an invariant under spin rotations. Defining two basis spinoperators L_± in 
Liouville space by (A is an arbitrary local operator) L^+ A=SA and
L^- A=-AS, the only two invariants are given by the identity and
L^+·L^-. Since the Liouvillian must also fulfil TrL(E)=0, we
find that the Liouvillian can be parametrized by

    L(E) = -i Γ(E) L^a   ,   L^a = 3 4 + L^+·L^-   .

Γ(E) is the energy dependent spin relaxation rate. The Liouvillian
has one zero eigenvalue with projector P_st=1-L^a and three degenerate
eigenvalues at -iΓ(E) with projector 1-P_st=L^a.
Therefore, by using (<ref>), we can write for the time evolution of the local
density matrix

    ρ(t) = ρ_st + i 2π ∫_-∞+i0^+^∞+i0^+dE 
    e^-iEt 1 E + iΓ(E) L^a ρ_t=0  ,

where ρ_st=(1-L^a)ρ_t=0 is the diagonal stationary density matrix with
equal probabilities for both spin directions (ρ_st)_ss=1 2. Using
TrSρ(t)=TrSL^aρ(t) and (L^a)^2=L^a, we find
that the spin relaxation rate Γ(E) determines the spin dynamics via

    ⟨S⟩(t) = i 2π ∫_-∞+i0^+^∞+i0^+dE 
    e^-iEt 1 E + iΓ(E) ⟨S⟩_t=0  .

As a consequence, the operator structure of the Liouvillian is no longer important, and we can
use all formulas derived in Section <ref> for the time evolution with the 
replacement L(E)→ -iΓ(E) or R(E)→1 E+iΓ(E). 
This means that all projectors P_k can be left out and Z'(E) can be used from the decomposition

    Γ(E) = Γ_Δ(E) + E Γ'(E)  ,   
    R(E) = 1 E + iΓ(E) = 1 E + iΓ̃_Δ(E) Z'(E)  ,

with 

    Γ̃_Δ(E) = Z'(E) Γ_Δ(E)  ,   
    Z'(E) = 1 1 + iΓ'(E)  .

Analysing the functions Γ̃_Δ(E) and Z'(E) in leading order from the 
RG equations (<ref>) and (<ref>), we will show below 
that the branching poles and branching points of the resolvent R(E) are given by

    z_0^p = -iΓ^*   ,   
    z_αα'^b = -iΓ^* + μ_α - μ_α'  for  αα'

where we used the notation of Eq. (<ref>) and  assumed
that all reservoirs have different chemical potentials (otherwise they can be taken together). 
z_0^p=iΓ̃_Δ(z_0^p) is the pole of the resolvent R(E) and
Γ^* is called the Korringa rate Korringa rate which, as will be shown below, is given by

    Γ^* = 2π J_V^2 ∑_αα' x_α x_α' |μ_α-μ_α'|
      ,   J_V = 1 2 lnV T_K  ,   T_K = D e^-1/(2J^(0))  ,

where we have used the notation of Eq. (<ref>). J_V is the renormalized exchange coupling 
at the scale V=max_αα'{|μ_α-μ_α'|} 
(which is the bias voltage for two reservoirs) and T_K is called the Kondo temperature
Kondo temperature. The result
for the Korringa rate holds in the weak-coupling case V≫ T_K which we will consider from now on.
In this case, we get J_V≪ 1 and Γ^*≪ V. 
Using (<ref>) and (<ref>), we get 
the following general form for the time evolution of the local spin 

    ⟨S⟩(t) = 
    F_0^p(t) e^-Γ^* t ⟨S⟩(t)_t=0 + ∑_αα' F_αα'^b(t) e^-Γ^*t e^-i(μ_α-μ_α')t ⟨S⟩(t)_t=0  .

For intermediate and long times t≳1 V, the pre-exponential functions can be calculated from
(<ref>) and (<ref>) as

    F^p_0(t)    = Z̅'(-iΓ^*-i/t)  ,
    
    F^b_αα'(t)    = -i 2π ∫_0^∞ dx e^-xt Z̅'(z_αα'^b-i/t)^2 δΓ(z_αα'^b-ix)
    (μ_α - μ_α' - ix)^2  ,

where we have neglected -iΓ^*+iΓ̅̃̅_Δ(z_αα'^b-ix)
in the dominator of the integrand of the last equation. This can be done for 
|μ_α-μ_α'|≫Γ^* for all αα', i.e. if the branch cuts 
are sufficiently apart from each other (other cases can be treated as well but need a special 
procedure <cit.>). For short times t≪1 V, we can use 
(<ref>) and get

    ⟨S⟩(t) = 
    Z'(1/t) e^-Γ̃_Δ(1/t) t ⟨S⟩(t)_t=0  .

To evaluate (<ref>), (<ref>) and (<ref>) explicitly,
we need the functions Γ̃_Δ(E), Z'(E) and δΓ(z_αα'-ix),
which we will derive in the following by considering the RG equations
(<ref>), (<ref>) and (<ref>),
together with the RG equation (<ref>) for the vertex.

In leading order it can be shown that the effective vertex G_11'(E) can be parametrized in the
same form as the initial vertex G^(0)_11', defined by (<ref>) and 
(<ref>). This gives the form

    G_11'(E) = - J_αα'(E) L^2·σ_σσ'  for  η=-η'=+  ,  L^2 = -1 2 (L^+  + L^-)   ,

together with G_11'(E)=-G_1'1(E) for η=-η'=-. Using this ansatz together with
the form (<ref>) for the Liouvillian in the RG equations and omitting the terms
∼ L^a in the RG equation for the vertex (which generate higher orders), one finds after
some straightforward algebra <cit.> the following RG equations

    ∂∂ EΓ_Δ(E)    = i χ_Δ(E,μ̂_αα') 
    R(Ê_αα') J_αα'(E) J_α'α(Ê_αα')  ,
    ∂∂ EZ'(E)    = Z'(E)^2 
    R(Ê_αα') J_αα'(E) J_α'α(Ê_αα')  ,
    ∂∂ EJ_αα'(E)    = -1 2 
    R(Ê_αα”) J_αα”(E) J_α”α'(Ê_αα”)
     - 1 2 
    R(Ê_α”α') J_α”α'(E) J_αα”(Ê_α”α)
      ,

where Ê_αα'=E+μ̂_αα', 
μ̂_αα'=μ_α-μ_α', 
χ_Δ(E,μ̂_αα')=μ̂_αα'
+iΓ̃_Δ(Ê_αα'), and Z'(E) and R(E) have been defined 
in (<ref>) and (<ref>). The initial conditions
at E=iD are Γ_Δ=0, Z'=1 and J_αα'=2√(x_α x_α')J^(0).

We first start with the analytic solution in the regime of high energies |E|≫ V. Neglecting 
μ̂_αα' everywhere gives the solution Γ_Δ(E)≈ 0 and
J_αα'(E)≈ 2√(x_α x_α')J(E), together with

    ∂∂ EZ'(E) = 4 E Z'(E)^2 J(E)^2
     = 4 E J(E)^2 + O(J^3)  ,
    ∂∂ EJ(E) = -2 E Z'(E) J(E)^2
     = -2 E J(E)^2 + O(J^3)  ,

where we have used Z'= 1+O(J) on the r.h.s. (see Eq. (<ref>)). We find that
-1 2 J(E)+ln(-iE)≡lnT_K is an invariant and ∂ Z'∂ J= -2, which
gives the solution

    Z'(E) = 1 - 2 (J(E) - J^(0)) → 1 - 2 J(E)   ,  
    J(E) = 1 2ln-iE T_K  ,

where the Kondo temperature T_K is defined in (<ref>) and we used the scaling limit 
in the first equation, defined by D→∞, J^(0)→ 0, such that
T_K remains a constant. We have chosen -iE in the argument of the logarithm to define a real value
for the Kondo temperature at E=iD and since we want the branch cut of the logarithm to point into the
direction of the negative imaginary axis. In the solution (<ref>) all
logarithmic terms ∼ (J^(0)lnD -iE)^k have been resummed.
J(E) is the poor man scaling solution of the Kondo model,
already introduced in the lecture B3 by T. Costi, but with the difference that -iE plays now the
role of the effective energy scale. Most importantly, the solution would diverge at E=iT_K when
extrapolated to small energies, indicating an increase of antiferromagnetic spin fluctuations at the
scale of the Kondo temperature. However, in this regime the solution can not be used since the 
RG flow becomes very different for |E|≲ V. The solution at high energies can be used to
evaluate the universal solution (<ref>) for short times t≪1 V

    ⟨S⟩(t) = (1 - 1|ln(T_K t)|) ⟨S⟩(t)_t=0  .

In this result all logarithmic terms ∼ (J^(0)ln(Dt))^k have been resummed, which can be
seen from 1/ln(T_K)=-2J^(0)/(1-2J^(0)ln(Dt)). Sub-leading terms are not included but
can be taken into account by truncating the RG equations at O(G^3) <cit.>.
For ferromagnetic Kondo models the universal short time behaviour has also been discussed 
in Ref. <cit.> using flow equation methods.


To find the solution at intermediate and small energies |E|≲ V but not exponentially close to the
singularites such that J(E)≪ 1 is still fulfilled (we state below what this precisely means),
we first set the initial value by expanding the solution (<ref>) at high energies
for the case when E starts to approach V such that |E|≫ V is still fulfilled 
but J_V|ln-iE V|≪ 1, where J_V=1/(2ln(V/T_K)) is the exchange coupling at high energies 
evaluated at the scale -iE=V, as introduced in (<ref>). This gives

    J_αα'(E) ≈ 2√(x_α x_α') J_V (1 - 2 J_V ln-iE V)  ,  
    Z'(E) ≈ 1 - 2 J_V + 4 J_V^2 ln-iE V  .

Using the first term of this expansion in the r.h.s. of the full RG equations
(<ref>), (<ref>) and (<ref>), we can 
easily integrate the equations up to O(J_V^2) by using
R(E)=∂∂ Eln-i(E+iΓ̃_Δ(E)) V + O(J_V^2)

    Γ̃_Δ(E)    = i 4 x_α x_α' J_V^2 
    (μ̂_αα' +iΓ̃_Δ(Ê_αα')) ln-i(Ê_αα'+iΓ̃_Δ(Ê_αα')) V  ,
    
    Z'(E)    = 1 - 2 J_V + 
    4 x_α x_α' J_V^2 ln-i(Ê_αα'+iΓ̃_Δ(Ê_αα')) V  ,
    
    J_αα'(E)    = 2√(x_α x_α') J_V {1 - 
    x_α” J_V ln-i(Ê_αα”+iΓ̃_Δ(Ê_αα”)) V.
       . - 
    x_α” J_V ln-i(Ê_α”α'+iΓ̃_Δ(Ê_α”α')) V}  .

The integration constants have been chosen such that for |E|≫ V and J_V|ln-iE V|≪ 1, 
the result (<ref>) at high energies is reproduced.
From the solution we can see that branch cuts appear starting at the singularities 
z_n=-iΓ^*+μ̂_αα', as stated in (<ref>). Furthermore, we
can see that the expansion is well-defined provided that J_V|lnV |E-z_n||≪ 1, which is the precise
condition that E should not be exponentially close to the branching points. This is the reason why the
scale V∼|z_n| has been chosen as reference scale in the logarithm to integrate the RG equations 
perturbatively for intermediate and small energies. In the solution all logarithmic terms
∼ (J^(0)lnD V)^k have been resummed in J_V, whereas a perturbative treatment has been
used for the logarithmic terms J_V|lnV |E-z_n||≪ 1.

Since Γ̃_Δ(E) is a weakly varying function for |E|≲|z_n|, we can replace 
Γ̃_Δ(Ê_αα')→Γ^* in the above equations and neglect
the term ∼ J_V^2Γ^* ∼ J_V^4 in (<ref>). Inserting 
E=-iΓ^*+δ in (<ref>) (where |δ|≪Γ^* is a small 
scale to exclude an exponentially small region around z_0^p=-iΓ^*), we find 
straightforwardly the result (<ref>) for Γ^*. Inserting the solution for Z'(E)
in (<ref>), we can calculate the pre-exponential function for the contribution from the
branching pole at E=z_0^p=-iΓ^*. For long times t≫1 V we obtain

    F^p_0(t) = 1 - 2 J_V - 4 ∑_α x_α^2 J_V^2 ln(Vt)  ,

whereas, for intermediate times t∼1 V, the contribution in O(J_V^2) is not 
logarithmic and unimportant (the precise coefficient is also influenced by other sub-leading terms).
Several interesting features appear in this result. The first term is the result from a Markov
approximation, where only the pole without the residuum is considered. We note that the
pole position is also influenced by non-Markovian contributions arising when 
(<ref>) is solved self-consistently for Γ^*. Here, this
is a very weak effect occuring in O(J_V^4). For quantum dot models such non-Markovian
contributions have been dicussed perturbatively in Ref. <cit.>.
All other terms of (<ref>) are of pure non-Markovian nature
arising from the term linear in E of the effective Liouvillian (leading to the Z'-factor). 
The second term linear in J_V can not be obtained from perturbation theory since this term
of the Z'-factor involves
the difference J_V-J^(0) (see Eq. (<ref>)), which reduces to J_V
only in the scaling limit. It arises from a resummation of a series of logarithmic terms
∼ (J^(0)lnD V)^k which starts at k=2, i.e. the k=1 term is absent.
The last term ∼ J_V^2 ln(Vt) is logarithmic in time and becomes
of O(1) for exponentially large times t∼1 Ve^1/J_V^2. In this regime the 
solution can no longer be used since it corresponds to the regime of energies E exponentially
close to z_0^p. In this regime, the full RG equation (<ref>) for the vertex
shows that J(E) does not stay small, i.e. a strong coupling problem arises and the truncation
scheme is no longer controlled. As a consequence we see that, concerning the long-time evolution
at exponentially large times, even in the regime V≫ T_K, a strong coupling method is needed
to calculate pre-exponential functions. On the other hand, the exponential decay e^-Γ^*t
leads to a very small contribution for exponentially large times, so that it is of no practical use to
know the pre-exponential function in this regime. However, for other problems with quantum critical
points, like e.g. multi-channel Kondo models or the sub-ohmic spin boson model, it happens that
the pole z_0^p=0 lies at the origin such that no exponential decay appears. For such models, it
is an interesting subject for the future to calculate the precise form of F_0^p(t) for exponentially
large times. E.g., for a multi-channel Kondo model with many channels N≫ 1, which turns out to
be a weak-coupling problem in the whole complex plane, it has been shown
in Ref. <cit.> that F_0^p(t)∼ (1 T_K t)^4/N.

Finally, to calculate the branch cut contribution (<ref>), we need also the
jump δΓ(z_αα'^b-ix) for x∼ 1/t≲ V. This is obtained from
the RG equation (<ref>), which reads with 
J_αα'(E)→ 2√(x_α x_α')J_V

    ∂∂ xδΓ(z_αα'^b-ix) = 
    - 8π i x_α x_α' J_V^2 θ(x)  .

This leads to δΓ(z_αα'^b-ix)=-8π i x_α x_α' J_V^2 xθ(x).
Inserting this result in (<ref>) and using Z'=1+O(J_V^2), we obtain for
long times t≫1|μ_α-μ_α'| (note that αα'
and we assumed that |μ_α-μ_α'|≫Γ^*)

    F_αα'^b(t) = - 4 x_α x_α' J_V^2 (1 (μ_α-μ_α')t)^2  .

Other time regimes t∼1|μ_α-μ_α'| can also be studied leading to
exponential integrals <cit.>. Inserting (<ref>) in 
(<ref>) we get an oscillating term 
∼ J_V^2(1 (μ_α-μ_α')t)^2 e^-Γ^* te^-i(μ_α-μ_α')t
for the time evolution of the local spin. It appears in second order in J_V and is again of
non-Markovian nature. In contrast to the Markov contribution it oscillates with a frequency set
by the differences of chemical potentials and the pre-exponential function decays as a power law ∼ 1/t^2 
for long times. This behaviour is quite generic for models with spin or orbital fluctuations. In 
higher orders the oscillation frequencies are set by the renormalized excitation energies of the system
associated with certain processes. E.g. a process where a particle is transferred from reservoir α
to reservoir α' involves an energy cost μ_α'-μ_α, which gives the oscillation
frequency. In the presence of a local magnetic field h^(0), the same process costs the energy
μ_α'-μ_α± h if the local spin is flipped, where h is the renormalized magnetic
field. As a consequence, these scales define further oscillation frequencies. In addition, each process
has its own decay rate, setting the scale of the exponential decay. These issues have been discussed in
detail in Ref. <cit.> for the case of the anisotropic Kondo model at finite magnetic field.




 §.§ Ohmic spin boson model


Here, we consider the ohmic spin boson model at zero bias ϵ=0 and zero temperature T=0.
We will follow Ref. <cit.> where the model has recently been solved for weak damping 
α≪ 1 by a systematic RG analysis using the E-RTRG method. In contrast
to the Kondo model it turns out that the effective vertex G(E)=G_1(E)_ω̅=0 at zero 
frequency stays small in the whole complex plane allowing for a full solution of the problem on all
time scales. We show here only the solution since the derivation is very similiar to the one for the Kondo
model, except that the algebra is more involved and the solution of the RG equations can also be derived
for E close to the branching points z_n. It turns out that the resolvent R(E) has four poles at 

    z_st = 0  ,   
    z_0 = -iΓ  ,   z_± = ±Δ̃ - iΓ/2   ,

where 

    Δ̃ = Δ (Δ̃ D)^α = Δ (Δ D)^α 1-α  ,  Γ = παΔ̃  .

Δ̃ is called the renormalized tunneling which is kept fixed in the scaling
limit D→∞ and α→ 0. In leading order truncation at O(α) it turns 
out that no branching poles appear, i.e. all poles are isolated. In addition, the eigenvalue
λ_0(E) has two branch cuts starting at z_± and the eigenvalues λ_±(E)
have a branch cut starting at z_0. Therefore, according to the general expression 
(<ref>) we get

    ρ(t) = ρ_st + ∑_k=0,± (F_k^p(t)  +  F_k^b(t)) 
    e^-i z_k t ρ_t=0  ,

i.e. all singularities can either act as a pole or as a branch cut.
Ordering the four possible states in Liouville space by ++,–,+-,-+, where ± are the
two local states, one can show that the stationary density matrix is given by 
ρ_st=1 2(1,1,Δ̃Δ,Δ̃Δ)
and the pre-exponential functions for long times Δ̃ t≫ 1 (note that this includes
the important regime Γ t∼ O(1) where the exponentials are of O(1)) are given by <cit.>

    F_0^p(t)    = Δ̃Δ ([    0    0;   -1 Δ̃/Δ ])
    ⊗τ_+   ,  
    F_±^p(t)  = 1 2 ([        1    ±Δ̃/Δ;    ±Δ̃/Δ (Δ̃/Δ)^2 ])
    ⊗τ_-   ,
    
    F_0^b(t)     =  - 2α 1(Δ̃ t)^2 ([ 1 0; 0 0 ])
    ⊗τ_-   ,  
    F_±^b(t)  =  - α s(t)(Δ t)^2 ([ 0 0; 0 1 ])
    ⊗τ_+   ,

where τ_±=1 2(1±σ_x), and, for two 2× 2-matrices A and B, 
we have defined the 4× 4-matrix 

    A⊗ B ≡ ([ A_11B A_12B; A_21B A_22B ])
      .

Furthermore, the logarithmic function s(t) is defined by 

    s(t) = (1 [1 + αln(Δ̃ t)] 
    [1 - ln(1 + αln(Δ̃ t))])^2
      .

In terms of the expectation values of the Pauli matrices 
⟨σ_i ⟩(t) =Trσ_iρ(t), these equations can also be written as

    ([         1; ⟨σ_x ⟩(t) ])
         = ([    1; Δ̃/Δ ])
     + Δ̃Δ ([    0    0;   -1 Δ̃/Δ ]) ([          1; ⟨σ_x ⟩_t=0 ])
     e^-i z_0 t
       
    - α s(t)(Δ t)^2 ([ 0 0; 0 1 ]) ([          1; ⟨σ_x ⟩_t=0 ])
     ∑_σ=± e^-i z_σ t  ,
    
    ([   ⟨σ_z ⟩(t); -i⟨σ_y ⟩(t) ])
         = 1 2 ∑_σ=± ([        1    σΔ̃/Δ;    σΔ̃/Δ (Δ̃/Δ)^2 ]) ([   ⟨σ_z ⟩_t=0; -i⟨σ_y ⟩_t=0 ])
     e^-i z_σ t  ,
       
    - 2α 1(Δ̃ t)^2 ([ 1 0; 0 0 ]) ([   ⟨σ_z ⟩_t=0; -i⟨σ_y ⟩_t=0 ])
     e^-i z_0 t  .

In this result all logarithmic terms at high energies ∼ (αlnDΔ)^k have
been resummed in the renormalized tunneling Δ̃, and all logarithmic terms at low energies 
(or large times) ∼(αln(Δ̃ t))^k are contained in s(t). For the pre-exponential
function F_0^b(t) it turns out that, in leading order, no logarithmic terms are present at large times.
This has to be contrasted to the solution within the noninteracting blip 
approximation (NIBA) NIBA <cit.>, where, for 
⟨σ_z ⟩_t=0=1 and ⟨σ_y ⟩_t=0=0, one obtains

    ⟨σ_z⟩(t)_NIBA = e^-Γ 2tcos(Δ̃ t) - 
    2α 1(Δ̃ t)^2-2α ,

whereas the correct result from (<ref>) reads

    ⟨σ_z⟩(t) = e^-Γ 2tcos(Δ̃ t) - 
    2α 1(Δ̃ t)^2 e^-Γ t .

Besides the missing exponential part in the second term, which has already been discussed at
the end of Section <ref>, the NIBA predicts a different exponent for
the pre-exponential power law. This shows that power-law exponents of pre-exponential functions
can only be calculated by resumming consistently all logarithmic terms for long times. The E-RTRG
method predicts that no such logarithmic terms are present for ⟨σ_y,z⟩(t) 
but they appear for ⟨σ_x⟩(t) within the logarithmic function s(t).
The leading power-law behaviour ∼ (1 t)^2 of the pre-exponential function is the same 
as for the Kondo model and can also be obtained from perturbative calculations <cit.>.
There are always two terms with different decay rates Γ and Γ/2 for the time evolution. 
If one transforms to the exact eigenbasis 
e_1/2=1√(2)(|+⟩± |-⟩) of the local system, the 
expectation values ⟨γ_i⟩ of the Pauli matrices in the new basis are related 
to the ones of the original basis by ⟨γ_x⟩=⟨σ_z⟩,
⟨γ_y⟩=-⟨σ_y⟩ and ⟨γ_z⟩=⟨σ_x⟩.
Thus the Markovian term ∼ e^-Γ t from the pole contribution describes the decay of
the diagonal matrix elements of the density matrix in the new basis, whereas the one
∼ e^-(Γ/2)te^±Δ̃ t corresponds to the decay of the nondiagonal matrix elements.
Therefore, Γ is called the relaxation rate, whereas Γ/2 is the decoherence rate, in
accordance with the general rule that, in the absence of pure dephasing, the relaxation rate is 
always twice as large as the decoherence rate. 


For large energies |E|≫Δ̃, one needs the function Z'(E) to determine the regime 
of short times t≪1Δ̃ from (<ref>) (the contribution from the 
exponential is a small correction and can be neglected). One obtains the result

    Z'(E) = ∑_σ=± ([      1      0;      0 Z_σ(E) ])⊗τ_σ  ,   Z_±(E) ≈ (-iE D)^2α

This gives rise to the universal short time behavior 

    ρ(t) = ([         1         0;         0 (1 Dt)^2α ])
    ⊗ 1 ρ_t=0  ,

or

    ⟨σ_x,y⟩(t) = (1 Dt)^2α ⟨σ_x,y⟩_t=0  ,  ⟨σ_z⟩(t) = ⟨σ_z⟩_t=0  .

This agrees with previous predictions and can also be obtained from the exact solution (<ref>)
at Δ=T=0 in the universal regime t≫1 D. Again we can see that all logarithmic terms
∼(αln(Dt))^k have been resummed in this result.







 §.§ Interacting resonant level model


Finally we discuss the IRLM for the special case of a single reservoir with chemical potential μ=0 
and zero level position ϵ=0 (i.e. in resonance with the reservoir). As discussed in 
Section <ref>, this model can be mapped to the ohmic spin boson model close to the
exactly solvable point α=1 2. In particular, we want to understand where the crossover
from coherent to incoherent time evolution by changing the sign of U=1-√(2α) comes from.
We follow Ref. <cit.>, where the IRLM has been studied by using
E-RTRG and functional RG. 

We concentrate on the time evolution of the occupation ⟨ n⟩(t) of the local level 
which is related via Eq. (<ref>) to the expectation 
⟨σ_z⟩(t) within the spin boson 
model by 2⟨ n ⟩(t)-1=⟨σ_z⟩(t). For ⟨σ_z⟩(t), one can
show that the result can be written in the form
⟨σ_z ⟩(t)=P(t)⟨σ_z⟩_t=0, with

    P(t) = i 2π ∫_-∞+i0^+^∞+i0^+ dE e^-iEt 1 E + iΓ_1(E)  .

For this special case there is no Z'-factor and Γ_1(E) is a slowly varying logarithmic
function describing the energy dependent charge relaxation rate. It is determined from the RG equations

    ∂∂ EΓ_1(E) = -g R_2(E) Γ_1(E)
      ,  ∂∂ EΓ_2(E) = -g R_1(E) Γ_1(E)
      ,

where g=2U-U^2=1-2α and the resolvents R_1/2(E) are defined by

    R_1(E)  = 1 E + iΓ_1(E)  ,  
    R_2(E)  = 1 E + iΓ_2(E)/2  .

The initial conditions are given by Γ_1/2(E=iD)=Γ^(0). 
Γ_2(E)/2 is also a slowly varying logarithmic function and describes the energy dependent
broadening of the local level corresponding to the decoherence mode for nondiagonal matrix elements
of the local density matrix w.r.t. the charge states (note, however, that such elements can not be prepared).
As we will see below the subtle coupling of the two RG equations for Γ_1/2(E) leads to the 
interesting effect that, for g>0, the resolvent R_1(E) can have poles with a finite real part 
although the local system has no finite excitation energy. 

We start by solving the RG equations at high energies E≫Γ_1/2(E). Neglecting Γ_1/2(E)
on the r.h.s. of the RG equations, we find the solution

    Γ_1/2(E) = Γ^(0) (D -iE)^g = Δ̃ (Δ̃ -iE)^g ,

where 

    Δ̃ = Γ^(0) (DΔ̃)^g = Γ^(0) (DΓ^(0))^g/(1+g)

is the renormalized tunneling which is kept fixed in the scaling limit D→∞
and α,Γ^(0)→ 0. Using the relation g=1-2α and 
Γ^(0)=Δ^2 D to the spin boson model, one can see that it is identical 
to the definition (<ref>) of the renormalized tunneling for the spin boson model.
As discussed in detail in Refs. <cit.>, the solution at
high energies contains all leading logarithmic terms ∼ (UlnD -iE)^k and 
all subleading ones ∼ U(UlnD -iE)^k. From the solution (<ref>) 
at high energies we can calculate with (<ref>) the time evolution 
for short times t≪ 1/Δ̃ as

    P(t) ≈ e^-Γ_1(1/t)t ≈ e^-(Δ̃ t)^1+g,

i.e. the relaxation rate in the exponent is cut off at the energy scale 1/t. In contrast
to the spin boson model at small α and the Kondo model, there is no Z'-factor and
therefore the exponential provides the leading order. Expanding the exponential we find
P(t)=1-(Δ̃ t)^1+g in agreement with previous results <cit.>.
Since (Δ̃ t)^1+g=(Dt)^g Γ^(0)t, we see again that all logarithmic
terms ∼ (gln(Dt))^k have been resummed for small times.


Next we study the analytic structure of the resolvent R_1(E) to find the time evolution for
intermediate and long times. As we will show below, for positive g>0, R_1(E) has two poles 
at z_± (followed by a branch cut with jump of O(g^2) which can be neglected) and one
branch cut starting at z_0 (with jump of O(g)), where the singularities z_n, n=0,±, are given by

    z_0 = - i Δ̃ 2  ,  
    z_± = ± Ω - i Δ̃  ,  Ω = π g Δ̃  .

For g<0, there is only a branch cut starting at z_0. Thereby, z_0 is the position
of the pole of the resolvent R_2(E), i.e. z_0 and z_± can be determined from the
equations

    z_± + iΓ_1(z_±) = 0   ,   z_0 + iΓ_2(z_0)/2 = 0   .

Note that, in contrast to the singularities (<ref>) for the spin boson model at
small α, for the IRLM (or the spin boson model at α=1 2) the renormalized 
tunneling determines the rate and not the oscillation frequency. Furthermore, we note that
the pole of R_1(E) describes the charge relaxation mode, whereas for the spin boson model
at small α it corresponds to the decoherence mode w.r.t. the exact eigenstates of the
local system. Therefore, z_0 corresponds to the decoherence mode for the IRLM and its 
imaginary part is half of the one of the relaxation poles z_±. To derive the result for the 
positions of the singularities we solve the RG equations for intermediate and small energies
|E|≲Δ̃ but glnΔ̃ |E-z_n|≪ 1, i.e. E should not
be exponentially close to the singularities. Expanding in the small parameter 
glnΔ̃ |E-z_n|≪ 1 and fixing the integration constants by comparing
with the solution (<ref>) at high energies in the usual way, we find

    Γ_1(E)/Δ̃ ≈ 1 - g ln-iE + Γ_2(E)/2Δ̃  ,  Γ_2(E)/Δ̃ ≈ 1 - g ln-iE + Γ_1(E)Δ̃ .

In contrast to the corresponding equation (<ref>) for the Kondo model,
there is a subtle coupling of the singularites of Γ_1(E) and Γ_2(E), which leads
to the new feature that z_± obtains a finite real part for g>0. We note that although
the equations can not be used for E exponentially close to the singularities, they can be
used for |E-z_n|∼ g^2 since gln(g)≪ 1 for g≪ 1. Therefore, the equations can be 
used to determine the positions of the branching points of Γ_1/2(E) up to O(g).
From the equations we can see that Γ_1(E) (Γ_2(E)) have a branch cut with jump 
of O(g) starting at the branching point of the logarithmic function where (<ref>)
is fulfilled, i.e. at z_0 (z_±). Thereby, the branch cut of Γ_2(E) starting at
z_± leads also to a branch cut for Γ_1(E) at the same position but this branch cut
has a jump of O(g^2) and can be neglected. Inserting the leading order results
Γ_1/2(E)≈Δ̃, z_0≈-iΔ̃/2 and z_±≈-iΔ̃
on the r.h.s. of (<ref>), we find for the position of the singularities
the result (<ref>)

    2i z_0/ Δ̃   = Γ_2(z_0)/Δ̃≈  1-gln(-iz_0/Δ̃+1) 
    ≈  1-gln( -1 2+1 )  ≈ 1   ,
    
    i z_±/Δ̃   = Γ_1(z_±)/Δ̃≈  1-gln( -iz_±/Δ̃+1 2) 
    ≈  1-gln( -1∓ iΩ/Δ̃+1 2) 
       ≈ 1± iπ g   .
 


Due to the analytic structure of the resolvent R_1(E) the time evolution can be written as

    P(t) = θ(g) ∑_σ=± e^-iz_σ t + F_0^b(t) e^-iz_0t = θ(g) 2 cos(Ω t) e^-Δ̃t + F_0^b(t) e^-(Δ̃/2)t ,

where the first term involves the contribution from the isolated poles (we have neglected
corrections of O(g) to the residuum) and the second term involves the analog of the branch cut integral
(<ref>), which can be written as

    F_0^b(t) = 1π Im ∫_0^∞ dx e^-xt 1 -i(z_0 + iΓ_1(z_0-i/t+0^+)) - x  .

For intermediate and long times t≳1Δ̃ but 
gln(Δ̃t)≪ 1, F_0^b(t) can be evaluated by using the result 
(<ref>), where we obtain 
-iΓ_1(z_0-i/t+ 0^+)=z_+(1 + O(gln(Δ̃t)).
In particular one has to consider the fact that x∼1 t
can not be neglected compared to the difference |z_0-z_±|∼Δ̃ for 
intermediate times t∼1Δ̃. This time regime is of particular interest here 
since the exponentials of the time evolution (<ref>)
decay on the time scale 1Δ̃. Therefore, the
integral (<ref>) has to be calculated more carefully in terms of
the exponential integral E_1(z)

    F_0^b(t) = - 1π Im {e^-i(z_0-z_+)t E_1(-i(z_0-z_+)t)}  .

This result has been used in Refs. <cit.> to discuss the
competition between the oscillating (i.e. coherent) and the purely decaying (i.e. incoherent) term of
the time evolution in Eq. (<ref>). Since the incoherent term decays on 
a longer time scale it turns out that it wins very rapidly such that the coherent term leads only
to a few number of oscillations, in contrast to the physics of a classical damped harmonic oscillator.
For long times t≫1Δ̃ but still gln(Δ̃t)≪ 1, the incoherent 
term dominates and, using the asymptotic expansion E_1(z)=e^z/z of the exponential integral, one obtains

    F_0^b(t) ≈ -4g 1Δ̃t  ,

i.e. a power law ∼ 1/t typical for models with charge fluctuations.

Finally, for exponentially large times gln(Δ̃t)∼ O(1), we need the solution for Γ_1(E)
for energies E exponentially close to the branching point z_0. In this regime, we can replace
Γ_2(E)/2→ iz_0 on the r.h.s. of the RG equation (<ref>) for Γ_1(E),
which gives the solution

    Γ_1(E) = Δ̃ (Δ̃ -i(E-z_0))^g  ,

where the integration constant has been fixed by comparison with the solution (<ref>)
at intermediate and small energies. Using this solution for the evaluation of the branch cut
integral (<ref>) for exponentially large times, we can neglect x  
in the denominator and find with 
Γ_1(z_0-i/t+0^+)=Δ̃(Δ̃t)^g+2π i gΔ̃(Δ̃t)^g 
the result (neglecting terms of O(g^2) in the denominator)

    F_0^b(t) ≈ -g 1 (1/2 - (Δ̃t)^g)^2 1 (Δ̃t)^1-g  .

This result holds for all times t≫1Δ̃. For long times with gln(Δ̃t)≪ 1 it 
reduces to the result (<ref>). However, for exponentially large times where
(Δ̃t)^g is some number of O(1), the result changes. In the extreme regime
(Δ̃t)^|g|≫ 1, it reduces to

    F_0^b(t) ≈ -g [1 + 3θ(-g)] 1 (Δ̃t)^1+|g|  .

This result agrees with the prediction of the NIBA <cit.> and its improved version 
<cit.> (where the exponential term e^-(Δ̃/2)t has also been obtained, 
see Eq. (<ref>)). However, as we have seen, it holds only for extremely 
long times and, for g>0, the prefactor is different from the
result (<ref>) for more realistically long times Δ̃t≫ 1 with 
gln(Δ̃t)≪ 1. Therefore, we see that the regime of long times is very subtle and
the result can change significantly by entering the regime of exponentially large times.

Finally, as already mentioned in Section <ref>, it has not yet been studied to a full extent
how the RG equation (<ref>) looks like in higher orders in the tunneling. There is some
evidence that all higher order terms in Γ are of the form

    U^n (Γ_iΔ̃)^k (lnE-z_nΔ̃)^l   for   n=1,2
      ,  
    U^n Γ_i E-z_n (Γ_iΔ̃)^k
    for   n>2  ,

i.e., after integration, either vanish in the limit E→ z_n or contribute to higher 
orders in U, but this is still under investigation. Furthermore, the results have been compared to 
functional RG in Refs. <cit.>, where all orders in the tunneling have been 
resummed keeping only the lowest order term in the Coulomb interaction. The numerical results of functional
RG agree quite nicely with the analytical result (<ref>) 
for intermediate and long times and, 
in particular for extremely long times, the result (<ref>) has been confirmed 
analytically by functional RG. Therefore, there is good evidence that also within E-RTRG 
higher orders in the tunneling will not change the results at least in leading order in U.

Acknowledgments. I am particularly thankful to M. Wegewijs for a thorough reading of the manuscript.



















99

projection_operator
H. Mori, Progress of Theoretical Physics 33, 423 (1965);
R. Zwanzig and R. D. Mountain, Journal of Chemical Physics 43, 4464 (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.

functional_integral
A. J. Leggett et al., Rev. Mod. Phys. 59, 1 (1987);
H. Grabert, P. Schramm, and G. L. Ingold, Phys. Rep. 168, 115 (1988);
U. Weiss, Quantum Dissipative Systems 
(World Scientific Publishing Company, Singapore, 2012).

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

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

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

kondo_E-RTRG
M. Pletyukhov and H. Schoeller, Phys. Rev. Lett. 108, 260601 (2012);
F. Reininghaus, M. Pletyukhov, and H. Schoeller, in preparation.

RTRG_irlm
C. Karrasch et al., Europhys. Lett. 90, 30003 (2010);
S. Andergassen et al., Phys. Rev. B 83, 205103 (2011).

kennes_PRL13
D. M. Kennes et al., Phys. Rev. Lett. 110, 100405 (2013).

kashuba_kennes_PRB13
O. Kashuba et al., Phys. Rev. B 88, 165133 (2013).

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

flow_equation
S. Kehrein, Phys. Rev. Lett. 95, 056602 (2005).
D. Lobaskin and S. Kehrein, Phys. Rev. B 71, 193303 (2005).

functional_RG
R. Gezzi, T. Pruschke, and V. Meden, Phys. Rev. B 75, 045324 (2007);
S. G. Jakobs, V. Meden, and H. Schoeller, Phys. Rev. Lett. 99, 150603 (2007);
W. Metzner et al., Rev. Mod. Phys. 84, 299 (2012);
D. M. Kennes et al., Phys. Rev. B 85, 085113 (2012).

TD-NRG
F. B. Anders and A. Schiller, Phys. Rev. Lett. 95,  196801  (2005);
F. B. Anders, R. Bulla, and M. Vojta, Phys. Rev. Lett. 98, 210402 (2007);
A. Hackl et al., Phys. Rev. Lett. 102, 219902 (2009).

TD-DMRG
A. Daley et al., J. Stat. Mech.: Theor. Exp. P04005 (2004);
S. R. White and A. Feiguin, Phys. Rev. Lett. 93, 076401 (2004);
P. Schmitteckert, Phys. Rev. B 70, 121302 (2004);
F. Heidrich-Meisner, A. E. Feiguin, and E. Dagotto, Phys. Rev. B 79, 235336 (2009).

ISPI
S. Weiss et al., Phys. Rev. B 77, 195316 (2008).

QMC
T. L. Schmidt et al., Phys. Rev. B 78, 235110 (2008).

exact
F. Lesage and H. Saleur, Phys. Rev. Lett. 80, 4370 (1998);
A. Schiller and S. Hershfield, Phys. Rev. B 62, R16271 (2000);
A. Komnik, Phys. Rev. B 79, 245102 (2009).

RG_noneq_D_cutoff
A. Rosch, J. Kroha and P. Wölfle, Phys. Rev. Lett. 87, 156802 (2001);
A. Rosch et al., Phys. Rev. Lett. 90, 076804 (2003);
L. I. Glazman and M. Pustilnik, in Nanophysics: Coherence and Transport
(H. Bouchiat et al., Elsevier, 2005), p. 427;
T.Korb et al., Phys. Rev. B 76, 165316 (2007);

saptsov_PRB12
R. B. Saptsov and M. R. Wegewijs, Phys. Rev. B 86, 235432 (2012). 

saptsov_arXiv
R. B. Saptsov and M. R. Wegewijs, arXiv:1311.1368.

correlation_function
D. Schuricht and H. Schoeller, Phys. Rev. B 80, 075120 (2009);
S. Y. Mueller et al., arXiv:1211.7072.

kashuba_EPL12
O. Kashuba, H. Schoeller, and J. Splettstoesser, EPL 98 (2012) 57003.

slutskin_EPL11
A. A. Slutskin et al., Europhys. Lett. 96, 40003 (2011).

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

pms
P. W. Anderson, J. Phys. C 3, 2436 (1970);
A. C. Hewson, The Kondo Problem to Heavy Fermions (Cambridge University Press, 1997).

goettel_preprint
S. Goettel, F. Reininghaus, and H. Schoeller, in preparation.

Hackl_09
A. Hackl et al., Phys. Rev. Lett. 102, 219902 (2009);
A. Hackl, M. Vojta, and S. Kehrein, Phys. Rev. B 80, 195117 (2009).

contreras_PRB12
L. D. Contreras-Pulido et al., Phys. Rev. B 85, 075301 (2012).

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


