Adiabatic Elimination and Sub-space Evolution of Open Quantum Systems

Efficient descriptions of open quantum systems can be obtained by performing an adiabatic elimination of the fast degrees of freedom and formulating effective operators for the slow degrees of freedom in reduced dimensions. Here, we perform the construction of effective operators in frequency space, and using the final value theorem or alternatively the Keldysh theorem, we provide a correction for the trace of the density matrix which takes into account the non trace-preserving character of the evolution. We illustrate our results with two different systems, ones where the eliminated fast subspace is constituted by a continuous set of states and ones with discrete states. Furthermore, we show that the two models converge for very large dissipation and at coherent population trapping points. Our results also provide an intuitive picture of the correction to the trace of the density matrix as a detailed balance equation.


I. INTRODUCTION
The adiabatic elimination method allows to reduce the dimensionality of a problem by discarding fast degrees of freedom and describing only the dynamics of the slow ones. Adiabatic elimination has played an important role in unifying dynamical patterns observed in very different phenomena, from laser and fluid dynamics to biological and chemical systems [1,2]. It has allowed to reduce these apparently very different problems to similar minimal sets of coupled differential equations. In quantum systems, adiabatic elimination dates back to the 1960s in atomic physics, with the development of a theory of the maser and laser which includes the quantum noise due to the spontaneous emission process [3]. It has also been essential to understand the mechanisms responsible for atom cooling [4].
While these first applications were concerned with dissipative systems, it seems that in the quantum arena, the adiabatic elimination procedure has been popularized mainly in the case of conservative Hamiltonian systems [5][6][7] and, in particular, in many-body systems [8,9], where it allows one to obtain effective Hamiltonians and open new perspectives for quantum simulations [9].
Meanwhile, the concept of quantum open systems has emerged and it is now taking over Hamiltonian systems as the elementary brick for the description of a quantum system. A quantum open system consists of subsystems interacting with its environment. Its state is described by the density operator, where the degrees of freedom of the bath have been traced out [10]. Among quantum open systems, the ones whose dynamics follows a one-parameter semigroup play a special role. Indeed, since the work of Lindblad, Gorini, Kossakowski, and Sudarshan [11,12], the form of its generator, the so-called Lindblad operator, is completely specified. Furthermore, this specific evolution is the one followed by a quantum subsystem interacting with a Markovian environment. The concept of open quantum system constitutes a first reduction by eliminating (tracing out) the bath degrees of freedom. Indeed, from a very high dimensional Hamiltonian dynamics, we end with a Lindblad dynamics in a Hilbert space of a smaller dimension. But even this reduced description can be cumbersome [13] and to get at least the steady states and the dynamics around these steady states can be very difficult and computationally intensive.
When this reduced-system Lindblad dynamics presents two different timescales, it should be useful to separate the fast-evolving degrees of freedom from the slow ones, that is, to perform an adiabatic elimination. In most cases, there is a unique steady state, and the adiabatic elimination consists in obtaining the dynamics in the proximity of the stationary state, where the fast part has already reached a stationary state while the slow part is still evolving to the steady state. In this way, the adiabatic approximation becomes a "long" time approximation, long with respect to the time needed for the fast part to reach a steady-state behavior. The main objective is then to be able to describe the dynamics of the slow part without the need to refer to the fast one.
To our knowledge the first work which addressed a general formalism to perform the adiabatic elimination with Lindbladian dynamics is the one by Mirrahimi et al. [14]. The main idea of this work and subsequent ones [15][16][17][18][19] from the QUANTIC group consists in preserving the Lindblad structure for the generator of the slow dynamics. To this end they built a bijective map from the exact density matrix to the couple of density matrices corresponding to fast and slow motions. Using singular perturbation theory [20][21][22], they are able, in principle, to obtain the slow motion at any given order of approximation. One of the main points is that the mapping is such that the dynamics of the slow density matrix is generated by an effective Lindblad operator. As a consequence, the dynamics of the slow density matrix is trace preserving.
With a completely different methodology, Reiter and Sørensen obtain an effective Lindblad operator which recovers the same result as obtained in [14] (up to an overall energy shift) for the case of a single excited state, but which can also be applied to more general systems where the energy level structure for the excited states takes into account arbitrary detunings [23].
We note that in these approaches the density matrix describing the slow part does not accurately describe the quantum state in the slow subspace when exchange of population between the fast and slow subspace cannot be neglected. Indeed, as the slow dynamics is described by a Lindblad operator, it is trace preserving and the initial population present in the slow subspace will remain in this subspace.
Adiabatic elimination for many-body systems, in particular for Rydberg atoms, has been addressed in [24,25] and relies mainly on perturbation methods applied to the Lindblad operator. In these works, the authors calculated the correction up to fourth order in the perturbation and concluded that the physical constraints of the solutions were only preserved to second order. Recently, Macieszczak et al. [26] recover a general formulation of long-time dynamics based on the eigenvalue decomposition of the Liouville operator and timedependent perturbation techniques, in order to describe a metastable manifold. A final application of adiabatic elimination techniques worth noting highlights its usefulness in finding conditions for dissipative-state preparation and noise suppression via interference effects. Recently an extension of Ref. [23] has presented an effective operator formulation including perturbations of the Hamiltonian and of the jump operators involved in the dissipative part of the Lindblad operators. They are able to show under very general terms how to understand and implement error correction strategies for steady-state subspaces of the Liouvillian [27]. Also, several publications have reported adiabatic eliminations in specific systems [28][29][30][31] but without a general recipe to make this approximation.
In this work, we follow an alternate route which consists in using Feshbach projectors P and Q = 1 − P [32] to develop a general strategy to approximate the evolution of Pρ(t ), the slow component of the quantum state ρ(t ) at time t. It is based on the projection PG(z)P of the resolvent G(z) = (z − L) −1 of the original Lindblad operator L in the slow subspace. We define L eff (z), a z-dependent operator defined on the slow subspace only, such that The operator L 0 = L eff (z = 0) is the analog of the effective Lindblad operator obtained previously by Mirrahimi [14] and Reiter [23]. Furthermore, we also show how to correct the trace-preserving evolution generated by L 0 to take into account possible population exchange between fast and slow subspace.
In this paper, we consider only the case where the projector P onto the space of operators themselves defined on H is built from a projector P onto the underlying Hilbert space H as Pρ = PρP, as in Refs. [14,23]. In others words, we assume that the fast/slow partition is linked to a partition of H in two complementary subspaces H = PH ⊕ QH. The application of our formalism to bipartite systems where the fast/slow partition is linked to a tensorial structure H = H slow ⊗ H fast will be the subject of a future publication. The partial trace over the fast subsystem can be formalized in terms of projector operators as in Ref. [33], for instance.
We apply our general result to several examples where the fast subspace is finite or infinite dimensional. In the last case, we consider that the Hamiltonian of the fast part has a continuous spectrum while the slow part has a discrete one. In other words, we address the problem of adiabatic elimination of the continuous set of states in dissipative Fano systems [34].
The generalization of Fano interferences from Hamiltonian to open quantum systems whose evolution is generated by a Lindblad operator has recently been the subject of great interest [35][36][37][38][39], in particular to describe mesoscopic systems or condensed matter systems. In the wideband approximation, corresponding to a "flat continuum," we are able to obtain the explicit expression for L eff (z) and therefore analyze in great detail the adiabatic approximation. In particular we show formally and numerically that in the limit where the fast dynamics reaches its steady state in a very short time, the Hamiltonian of the fast part can be approximated by a flat continuous spectrum.
The paper is organized as follows: In Sec. II the general formalism is developed and in Sec. III our general results are illustrated with several examples with excited-state continua or excited-state discrete levels. In Sec. IV we draw the connection between the effective Liouvillians of these two classes of systems, and conclude in Sec. V.

II. THEORY
The Hilbert space H of the system is partitioned into two subspaces with the help of two orthogonal projectors P and Q = 1 H − P, where 1 H is the identity operator on H. The QH subspace represents the fast degrees of freedom which reach a stationary regime in a short time. Our goal is to describe the slow motion in the subspace PH only, after the QH has reached its stationary state.
We suppose that the system is coupled to a bath that opens dissipation channels between QH and PH, or within PH and QH. Hamiltonian couplings (PHQ, or QHP) can also open transitions between PH and QH. Associated with P and Q, we define superprojector operators P and Q such that where 1 is the identity superoperator on the space of operator on H, and ρ is an operator on H.
We assume that the bath is Markovian so that the density matrix evolves according to a Lindblad equation [11,12]. For convenience, we will use the operator-vector isomorphism [40], which maps the operator |a b| in the Hilbert space H onto the vector |b ⊗ |a in the H ⊗ H Hilbert space, or equivalently maps any n × n density matrix ρ to a column vector ρ with n 2 elements, by stacking the columns of the ρ matrix. Under this isomorphism, the operation AρB † is mapped to B ⊗ A ρ, where A and B are operators on H and B denotes the complex conjugate of B; that is, where B † is the adjoint and B T is the transpose of B (see Appendix A). From now on, we drop the arrow in ρ as we assume that ρ is in vector form. The only exception is when a density matrix ρ is inside a bracket like in tr[ρ]. With this notation, the superprojector operators P and Q read Also, the general form of the Lindblad operator L, generator of the evolution,ρ = Lρ, can be written as [41] We start by expressing the density matrix evolution in an integral form through the Laplace transform: where G(z) = (z − L) −1 is the resolvent of L, and the integral on the complex plane is performed on a straight line D = {z ∈ C; Re z = a > 0}. Projecting Eq. (4) using P and Q gives In the remainder of the text, we make the assumption that at time t = 0 the population is entirely in the slow subspace PH so that Qρ(0) = 0. Hence the evolution in the PH subspace is simply given by We define the operator L eff (z), a z-dependent operator defined on PH, such that PG(z)P = [z − L eff (z)] −1 . Using the definition of the resolvent and the orthogonality of the P and Q projectors, we have where QG 0 (z)Q = [z − QLQ] −1 is the resolvent of QLQ. Equation (6) with Eq. (7) is an exact description of the dynamics (restricted to PH subspace) of a system coupled to a Markovian bath, and so is a completely positive map; however it is not trace preserving because the PH and QH partitions can exchange population during the evolution.
Generator of the slow dynamics. We notice that L 0 = L eff (z = 0) is the generator of the slow time dynamics. Indeed, projecting the Lindblad equationρ = Lρ on PH and QH we have To obtain the approximate slow time dynamics in the subspace PH, we assume that Qρ has reached a stationary regime, Qρ = 0. Using Eq. (9) to express Qρ as a function of Pρ, and inserting the result in Eq. (8), we obtain In Appendix E, we show a sufficient condition for L 0 to be the generator of a trace-preserving evolution. In all the examples we will present below, this condition is fulfilled. In addition, we have found, explicitly or numerically, that the operator L 0 is of Lindblad form. But we know that we are looking for a non-trace-preserving evolution as the total initial population may be distributed on Pρ and Qρ. We must then correct this evolution to take into account the possible variation of the trace of Pρ. To this end, we look for the exact final state, Mapping to the final state. By Eq. (10), we know that the final state ρ f , in PH subspace, is in the kernel of L 0 = L eff (z = 0). We assume that the kernel is one dimensional and define ρ its unique element with tr[ρ] = 1. Then ρ f = αρ, and we are left to determine α = tr[ρ f ]. The final stationary state ρ f can be obtained taking the limit of Eq. (6) when t → ∞. This limit can be obtained using the final value theorem: As we show in Appendix G, this limit can be calculated explicitly as where L 1 = dL eff (z) dz | z=0 . We notice that α given by Eq. (12) does not depend on the initial state ρ 0 . This is a consequence of assuming that the kernel of L 0 is one dimensional. The generalization to the case where the kernel is multidimensional will be reserved for future work. In this paper we focus on the generic case where the dynamics has only one stationary state. The mapping ρ 0 → ρ f given by Eq. (12) is exact and only requires obtaining L 1 , L 0 , and its kernel. Using the definition of L eff [see Eq. (7)], both operators L 0 = L eff (z = 0) and L 1 = dL eff (z) dz | z=0 can be written in terms of the original Lindblad operator L: Slow time non-trace-preserving evolution. We finally correct the evolution given by Eq. (10) by normalizing the state by α = tr[ρ f ] given by Eq. (12). We will drop the projector P in Pρ(t ) for clarity, given that we are only interested in the dynamics in the reduced subspace P, and we further assume that Qρ(0) = 0. Thus, the corrected evolution reads Equation (15), along with Eqs. (13) and (14) defining L 0 and L 1 , is one of the main results of the paper. It consists of a dynamical evolution rescaled by the population loss to the Q partition in the steady state. The difficult part in the calculation of L 0 and L 1 given by Eqs. (13) and (14) consists in the computation of the inverse of QLQ. As we will see in the next section, this inversion can be obtained explicitly only in specific cases. In general, a numerical inversion can be attempted but can be cumbersome, for instance when ran[Q] is an infinite-dimensional space. In that case, the inverse can be computed using perturbation theory. Indeed, QLQ can be written as QLQ = L D + W , where the matrix representation of L D is diagonal in the basis formed by the eigenvectors of PHP and QHQ, and W is nondiagonal. The inversion of QLQ can be written as As we show in Appendix B, in all cases where the relaxation processes inside the ran[Q] subspace can be neglected, W will depend only on the Hamiltonian couplings PHQ and QHP, and does not depend on the dissipative part. The fast dissipation of the QH part is involved in L D only. Therefore when the adiabatic elimination is a good approximation it is justified to consider that L D W . In most cases, retaining only the second-order terms (n = 2) at most, in the sum of Eq. (16), is enough to obtain a good approximation of the dynamics. Indeed, the level shift operator [L eff (z) − L 0 ] of Eq. (7) involves the operators PLQ and QLP which can each be first or zeroth order in the Hamiltonian coupling QHP or PHQ, so that only terms n = 0, 1, 2 for (QLQ) −1 are needed.
In the next section, we will illustrate in several examples how our result gives a very good approximation to the true dynamics.

III. EXAMPLES
We examine the evolution generated by the effective operator L 0 derived in the previous section with the correction given by Eq. (15), for a few specific cases when the excited states which are eliminated are (i) continuous manifolds and (ii) discrete states. We use continuous manifolds because they are part of fundamental toy models for both basic quantum evolution and spectroscopy, and also because they allow simplifications in the wideband approximation. In such an approximation, analytical expression of L eff (z) can be obtained. In general, using a continuous set of states in the wideband approximation, instead of a set of discrete levels, gives a zero real part of the level-shift operator (also called self-energy) leaving only the imaginary dissipative contribution. We then investigate systems with discrete excited states since they are more prevalent. We finally show that in the limit of large dissipation, the adiabatic evolutions where continuous and discrete excited-state manifolds are eliminated coincide. We only consider time-independent Hamiltonians; however the expressions derived in this work can describe the case where coherent radiation couples ground and excited states as long as the rotating wave approximation can be applied (which converts the problem to a time-independent one). The developed formalism is valid for arbitrary temperatures, which are captured in the Lindblad equation by the relative strength of the incoherent transitions from Q to P and vice versa. In the limit of zero temperature only incoherent transitions from Q to P are nonzero.

A. Elimination of continua excited states
Hamiltonians with a continuous spectrum have been part of the spectroscopist toolbox for several decades to describe atomic, molecular, and condensed matter systems [34][35][36][42][43][44][45][46][47][48][49][50][51]. Their distinctive property is that they result in an asymmetric line shape in the photodissociation and photofragmentation cross section as a function of the frequency of the impinging light, arising from interference processes [34]. The Hamiltonian structure as well as dissipative transitions are shown in Fig. 1. A set of N g ground states |g i are coupled among themselves by Hamiltonian couplings V i j (i, j = 1, . . . , N g ) as well as to N e continuous sets of excited states |k j ( j = 1, 2, . . . , N e ). Continua are not coupled among themselves (any coupling between continua can be removed by a unitary transformation which redefines all the other couplings); they are coupled to the ground states through Hamiltonian couplings V ( j) i and through dissipation at rates ( j) i (i = 1, 2, . . . , N g and j = 1, 2, . . . , N e ). In the following we adopt the wideband approximation where the couplings V ( j) i , the rates ( j) i , and the density of states n ( j) = dk j dE per unit of energy E are considered to be independent of k j . The continuum states are normalized so that k|k = δ(k − k ). The general problem with N g ground states coupled to N e excited states is considered in Appendix C, while in the following we examine in detail the case of one continuum coupled to either one or two ground states.
Single ground-state level coupled to a single continuum. We first consider a single discrete level coupled to a continuum of states via a Hamiltonian coupling V (1) 1 . The continuum can dissipate back to the ground state with a rate (1) 1 (Fig. 2). The Liouvillian for this system is Ground-state population of a single discrete system coupled to a continuum as a function of the rescaled time γ (1) (19)]. The dotted line corresponds to the limit (1) 1 → ∞ which in this case leaves the population unchanged in the discrete state, and the dashed line corresponds to the limit (1) 1 = 0, which is the limit of a discrete level unitarily coupled to a continuum of states that corresponds to a description for particle decay.
The effective operator L eff (z) to describe the ground-state dynamics after elimination of the continuous set of excited states can be obtained explicitly using Eq. (7) (see Appendix C): with F eff = γ (1) 1 |g 1 g 1 |, and where γ (1) 1 = 2n (1) π (V (1) 1 ) 2 ; it represents the injection rate from discrete to continuum due to the Hamiltonian coupling. The operator can be expanded in powers of z as L eff (z) ≈ zL 1 where L 1 = −F eff ⊗F eff (1) 1 and where the z-independent term L 0 is zero. This means that the approximate dynamics given by e L 0 t = 1 [see Eq. (15)] has no dynamics. The correction to the ground state is then We can readily solve the exact dynamics of the ground state in terms of the dimensionless constant β and the rescaled time τ = γ (1) 1 t: As expected, we see that the correction 1 − L 1 −1 gives the exact population lost β/(β + 1) for the steady state. We can already see from this simple example that this correction is nothing else than the detailed balance obtained from a kinetic equation between two sites P and Q in the steady state. Indeed, considering temporarily that P and Q are sites connected by classical rates, and taking n P and n Q to be the populations of the two sites and k Q→P , k P→Q the transition rates, we can writė which readily yield the steady-state population in P as It is expected that the site populations can be written in the form of a detailed balance equation since the Kolmogorov criterion is trivially satisfied. However, we remark that the relevant transition rate from Q to P is the dissipative relaxation rate (1) 1 while the relevant transition from P to Q is the Hamiltonian rate γ (1) 1 . This identification will be recovered in the more complicated case of a two-level system coupled to a continuum and then in a different form in the case of a system.
It is also illustrative to look at the exact solution given by Eq. (19) in the two limits of absent ( (1) 1 = 0) and very large dissipation ( (1) 1 γ (1) 1 ) from continuum to the ground state. As the dissipation rate (1) 1 goes to zero, we have a discrete level coupled to a continuum through Hamiltonian couplings only. This is the standard model for particle decay or injection into a band [52,53]. The evolution of the discrete state only can be fully described by a non-Hermitian Hamiltonian alone, entirely in Hilbert space, without the need for a Lindblad operator. In this case, the final state has zero population in the discrete ground state as all the population has been lost in the continuum. The opposite limit of infinitely high dissipation results in no dynamics whatsoever with the single discrete level being always populated. Because both cases are expressed in superoperator space as the limits of a continuous function of (1) 1 , we provide a rigorous connection between non-Hermitian Hamiltonian decay dynamics ( (1) 1 /γ (1) 1 → 0) and fully trace-preserving dissipative dynamics ( (1) 1 /γ (1) 1 → ∞) thanks to the nonlinear term of the form − z z+ . This connection is not restricted to the single discrete level system but is a general feature of discrete levels coupled to a manifold of continua where the evolution presents a transition from non-Hermitian decay Hamiltonians to trace-preserving generators, when the dissipation rate from the continuum is varied, and which could provide insight into comparisons of both approaches [54,55].
Two discrete states coupled to a single continuum. The model of a two-level system coupled to a continuous set of states is the standard Fano model invoked so often in spectroscopy [35]. Once more, the Liouvillian is written as and the quantum jump operators are Using Eq. (7) for the effective operator L eff (z) (see Appendix C), we obtain and = (1) 1 + (1) 2 . The effective Liouvillian can be expanded in powers of z as where and the correction coefficient is α = 1 − L 1 −1 as in Eq. (12). We calculate the time evolution with and without the correction α to the trace of the density matrix. In Fig. 3, we compare the exact evolution (solid line), the evolution with the effective Liouvillian ρ(t ) = e L 0 t ρ(0) (dash-dotted line), and the corrected evolution with ρ(t ) = αe L 0 t ρ(0) (dashed line). For each case, we show the expectation tr[ρ(t )σ k ] of the Pauli matrices σ k (k = x, y, z). The initial condition is ρ(0) = |g 1 g 1 |. For large values of the dissipation, all three evolutions coincide as expected since there is a negligible amount of population in the excited state. For small values of the dissipation, there is a fraction of the population that remains in the excited state so that evolution without the correction factor no longer asymptotically reaches the exact steady state. In addition to the evolution, we show the eigenvalues of L 0 and the nonlinear eigenvalues of L eff (z) [56]. We see that the first eigenvalues of L 0 are in good agreement with those of L eff (z). As a consequence of the Keldysh theorem [57][58][59] (see Appendix F), the nonlinear eigenvalues and eigenvectors of L eff (z) completely determine the timescales of the dynamics. In particular the gap of L eff (z), that is, the largest and nonzero real part of the nonlinear eigenvalues of L eff (z), determines the 4. Energy levels and transitions of a three-level system with dissipation. Hamiltonian couplings are indicated by straight arrows, dissipative processes by twisted arrows. Only decay from |e 1 to states |g i is represented, but the complete model includes also the reverse processes, that is, incoherent pumping. typical timescale to reach the stationary state. We see that the gap of L eff (z) is well reproduced by the gap of L 0 .
The correction factor α = 1 + i J i −1 can also be interpreted as arising from a detailed balance problem between P and Q. To make this more transparent, we recognize that L 0 = 0 (where the trace is taken over the steady-state density matrix) so that we may write α = 1 − N −1 where N = L 0 − i=1,2 J i is the non-Hermitian Hamiltonian superoperator that describes the decay of a two-level system into a continuum. Indeed, j . Therefore, the correction factor α can be interpreted again as the detailed balance factor arising from two sites P and Q equilibrating with rates k P→Q , corresponding to that of a non-Hermitian Hamiltonian decaying into a continuum, and k Q→P corresponding to a purely incoherent transition equal to the sum of decay rates from continuum to the discrete manifold. The dynamics themselves however cannot be described by a simple detailed balance problem.

B. Elimination of excited discrete states
The system. The system is one of the most used model systems in adiabatic elimination [23]. Its usefulness lies in that it sustains most of the useful features for applications in metrology, quantum computing, and thermometry, in particular in cold-ion traps [60][61][62][63][64][65][66][67][68][69][70][71][72]. The Liouvillian is Fig. 4) and the jump operators are The operators F j i and F j i represent incoherent channels going from the excited to the ground-state manifold, and from the ground-state manifold to the excited state, respectively.
The effective operator L eff (z) = PLP + PLQG 0 (z)QLP can be written in the perturbative limit up to order O((V 1 i ) 2 / ) for i = 1, 2, and in the zero-temperature limit, as where we have used the notation where The operator σ is defined as σ = |g 2 g 1 |.
After some algebra we get The form of the operators in the finite-temperature limit (with incoherent pumping from ground to excited state) is given in Appendix D. As in Fig. 3, in Figs. 5 and 6, we show the evolution of the expectation of the Pauli matrices as a function of time t for the same initial state. In Fig. 5 zero temperature is considered where only dissipation from excited to discrete states takes place ( 1 i = 0). On the contrary, in Fig. 6 the temperature is taken as infinite with equal rates for the dissipation from excited to ground and from ground to excited states ( 1 i = 1 i ). In the zero-temperature case, there is a negligible amount of population in the excited state [for the perturbative calculation of (QLQ) −1 to remain valid], and both the evolution with U (t ) = e L 0 t and that with U (t ) = αe L 0 t work well. We notice that the gap of L eff (z) in the infinite-temperature case is very well reproduced by the one of L 0 even though there is a significant population in the excited state (lower right panel of Fig. 6).
In the case of infinite temperature, the weak-field approximation is valid (so we can calculate the inverse of QLQ perturbatively) but there is a non-negligible population in the excited state. In this case, the correction α introduced in this article works very well in reproducing the final dynamics, while using a trace-preserving map does not. Writing the density matrix as a linear combination of Pauli matrices and the identity operator makes evident that the dynamics is well reproduced by L 0 (the Pauli matrices' evolutions with all operators are very close) as long as we use the correct normalization. Fig. 3 for a system in the zero-temperature limit. We show the trace of the density matrix and all three expectation values of the Pauli matrices. The parameters are E 1 = 7,

FIG. 5. Same as
All energy values are in units of V 1 1 , so that time is in units ofh/V 1 1 . The starkest difference between the exact evolution and the evolution with U (t ) = e L 0 t can be seen in the trace of the subsystem. As such this approximation sometimes fails to faithfully describe the population dynamics, which are recovered with the rescaled operator. The calculation of the resolvent (QLQ) −1 is done perturbatively in the field to the lowest order. A simulation with the exact resolvent is presented in Appendix H. The eigenvalues of L 0 are obtained by solving the linear eigenvalue problem (z − L 0 )v = 0 and in this case are restricted to four, while those of L eff (z) correspond to the nonlinear eigenvalue problem [z − L eff (z)]v = 0, which can have in principle up to 9 distinct eigenvalues.

IV. CONNECTION BETWEEN MODELS WITH ELIMINATION OF CONTINUOUS AND DISCRETE STATES
In this section we consider the connection between models where the states to be eliminated belong to a continuous set and models where theses states are discrete. Although Hamiltonians with continuous spectra represent a myriad of physical systems in their own right, they can also be viewed as useful ancillary mathematical structures that make the physics behind the more complicated Hamiltonians with a discrete spectrum more transparent. The reason for this is that Lamb shifts (or the conservative part of the level-shift operator) are absent in the case of a flat continuum (in the wideband approximation). The question we ask is, When does it matter whether we describe the excited states (which we would like to eliminate) as discrete states or as approximate continua?
Intuitively, both classes of models should coincide when the population of the excited states is negligible. We will show that this happens in two cases: (i) as the dissipation rate increases, the population of the excited state asymptotically vanishes, and (ii) at the points of coherence population FIG. 6. Same as Fig. 5 with the same parameters and in the infinite-temperature limit which opens incoherent transitions from ground to excited state with the same rate as the dissipation from the excited to the ground state. We show the trace of the density matrix and all three expectation values of the Pauli matrices. The starkest difference between the exact evolution and the evolution with U (t ) = e L 0 t can be seen in the trace of the subsystem. As such this approximation sometimes fails to faithfully describe the population dynamics, which are recovered with the rescaled operator. The calculation of the resolvent (QLQ) −1 is done perturbatively in the field to the lowest order. A simulation with the exact resolvent is presented in Appendix H. The eigenvalues of L 0 are obtained by solving the linear eigenvalue problem (z − L 0 )v = 0 and in this case are restricted to four, while those of L eff (z) correspond to the nonlinear eigenvalue problem [z − L eff (z)]v = 0, which can have in principle up to 9 distinct eigenvalues. trapping (CPT), the transition probability amplitudes to the excited state interfere destructively and the population of the excited state exactly vanishes [60][61][62][63][64]73,74].
Coincidence for large values of the dissipation. We calculate the limit of the effective operators L 0 and L 1 , as /ω j → ∞. For this we recast them in terms of the smallness parameters δ j = −iω j /( /2), for j ∈ {e 1 g 1 , g 1 e 1 , e 1 g 2 , g 2 e 1 }. As we take the limit of large dissipation lim →∞ δ j = 0, we get for the effective operators where the labels (3LS) and (cont) mean 3-level system and continuum models, respectively. We find that in the limit of large dissipation, the continuum and discrete effective operators for L 0 are the same as long as we set the density of states in the continuum model as n = 2/(π ), while they differ for L 1 . We can understand this convergence of operators as follows. The level-shift operator for the discrete excited states consists of a real part related to the dissipation and an imaginary part related to the Lamb shift. That of a flat continuum only has the real dissipative part. As the dissipation rate increases, the Lamb shift part of the operator becomes negligibly small and a discrete excited state becomes analogous to a continuum manifold as far as the evolution of the ground states is involved. In Fig. 7, we show the convergence of these models toward the exact solution of a system. For this we plot the steady-state population in the ground states and the steady-state fidelity as a function of dissipation rate from excited states to ground states. We rescale the fidelity F = Tr( √ ρ exact ρ √ ρ exact ) 2 /[Tr(ρ)Tr(ρ exact )] by a factor 1 − [Tr(ρ) − Tr(ρ exact )] 2 which penalizes evolution operators that do not have the correct asymptotic trace. We clearly see that the rescaled steady state for a three-level system performs best, and that the rescaled steady state for an equivalent continuum and the unscaled steady state for the three-level system approach the correct solution for similar values of the dissipation rate.
Coincidence at the coherent population trapping points. It can be shown that as long as we do not have dissipation within the ground-state manifold, or from the ground-state manifold to the excited state, that there will be points of coherent population trapping as long as the steady state ρ fulfills the following conditions [74]: FIG. 9. Steady-state population in the ground-state manifold ρ g 1 g 1 and ρ g 2 g 2 , and in the excited-state manifold ρ e 1 e 1 , with a threelevel system, a Fano model, and an effective Liouvillian evolution with L 0 for the three-level system. Parameters are V (1) 1 = 1.7, V (1) 2 = 1.0. The detuning between g 2 and e 1 is set to zero while the detuning between g 1 and e 1 is varied. All energies are given in units of V (1) 2 .
Remarkably, this condition is independent of the value of the dissipation rate (as long as ρ is a dark state), so that we are free to choose an arbitrarily large value and still retain the property of CPT where the population is restricted to the ground-state manifold. Accordingly, it follows from the previous paragraph that if we scale n (1) = 2/(π ) then the effective operators will be the same. It also follows that since α = 1, then L 1 = 0.
We illustrate the effect of coherence population trapping points on our models in Figs. 8 and 9. By plotting the groundstate population and fidelity of the three models at the CPT condition we see that all four models coincide (Fig. 8). To further stress the equivalence of the models around CPT, we plot the steady state of a three-level system, of a Fano model, and of the effective Liouvillian L 0 for a three-level system, as a function of the detuning of the ground states E g 1 − E g 2 (Fig. 9). We observe the CPT point at zero detuning where all three models coincide. The unscaled L 0 only agrees at the CPT condition since it preserves the population in the ground-state manifold while both the continuum and the exact system agree around a neighborhood of the CPT point.
We have shown that replacing discrete excited states by continua corresponds to taking the limit of large dissipation, or alternatively finding the CPT points. This is important since the effective operator with a continuum is much more straightforward to calculate exactly than that of a discrete level. Thus calculations that fulfill these conditions, if carried out using these simplified operators, can be more easily solved analytically.

V. CONCLUSION
We have derived expressions for the adiabatic elimination of a fast manifold in frequency space. This has allowed us to correct for particle density loss to the fast manifold and rescale the evolution operator. We have illustrated this with examples spanning discrete and excited state continua which show the advantages of the correction factor as well as its physical meaning. We have provided an equivalence between the discrete and continuum models at the CPT condition and in the limit of large dissipation, giving insight into commonly used adiabatic elimination approaches. For an N-level system, the underlying Hilbert space H is of dimension N and the states of the quantum systems are described by positive operators acting on H that can be represented by N × N density matrices. The superoperators such as the Lindblad operator L or its resolvent G are linear operators acting on operators themselves acting on H.
To describe an open quantum system we need to know the evolution of the density matrix using the Lindblad equation, which in Hilbert H space is written aṡ A disadvantage of this form is that neither the exponential map nor the resolvent can be straightforwardly expressed or calculated numerically. It is therefore convenient to represent the density matrix as a vector ρ with N 2 components, obtained from the column-stretched form of the N × N density matrix. This representation is obtained by considering the density matrix as an element ρ of the Hilbert space H ⊗ H [40].
In that way, superoperators are linear operators acting on H ⊗ H and they can be represented by N 2 × N 2 matrices. The linear superoperator acting on ρ, built from 2 arbitrary operators S n and S m on H and acting on ρ as S n ρS † m , is given by the mapping S n ρS † m → (S m ⊗ S n ρ ). With the help of this mapping, the Lindblad operator operating on the vector form of the density matrix as d dt

APPENDIX B: PERTURBATIVE INVERSION OF QLQ
We consider a generic system, with a Hamiltonian H written as H = H 0 + V , where H 0 = PHP + QHQ and V = PHQ + QHP. Let |i; p (| j, q ) be the eigenstates of PHP (QHQ), with i = 1, 2, . . . , N p ( j = 1, 2, . . . , N q ). For the dissipation processes, we consider relaxation from the fast subspace ran[Q] to the slow subspace ran[P], described by jump operators F i j = i j |i; p j; q|, relaxation from the slow subspace to the fast subspace described by jump operators J ji = √ γ i j | j; q i; p|, and finally we also consider relaxation inside ran[P], described by jump operators N m which we do not specify as they do not intervene in QLQ. We neglect all the dissipation processes between states belonging to ran [Q].
It is convenient to define a non-Hermitian Hamiltonian has a diagonal matrix representation in the basis {|i; p , | j, q }.
We can then rewrite the Lindblad operator L as [see Eq. (A3)]

042102-10
Using the expression of Q given by Eq.
(2), we notice that Therefore QLQ can be written as where L D = −ıQ(1 ⊗ K 0 −K 0 ⊗ 1)Q has a diagonal matrix representation in the basis {|i; p ⊗ | j, q } and W = −ıQ(1 ⊗ V − V ⊗ 1)Q has a nondiagonal matrix representation in the same basis. The nondiagonal part of QLQ depends only on the Hamiltonian coupling V which can be considered as a small perturbation with respect to the diagonal part when the relaxation of the fast space ran[Q] is fast ( i j V i j ).

APPENDIX C: GENERAL CASE OF N g GROUND STATES COUPLED TO N e CONTINUA
We provide here the general expressions to calculate the effective Liouvillian for N g discrete ground states coupled to N e continua, from which the more specific examples detailed in the main text can be derived. The complete Liouvillian for such a system is i * |k a i| , With the dissipative part the Liouvillian is The effective operators are obtained in a calculation similar to what we have done previously [38] but keeping the z dependence of the operators. Briefly, we define the projection operators for the continuous part (Q) and the discrete part (P). The effective Liouvillian [see Eq. (7)] hinges on the resolvent operator in Q. This operator can be expanded in a Lippman-Schwinger series that is exactly resummed for the wideband approximation, where the parameters of the continuum do not depend on the continuum energy. We obtain In order for L 0 to be a generator of a trace-preserving map, the maximally mixed state [ρ] = 1 N 1 H must be a left eigenvector for L 0 with eigenvalue 0, where N is the dimension of H. In vector form, we can associate with ρ the maximally entangled state |1 H ∈ H ⊗ H. Therefore, the trace-preserving condition can be written as where we have used the fact that L is a Lindblad operator, hence generating a trace-preserving dynamics. If we define the set X as then a sufficient condition for L 0 to generate a trace-preserving dynamics is that |1 H is orthogonal to the set X . Equivalently, since X ⊆ ran[Q], we can write this condition as where we have defined the state | Q = i q |i q ⊗ |i q , with the index i q enumerating the left eigenvectors of Q corresponding to eigenvalue 1.
In all the examples considered in this article, we had ran[QLP] ⊆ ran[QLQ] = ran [Q] implying that X is an empty set and, therefore, the fulfillment of the above condition.

APPENDIX F: KELDYSH THEOREM
For the sake of completeness, we recall here the Keldysh theorem. We consider only the case where the nonlinear eigenvalues are simple. This section is based on the material of Ref. [59]. To connect our notation with the usual statement of the theorem, we define T (z) such that T (z) = z1 − L eff (z); therefore PG(z)P = [T (z)] −1 .
First we recall the definition of a nonlinear eigenvalue λ of T (z): λ is an eigenvalue of T (z) if T (λ)v = 0 for some nonzero vector v. The vector v is the right eigenvector of T . The eigenvalue is called simple if in addition ker[T (λ)] = span{v}, v = 0, T (λ)v / ∈ ran[T (λ)].
In this case the adjoint T † of T satisfies ker[T † (λ)] = span{w} for some nonzero vector w, and furthermore, w † T (λ)v = 0. Without loss of generality we can choose where T (λ) is the value of the derivative of T (z) with respect to z, taken at z = λ. The Keldysh theorem states the following: Let D be a compact subset that contains only simple eigenvalues λ n , n = 1, . . . , N, with right and left eigenvectors v n and w n , respectively; then there is a neighborhood U of D and a holomorphic function R(z) such that Now, if we assume that all the eigenvalues of T are simple, then we can use Eq. (F2) to calculate Pρ(t ), performing the integration of Eq. (6), and we obtain where we recall that we use the normalization w † n T (λ)v n = 1. In particular for λ = 0, we have T (0) = 1 − L 1 ; thus w † 0 (1 − L 1 )ρ = 1, where ρ = v 0 is the kernel of L(z = 0) = L 0 .