Dynamics of Open Quantum Systems I, Oscillation and Decay

We develop a framework to analyze the dynamics of a finite-dimensional quantum system $\rm S$ in contact with a reservoir $\rm R$. The full, interacting $\rm SR$ dynamics is unitary. The reservoir has a stationary state but otherwise dissipative dynamics. We identify a main part of the full dynamics, which approximates it for small values of the $\rm SR$ coupling constant, uniformly for all times $t\ge 0$. The main part consists of explicit oscillating and decaying parts. We show that the reduced system evolution is Markovian for all times. The technical novelty is a detailed analysis of the link between the dynamics and the spectral properties of the generator of the $\rm SR$ dynamics, based on Mourre theory. We allow for $\rm SR$ interactions with little regularity, meaning that the decay of the reservoir correlation function only needs to be polynomial in time, improving on the previously required exponential decay. In this work we distill the structural and technical ingredients causing the characteristic features of oscillation and decay of the $\rm SR$ dynamics. In the companion paper [27] we apply the formalism to the concrete case of an $N$-level system linearly coupled to a spatially infinitely extended thermal bath of non-interacting Bosons.


Introduction
The fundamental evolution equation of quantum theory is the Schrödinger equation, which governs the dynamics of closed quantum systems, isolated from their surroundings. In contrast, when a system is interacting with its surroundings, then the evolution equation for the degrees of freedom of the system alone have to be deduced from other fundamental principles. One approach is to view the system plus its surroundings together as a closed system, meaning that they evolve according to the Schrödinger equation involving both, the degrees of freedom of the system and those of the surroundings. This amounts to a huge number of variables, even if the system itself is small. The evolution equation for the system is obtained by "tracing out" the degrees of freedom of the environment, and its solution is called the reduced system dynamics. Generically, this situation is described as follows. The Hilbert space of states is thermal reservoir if one has in mind a big system in thermal equilibrium). The total system (S + R) being closed, the dynamics of an initial (time t = 0) SR density matrix ρ(0) is given by in accordance with the Schrödinger equation -also called the von-Neumann equation in the case of density matrices. Here, H is the total Hamiltonian, which is the sum of individual system and reservoir Hamiltonians plus an interaction term λV , where λ is a coupling constant. The operators H S and H R act non-trivially only on the system and reservoir factors of H, while V mixes up those factors. Typically V = A ⊗ B (or a sum of such terms), where A and B are system and reservoir operators, respectively. The reduced system density matrix at time t is given by

3)
where tr R denotes the partial trace over the reservoir part of H. For λ = 0, the dynamics is uncoupled, e −itH = e −itH S ⊗ e −itH R and each factor (subsystem) evolves individually and independently. The dynamics is usually known explicitly in the uncoupled case, and the challenge is to find out what happens when λ = 0. The idea then is, usually, to carry out perturbation theory for small λ. If the initial state is not correlated, ρ(0) = ρ S (0) ⊗ ρ R (0), then the map is well defined. It is the flow, mapping an arbitrary system initial state to its state at time t ≥ 0. For λ = 0, we have V t = e tL S , where L S ρ = −i[H S , ρ]. However, for λ = 0, generically, S and R become entangled during the evolution, and this prevents V t from having the above simple structure. Namely, V t cannot be written as e tL for any L; this means that V t is not a semigroup in t any longer (we call e tL a semigroup instead of a group, if we restrict the values of t to t ≥ 0). Apart from encoding the Markovian characteristics, any dynamics of the form e tL has the practical advantage that the eigenvalues and eigenvectors of the operator L encode the dynamical behaviour. Spectral theory thus becoms a main tool in the analysis of such a Markovian dynamics.
Markovian approximation. Even though the flow V t is not of exponential (semigroup) form, one expects that under suitable conditions, V t can be approximated by a semigroup, with L(0) = L S . This is called the Markovian approximation [5,36]. The approximate evolution equation, which reads in differential forṁ ρ M (t) = L(λ)ρ M (t), (1.4) is the ubiquitous Markovian master equation. On physical grounds, its validity is plausible under two main assumptions: Firstly, the reservoir loses its memory quickly and secondly, the reservoir is 'large' and the interaction is weak, so as to guarantee that the reservoir dynamics is not affected much by the inteaction with the system. These assumptions, also called the Markov-and Born approximations, respectively, are quantified by the decay of the reservoir correlation function (quick memory loss) and a smallness condition on λ (weak coupling). The main two issues are (a) How to construct L(λ) and (b) How to estimate the error made by the Markovian approximation.
The generator L(λ) has to be constructed starting from the full SR dynamics (interacting Hamiltonian H as above) and reducing the evolution to the system component alone by tracing out the reservoir degrees of freedom. The literature on this topic is enourmous, and there have been many proposals for L(λ). Eventually, the so-called Davies generator [9,10] L = L S + λ 2 K (1.5) emerged as the 'right' generator. Here, K is an operator describing the effect of R on S, which is explicitly determined by the interaction operator V and the reservoir correlation function. The approximate Markovian dynamics e tL is a so-called CPTP semigroup [1,6] (completely positive, trace preserving), which guarantees, in particular, that the approximate solution ρ M (t) is a (positive definite, normalized) density matrix for all t. For other generators proposed in the literature, this is not the case [13]. Of course, the Davies generator is not the remedy for all ailments. There are plenty of interesting situations where the usual Davies generator (1.5) does not lead to a good approximation. An important class of examples are systems with small energy level spacings, for which a different generator has to be used [32,22,40,33,35,29,30].

Goals
We have two goals: G1 Identify a general framework for the perturbation theory of the dynamics of a unitarily evolving total system-reservoir complex. Isolate the dominant term of the total dynamics for small coupling. Derive from it the open system dynamics. Do this with weak regularity requirements on the interaction between the parts.
G2 Show how a widely used model of open quantum systems, an N -level system coupled linearly to a spatially infinitely extended system of non-interacting Bosons in thermal equilibrium, fits the general framework. Show the validity of the Markovian master equation for this model, under weak regularity conditions.
We achieve the first goal in the current manuscript, while the second one is addressed in the companion paper [27]. The setup for G1 is more general than the usual setup found in the open quantum systems literature. This is motivated by us trying to distill the essential ingredients that imply characteristic oscillation and decay features of the total dynamics, as presented in our main result, Theorem 3.1. The companion paper [27] can be read independently of the current work, and the concrete model considered there will be more familiar to many researchers in the field. We hope that a reader finding the current manuscript somewhat abstract may still consult [27] for the application part, which presents new results on this classical family of models.
We consider a bipartite Hilbert space H, of the form (1.1), carrying a unitary dynamics e itL λ , with a self-adjoint generator (compare to (1.2)) were L S and L R are called the system and reservoir Liouville operators, acting non-trivially only on the respective factors of H, λ is a coupling constant and I is an interaction operator. Our assuptions are expressed mathematically precisely in (A1)-(A5), Section 2. The main two assumptions are described as follows.
• We assume that the reservoir dynamics has a unique stationary state. In other words, the reservoir Liouvillean L R has a unique eigenvalue. This eigenvalue is zero and nondegenerate. The unique stationary state (eigenvector) is denoted by Ω R .
• We assume that away from the reservoir stationary state, the total dynamics is dissipative. By this we mean that for vectors φ, ψ ∈ H satisfying P R φ = P R ψ = 0, where P R = 1l H S ⊗ |Ω R Ω R | is the projection onto the reservoir stationary state, we have | φ, e itL λ ψ | ∼ t −k for large t, for some k ≥ 2. (1.6) Here,L λ denotes the operator L λ restricted to the range of P ⊥ R = 1l − P R (the space of vectors ψ that have a vanishing overlap with the reservoir stationary state Ω R ).
How do we phrase (1.6) precisely? This is done by representing e itL λ by the Fourier-Laplace transform of the resolvent operator (L λ − z) −1 , e itL λ = − 1 2πi R e itx (L λ − x − i0 + ) −1 dx. Using e itx = 1 it d dx e itx and integrating by parts with respect to x gives Stronger decay in t involves the Fourier-Laplace transform of higher powers (L λ − x − i0 + ) −1 (repeat the integration by parts). To estimate such integrals, we require some powers of the resolvent (L λ − z) −1 to be bounded in z ∈ C − = {z : Imz < 0} (weakly on a suitable set of vectors). We give a rigorous discussion of this point below, after (2.26). Our method requires decay ∼ t −2 , which results in the assumption that ∂ j z (L λ − z) −1 is bounded in z ∈ C − , for j = 0, 1, 2, see (2.6).
• Further conditions: We assume some non-degeneracy of the eigenvalues of L λ and some suitable relative boundedness conditions on the interaction operator I (see (A1)-(A5)).
These further conditions, given in detail in Section 2, are less fundamental (and can be weakened, or entirely removed).
We now explain what we mean by weak regularity requirement on the interaction, stated in the Goals G1, G2 above. A result similar to Theorem 3.1 (which implies the validity of the Markovian approximation) was proven in [21,26] under an analyticity condition on the reduced resolvent. Instead of requiring the first three derivatives of (L λ − z) −1 to be (weakly) bounded in z ∈ C − , it was assumed in those works that (L λ −z) −1 has an analytic extension as z moves from the lower complex plane across the real line into the upper complex plane. The difference between the two cases is not merely a technical issue, as the regularity determines physical features. For instance, the stronger analyticity condition implies that overlaps (1.6), as well as the reservoir correlation function, decay exponentially quickly in t. Improving the 'analytic deformation theory' used in [20,21,26] to the weaker regularity regime considered in the current paper and in [27], is not totally easy. It is done using 'Mourre theory', which is more delicate and technically trickier than the analytic theory. This improvement is worth pursuing as it sheds light on the question of how quickly the reservoir has to lose its memory in order for the Markovian approximation still to hold. The current work and [27] show that polynomial decay of the memory (reservoir correlation function) suffices. Incidentally, this polynomial memory loss still drives an exponentially quick approach of the system to its final state (c.f. Theorem 3.1). It is well known that if the memory loss is not quick enough, then the Markovian approximation is not valid [35,37,23]. We plan on extending our methods to this case in the future.
Remarks in view of applications.
(1) In applications [27,28], we need, in general, to describe mixed states (such as equilibrium states). The Hilbert space H considered in the current paper does not, in general, describe pure states. Rather, H here is the space of purification (or, vectorization) of mixed states in question. Similarly, L λ is the Liouville operator, describing the dynamics expressed in the purification Hilbert space. This setup allows us to treat mixed and pure states on the same footing. Example. The system equilibrium state of a qubit with Hamiltonian H S = ω 0 2 σ z = ω 0 2 | ↑ − ω 0 2 | ↓ is the mixed state ρ S,β ∝ e −βH S . Its purification is given by Ψ β,S ∝ e −βω 0 /4 | ↑ ⊗ | ↑ + e βω 0 /4 | ↓ ⊗ | ↓ , a normalized vector in the purification space The key link between ρ S,β and its purification Ψ S,β is that for all observables A and times t, the expectations can be expressed as, Given ρ S,β and H S , neither Ψ S,β nor L S are unique -but there are standard choices for them [4,26,3,18,12]. The reservoir we consider in applications [27,28] consists of non-interacting Bose particles in infinite position space R 3 . They are in a state of thermal equilibrium at temperature 1/β > 0, determined by Planck's black body momentum density distribution: there are (e −β|k| − 1) −1 d 3 k particles having momentum k ∈ R 3 in a volume d 3 k centered around k, for each unit volume in position space. Usual Bosonic Fock space F(L 2 (R 3 , d 3 k)) does not accommodate this state, as any state ψ ∈ F(L 2 (R 3 , d 3 k)) (or denisty matrix) describes finitely many particles, which results in a vanishing particle density at infinite volume! However, one can find a Hilbert space H R (turning out to be the tensor product of two copies of Fock space) and a vector Ω R (turning out to be the tensor product of two vacua) which do represent the equilibrium state at positive density. The details of this construction are not the topic here, they can be found in textbook style in [31] (see also [26], or the original paper [2]). In the mathematical literature, representing a mixed state as a normalized vector in a Hilbert space is known as the Gelfand-Naimark-Segal construction (GNS); in the physics and chemistry literature it is also sometimes called the thermofield method.
(2) We assume that dim H S < ∞. However, the decay (1.6) necessitates thatL λ does not have any eigenvalues, for otherwise clearly the overlap in (1.6) would be independent of time on eigenvectors. This means that the spectrum ofL λ must be continuous. In particular, the Hilbert space H R has to be infinite dimensional.
Example. In the case of a reservoir of infinitely extended Bose particles, continuous energy spectrum arises due to the infinite volume limit. We prove in Theorem 3.1 of [27] that (1.6) holds for this reservoir. The proof is based on the fact that the estimate holds forL λ replaced byL 0 , together with a suitable perturbation theory in λ.

Explanation of the main result
Our main result is Theorem 3.1. It states that for coupling constants λ small enough, we have the expansion valid for all t ≥ 0. Here, P R = 1l H S ⊗ |Ω R Ω R | ≡ |Ω R Ω R | is the projection onto the reservoir stationary state Ω R and P ⊥ R = 1l−P R . The operator on the left side, e itL λ , is the propagator of the full, unitary SR evolution. The first term on the right side describes a non-trivial evolution of S, generated by an operator M (λ), while R is projected onto the stationary state Ω R . This part of the dynamics is Markovian. The second term on the right side describes the dynamics of states (vectors) in the range of P ⊥ R , so states with vanishing overlaps with the reservoir stationary state Ω R . Finally, R(λ, t) is a remainder term. The equality (1.7) is understood in the weak sense, that is, when φ, · ψ is applied to both sides, for suitable vectors ψ, φ ∈ H belonging to a dense set. This set of vectors includes all uncorrelated SR states in which the system is in any state and the reservoir is in the stationary state Ω R . However, the set also includes a wide class of correlated initial SR states; this aspect is exploited in [28] to show the validity of the Markovian approximation even for correlated initial states. In Theorem 3.1, • We find the detailed structure of e itM (λ) . It consists of explicit oscillating terms and explicit decaying terms with decay rates ∝ 1/λ 2 .
We apply our main result (1.7) to the concrete class of open quantum systems explained in goal G2 above in [27,28].

New results in the theory of open quantum systems
We verify in [27] that the assumptions (A1)-(A5) of Section 2 below are satisfied, and hence that Theorem 3.1 holds, for a standard class of open quantum systems described in goal G2 above. For this model, the generator M (λ) is just the Davies generator L, (1.5), as shown in [26], and Theorem 3.1 yields the following new results.
1. Proof of the validity of the Markov approximation for all times [27].
Davies [9,10] showed the validity of the approximation in the ultra-weak coupling limit. 1 Namely, for any a > 0, lim where · is the norm of super operators and L is the Davies generator (1.5). This means that e tL is a good approximation of the true system dynamics V t for small values of λ, but only for times up to t ∼ λ −2 . In particular, Davies' result does not guarantee that the asymptotic system dynamics is approximated by the Markovian master equation.
Our results overcomes this defect. Namely, we show in [27] that there is a λ 0 > 0 such that if |λ| < λ 0 , then sup for a constant C independent of λ, t.
We can rephrase the result (1.8). Let ρ S (t) be the system density matrix (1.3) and let ρ M (t) be the solution of the Markovian master equation (1.4), with L(λ) = L the Davies generator, and with equal initial conditions, ρ S (0) = ρ M (0). Then (1.8) asserts that where σ 1 = tr S (|σ|) is the trace norm and C is a constant independent of the initial condition.
2. Proof of the validity of the Born approximation [28].
Our main result (1.7) gives an expansion of the full, interacting SR dynamics. When reduced to the system dynamics alone, it results in (1.8). However, when analyzed in full, (1.7) is shown in [28] to yield the following results.
-For uncorrelated initial states ρ(0) = ρ S (0) ⊗ |Ω R Ω R |, where S is in any state and R is in equilibrium, the reservoir stays in equilibrium during the coupled evolution up to an error of O(|λ| 1/4 ), for all times t ≥ 0. This is a proof of the validity of the Born approximation, for all times.
-For a large class of correlated initial SR states, the correlations decay polynomially in time. After this decay has happened, the reservoir is in thermal equilibrium and the system evolves according to the Markovian dynamics generated by the Davies generator, up to errors O(|λ| 1/4 ), uniformly in times t ≥ 0.

Relation to earlier work
It is not our aim to present a detailed discussion of the huge literature on the dynamics of open quantum systems, as the goal of the current manuscript is the construction of a general mathematical framework. Let us rather discuss some mathematically rigorous works related to ours. The spectral methods developed here have a certain similarity with the theory of metastable states in many-body quantum theory. There, the Hilbert space does not have the structure (2.1), but rather, L 0 is the kinetic energy operator of N particles and I is an (interaction) potential. In this Schrödinger operator setup, initial states close to an eigenstate of L 0 stay bound (spatially localized) for a long time under the evolution e −itL λ , but eventually decay for large times [38,16,17]. It is intuitively plausible, quite generally, that eigenvectors of L 0 associated to an unstable eigenvalue e, describing a bound states of a quantum system for λ = 0, turn into 'almost-bound' states for small λ = 0. This phenomenon, and some related mathematical tools (complex spectral deformation and Mourre theory, resolvent representation of the propagator), are common to the many-body Schrödinger and the open system setups.
The dynamics of an N -level system coupled linearly to a spatially infinitely extended reservoir of non-interacting Bose particles (as we will consider as an application of the current results in [27,28]) has been investigated in detail before. It is based on the representation of the infinite volume equilibrium state as a vector in a purification Hilbert space, which was first constructed in the early sixties in [2]. However, it was not exploited to analyze the dynamics of open systems until the pioneering papers [18,3] developed the spectral approach. In these works, the phenomenon of Return to Equilibrium (RtE) was proven: States 'close' to the (coupled) SR equilibrium state are driven to equilibrium in the long time limit. In contrast with usual open system dynamics results, RtE is a result about the full system-reservoir dynamics, not only about the reduced system dynamics. Averages of reservoir observables also converge to the equilibrium values. While [18,3] based their analysis on analytic deformation methods, a Mourre theory version for RtE was developed in [15,25]. A bit later, in [24], the formalism used to show RtE was refined and a detailed description of the dynamics of each system density matrix element was obtained. Then in [20,21,26] is was shown that the main term of the reduced system dynamics is a completely positive trace preserving semigroup. Under general assumptions, the generator of the semigroup is the Davies generator, as is shown in [26].
In terms of technique, the paper [20], where Mourre theory is used to analyze the dynamics, is closest to the present work. However, as explained in discussion point (v) after Theorem 3.1, the approach of [20] is not suited to show the Markovian approximation. In the present work, we make substantial changes to the method of [20], changes which we explain in Section 4.
Further rigorous work. Establishing the validity of Markovian master equations is an active research field. In [32,22], the authors develop a so-called coarse-grained Master equation (CGME) and they compare its validity to that of the Davies and the Redfield master equations. They show that the CGME combines the advantages of the other two, but without incorporating their disadvantages. Namely, the CGME is a good approximation of the true system dynamics regardless of the system level spacing (Davies cannot do this) and the equation is CPTP (Redfield is not). However, the error bounds for neither of the three Master equations are uniform in time; the error of the CGME is ∝ O(|λ|e cλ 2 t ) for some constant c (linked explicitly to the fastest system decoherence time scale). In a similar vein, [40,33] develop a Master equation, again a CPTP equation, which describes the approximation of the true system dynamics for arbitrarily small system level spacings, but again with remainders guaranteed to be small only for finite times. In [33] the authors also provide a bound on the error generated by the Born approximation, for finite times. We have applied our resonance theory to systems with almost degenerate system levels in [29,30] and shown that almost degeneracy leads to a separation of time scales in the system dynamics. The papers [29,30] assume the stronger regularity (analyticity), but we do not see any obstacle to adapting the general theory for weak regularity developed here to the case of degenerate or almost degenerate system levels.

Approach and assumptions
The purpose of this section is to define the setup and state precisely the assumptions we make. The core of our overall strategy is to link the dynamics to the spectral data of the generator L λ = L 0 + λI. A basic mechanism producing irreversible dynamics, or decay, is the disappearence (instability) of eigenvalues of L 0 under the perturbation λI. The instability occurs because the eigenvalues of L 0 are embedded in continuous spectrum. Their fate after perturbation is analyzed using 'singular perturbation theory'. We explain the main ideas of it in Section 2.2. In Section 2.3 we explain how stable and unstable eigenvalues of L λ affect the propagator e itL λ .
Instead of making a list of assumptions, we are trying to proceed in an inductive manner: We explain the main strategy of our approach, and while doing so, the assumptions should emerge naturally.

Basic properties of the Liouville operator L λ
We consider a bipartite Hilbert space where dim H S < ∞, and a family of self-adjoint operators where both L 0 and I are self-adjoint, λ ∈ R and L 0 is of the form This is the setup describing the composition of a system (S) plus reservoir (R) arrangement, in which L S and L R generate the free (non interacting) dynamics of the individual components, I is an interaction operator and λ is a coupling constant. We suppose that L 0 has finitely many eigenvalues e, of finite multiplicity m e , possibly embedded in continuous spectrum, which can cover parts or all of R. The set of eigenvalues e of L 0 is denoted by E 0 and the associated orthogonal eigenprojection is P e . We are going to impose a regularity condition (c.f. (A1) below) which implies the following picture for small λ: All eigenvalues of L λ lie inside an O(λ) neighbourhood of E 0 . Moreover, within such a neighbourhood around any given e ∈ E 0 , either L λ does not have any eigenvalues for λ = 0 (we say e is unstable), or L λ does have some eigenvalues, with summed multiplicity m e not exceeding m e (we say e is stable if m e = m e , and partially stable if 0 < m e < m e ). In the analytic perturbation theory of isolated eigenvalues, the summed multiplicity m e would always equal m e , but for embedded eigenvalues, it is generically strictly less than m e . Fig. 1 gives a graphical depiction of the situation.
Let us now introduce the assumptions and discuss their meaning. We start with an assumption to simplify the bookkeeping. It is not essential for our method to work and could be removed at the cost of a more cumbersome presentation.
We call these eigenvalues E The next assumption is a key characteristic for the physical situation we want to describe. We assume that the reservoir dynamics has a single stationary state and that the coupled dynamics is dissipative on the orthogonal complement of this state.

(A2)
The reservoir dynamics has a unique stationary state Ω R ∈ H R , that is, On the orthogonal complement the full, coupled dynamics is dissipative in the sense that there exists a λ * > 0 and a dense set D ⊂ H such that ∀λ with 0 ≤ |λ| < λ * and ∀φ, ψ ∈ D, Here C − = {z ∈ C : Imz < 0} is the open lower complex half plane and is the reduced resolvent. Here, RanP ⊥ R denotes the restriction of an operator to the subspace RanP ⊥ R . In (2.6), (2.7), C 1 is well defined (finite) on D × D. Discussion of Assumption (A2).
(i) The estimate (2.7) is of technical nature, but the estimate (2.6) is key as it implies that the dynamics generated byL λ ≡ P ⊥ We prove (2.8) at the end of this section.
(ii) We argue that (2.6) is a natural assumption. Namely, (2.6) for λ = 0 is implied by . This is readily seen by writing where P S,e is the spectral projection of L S associated to e, so that e∈E 0 P S,e = 1l S . In turn, (2.9) is a natural assumption on the dynamics of a reservoir, since that dynamics should be dissipative away from the stationary state. In concrete applications, one starts with (2.6) for λ = 0 and then proves its validity for small λ = 0 by perturbation theory [27].
(iii) In some recent works on the dynamics on open quantum systems [21,26], the assumption (2.6) is replaced by the condition that z → φ, R P R z (λ)ψ have a meromorphic continuation from z in the lower complex plane across the real axis into the upper plane. This is a much stronger condition than (2.6). In applications to open quantum systems, this difference means that the reservoir correlation function has to decay exponentially quickly in time for the meromorphic situation, while under the present assumption, the decay only needs to be polynomial, c.f. [27].
(iv) Without loss of generality, we may assume that As is well known (see e.g. Proposition 4.1 of [8]) if A is a self-adjoint operator and for each vector φ in some dense set, there exists a constant C(φ) such that then the spectrum of A in the interval (a, b) is purely absolutely continuous. Thus the estimate (2.6) with j = 0 implies that the spectrum of P ⊥ R L λ P ⊥ R acting on RanP ⊥ R is purely absolutely continuous and for λ = 0, this implies that the spectrum of L 0 reduced to RanP ⊥ R is purely absolutely continuous. On the finite dimensional part RanP R ∼ = H S , the operator L 0 is (identified with) L S which has pure point spectrum E 0 . Therefore L 0 has absolutely continuous spectrum except for the eigenvalues E 0 , the same as those of L S and the eigenprojection P e of L 0 associated to e is given by where P S,e is the eigenprojection associated to e as an eigenvalue of L S .

Instability of embedded eigenvalues and how to track them
If e is an isolated eigenvalue of L 0 then by standard analytic perturbation theory [19] L λ has eigenvalues E (s) e (λ) close to e, for λ small. Those eigenvalues coincide with the eigenvalues of the operator .11) is not defined as a bounded operator. However, its regularization (L 0 − e + i ) −1 P ⊥ e certainly is, for any > 0. We may hope that in a sense, the perturbation expansion (2.11) stays valid also for embedded eigenvalues, upon regularizing the resolvent and taking → 0 + . But due to the regularization, the operator (L 0 − e + i ) −1 P ⊥ e is not self-adjoint any longer! So according to (2.11), the second order corrections to the embedded eigenvalue would become complex numbers. This is, however, not compatible with L λ being self-adjoint. We may then intuit that for small nonzero λ, the number m e of eigenvalues of L λ close to e might be strictly reduced, m e < m e , and that this reduction is accounted for by the existence of complex eigenvalues of the so-called level shift operator Here, (L 0 − e + i0 + ) −1 is the limit of (L 0 − e + i ) −1 as → 0 + , taken in the sense of the operator norm in the expression (2.12). This mechanism has the following precise formulation. Let Q be an orthogonal projection on H and set Q ⊥ = 1l − Q. The Feshbach map applied to an operator A on H is defined by where it is assumed that Q ⊥ AQ ⊥ RanQ ⊥ is invertible. The Feshbach map satisfies the following isospectrality property: Let a ∈ C be in the resolvent set of the operator Q ⊥ AQ ⊥ RanQ ⊥ , so that F(A − a; Q) is well defined. Then the isospectrality property says that a is an eigenvalue of A if and only if zero is an eigenvalue of F(A − a; Q). The redeeming quality of this mapping is that F(A − a; Q) acts on RanQ, a space smaller than H (reduction in dimension). Consider now, for z ∈ C − ,

14)
where R Pe RanP ⊥ e is the reduced resolvent. The isospectrality is not of any good use to analyze the spectrum of L λ directly, since R Pe z (λ) is not defined for real z. However, one can show (see Theorem A1 of [20] and also Theorem 4.1 below) that the condition (A2), together with the assumption 2 that (A3) IP e is a bounded operator and RanIP e ⊂ D, (2.15) imply that the derivatives of order up to two of denotes the minimal gap of the eigenvalues of L 0 (which is the same as that of L S ). The Feshbach map (2.14) is then well defined (by continuity) for z ∈ R, |z − e| ≤ g/2. Now the isospectrality property can be extended to real values of z. Namely, one can show (see Theorem 6.1 below and also [11,3,20]) that any E ∈ R is an eigenvalue of L λ if and only if zero is an eigenvalue of F(L λ − E; P e ), and that Due to (2.14), eigenvalues of F(L λ − E; P e ) are located in an O(λ) neighbourhood of E 0 , and hence by isospectrality, so are those of L λ . The multiplicity of E as an eigenvalue of L λ is controlled by (2.18). We assume now as this condition does not alter the emergence of complex eigenvalues (because P e IP e is selfadjoint). We could dispense with the condition (2.19) by a simple modification of our arguments. In applications [26,27,28], this assumption is naturally satisfied. According to (2.12), (2.14), (2.19), As it acts on the m e -dimensional space RanP e , the operator Λ e has m e generally complex eigenvalues. Since Λ e it is a dissipative operator, meaning that However, Λ e may have real eigenvalues without L λ having any eigenvalues close to e (for λ small, nonzero). This is so since the O(λ 3 ) term in (2.20) may cause the spectrum of (2.20) to be non-real. In this case, L λ does not have any eigenvalues close to e, according to (2.18).
To simplify the analysis, we do not consider this higher order effect. Instead, we assume that the real eigenvalues of Λ e are in bijection with the eigenvalues E (s) e (λ) of L λ close to e. One way to ensure this is to impose the condition The eigenvalues of Λ e are simple and Λ e has exactly m e real eigenvalues.
We recall that m e , defined before (2.4), is the number of (distinct and simple) eigenvalues of L λ close to e, for small λ = 0. Assuming Λ e to have purely simple spectrum is done for convenience of the presentation. This restriction can be removed easily and our approach still works, as long as Λ e is diagonalizable. In some models, it can happen though that Λ e is not diagonalizable at so-called exceptional points of parameters; one then expects a qualitatively different behaviour of the dynamics (Jordan blocks of Λ e cause polynomial corrections to exponential decay in time). We do not further explore this interesting aspect here. The condition (A5) without the simplicity assumption is also called the Fermi Golden Rule Condition. It ensures that stability or instability of eigenvalues of L λ is detected at the lowest order (λ 2 ) in the perturbation. In terms of dynamical properties, the condition (A5) means that metastable states have life-times of O(λ −2 ), as shown in Theorem 3.1.

How to link stable and unstable eigenvalues to the dynamics
We have seen above that the (partially) stable eigenvalues of L λ are the real eigenvalues of the level shift operators Λ e and that the Λ e may have eigenvalues a (s) e with strictly positive imaginary part. We introduced the level shift operators Λ e as the second order (λ 2 ) contributions in the Feshbach map, (2.20). The isospectrality property of the Feshbach map then linked Λ e to the spectrum of L λ . As it turns out, the Feshbach map is also an ingredient in the block decomposition of an operator H acting on H = RanQ ⊕ RanQ ⊥ , where Q is an orthogonal projection. More precisely, one can readily verify that for any operator H such that F(H; Q) exists, the following identity holds (see Section 6, in particular (6.5)), (2.23) The first component in the 2 × 2 decomposition of (2.23) is that of RanQ, the second one that of RanQ ⊥ . Choosing in (2.23) Q = P e and H = L λ − z for z ∈ C\R gives the decomposition The first term on the right side is the inverse of the Feshbach map, an operator acting on RanP e and the last term acts on RanP ⊥ e . The resolvent, (2.24), is linked to the propagator via the Fourier-Laplace transform, for w > 0, Using (2.24) in (2.25) provides a link between the Feshbach map and the dynamics. We then expand the (inverse) of the first term on the right side of (2.24), e at the eigenvalues of e + λ 2 Λ e . Upon "integration around" those poles, in accordance with (2.25), one can extract the dynamical factors e it(e+λ 2 a (s) e ) . This procedure works locally, that is, for z close to a fixed e. So we will subdivide the integration contour in (2.25) into regions close to e ∈ E 0 , for each e, and apply the Feshbach map with the appropriate P e . On the complement, where z is away from all eigenvalues of L 0 , a similar analysis is done using the Feshbach map with projection P R . This is the outline of the idea, and we refer to Section 4 for the detailed analysis.
Proof of (2.8). Let φ, ψ ∈ RanP ⊥ R and use the resolvent representation (Cauchy formula or Fourier-Laplace transform), where w > 0 is arbitrary. We have for φ, ψ ∈ Dom(L λ ). We now use e itz = (it) −1 ∂ z e itz and integrate in (2.28) by parts k times, It now follows from (2.29) and (2. In this way, smoothness of the resolvent gives rise to dissipation in the dynamics. This proves (2.8).

Main result
Our main result is Theorem 3.1 below. It involves a reduced system dynamics e itM (λ) , and we explain the generator M (λ) now. The simplicity of the spectrum assumed in (A5) implies that the level shift operators Λ e is diagonalizable, The a Given an eigenvalue e ∈ E 0 of L 0 , we partition the indices s = 1, . . . , m e into the oscillating and decaying classes We define the operators e we view each M e as an operator on RanP S,e , the eigenspace of L S associated to the eigenvalue e. The operator M (λ) is then also an operator acting on the system Hilbert space H S alone. By construction, M (λ) and L S commute, and M (0) = L S .

5)
for constants C and K independent of t ≥ 0.
If interested in φ, e −itL λ ψ , t ≥ 0, then it suffices to take the adjoint of (3.4). The choice of the sign in the exponent on the left side of (3.4) is motivated by the application of our result to the reduced system dynamics, [27]. The generator M (λ) is defined in its diagonalized form (3.3) and so e itM (λ) is easily obtained from the functional calculus: The equality (3.4) means that for all φ, ψ ∈ D, . It is clear that we need to include the exact phases (eigenvalues E (s) e (λ) of L λ to all orders in λ) in the approximate dynamics on the right side of (3.4), for otherwise we have diverging phase differences, and we cannot get an approximation small in λ for all times. However, in the decaying terms, we can truncate the exponents at the second order in λ, keeping only e + λ 2 a (s) e , and the time decay still allows us a control uniform in t.
(ii) On the finite-dimensional part, the dynamics e itM (λ) ⊗P R on the right side of (3.4), there are terms oscillating in time, ∝ e itE (s) e (λ) , and terms decaying in time, ∝ e it(e+λ 2 a (s) e ) . The decay is exponential, with rates λ 2 Ima (s) On the infinite-dimensional part (i.e., on RanP ⊥ R ), the dynamics is also decaying, but only at a polynomial rate 1/t 2 , as per (2.8).
(iii) The asymptotic dynamics for times larger than all the life times (λ 2 Ima (s) e ) −1 and large enough so that the dissipative term φ, ) the quasi-static and dissipative parts, respectively. The decomposition (3.4) is reminiscent of the structure of solutions of dispersive partial differential equations, which split into a dispersive wave plus a part converging towards an invariant manifold [34,39]. In our case, after the decay of the dissipative term, the dynamics first approaches the quasi-static manifold spanned by the Q (v) The approach we take to prove Theorem 3.1 is similar to [20]. In that paper, an expansion of the dynamics was given as a sum of a main part plus a remainder, similar to (3.4).
The remainder of [20] converges to zero in the limit of large times, but it is not shown to be small in λ. In the current work, the remainder is small in λ for all times (3.5), but it does not converge to zero for large times. This means, incidentally, that the main term of the current work has a simpler form than the one of [20]. Indeed, the current main term is a second order approximation along the decaying terms, with s ∈ S dec , of that of [20]. The result of [20] is not suitable for proving the validity of the Markovian approximation, but the current result here is.

Uncorrelated initial states. For vectors of the form
We then obtain the following result directly from (3.4).
This corollary is the starting point for a proof that the Markovian approximation is valid for all times, which we give in [27]. We show that M (λ) gives rise to a completely positive trace-preserving semigroup acting on system density matrices, which coincides with the wellknown Gorini-Kossakovski-Sudarshan-Lindblad semigroup [1,6,26]. We point out, though, that correlated intitial states are also treatable with our method -indeed, the main result, Theorem 3.1, does not make an assumption on intital SR correlations. We explore this in [28] (see also Section 1.4).
Parameter dependence. We take some care in estimating c 0 , C and K appearing in Theorem 3.1 and Corollary 3.2 in terms of various parameters, as we explain now. The analysis is based on a perturbation theory (λ small) and the corrections of the spectrum are governed by the level shift operators Λ e , whose spectral decomposition is given in (3.1). To quantify the perturbation theory, we introduce the quantities a = min We also define is an orthonormal eigenbasis of L S . The parameter κ 1 is well defined since due to (2.15), the vectors Iϕ m ⊗ Ω R are in D. We combine all these constants and the spectral gap g of L 0 (see (2.17)) into the following two effective constants, κ 0 = max 1, 1/g, κ/a, ακ, κ 1 (1 + κ 1 )(1 + g + 1/g), (3.14) Furthermore, for vectors φ, ψ ∈ D such thatL λ φ, In (3.15) we use the notation for any function E of φ and ψ. and C = C κ 0 , where c , C do not depend on λ, g, a, α, κ, δ and λ 0 , κ 0 are given in (3.14) and (3.13), respectively. The constant K in (3.5), which is a function on D × D, is given by (3.15).
Remark. The parameter κ 1 , (3.12), is defined in terms of the function C 1 , (2.6). The latter describes the dissipative nature of the interacting dynamics and, a priori, depends on λ. However, the dissipation is due to the nature of the environment alone and so in applications, one can bound κ 1 uniformly in λ, for small enough λ (see [27]). See also the point (ii) in the discussion of condition (A2) above. In this sense, κ 1 and thus λ 0 , are independent of λ.

Proof of Theorem 3.1
We are presenting the core strategy in Section 4.1. It involves estimates which we derive in Sections 4.3 and 4.4. We then collect the estimates and implement the outlined strategy in Section 4.5. As we explained after Theorem 3.1, our approach is similar to that of [20]. However, in order to be able to show that the remainder is small in λ for all times, (3.5), we need to substantially modify the analysis of [20]. The main alteration is the introduction of two new energy scales η and ϑ which allow for a more detailed estimation of the resolvent. See also Fig. 2 and Fig. 4.

Notational convention.
For X, Y ≥ 0, we write X ≺ Y to mean that there is a constant C independent of the coupling constant λ as well as the parameters g, a, α, δ, κ, such that X ≤ CY .

Strategy
Given an orthogonal projection Q on H we set Q ⊥ = 1l − Q and we denote the resolvent and reduced resolvent operators by The resolvent representation of the propagator is

2)
valid for any w > 0 if either of ψ or φ belong to Dom(L) (see [14]). Given an orthogonal projection Q, the resolvent has a decomposition into a sum of three parts: 4 a part acting on RanQ, one acting on RanQ ⊥ and one part mixing these subspaces. This decomposition reads, for z ∈ C\R, Here, F is the Feshbach map (2.13), The operator B(z) in (4.3) is given by (4.5) Here and in what follows, it is convenient to simply write R Q z for R Q z (λ). The existence of F(z) −1 is automatic, this operator equals QR z (λ)Q, see (4.3). It is elementary to establish the relation (4.3), see Section 6 for more detail.
Our strategy is to partition the integration contour R−iw in (4.2) into segments, see Fig. 2, and apply (4.3) with suitable projections Q on each segment. To describe the partition, fix an η > 0. It will be chosen as a certain positive power of |λ| later on, so η is a small parameter (relative to the eigenvalue gap g of L 0 , (2.17)). For every eigenvalue e of L 0 define the segment (see Fig. 2) where w > 0 is the parameter in (4.2) and E 0 is the set of all eigenvalues of L 0 . We have the disjoint decomposition On G e we will apply the Feshbach map with projection Q = P e , (2.10), and on G ∞ we use the projection Q = P R , (2.5). Accordingly, it is important that we have dissipative bounds on the reduced resolvents, i.e., the following local limiting absorption principles, which follow from the global ones given in (2.6), (2.7).

Cheat sheet
We keep track of several constants during the estimates to follow. The following cheat sheet is presented for the convenience of the reader.

Analysis of J e (t)
Fix an eigenvalue e ∈ E 0 of L S and choose Q = P e in the Feshbach decompositon (4.3), where P e is the eigenprojection (2.10). The Feshbach operator (4.4) reads

.16)
For z = e and λ = 0, A e (z, λ) reduces to the level shift operator Λ e . We show in Lemma 5.1 below that for z ∈ C − with |z − e| not too large and λ not too large, A e (z, λ) has simple eigenvalues (a property inherited from Λ e ) and we denote the spectral representation by Here a The main result of this section is Then where S e ≺ (4.21) with the various constants defined in Section 4.2.
Proof of Proposition 4.2. We show (4.20) by estimating separately the three contributions coming from the decomposition (4.3) of the resolvent R z , The corresponding bounds are presented in (4.65), (4.83) and (4.84) below.
Estimate of D 1 in (4.23). Our main estimate for D 1 is given in (4.65) below. All estimates in this section are controlled in the operator norm sense for operators on RanP e . We have where the properties of A e (z, λ) = −P e IR Pe z IP e , are discussed in Section 5.2. It follows from (5.15) that . (4.25) • Consider (4.25) for an s ∈ S dec e . The case s ∈ S osc e is dealt with below. Let z ∈ G e . We show that Q where (recall the definitions (3.8) and (3.9) of a, α) To show (4.26), we start with the expression for T (z, λ), (4.28) Recall that a (s) e (e, 0) = a (s) e is an eigenvalue of Λ e . We have where we have used (5.17) (see also Corollary 5.2). Now |e − z + λ 2 a (s) 2 η (we will take w → 0 at fixed η, so w < η without loss of generality)). We get  The combination of (4.32) and (4.34) gives (4.27).
• Next we consider a term in (4.25) with s ∈ S osc e . We have to modify the above argument, since it used that the imaginary part of a (s) e was strictly positive, which is not the case any more now. By the isospectrality of the Feshbach map, there is exactly one eigenvalue E (s) Moreover, the projection Q Using (4.35), the second term on the right side of (4.36) becomes (4.37) We estimate the last fraction in (4.37) by using (5.17) as where the last relations holds since λ 2 κ 4 ≺ 1 (see (4.19)). Estimates  Consider now the last term on the right side of (4.36). Using (4.35), the fraction reads (4.40) The second fraction in (4.40) is ≺ 1 by (4.38). Next, we obtain from point 3. of Lemma 5.1   e (λ), λ)| ≺ κ 4 and so the second fraction in (4.43) is bounded ≺ 1 since λ 2 κ 4 ≺ 1 (see (4.19)).
where the remainder satisfies We now analyze D 1 , see (4.23). From (4.48), where the remainder S 1 is due to the integration of T (z, λ), estimated by (take into account that |G e | = 2η) We used (4.19) in the second estimate. Now where A e is given in (4.18). Note that by (4.35) and (5.17), |E (4.19)), we can solve for (see (3.9)): In the last estimate, we used that (see (4.19)). Thus (4.51) gives where A e is the contour (4.18). For z ∈ A e we have |z − e| ≥ η and so λ 2 |a (s) e (e, λ)/(e − z)| < 1 (because we have λ 2 (α + |λ|κ 4 )/η < 1). Thus the geometric series converges, In the same way we obtain Next we replace the projections in the remaining integrals in (4.61) by their values at λ = 0: (4.19). Upon replacing the projections in the integrals of (4.61) by their values at λ = 0, we thus make an error ≺ (w+η)e wt η κ 3 (|λ| + η). Once this replacement is made, we use the fact that where We now further analyze the integral in (4.62). Denoting P ⊥ S,e = 1l S − P S,e , where P S,e is the spectral projection of L S onto the eigenspace associated to e (see (2.10)), we have (as η < g/2) where g is the spectral gap of L S , (2.17). So we can replace P e = P S,e ⊗ P R in (4.62) by simply P R , incurring an error of size of the right hand side of (4.64), Estimate of D 2 in (4.23). Our main estimate is given in (4.83) below. We have from (4.5) We begin by analyzing the first term on the right side of (4.67). Using (5.15), (4.68) • For s ∈ S dec e we have from (4.30) To get the last bound, we took into account that P e = P e P Ω , that IP e ⊂ D, so that from (5.30) and (4.8), we obtain • For s ∈ S osc e we use (4.35) to get (4.72) The first two terms of the right side of (4.72), when used in (4.71), are estimated by taking into account that the second factor on the right side of (4.71) is ≺ 1 due to (4.38). However, the last term on the right side of (4.72) has a singularity in z which must be removed by integrating over G e . To do this, we write The denominator of the second summand on the right side is 1 and its numerator is ≺ (4.19)) and so by (4.73) and (4.54), we have ψ | ≺ φ max j C 2 (Iϕ j ⊗ Ω, ψ). Next, by (5.16), (4.8) and again by the argument of (4.70), the last two terms on the right hand side of (4.75) has an upper bound ≺ |λ|η φ (1 + κ 3 ) max j C 2 (Iϕ j ⊗ Ω, ψ). We thus obtain from give The integral Ge e itz φ, λR Pe z IP e F(z) −1 ψ dz in (4.67) has the same upper bound (4.77), but with φ and ψ interchanged.
• We now estimate the term in (4.67) involving λ 2 R Pe z IF(z) −1 IR Pe z . We use again (4.25). Proceeding as above, we obtain from (4.30) that for s ∈ S dec e , and and similarly, P e IR Pe z ψ ≺ max j C 2 (Iϕ j ⊗ Ω, ψ). Furthermore by (5.30), Q (s) e (z, λ) ≺ κ. Combining these estimates with (4.78) yields that for all s ∈ S dec e , Next we treat the case s ∈ S osc e , for which we show the bound (recall that e = E (s) as follows. We proceed as in (4.71) and (4.72) and get terms with [Q We find that they are ≺ λ 2 (η + w)e wt (1 + κ 3 ) max j C 2 (Iϕ j ⊗ Ω, φ) max j C 2 (Iϕ j ⊗ Ω, ψ).

The only remaining term to consider is
where we have used (4.74). This yields (4.81).
We collect the estimates (4.77), (4.81) and use them in (4.67) to get

Analysis of J ∞ (t)
We introduce a new parameter ϑ > 0 which will be chosen suitably small (a power of |λ|) below, and which is depicted in Fig. 4. The goal of this section is to show the following result. (4.85) In the last line, we use the notationL λ ≡ P ⊥ Proof of Proposition 4.3. Consider two adjacent eigenvalues e < e of L S and set We introduce a new small parameter ϑ > 0 and set D 1 = {x − iw : e + η ≤ x ≤ e + η + ϑ} and D 2 = {x − iw : e + η + ϑ ≤ x ≤ e + e −e 2 }, as depicted in Fig. 4. We use the Feshbach decomposition (4.3) with Q = P R = |Ω R Ω R | and accordingly, we get three contributions to the resolvent as in (4.23), Estimating D 1 of (4.87). The main bound we prove is given in (4.95). We have, in the sense of operators acting on RanP R , For z ∈ G e ∞ , (L S − z) −1 ≤ 1/η and by (5.2) we get λ 2 P R IR P R z IP R (L S − z) −1 ≺ κ 1 λ 2 η and so we get from (4.88), since κ 1 λ 2 η ≺ 1, (4.89) We consider now z ∈ D 1 ∪ D 2 (the region closer to e than to e , see Fig. 4). By inserting 1l S = P S,e + P ⊥ S,e , using that P e = P S,e ⊗ P R and that P ⊥ S,e (L S − z) −1 ≤ 2/g (see (2.17) for the gap g), we estimate the term in (4.89) involving the resolvent R P R z as follows, By (5.2) we have P e IR P R z IP e − P e I(R P R e | λ=0 )IP e ≺ κ 1 (|z − e| + |λ|). Thus, taking into account the definition (2.12) of the level shift operator Λ e and using that P ⊥ R IP e = P ⊥ e IP e , we have The reason for introducing the new small parameter ϑ (see Fig. 4), is that for z ∈ D 1 , the right hand side of (4.91) is bounded above by κ 1 (η + ϑ + |λ|). It then follows from (4.91) and (4.90) that ∀z ∈ D 1 , and hence, as Λ e ≺ ακ, we get Next we find a bound for the integral on the left side of (4.92) with D 1 replaced by D 2 . For z ∈ D 2 , we have (L S − z) −1 ≤ (η + ϑ) −1 and so it follows from We now combine (4.92) and (4.93) with (4.89), For z ∈ G e ∞ lying closer to e than e (i.e., z ∈ G e ∞ \(D 1 ∪ D 2 )) the analysis is the same, as only the distance from x to the nearest eigenvalue of L S plays a role in the estimates. We conclude that the bound (4.94) holds with D 1 ∪ D 2 replaced by G e ∞ in the integrals on the left side: Estimating D 2 of (4.87). Our main etimate is given in (4.102). According to (4.5), the term B(z) gives three contributions. The ones involving only one resolvent R P R z are all estimated in the same way as follows. We have for all z ∈ C − (use (2.6)), According to (4.89), the main term of F(z) −1 is (L S − z) −1 the norm of which is bounded above by 1/η for z ∈ D 1 and by 1/(η + ϑ) for z ∈ D 2 . Then, The terms in (4.89), which are of order two and higher in λ, are estimated as for z ∈ G e ∞ . It follows from (4.89), (4.97) and (4.98) that (4.99) The same upper bound is achieved for the term involving R P R z λIP R F(z) −1 in (4.5). To deal with the term in (4.5) involving the resolvent twice, we use the bound From (4.89), F(z) −1 ≺ 1 η [1 + λ 2 κ 1 /η] for z ∈ G e ∞ , which gives, combined with (4.100) and (4.99) the result Estimating D 3 of (4.87). Recalling the definition (4.86) of G e ∞ and using (2.6), we extend the integration by an amount of 2η to e ≤ x ≤ e , Estimates on the infinite parts of G ∞ . So far in this section, we have dealt with all z between any two eigenvalues of L S , see Figs. f1 and f3. Let e + and e − be the largest and smallest eigenvalues of L S , respectively and set We choose Q = P R for the projection in the Feshbach decomposition (4.3). We use (4.88) to expand which converges for λ 2 κ 1 (L S − z) −1 ≺ 1. For z ∈ G + ∞ , we have (L S − z) −1 ≤ 1/η and so, since λ 2 κ 1 /η ≺ 1, the series (4.106) converges uniformly in z ∈ G + ∞ an uniformly in w. Splitting off the first term in the series (4.106) gives Cκ 1 λ 2 η n ≺ e wt κ 1 λ 2 η , (4.108) since κ 1 λ 2 η ≺ 1. The analysis and estimates for z ∈ G − ∞ are the same. We conclude that Next we need to estimate G + is the sum of three terms according to (4.5). In each term, we use (4.106) to split off the main part. Accordingly, we have which defines the quantity B 2 (z). Proceeding as in the derivation of (4.108), we get Now we analyze the integral associated to the first term in (4.110), which is given by −λ G + ∞ e itz φ, (L S − z) −1 P R IR P R z ψ dz. We split the integration domain into e + + η ≤ x ≤ e + + η + ϑ and x ≥ e + + η +ϑ. On the compact domain, the integral has the bound ≺ e wt |λ|ϑ η φ max j C 1 (Iϕ j ⊗Ω R , ψ), where the factor 1/η is due to the resolvent (L S − z) −1 and ϑ is the length of the interval of integration. On the infinite integration domain, we need to control the integrability for large Of course, we get the same upper bound for −λ G + ∞ e itz φ, R P R z IP R (L S − z) −1 ψ dz (see (4.110)), with φ and ψ exchanged on the right side of (4.114). Finally, we estimate The last integral is ≺ 1/η. Collecting the bounds (4.111), (4.114) and (4.115) and using them in (4.110), we obtain (4.116) Next, extending the integration domain by 2η similar to (4.103), we estimate where the contour Γ appearing in (4.120) is given by Γ = e A e ∪ G ∞ . It is represented in Fig. 5. The error term S(w, φ, ψ) in (4.120) is the sum of the error terms in (4.20) and (4.85), satisfying .

All eigenvalues of
3. On the domain of z, λ determined by the constraint stated at the beginning of the lemma, we have The previous result leads readily to the fact that one can extend the eigenvalues and eigenprojections as functions of z continuously to the real axis:  Proof of Corollary 5.2. Let z n be a sequence in C − with |z n − e| < δ κ 2 κ min{1, 1/κ} and z n → x ∈ R. Then (5.16) shows that Q (s) e (z n , λ) is a Cauchy sequence, so it converges. Next, again by (5.16), Q e (z n , λ) exists and is ≺ κ 3 |z − x|. If z = x ∈ R then approximate it by z n ∈ C − and the above argument works the same. The argument is also the same for the a since |z − e| + |λ| ≺ δ/(κ 2 κ 2 ). By standard perturbation theory [19], the ranks of Q Taking the z derivative twice (or the λ derivative) and proceeding as in (5.26) yields the estimates We combine (5.26) and (5.27) to obtain e (z, λ)} ≺ κ 2 κ(1 + κ 2 κ/δ). Proceeding in the same way we find |∂ λ a (s) e (z, λ)| ≺ κ 2 κ(1 + κ 2 κ/δ). By integrating these bounds similarly to what we did in (5.28), we get (5.17).
We finally prove point 3. The relation (5.29) yields In the case O = H − z, for a self-adjoint H (like H = L λ ) and z ∈ R, we have indeed that D = Q ⊥ (H − z)Q ⊥ RanQ ⊥ is invertible. Also, A − BD −1 C = F(H − z; Q) and hence the sum of X, Y, Z and W give the right hand side of (4.3).
We now show that χ is Cauchy. This part of the analysis is inspired by [DJ], Theorem 3.8. Denote by p the spectral projection of F(H − E; Q) associated to the eigenvalue 0. In particular, pϕ = ϕ. Since F(H − E; Q) is a dissipative operator and 0 is on the boundary of its numerical range, 0 is a semisimple eigenvalue and p is the orthogonal projection onto the kernel of F(H − E; Q), in particular, p * = p. First we show that z → pHR Q z Hp is a continuously differentiable function of z ∈ C − ∪ C + ∪ {E}. (Note: QHR Q z HQ, without the restriction to the range of p, is not even continuous at E, because the limits coming from C + and C − do not coincide.) We point out that z → QHQ ⊥ R Q z Q ⊥ HQ is continuously differentiable on z ∈ C + ∪ {E}, which follows by taking the adjoint in the assumption of the theorem. Next, we have pF(H − E; Q)p = 0 and hence p(H − E)p = pHR Q E−i0 + Hp, and taking the adjoint gives p(H − E)p = pHR Q E+i0 + Hp. So we have which shows that z → pHR Q z Hp is continuous at z = E and its value at z = E is selfadjoint. Next, let γ(τ ), τ ∈ [−1, 1], be a smooth curve in C − ∪ {E} which is tangent to the point E ∈ C, satisfying γ(0) = E and γ (0) = 1. Then τ → pHR Q γ(τ ) Hp is continuously differentiable at τ = 0 and d dτ | τ =0 Im pHR Q γ(τ ) Hp = Im pH(R Q E−i0 + ) 2 Hp. On the other hand, this last derivative has to be equal to zero for the following reason. For all τ we have Im pHR Q γ(τ ) Hp ≤ 0 and from (6.7), Im pHR Q γ(τ ) Hp = 0 at τ = 0. If d dτ | τ =0 Im pHR Q γ(τ ) Hp was > 0 or < 0 this would contradict the fact that Im pHR Q γ(τ ) Hp ≤ 0. Since this derivative vanishes, we get Im pH(R Q E−i0 + ) 2 Hp = 0, so pH(R Q E−i0 + ) 2 Hp is self-adjoint, which means that pH(R Q E−i0 + ) 2 Hp = pH(R Q E−i0 − ) 2 Hp. So z → pHR Q z Hp is continuously differentiable on C − ∪ C + ∪ {E}.
We now use this regularity property to show that χ converges. Let , > 0. We have From the resolvent identity, R Q E−i R Q E+i = 1 i( + ) (R Q E−i − R Q E+i ) and similarly for the other three terms on the right side of (6.8). Therefore, (also using that pQ = Qp = p), Hpϕ − ϕ, pH Hpϕ + ϕ, pH Hpϕ . (6.10) Since z → pHR Q z Hp is a continuously differentiable function of z ∈ C − ∪ C + ∪ {E}, the right hand side of (6.10) converges to zero as , → 0. Thus χ converges.
was supported by a Discovery Grant from the National Sciences and Engineering Research Council of Canada (NSERC).

Conclusion
In this paper we consider the evolution of a finite-dimensional system coupled to a reservoir. The coupling is small, but fixed. The reservoir is assumed to have a unique stationary state and otherwise dissipative dynamics. We develop a method to separate the total, unitary evolution of a SR complex into a Markovian part, a dissipative part, plus a remainder. The Markovian part describes the system dynamics with the reservoir in its stationary state and contains explicit oscillating (in time) and exponentially decaying terms. The dissipative term governs the dynamics on SR states which do not have any overlap with the system stationary state. This term decays polynomially in time. The remainder term is small in the coupling, uniformly for all t ≥ 0. The strengths of our results are: -The decomposition is derived mathematically rigorously, with a controlled remainder which is small in the SR coupling, for all times t ≥ 0.
-We can treat SR interactions with significantly less regularity than previously needed. In terms of the reservoir correlation function, this means that only polynomial decay is required (versus exponential decay in previous rigorous works).
-We analyze the whole SR dynamics, not only the reduced system dynamics and moreover, the initial SR states do not need to be of product form (they can be classically correlated or entangled).
-Our method applies to a widely used class of open systems: An N -level system coupled linearly to a spatially infinitely extended reservoir of thermal, non-interacting Bose particles. Detailed results on this model are presented in [27,28]. In particular, it is shown there that the system dynamics is well approximated by the Markovian master equation for all times t ≥ 0, even for initially correlated SR states. It is also shown there that the Born approximation is valid for all times.