# Adiabatic elimination and subspace evolution of open quantum systems

Daniel Finkelstein-Shapiro, David Viennot, Ibrahim Saideh, Thorsten Hansen,

Tönu Pullerits, Arne Keller

► To cite this version:

Daniel Finkelstein-Shapiro, David Viennot, Ibrahim Saideh, Thorsten Hansen, Tönu Pullerits, et al. Adiabatic elimination and subspace evolution of open quantum systems. Physical Review A, 2020, 101 (4), 10.1103/PhysRevA.101.042102. [https://arxiv.org/abs/1908.09476](https://arxiv.org/abs/1908.09476)

HAL Id: hal-02564106

https://hal.science/hal-02564106

Submitted on 29 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Adiabatic elimination and subspace evolution of open quantum systems

Daniel Finkelstein-Shapiro<sup>1,*</sup> David Viennot,<sup>2</sup> Ibrahim Saideh<sup>3,4</sup> Thorsten Hansen<sup>5</sup>, Tönu Pullerits,<sup>1</sup> and Arne Keller<sup>3,4,†</sup>

<sup>1</sup>Division of Chemical Physics, Lund University, Box 124, 221 00 Lund, Sweden

<sup>2</sup>Institut UTINAM, CNRS UMR 6213, Université de Bourgogne-Franche-Comté, Observatoire de Besançon, 25010 Besançon, France

<sup>3</sup>Laboratoire Matériaux et Phénomènes Quantiques, Université Paris Diderot, CNRS UMR 7162, 75013 Paris, France

<sup>4</sup>Université Paris-Saclay, 91405 Orsay, France

<sup>5</sup>Department of Chemistry, University of Copenhagen, DK 2100 Copenhagen, Denmark

<figure><img src="image_2.png" /></figure>

(Received 28 September 2019; accepted 31 January 2020; published 1 April 2020)

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, one 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.

DOI: 10.1103/PhysRevA.101.042102

## 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–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.

<footnote>*daniel.finkelstein_shapiro@chemphys.lu.se</footnote>

<footnote>†arne.keller@u-psud.fr</footnote>
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–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–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 time-dependent 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–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\rho(t)$, the slow component of the quantum state $\rho(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_{\text{eff}}(z)$, a z-dependent operator defined on the slow subspace only, such that $PG(z)P = [z - L_{\text{eff}}(z)]^{-1}$.

The operator $L_0 = L_{\text{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 $\mathcal{H}$ is built from a projector $P$ onto the underlying Hilbert space $\mathcal{H}$ as $P\rho = P\rho P$, as in Refs. [14,23]. In others words, we assume that the fast/slow partition is linked to a partition of $\mathcal{H}$ in two complementary subspaces $\mathcal{H} = PH \oplus QH$. The application of our formalism to bipartite systems where the fast/slow partition is linked to a tensorial structure $\mathcal{H} = \mathcal{H}_{\text{slow}} \otimes \mathcal{H}_{\text{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–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_{\text{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 $\mathcal{H}$ of the system is partitioned into two subspaces with the help of two orthogonal projectors $P$ and $Q = 1 - P$, where $1_{\mathcal{H}}$ is the identity operator on $\mathcal{H}$. The $Q\mathcal{H}$ 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 $Q\mathcal{H}$ has reached its stationary state.

We suppose that the system is coupled to a bath that opens dissipation channels between $Q\mathcal{H}$ and $PH$, or within $PH$ and $Q\mathcal{H}$. Hamiltonian couplings ($PHQ$, or $QHP$) can also open transitions between $PH$ and $Q\mathcal{H}$. Associated with $P$ and $Q$, we define superprojector operators $P$ and $Q$ such that

$$P\rho = P\rho P, \quad Q = 1 - P, \qquad (1)$$

where $1$ is the identity superoperator on the space of operator on $\mathcal{H}$, and $\rho$ is an operator on $\mathcal{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\rangle\langle b|$ in the Hilbert space $\mathcal{H}$ onto the vector $|\bar{b}\rangle \otimes |a\rangle$ in the $\mathcal{H} \otimes \mathcal{H}$ Hilbert space, or equivalently maps any $n \times n$ density matrix $\rho$ to a column vector $\bar{\rho}$ with $n^2$ elements, by stacking the columns of the $\rho$ matrix. Under this isomorphism, the operation $A\rho B^\dagger$ is mapped to $\bar{B} \otimes A\bar{\rho}$, where A and B are operators on $\mathcal{H}$ and $\bar{B}$ denotes the complex conjugate of B; that is, $\bar{B} = (B^\dagger)^T$, where $B^\dagger$ is the adjoint and $B^T$ is the transpose of B (see Appendix A). From now on, we drop the arrow in $\bar{\rho}$ as we assume that $\rho$ is in vector form. The only exception is when a density matrix $\rho$ is inside a bracket like in tr[$\rho$].

With this notation, the superprojector operators $\mathcal{P}$ and $\mathcal{Q}$ read

$$ \mathcal{P} = \mathcal{P} \otimes \mathcal{P}, \quad \mathcal{Q} = \mathcal{Q} \otimes \mathcal{Q} + \mathcal{P} \otimes \mathcal{Q} + \mathcal{Q} \otimes \mathcal{P}. \quad (2) $$

Also, the general form of the Lindblad operator L, generator of the evolution, $\dot{\rho} = L\rho$, can be written as [41]

$$ L = -i[\mathbb{1} \otimes H - \bar{H} \otimes \mathbb{1}] + \sum_i \mathcal{D}(F_i), $$

where

$$ \mathcal{D}(F) = \bar{F} \otimes F - \frac{1}{2}[\mathbb{1} \otimes F^\dagger F + (F^\dagger F)^T \otimes \mathbb{1}]. \quad (3) $$

We start by expressing the density matrix evolution in an integral form through the Laplace transform:

$$ \rho(t) = \frac{1}{2\pi i} \int_D e^{zt} G(z) \rho_0 dz, \quad (4) $$

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 \in \mathbb{C}; \text{Re } z = a > 0\}$. Projecting Eq. (4) using $\mathcal{P}$ and $\mathcal{Q}$ gives

$$ \begin{aligned} \mathcal{P}\rho(t) &= \frac{1}{2\pi i} \int_D dz e^{zt} [\mathcal{P}G(z)\mathcal{P}\rho(0) + \mathcal{P}G(z)\mathcal{Q}\rho(0)], \\ \mathcal{Q}\rho(t) &= \frac{1}{2\pi i} \int_D dz e^{zt} [\mathcal{Q}G(z)\mathcal{P}\rho(0) + \mathcal{Q}G(z)\mathcal{Q}\rho(0)]. \end{aligned} \quad (5) $$

In the remainder of the text, we make the assumption that at time $t = 0$ the population is entirely in the slow subspace $\mathcal{PH}$ so that $\mathcal{Q}\rho(0) = 0$. Hence the evolution in the $\mathcal{PH}$ subspace is simply given by

$$ \mathcal{P}\rho(t) = \frac{1}{2\pi i} \int dz e^{zt} \mathcal{P}G(z)\mathcal{P}\rho(0). \quad (6) $$

We define the operator $L_{\text{eff}}(z)$, a z-dependent operator defined on $\mathcal{PH}$, such that $\mathcal{P}G(z)\mathcal{P} = [z - L_{\text{eff}}(z)]^{-1}$. Using the definition of the resolvent and the orthogonality of the $\mathcal{P}$ and $\mathcal{Q}$ projectors, we have

$$ L_{\text{eff}}(z) = \mathcal{P}\mathcal{L}\mathcal{P} + \mathcal{P}\mathcal{L}\mathcal{Q}G_0(z)\mathcal{Q}\mathcal{L}\mathcal{P}, \quad (7) $$

where $\mathcal{Q}G_0(z)\mathcal{Q} = [z - \mathcal{Q}\mathcal{L}\mathcal{Q}]^{-1}$ is the resolvent of $\mathcal{QLQ}$. Equation (6) with Eq. (7) is an exact description of the dynamics (restricted to $\mathcal{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 $\mathcal{PH}$ and $\mathcal{QH}$ partitions can exchange population during the evolution.

**Generator of the slow dynamics.** We notice that $L_0 = L_{\text{eff}}(z = 0)$ is the generator of the slow time dynamics. Indeed, projecting the Lindblad equation $\dot{\rho} = L\rho$ on $\mathcal{PH}$ and $\mathcal{QH}$ we have

$$ \mathcal{P}\dot{\rho}(t) = \mathcal{P}\mathcal{L}\mathcal{P}\rho(t) + \mathcal{P}\mathcal{L}\mathcal{Q}\rho(t), \quad (8) $$

$$ \mathcal{Q}\dot{\rho}(t) = \mathcal{Q}\mathcal{L}\mathcal{Q}\rho(t) + \mathcal{Q}\mathcal{L}\mathcal{P}\rho(t). \quad (9) $$

To obtain the approximate slow time dynamics in the subspace $\mathcal{PH}$, we assume that $\mathcal{Q}\rho$ has reached a stationary regime, $\mathcal{Q}\dot{\rho} = 0$. Using Eq. (9) to express $\mathcal{Q}\rho$ as a function of $\mathcal{P}\rho$, and inserting the result in Eq. (8), we obtain

$$ \mathcal{Q}\dot{\rho} = 0 \Rightarrow \mathcal{P}\dot{\rho} = L_0\mathcal{P}\rho. \quad (10) $$

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 $\mathcal{P}\rho$ and $\mathcal{Q}\rho$. We must then correct this evolution to take into account the possible variation of the trace of $\mathcal{P}\rho$. To this end, we look for the exact final state, reached in $\mathcal{PH}$ subspace, $\rho_f = \lim_{t \to \infty} \mathcal{P}\rho(t)$, from a given initial state $\rho_0 = \rho(t = 0)$.

**Mapping to the final state.** By Eq. (10), we know that the final state $\rho_f$, in $\mathcal{PH}$ subspace, is in the kernel of $L_0 = L_{\text{eff}}(z = 0)$. We assume that the kernel is one dimensional and define $\bar{\rho}$ its unique element with tr[$\bar{\rho}$] = 1. Then $\rho_f = \alpha\bar{\rho}$, and we are left to determine $\alpha = \text{tr}[\rho_f]$. The final stationary state $\rho_f$ can be obtained taking the limit of Eq. (6) when $t \to \infty$. This limit can be obtained using the final value theorem:

$$ \rho_f = \lim_{z \to 0} z \mathcal{P}G(z)\mathcal{P}\rho(0) = \lim_{z \to 0} z[z - L_{\text{eff}}(z)]^{-1}\rho(0). \quad (11) $$

As we show in Appendix G, this limit can be calculated explicitly as

$$ \rho_f = \alpha \bar{\rho}, \quad \text{with } \alpha = \frac{1}{\text{tr}[(\mathbb{1} - L_1)\bar{\rho}]} = \frac{1}{1 - \langle L_1 \rangle}, \quad (12) $$

where $L_1 = \frac{dL_{\text{eff}}(z)}{dz}|_{z=0}$. We notice that $\alpha$ given by Eq. (12) does not depend on the initial state $\rho_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 $\rho_0 \to \rho_f$ given by Eq. (12) is exact and only requires obtaining $L_1$, $L_0$, and its kernel. Using the definition of $L_{\text{eff}}$ [see Eq. (7)], both operators $L_0 = L_{\text{eff}}(z = 0)$ and $L_1 = \frac{dL_{\text{eff}}(z)}{dz}|_{z=0}$ can be written in terms of the original Lindblad operator L:

$$ L_0 = \mathcal{P}\mathcal{L}\mathcal{P} - \mathcal{P}\mathcal{L}\mathcal{Q}(\mathcal{Q}\mathcal{L}\mathcal{Q})^{-1}\mathcal{Q}\mathcal{L}\mathcal{P}, \quad (13) $$

$$ L_1 = -\mathcal{P}\mathcal{L}\mathcal{Q}(\mathcal{Q}\mathcal{L}\mathcal{Q})^{-2}\mathcal{Q}\mathcal{L}\mathcal{P}. \quad (14) $$

**Slow time non-trace-preserving evolution.** We finally correct the evolution given by Eq. (10) by normalizing the state by $\alpha = \text{tr}[\rho_f]$ given by Eq. (12). We will drop the projector $\mathcal{P}$ in $\mathcal{P}\rho(t)$ for clarity, given that we are only interested in the
dynamics in the reduced subspace P, and we further assume that $Q\rho(0) = 0$. Thus, the corrected evolution reads

$$ \rho(t) = \frac{1}{1 - \langle L_1 \rangle} e^{L_0 t} \rho_0. \qquad (15) $$

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

$$ (QLQ)^{-1} = L_D^{-1} \sum_{n=0}^{\infty} (WL_D^{-1})^n. \qquad (16) $$

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 \gg 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_{\text{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_{\text{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

<figure><img src="image_3.png" /><figcaption>FIG. 1. Energy levels and transitions of a Fano-type model with dissipation. Hamiltonian couplings are indicated by straight arrows, dissipative processes by twisted arrows.</figcaption></figure>

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–36,42–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\rangle$ are coupled among themselves by Hamiltonian couplings $V_{ij}$ ($i, j = 1, \dots, N_g$) as well as to $N_e$ continuous sets of excited states $|k_j\rangle$ ($j = 1, 2, \dots, 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_i^{(j)}$ and through dissipation at rates $\Gamma_i^{(j)}$ ($i = 1, 2, \dots, N_g$ and $j = 1, 2, \dots, N_e$). In the following we adopt the wideband approximation where the couplings $V_i^{(j)}$, the rates $\Gamma_i^{(j)}$, and the density of states $n^{(j)} = \frac{dk_j}{dE}$ per unit of energy $E$ are considered to be independent of $k_j$. The continuum states are normalized so that $\langle k|k' \rangle = \delta(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 $\Gamma_1^{(1)}$ (Fig. 2). The Liouvillian for this system is $L = -i(1 \otimes H - \bar{H} \otimes 1) +$
<figure><img src="image_4.png" /><figcaption>FIG. 2. Ground-state population of a single discrete system coupled to a continuum as a function of the rescaled time $\gamma_1^{(1)}t$ with $\gamma_1^{(1)} = n\pi(V_1^{(1)})^2$, for different values of the dimensionless constant $\beta = \Gamma_1^{(1)}/\gamma_1^{(1)}$ [given by Eq. (19)]. The dotted line corresponds to the limit $\Gamma_1^{(1)} \to \infty$ which in this case leaves the population unchanged in the discrete state, and the dashed line corresponds to the limit $\Gamma_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.</figcaption></figure>

$\mathcal{D}(F_1^{(1)})$ [see Eq. (3)], where

$$H = \int dk_1 V_1^{(1)} |g_1\rangle\langle k_1| + \text{c.c.}, \\ F_1^{(1)} = \sqrt{\Gamma_1^{(1)}} \int dk_1 |g_1\rangle\langle k_1|. \qquad (17)$$

The effective operator $L_{\text{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):

$$L_{\text{eff}}(z) = \left( \frac{-z}{z + \Gamma_1^{(1)}} \right) \bar{F}_{\text{eff}} \otimes F_{\text{eff}}, \qquad (18)$$

with $F_{\text{eff}} = \sqrt{\gamma_1^{(1)}}|g_1\rangle\langle g_1|$, and where $\gamma_1^{(1)} = 2n^{(1)}\pi(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_{\text{eff}}(z) \approx zL_1$ where $L_1 = -\frac{\bar{F}_{\text{eff}}\otimes F_{\text{eff}}}{\Gamma_1^{(1)}}$ and where the $z$-independent term $L_0$ is zero. This means that the approximate dynamics given by $e^{L_0t} = 1$ [see Eq. (15)]

has no dynamics. The correction to the ground state is then $(1 - L_1)^{-1} = \frac{1}{1 + \beta^{-1}}$ where $\beta = \Gamma_1^{(1)}/\gamma_1^{(1)}$.

We can readily solve the exact dynamics of the ground state in terms of the dimensionless constant $\beta$ and the rescaled time $\tau = \gamma_1^{(1)}t$:

$$\rho(\tau) = \frac{\beta + e^{-(\beta+1)\tau}}{\beta + 1} |g_1\rangle\langle g_1|. \qquad (19)$$

As expected, we see that the correction $(1 - L_1)^{-1}$ gives the exact population lost $\beta/(\beta + 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 \to P}$, $k_{P \to Q}$ the transition rates, we can write

$$\dot{n}_P = k_{Q \to P} n_Q - k_{P \to Q} n_P, \\ \dot{n}_Q = -k_{Q \to P} n_Q + k_{P \to Q} n_P, \qquad (20)$$

which readily yield the steady-state population in $P$ as $n_P = \frac{1}{1 + \frac{k_{P \to Q}}{k_{Q \to P}}}$ . We thus identify $k_{Q \to P} \equiv \Gamma_1^{(1)}$ and $k_{P \to Q} \equiv \gamma_1^{(1)} = n^{(1)}\pi(V_1^{(1)})^2$. 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 $\Gamma_1^{(1)}$ while the relevant transition from $P$ to $Q$ is the Hamiltonian rate $\gamma_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 $\Lambda$ system.

It is also illustrative to look at the exact solution given by Eq. (19) in the two limits of absent $(\Gamma_1^{(1)} = 0)$ and very large dissipation $(\Gamma_1^{(1)} \gg \gamma_1^{(1)})$ from continuum to the ground state. As the dissipation rate $\Gamma_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 $\Gamma_1^{(1)}$, we provide a rigorous connection between non-Hermitian Hamiltonian decay dynamics $(\Gamma_1^{(1)}/\gamma_1^{(1)} \to 0)$ and fully trace-preserving dissipative dynamics $(\Gamma_1^{(1)}/\gamma_1^{(1)} \to \infty)$ thanks to the nonlinear term of the form $-\frac{z}{z+\Gamma}$. 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 $L = -i(1 \otimes H - \tilde{H} \otimes 1) + \sum_{i=1,2} \mathcal{D}(F_i^{(1)})$, where the Hamiltonian $H = H_0 + H_V$ is

$$ 
\begin{aligned}
H_0 &= E_1|g_1\rangle\langle g_1| + E_2|g_2\rangle\langle g_2| + \int dk \epsilon_{k_1} |k_1\rangle\langle k_1|, \\
H_V &= V_{12}|g_1\rangle\langle g_2| + V_{12}^*|g_2\rangle\langle g_1| \\
&\quad + \int dk_1 [V_1^{(1)}|g_1\rangle\langle k_1| + V_1^{(1)*}|k_1\rangle\langle g_1|] \\
&\quad + \int dk_1 [V_2^{(1)}|g_2\rangle\langle k_1| + V_2^{(1)*}|k_1\rangle\langle g_2|],
\end{aligned} 
\quad (21) $$

and the quantum jump operators are

$$ 
\begin{aligned}
F_1^{(1)} &= \int dk_1 \sqrt{\Gamma_1^{(1)}} |g_1\rangle\langle k_1|, \\
F_2^{(1)} &= \int dk_1 \sqrt{\Gamma_2^{(1)}} |g_2\rangle\langle k_1|.
\end{aligned} 
\quad (22) $$

Using Eq. (7) for the effective operator $L_{\text{eff}}(z)$ (see Appendix C), we obtain

$$ L_{\text{eff}}(z) = L_0 + \Delta(z) \sum_{i=1,2} J_i, \quad (23) $$

where $L_0 = -i(1 \otimes \text{PHP} - \tilde{P}\text{HP} \otimes 1) + \sum_{i=1,2} \mathcal{D}(F_{\text{eff},i})$ and

$$ 
\begin{aligned}
F_{\text{eff},i} &= \sum_{j=1,2} \sqrt{\frac{\Gamma_i^{(1)}}{\Gamma}} V_j^{(1)} |g_i\rangle\langle g_j|, \\
J_i &= \tilde{F}_{\text{eff},i}^{(1)} \otimes F_{\text{eff},i}^{(1)}, \\
\Delta(z) &= -\frac{z}{(z+\Gamma)},
\end{aligned} 
\quad (24) $$

and $\Gamma = \Gamma_1^{(1)} + \Gamma_2^{(1)}$. The effective Liouvillian can be expanded in powers of $z$ as

$$ L_{\text{eff}}(z) = L_0 + zL_1 + \dots, \quad (25) $$

where

$$ L_1 = - \sum_{i=1,2} \frac{J_i}{\Gamma}, \quad (26) $$

and the correction coefficient is $\alpha = \langle 1 - L_1 \rangle^{-1}$ as in Eq. (12).

We calculate the time evolution with and without the correction $\alpha$ to the trace of the density matrix. In Fig. 3, we compare the exact evolution (solid line), the evolution with the effective Liouvillian $\rho(t) = e^{L_0 t} \rho(0)$ (dash-dotted line), and the corrected evolution with $\rho(t) = \alpha e^{L_0 t} \rho(0)$ (dashed line). For each case, we show the expectation tr[$\rho(t)\sigma_k$] of the Pauli matrices $\sigma_k$ ($k = x, y, z$). The initial condition is $\rho(0) = |g_1\rangle\langle 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.

<figure><img src="image_5.png" /></figure>

$\Gamma_1^{(1)} = \Gamma_2^{(2)} = 0.01$

$\Gamma_1^{(1)} = \Gamma_2^{(2)} = 10.0$

<figure><img src="image_6.png" /><figcaption>FIG. 3. Fano model in the zero-temperature limit. The energies are given in units of $V_1^{(1)}$ and times in units of $1/V_1^{(1)}$. Left column: Small values of the dissipation rates. Right column: High values of the dissipation rates. First row: Evolution of the expectation value of the Pauli matrices $\sigma_i$ ($i = x, y, z$) as indicated in the inset. Dash-dotted line: Trace-preserving evolution with $e^{L_0 t}$. Dashed line: Trace rescaled evolution $\alpha e^{L_0 t}$. Second row, red circle: Nonlinear eigenvalues of $L_{\text{eff}}(z)$; black cross: linear eigenvalues of $L_0$. Parameters for the simulations are $E_1 = 0.0$, $E_2 = 0.9$, $V_{12} = 0.0$, $V_1^{(1)} = 1.0$, $V_2^{(1)} = 0.2$, $n = 1$. All energy values are in units of $V_1^{(1)}$, so that time is in units of $\hbar/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 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_{\text{eff}}(z)$ correspond to the nonlinear eigenvalue problem [$z - L_{\text{eff}}(z)]v = 0$, which can have in principle up to 6 distinct eigenvalues.</figcaption></figure>

In addition to the evolution, we show the eigenvalues of $L_0$ and the nonlinear eigenvalues of $L_{\text{eff}}(z)$ [56]. We see that the first eigenvalues of $L_0$ are in good agreement with those of $L_{\text{eff}}(z)$. As a consequence of the Keldysh theorem [57–59] (see Appendix F), the nonlinear eigenvalues and eigenvectors of $L_{\text{eff}}(z)$ completely determine the timescales of the dynamics. In particular the gap of $L_{\text{eff}}(z)$, that is, the largest and nonzero real part of the nonlinear eigenvalues of $L_{\text{eff}}(z)$, determines the
<figure><img src="image_7.png" /><figcaption>FIG. 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\rangle$ to states $|g_i\rangle$ is represented, but the complete model includes also the reverse processes, that is, incoherent pumping.</figcaption></figure>

typical timescale to reach the stationary state. We see that the gap of $L_{\text{eff}}(z)$ is well reproduced by the gap of $L_0$.

The correction factor $\alpha = \langle 1 + \sum_i \frac{J_i}{\Gamma} \rangle^{-1}$ can also be interpreted as arising from a detailed balance problem between P and Q. To make this more transparent, we recognize that $\langle L_0 \rangle = 0$ (where the trace is taken over the steady-state density matrix) so that we may write $\alpha = \langle 1 - \frac{N}{\Gamma} \rangle^{-1}$ where $N = L_0 - \sum_{i=1,2} J_i$ is the non-Hermitian Hamiltonian superoperator that describes the decay of a two-level system into a continuum. Indeed, $N = -i(1 \otimes H_D - \hat{H}_D \otimes 1)$ with $H_D = \text{PHP} - i \sum_{j=1,2}[F_j^{(1)}]^\dagger F_j^{(1)}$. Therefore, the correction factor $\alpha$ can be interpreted again as the detailed balance factor arising from two sites P and Q equilibrating with rates $k_{P \to Q}$, corresponding to that of a non-Hermitian Hamiltonian decaying into a continuum, and $k_{Q \to 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 $\Lambda$ system. The $\Lambda$ 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–72]. The Liouvillian is $L = -i(1 \otimes H - \hat{H} \otimes 1) + \sum \mathcal{D}(F_i)$, where (see Fig. 4)

$$ 
\begin{aligned}
H_0 &= E_1|g_1\rangle\langle g_1| + E_2|g_2\rangle\langle g_2| + E_3|e_1\rangle\langle e_1|, \\
H_V &= V_{12}|g_1\rangle\langle g_2| + V_{12}^*|g_2\rangle\langle g_1| + V_1^*|g_1\rangle\langle e_1| + V_1^*|e_1\rangle\langle g_1| \\
&\quad + V_2^*|g_2\rangle\langle e_1| + V_2^*|e_1\rangle\langle g_2|,
\end{aligned} 
\quad (27) 
$$

and the jump operators are

$$ 
\begin{aligned}
F_1^1 &= \sqrt{\Gamma_1^1}|g_1\rangle\langle e_1|, \\
F_2^1 &= \sqrt{\Gamma_2^1}|g_2\rangle\langle e_1|, \\
F_1^{1'} &= \sqrt{\Gamma_1^{1'}}|e_1\rangle\langle g_1|, \\
F_2^{1'} &= \sqrt{\Gamma_2^{1'}}|e_1\rangle\langle g_2|.
\end{aligned} 
\quad (28) 
$$

The operators $F_i^j$ and $F_i^{j'}$ 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_{\text{eff}}(z) = \text{PLP} + \text{PLQG}_0(z)\text{QLP}$ can be written in the perturbative limit up to order $O((V_i^1)^2/\Gamma)$ for $i = 1, 2$, and in the zero-temperature limit, as

$$ L_{\text{eff}}(z) = \sum_{j=1}^{4} \left(1 + \frac{M\Gamma}{z + \Gamma}\right) \frac{h_j}{z - \xi_j}, \quad (29) $$

where we have used the notation

$$ 
\begin{aligned}
h_1 &= \sigma\sigma^\dagger \otimes S_0^\dagger S_0, & \xi_1 &= -i\omega_{e_1 g_1} - \Gamma/2, \\
h_2 &= S_0^\dagger S_0 \otimes \sigma\sigma^\dagger, & \xi_2 &= i\omega_{e_1 g_1} - \Gamma/2, \\
h_3 &= \sigma^\dagger\sigma \otimes S_0^\dagger S_0, & \xi_3 &= -i\omega_{e_1 g_2} - \Gamma/2, \\
h_4 &= S_0^\dagger S_0 \otimes \sigma^\dagger\sigma, & \xi_4 &= i\omega_{e_1 g_2} - \Gamma/2.
\end{aligned} 
\quad (30) 
$$

where $\omega_{e_1 g_i} = (E_{e_1} - E_{g_i})/\hbar$, $S_0 = \sum_i V_i^1|e_1\rangle\langle g_i|$, and $M = -\sum_{i=1}^2 \frac{M_i \Gamma_i^1}{\Gamma}$, $M_1 = \sigma^\dagger\sigma \otimes \sigma^\dagger\sigma + \sigma^\dagger\sigma \otimes \sigma^\dagger\sigma$, $M_2 = \sigma\sigma^\dagger \otimes \sigma\sigma^\dagger + \sigma \otimes \sigma$. The operator $\sigma$ is defined as $\sigma = |g_2\rangle\langle g_1|$. After some algebra we get

$$ 
\begin{aligned}
L_{\text{eff}}(z) &= L_0 + zL_1 + O(z^2), \\
L_0 &= \sum_{i=1}^{4} (1 + M) \frac{h_i}{-\xi_i}, \\
L_1 &= \frac{M}{\Gamma} \sum_{i=1}^{4} \frac{h_i}{\xi_i} + (1 + M) \sum_{i=1}^{4} \frac{h_i}{-\xi_i^2}.
\end{aligned} 
\quad (31) 
$$

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 ($\Gamma_i^{1'} = 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 ($\Gamma_i^1 = \Gamma_i^{1'}$). 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) = \alpha e^{L_0 t}$ work well. We notice that the gap of $L_{\text{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 $\alpha$ 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.
<figure><img src="image_8.png" /><figcaption>FIG. 5. Same as 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$, $E_2 = 9$, $E_3 = 6$, $V_{12} = 0$, $V_1^1 = 1.0$, $V_2^1 = 0.7$. All energy values are in units of $V_1^1$, so that time is in units of $\hbar/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_{\text{eff}}(z)$ correspond to the nonlinear eigenvalue problem $[z - L_{\text{eff}}(z)]v = 0$, which can have in principle up to 9 distinct eigenvalues.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>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_{\text{eff}}(z)$ correspond to the nonlinear eigenvalue problem $[z - L_{\text{eff}}(z)]v = 0$, which can have in principle up to 9 distinct eigenvalues.</figcaption></figure>

## 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 these 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

trapping (CPT), the transition probability amplitudes to the excited state interfere destructively and the population of the excited state exactly vanishes [60–64,73,74].

*Coincidence for large values of the dissipation.* We calculate the limit of the effective operators $L_0$ and $L_1$, as $\Gamma/\omega_j \to \infty$. For this we recast them in terms of the smallness parameters $\delta_j = -i\omega_j/(\Gamma/2)$, for $j \in \{e_1g_1, g_1e_1, e_1g_2, g_2e_1\}$. As we take the limit of large dissipation $\lim_{\Gamma \to \infty} \delta_j = 0$, we get for the effective operators

$$ 
\begin{aligned}
L_{\text{eff}}^{(3\text{LS})}(z) &= L_0^{(3\text{LS})} + zL_1^{(3\text{LS})} + O(z^2), \\
L_0^{(3\text{LS})} &= \sum_j (1+M) \frac{h_j}{(1-\delta_j)\Gamma/2} \\
&\approx \sum_j (1+M) \frac{h_j}{\Gamma/2} \\
&= \frac{2}{n^{(1)}\pi\Gamma} L_0^{(\text{cont})},
\end{aligned}
$$
<figure><img src="image_10.png" /><figcaption>FIG. 7. Ground-state population and fidelity of the steady-state density matrix evolved according to the exact and several approximate operators as a function of the dissipation rate. The parameters of the calculation are $E_1 = 0.5$, $E_2 = -0.1$, $E_3 = 0.01$, $V_{12} = 0$, $V_1^{(1)} = 0.1$, $V_2^{(1)} = 0.15$. The calculation of $(QLQ)^{-1}$ is done perturbatively to lowest order in the Hamiltonian coupling.</figcaption></figure>

$$L_1^{(3LS)} = -\frac{M}{\Gamma} \sum_j \frac{h_j}{(1-\delta_j)\Gamma/2} + (1+M) \sum_j \frac{h_j}{-(1-\delta_j)^2(\Gamma/2)^2} \\ \approx \frac{2}{n^{(1)}\pi\Gamma} L_1^{(\text{cont})} - \frac{2}{n^{(1)}\pi\Gamma} L_0^{(\text{cont})}, \quad (32)$$

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/(\pi\Gamma)$, 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 $\Lambda$ 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 = \text{Tr}(\sqrt{\rho_{\text{exact}}\rho\sqrt{\rho_{\text{exact}}}})^2/[\text{Tr}(\rho)\text{Tr}(\rho_{\text{exact}})]$ by a factor $1 - [\text{Tr}(\rho) - \text{Tr}(\rho_{\text{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 $\rho$ fulfills the

<figure><img src="image_11.png" /><figcaption>FIG. 8. Ground-state population and fidelity of the steady-state density matrix evolved according to the exact and several approximate operators as a function of the dissipation rate, at the CPT condition. All traces overlap.</figcaption></figure>

following conditions [74]:

$$[PHP, \rho] = 0, \\ QHP\rho = 0. \quad (33)$$

<figure><img src="image_12.png" /><figcaption>FIG. 9. Steady-state population in the ground-state manifold $\rho_{g1g1}$ and $\rho_{g2g2}$, and in the excited-state manifold $\rho_{e1e1}$, with a three-level 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_2^{(1)} = 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_2^{(1)}$.</figcaption></figure>
Remarkably, this condition is independent of the value of the dissipation rate (as long as $\rho$ 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/(\pi\Gamma)$ then the effective operators will be the same. It also follows that since $\alpha = 1$, then $\langle L_1 \rangle = 0$.

We illustrate the effect of coherence population trapping points on our models in Figs. 8 and 9. By plotting the ground-state 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 $\Lambda$ 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.

## ACKNOWLEDGMENT

D.F.S. acknowledges support from the European Union through Marie Sklodowska-Curie Grant Agreement No. 590 702694.

## APPENDIX A: OPERATOR VERSUS SUPEROPERATOR NOTATION

For an $N$-level system, the underlying Hilbert space $\mathcal{H}$ is of dimension $N$ and the states of the quantum systems are described by positive operators acting on $\mathcal{H}$ that can be represented by $N \times N$ density matrices. The superoperators such as the Lindblad operator $L$ or its resolvent $G$ are linear operators acting on operators themselves acting on $\mathcal{H}$.

To describe an open quantum system we need to know the evolution of the density matrix using the Lindblad equation, which in Hilbert $\mathcal{H}$ space is written as

$$ \dot{\rho}(t) = L\rho = -i[H, \rho(t)] + \sum_i \left[ F_i \rho F_i^\dagger - \frac{1}{2} \{F_i^\dagger F_i, \rho(t)\} \right] \quad (A1) $$

or as

$$ \rho(t) = e^{Lt} \rho(0) = \frac{1}{2\pi i} \int dz \frac{e^{zt}}{z-L} \rho(0). \quad (A2) $$

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 $\vec{\rho}$ with $N^2$ components, obtained from the column-stretched form of the $N \times N$ density matrix. This representation is obtained by considering the density matrix as an element $\vec{\rho}$ of the Hilbert space $\mathcal{H} \otimes \mathcal{H}$ [40]. In that way, superoperators are linear operators acting on $\mathcal{H} \otimes \mathcal{H}$ and they can be represented by $N^2 \times N^2$ matrices. The linear superoperator acting on $\vec{\rho}$, built from 2 arbitrary operators $S_n$ and $S_m$ on $\mathcal{H}$ and acting on $\rho$ as $S_n \rho S_m^\dagger$, is given by the mapping $S_n \rho S_m^\dagger \to (\bar{S}_m \otimes S_n \rho)$. With the help of this mapping, the Lindblad operator operating on the vector form of the density matrix as $\frac{d}{dt}\vec{\rho}(t) = L\vec{\rho}(0)$ is

$$ L = -i(\mathbb{1} \otimes H - \bar{H} \otimes \mathbb{1}) + \sum_i \left\{ \bar{F}_i \otimes F_i - \frac{1}{2} [\mathbb{1} \otimes F_i^\dagger F_i + (F_i^\dagger F_i)^T \otimes \mathbb{1}] \right\}. \quad (A3) $$

## 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\rangle$ ($|j, q\rangle$) be the eigenstates of $PHP$ ($QHQ$), with $i = 1, 2, \dots, N_p$ ($j = 1, 2, \dots, 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_{ij} = \sqrt{\Gamma_{ij}}|i; p\rangle\langle j; q|$, relaxation from the slow subspace to the fast subspace described by jump operators $J_{ji} = \sqrt{\gamma_{ji}}|j; q\rangle\langle 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 operator $K = K_0 + V$ where

$$ K_0 = H_0 - \frac{i}{2} \left( \sum_{ij} F_{ij}^{\dagger} F_{ij} + \sum_{ij} J_{ji}^{\dagger} J_{ji} + \sum_{n} N_{m}^{\dagger} N_{m} \right) $$

has a diagonal matrix representation in the basis $\{|i; p\rangle, |j, q\rangle\}$. We can then rewrite the Lindblad operator $L$ as [see Eq. (A3)]

$$ L = -i(\mathbb{1} \otimes K - \bar{K} \otimes \mathbb{1}) + \sum_{ij} (\bar{F}_{ij} \otimes F_{ij} + \bar{J}_{ji} \otimes J_{ji}) + \sum_{n} \bar{N}_{n} \otimes N_{n}. $$
Using the expression of $\mathcal{Q}$ given by Eq. (2), we notice that $\bar{F}_{ij} \otimes F_{ij} = \mathcal{P}\bar{F}_{ij} \otimes F_{ij}\mathcal{Q}$, $\bar{J}_{ji} \otimes J_{ji} = \mathcal{Q}\bar{J}_{ji} \otimes J_{ji}\mathcal{P}$, and $\bar{N}_n \otimes N_n = \mathcal{P}\bar{N}_n \otimes N_n\mathcal{P}$.

Therefore $\mathcal{QLQ}$ can be written as

$$ \mathcal{QLQ} = L_D + W, $$

where $L_D = -i\mathcal{Q}(1 \otimes K_0 - \bar{K}_0 \otimes 1)\mathcal{Q}$ has a diagonal matrix representation in the basis $\{|i; p\rangle \otimes |j, q\rangle\}$ and $W = -i\mathcal{Q}(1 \otimes V - V \otimes 1)\mathcal{Q}$ has a nondiagonal matrix representation in the same basis.

The nondiagonal part of $\mathcal{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 $\text{ran}[\mathcal{Q}]$ is fast ($\Gamma_{ij} \gg V_{ij}$).

## 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

$$ H = H_0 + H_V + H_{V_p}, \\ H_0 = E_i|i\rangle\langle i| + \sum_a \int dk_a E_{k_a} |k_a\rangle\langle k_a|, \quad (C1) $$

$$ H_V = \sum_a \int dk_a [V_i^{(a)} |i\rangle\langle k_a| + (V_i^{(a)})^* |k_a\rangle\langle i|], \\ H_{V_p} = \sum_{i,j} [V_{ij} |i\rangle\langle j| + V_{ij}^* |j\rangle\langle i|]. \quad (C2) $$

With the dissipative part the Liouvillian is

$$ L = -i(1 \otimes H - \bar{H} \otimes 1) + \sum_{i,a} \mathcal{D}(F_i^{(a)}), \quad (C3) $$

$$ F_i^{(a)} = \sqrt{\Gamma_i^{(a)}} \int dk_a |g_i\rangle\langle k_a|. \quad (C4) $$

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 ($\mathcal{Q}$) and the discrete part ($\mathcal{P}$). The effective Liouvillian [see Eq. (7)] hinges on the resolvent operator in $\mathcal{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

$$ L_{\text{eff}}(z) = -i(1 \otimes H - \bar{H} \otimes 1) + \sum_i \sum_a^{N_g} \sum_{N_e} \mathcal{D}(F_{\text{eff},i}^{(a)}) \\ + \sum_i \sum_a^{N_g} \sum_{N_e} f^{(a)}(z) \bar{F}_{\text{eff},i}^{(a)} \otimes F_{\text{eff},i}^{(a)}, \quad (C5) $$

where $F_{\text{eff},i}^{(a)} = \sum_j \sqrt{\frac{\Gamma_j^{(a)}}{\sum_i \Gamma_i^{(a)}}} n^{(a)} \pi V_j^{(a)} |i\rangle\langle j|$ and $f^{(a)}(z) = -\frac{z}{z+\sum_i \Gamma_i^{(a)}}$ clearly vanishes when $z=0$. We recognize that the nonlinear operator can be written as a $z$-independent part in Lindblad form and a $z$-dependent part which involves only the incoherent jump $\bar{F}_{\text{eff},i}^{(a)} \otimes F_{\text{eff},i}^{(a)}$ from $|k_a k_a\rangle$ to $|g_i g_i\rangle$ that restores population to the states $g_i$. From the above expressions the specific cases in the examples section (Sec. III) can be straightforwardly derived.

## APPENDIX D: $L_{\text{eff}}(z)$ FOR A $\Lambda$ SYSTEM AT FINITE TEMPERATURE

We give the general expression for the effective operator of a $\Lambda$ system with incoherent transitions from the ground-state manifold to the excited states. The generalization of the operator presented in the main text is

$$ L_{\text{eff}}(z) = \left(1 + \frac{M\Gamma}{z + \Gamma}\right) \sum_{j=1}^{4} \frac{h_j}{z - \xi_j} \left(1 + \frac{M'\Gamma'}{z + \Gamma}\right) \\ + \frac{1}{2} \frac{MM'\Gamma\Gamma'}{z + \Gamma}, \quad (D1) $$

where we have used the notation $\xi_j = \{\pm i\omega_{e_1 g_1} - \Gamma/2 - \Gamma_1^{\prime 1}/2, \pm i\omega_{e_1 g_2} - \Gamma/2 - \Gamma_2^{\prime 1}/2\}$ and $\Gamma = \Gamma_1^{\prime 1} + \Gamma_2^{\prime 1}$. The $h_j$ matrices are defined as follows:

$$ \begin{aligned} h_1 &= \sigma \sigma^\dagger \otimes S_0^\dagger S_0, & \xi_1 &= -i\omega_{e_1 g_1} - \Gamma/2 - \Gamma_1^{\prime 1}/2, \\ h_2 &= S_0^\dagger S_0 \otimes \sigma \sigma^\dagger, & \xi_2 &= i\omega_{e_1 g_1} - \Gamma/2 - \Gamma_1^{\prime 1}/2, \\ h_3 &= \sigma^\dagger \sigma \otimes S_0^\dagger S_0, & \xi_3 &= -i\omega_{e_1 g_2} - \Gamma/2 - \Gamma_2^{\prime 1}/2, \\ h_4 &= S_0^\dagger S_0 \otimes \sigma^\dagger \sigma, & \xi_4 &= i\omega_{e_1 g_2} - \Gamma/2 - \Gamma_2^{\prime 1}/2, \end{aligned} \quad (D2) $$

where $S_0 = \sum_i V_i^\dagger |e_1\rangle\langle g_i|$, $M = -\sum_{i=1}^2 \frac{M_i \Gamma_i^{\prime 1}}{\Gamma}$, $M_1 = \sigma^\dagger \sigma \otimes \sigma^\dagger \sigma + \sigma^\dagger \sigma \otimes \sigma^\dagger \sigma$, $M_2 = \sigma \sigma^\dagger \otimes \sigma \sigma^\dagger + \sigma^\dagger \sigma \otimes \sigma^\dagger \sigma$, $M'_1 = \sigma^\dagger \sigma \otimes \sigma^\dagger \sigma + \sigma \otimes \sigma$, $M'_2 = \sigma \sigma^\dagger \otimes \sigma \sigma^\dagger + \sigma^\dagger \sigma \otimes \sigma^\dagger \sigma$. The operator $\sigma$ is defined as $\sigma = |g_2\rangle\langle g_1|$.

## APPENDIX E: $L_0$ GENERATOR OF A TRACE-PRESERVING DYNAMICS

Let us recall the expression for $L_0$:

$$ L_0 = \mathcal{P}\mathcal{L}\mathcal{P} - \mathcal{P}\mathcal{L}\mathcal{Q}(\mathcal{Q}\mathcal{L}\mathcal{Q})^{-1}\mathcal{Q}\mathcal{L}\mathcal{P}. \quad (E1) $$

We can rewrite this equation as

$$ L_0 = L\mathcal{A} - \mathcal{G}, \quad (E2) $$

where we have defined the operators $\mathcal{A}$ and $\mathcal{G}$ as follows:

$$ \mathcal{A} = \mathcal{P} - (\mathcal{Q}\mathcal{L}\mathcal{Q})^{-1}\mathcal{Q}\mathcal{L}\mathcal{P}, \quad (E3) $$

$$ \mathcal{G} = [1 - \mathcal{Q}\mathcal{L}\mathcal{Q}(\mathcal{Q}\mathcal{L}\mathcal{Q})^{-1}]\mathcal{Q}\mathcal{L}\mathcal{P} = \mathcal{Q}'\mathcal{Q}\mathcal{L}\mathcal{P}. \quad (E4) $$

In all the above equations, $A^{-1}$ signifies the Moore-Penrose inverse of $A$ [75] which coincides with the matrix inverse when $A$ is invertible. Finally, we have defined $\mathcal{Q}' = 1 - \mathcal{Q}\mathcal{L}\mathcal{Q}(\mathcal{Q}\mathcal{L}\mathcal{Q})^{-1}$, which is an orthogonal projector [76] to ker $[(\mathcal{Q}\mathcal{L}\mathcal{Q})^\dagger] = \text{ran}[\mathcal{Q}\mathcal{L}\mathcal{Q}]^\perp$, where $\perp$ stands for the orthogonal complement.
In order for $L_0$ to be a generator of a trace-preserving map, the maximally mixed state $[\rho] = \frac{1}{N}1_H$ must be a left eigenvector for $L_0$ with eigenvalue 0, where $N$ is the dimension of $\mathcal{H}$. In vector form, we can associate with $\bar{\rho}$ the maximally entangled state $|1_H\rangle \in \mathcal{H} \otimes \mathcal{H}$. Therefore, the trace-preserving condition can be written as

$$L_0^\dagger |1_H\rangle = \mathcal{G}^\dagger |1_H\rangle = 0, \qquad (E5)$$

where we have used the fact that $L$ is a Lindblad operator, hence generating a trace-preserving dynamics. If we define the set $\mathcal{X}$ as

$$\mathcal{X} = \text{ran}[QLQ]^\perp \cap \text{ran}[QLP], \qquad (E6)$$

then a sufficient condition for $L_0$ to generate a trace-preserving dynamics is that $|1_H\rangle$ is orthogonal to the set $\mathcal{X}$. Equivalently, since $\mathcal{X} \subseteq \text{ran}[Q]$, we can write this condition as

$$|\Psi_Q\rangle \perp \mathcal{X}, \qquad (E7)$$

where we have defined the state $|\Psi_Q\rangle = \sum_{i_q} |i_q\rangle \otimes |i_q\rangle$, 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 $\text{ran}[QLP] \subseteq \text{ran}[QLQ] = \text{ran}[Q]$ implying that $\mathcal{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_{\text{eff}}(z)$; therefore $\mathcal{P}G(z)\mathcal{P} = [T(z)]^{-1}$.

First we recall the definition of a nonlinear eigenvalue $\lambda$ of $T(z)$: $\lambda$ is an eigenvalue of $T(z)$ if $T(\lambda)v = 0$ for some nonzero vector $v$. The vector $v$ is the right eigenvector of $T$. The eigenvalue is called simple if in addition

$$\text{ker}[T(\lambda)] = \text{span}\{v\}, \quad v \neq 0, \quad T'(\lambda)v \notin \text{ran}[T(\lambda)].$$

In this case the adjoint $T^\dagger$ of $T$ satisfies

$$\text{ker}[T^\dagger(\lambda)] = \text{span}\{w\}$$

for some nonzero vector $w$, and furthermore, $w^\dagger T'(\lambda)v \neq 0$. Without loss of generality we can choose

$$w^\dagger T'(\lambda)v = 1, \qquad (F1)$$

where $T'(\lambda)$ is the value of the derivative of $T(z)$ with respect to $z$, taken at $z = \lambda$.

The Keldysh theorem states the following: Let $D$ be a compact subset that contains only simple eigenvalues $\lambda_n$, $n = 1, \dots, 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

$$T(z)^{-1} = \sum_{n=1}^{N} \frac{1}{z - \lambda_n} v_n w_n^\dagger + R(z). \qquad (F2)$$

Now, if we assume that all the eigenvalues of $T$ are simple, then we can use Eq. (F2) to calculate $\mathcal{P}\rho(t)$, performing the

integration of Eq. (6), and we obtain

$$\mathcal{P}\rho(t) = \sum_{n=1}^{N} e^{\lambda_n t} v_n w_n^\dagger \rho(0), \qquad (F3)$$

where we recall that we use the normalization $w_n^\dagger T'(\lambda) v_n = 1$. In particular for $\lambda = 0$, we have $T'(0) = 1 - L_1$; thus $w_0^\dagger(1 - L_1)\bar{\rho} = 1$, where $\bar{\rho} = v_0$ is the kernel of $L(z = 0) = L_0$.

## APPENDIX G: CORRECTION TO THE STEADY-STATE TRACE

We know that the steady state $\rho_f$ of the system will be in the kernel of $L_0$; that is, $\rho_f = \alpha\bar{\rho}$, where $L_0\bar{\rho} = 0$ and $\text{tr}[\bar{\rho}] = 1$. We are only left with determining the constant $\alpha$. This can be done from the final value theorem:

$$\rho_f = \lim_{z \to 0} z \mathcal{P}G(z) \mathcal{P}\rho(0) = \lim_{z \to 0} z[z - L_{\text{eff}}(z)]^{-1} \rho(0). \qquad (G1)$$

We expand $L_{\text{eff}}(z)$ as $L_{\text{eff}}(z) = L_0 + zL_1 + z^2L_2 + \dots$ and get

$$\begin{aligned} \rho_f &= \lim_{z \to 0} z[z - L_{\text{eff}}(z)]^{-1} \rho(0) \\ &= \lim_{z \to 0} \frac{z}{z(1 - L_1) - L_0 + O(z^2)} \rho(0). \end{aligned} \qquad (G2)$$

Multiplying by $(1 - L_1)$, and taking the limit, we obtain

$$[1 - L_1]\rho_f = [1 - L_1]\alpha\bar{\rho} = \lim_{z \to 0} \left[1 - \frac{1}{z}L_0(1 - L_1)^{-1}\right]^{-1} \rho(0).$$

Taking the trace of both sides, we obtain

$$\alpha = \frac{1}{\text{tr}[(1 - L_1)\bar{\rho}]}, \qquad (G3)$$

where we have used the fact that the dynamics generated by $L_0$ is trace preserving, implying that $\text{tr}[L_0\rho] = 0$ for all operators $\rho$ and where we have considered that $\text{tr}[\rho(0)] = 1$.

The same result can be obtained using the Keldysh theorem. Indeed, taking the limit $t \to \infty$ of Eq. (F3), we get $\rho_f = \bar{\rho}w_0^\dagger\rho(0) = \alpha\bar{\rho}$ and thus $\alpha = w_0^\dagger\rho(0)$, which in matrix form means

$$\alpha = \text{tr}[w_0^\dagger \rho(0)],$$

where $w_0$ is such that $L_0^\dagger w_0 = 0$, and $w_0^\dagger T'(0)\bar{\rho} = 1$; that is,

$$w_0^\dagger(1 - L_1)\bar{\rho} = 1. \qquad (G4)$$

But in matrix form, $[w_0]$ is proportional to the identity, $[w_0] = \beta1_H$. This is a consequence of the trace-preserving dynamics induced by $L_0$. Therefore, Eq. (G4) gives $\beta^*\text{tr}[(1 - L_1)\bar{\rho}] = 1$, and $\alpha = \beta^*\text{tr}[\rho(0)]$. Considering that $\text{tr}[\rho(0)] = 1$, we obtain the same result as in Eq. (G3).

## APPENDIX H: SIMULATIONS FOR THE A SYSTEM USING AN EXACT RESOLVENT

Figures 10 and 11 show the same dynamics as Figs. 5 and 6, but using an exact form of the resolvent $(QLQ)^{-1}$ valid to all orders of the Hamiltonian couplings between $P$ and $Q$.
<figure><img src="image_13.png" /><figcaption>FIG. 10. Same as Fig. 5, but using an exact form of the resolvent $(QLQ)^{-1}$ valid to all orders of the Hamiltonian couplings between $P$ and $Q$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 11. Same as Fig. 6, but using an exact form of the resolvent $(QLQ)^{-1}$ valid to all orders of the Hamiltonian couplings between $P$ and $Q$.</figcaption></figure>

[1] H. Haken, Z. Phys. B **20**, 413 (1975).

[2] H. Haken, *Synergetics: An Introduction* (Springer, Berlin, 1977).

[3] M. Lax, Phys. Rev. **157**, 213 (1967).

[4] C. Cohen-Tannoudji, Phys. Rep. **219**, 153 (1992).

[5] V. Paulisch, H. Rui, H. K. Ng, and B.-G. Englert, Eur. Phys. J. Plus **129**, 12 (2014).

[6] E. Brion, L. H. Pedersen, and K. Mølmer, J. Phys. A: Math. Theor. **40**, 1033 (2007).

[7] L. You, X. X. Yi, and X. H. Su, Phys. Rev. A **67**, 032308 (2003).

[8] D. Nagy, G. Kónya, G. Szirmai, and P. Domokos, Phys. Rev. Lett. **104**, 130401 (2010).

[9] J. S. Douglas, H. Habibian, C.-L. Hung, A. V. Gorshkov, H. J. Kimble, and D. E. Chang, Nat. Photon. **9**, 326 (2015).

[10] H.-P. Breuer, Phys. Rev. A **69**, 022115 (2004).

[11] G. Lindblad, Commun. Math. Phys. **48**, 119 (1976).

[12] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J. Math. Phys. **17**, 821 (1976).

[13] F. Minganti, A. Biella, N. Bartolo, and C. Ciuti, Phys. Rev. A **98**, 042118 (2018).

[14] M. Mirrahimi and P. Rouchon, IEEE Trans. Automatic Control **54**, 1325 (2009).

[15] R. Azouit, A. Sarlette, and P. Rouchon, arXiv:1603.04630.

[16] R. Azouit, F. Chittaro, A. Sarlette, and P. Rouchon, Quantum Sci. Technol. **2**, 044011 (2017).

[17] R. Azouit, F. Chittaro, A. Sarlette, and P. Rouchon, IFAC-PapersOnLine **50**, 13026 (2017).

[18] R. Azouit, Ph.D. thesis, PSL Research University, 2017.

[19] P. Forni, A. Sarlette, T. Capelle, E. Flurin, S. Deléglise, and P. Rouchon, arXiv:1803.07810.

[20] A. N. Tikhonov, Mat. Sb. (N.S.) **31**, 575 (1952).

[21] N. Fenichel, J. Diff. Equ. **31**, 53 (1979).

[22] L. Noethen and S. Walcher, Discr. Contin. Dyn. Syst. Ser. B **16**, 945 (2011).

[23] F. Reiter and A. S. Sørensen, Phys. Rev. A **85**, 032111 (2012).

[24] I. Lesanovsky and J. P. Garrahan, Phys. Rev. Lett. **111**, 215305 (2013).

[25] M. Marcuzzi, J. Schick, B. Olmos, and I. Lesanovsky, J. Phys. A: Math. Theor. **47**, 482001 (2014).

[26] K. Macieszczak, M. Guta, I. Lesanovsky, and J. P. Garrahan, Phys. Rev. Lett. **116**, 240404 (2016).

[27] V. V. Albert, K. Noh, and F. Reiterr, arXiv:1809.07324.

[28] N. Lutkenhaus, J. I. Cirac, and P. Zoller, Phys. Rev. A **57**, 548 (1998).

[29] F. Damanet, A. J. Daley, and J. Keeling, Phys. Rev. A **99**, 033845 (2019).

[30] P. Warszawski and H. M. Wiseman, Phys. Rev. A **63**, 013803 (2000).

[31] D. Burgarth, P. Facchi, H. Nakazato, S. Pascazio, and K. Yuasa, Quantum **3**, 152 (2019).

[32] H. Feshbach, Ann. Phys. **19**, 287 (1962).

[33] P. Degenfeld-Schonburg and M. J. Hartmann, Phys. Rev. B **89**, 245108 (2014).

[34] U. Fano, Phys. Rev. **124**, 1866 (1961).

[35] A. E. Miroshnichenko, S. Flach, and Y. S. Kivshar, Rev. Mod. Phys. **82**, 2257 (2010).

[36] B. Lukyanchuk, N. I. Zheludev, S. A. Maier, N. J. Halas, P. Nordlander, H. Giessen, and C. T. Chong, Nat. Mater. **9**, 707 (2010).

[37] D. Finkelstein-Shapiro, I. Urdaneta, M. Calatayud, O. Atabek, V. Mujica, and A. Keller, Phys. Rev. Lett. **115**, 113006 (2015).

[38] D. Finkelstein-Shapiro and A. Keller, Phys. Rev. A **97**, 023411 (2018).
[39] D. Finkelstein-Shapiro, M. Calatayud, O. Atabek, V. Mujica, and A. Keller, Phys. Rev. A **93**, 063414 (2016).

[40] T. F. Havel, *J. Math. Phys.* **44**, 534 (2003).

[41] In the case where the system Hamiltonian has a continuous spectrum, the discrete sum can be replaced by an integral over the generalized Hamiltonian eigenstates.

[42] U. Fano, *Nuovo Cimento* **12**, 154 (1935).

[43] G. Baldini, Phys. Rev. **128**, 1562 (1962).

[44] K. P. Jain, Phys. Rev. **139**, A544 (1965).

[45] S. Glutsch, U. Siegner, M.-A. Mycek, and D. S. Chemla, Phys. Rev. B **50**, 17009 (1994).

[46] U. Siegner, M. A. Mycek, S. Glutsch, and D. S. Chemla, Phys. Rev. Lett. **74**, 470 (1995).

[47] U. Siegner, M.-A. Mycek, S. Glutsch, and D. S. Chemla, Phys. Rev. B **51**, 4953 (1995).

[48] U. Siegner, S. Bar-Ad, and D. Chemla, Chem. Phys. **210**, 155 (1996).

[49] R. P. Seisyan, Phys. Solid State **58**, 859 (2016).

[50] C. P. Holfeld, F. Löser, M. Sudzius, K. Leo, D. M. Whittaker, and K. Köhler, Phys. Rev. Lett. **81**, 874 (1998).

[51] S. Yoshino, G. Oohata, and K. Mizoguchi, Phys. Rev. Lett. **115**, 157402 (2015).

[52] V. May and O. Kühn, *Charge and Energy Transfer Dynamics in Molecular Systems* (Wiley-VCH, Weinheim, Germany, 2011).

[53] G. Schatz and M. Ratner, *Quantum Mechanics in Chemistry*, Dover Books on Chemistry (Dover Publications, Mineola, New York, 2002).

[54] K. G. Zloshchastiev and A. Sergi, J. Mod. Opt. **61**, 1298 (2014).

[55] S. Echeverri-Arteaga, H. Vinck-Posada, and E. A. Gómez, Optik **180**, 505 (2019).

[56] The nonlinear eigenvalues of $L_{\text{eff}}(z)$ are the complex numbers $\lambda$ satisfying $[L_{\text{eff}}(\lambda) - \lambda I]\rho = 0$, for some nonzero $\rho$.

[57] M. V. Keldysh, Dokl. Akad. Nauk SSSR (N.S.) **77**, 11 (1951).

[58] M. V. Keldysh, Usp. Mat. Nauk **26** (4), 15 (1971).

[59] W.-J. Beyn, *Lin. Alg. Its Applic.* **436**, 3839 (2012).

[60] P. M. Radmore and P. L. Knight, *J. Phys. B: At. Mol. Phys.* **15**, 561 (1982).

[61] K. Bergmann, H. Theuer, and B. W. Shore, Rev. Mod. Phys. **70**, 1003 (1998).

[62] K.-J. Boller, A. Imamoğlu, and S. E. Harris, Phys. Rev. Lett. **66**, 2593 (1991).

[63] M. Fleischhauer, A. Imamoğlu, and J. P. Marangos, Rev. Mod. Phys. **77**, 633 (2005).

[64] N. V. Vitanov, A. A. Rangelov, B. W. Shore, and K. Bergmann, Rev. Mod. Phys. **89**, 015006 (2017).

[65] J. Vanier, A. Godone, and F. Levi, Phys. Rev. A **58**, 2345 (1998).

[66] S. Sevinçli, C. Ates, T. Pohl, H. Schempp, C. S. Hofmann, G. Günter, T. Amthor, M. Weidemüller, J. D. Pritchard, D. Maxwell, A. Gauguet, K. J. Weatherill, M. P. A. Jones, and C. S. Adams, J. Phys. B **44**, 184018 (2011).

[67] M. Kasevich and S. Chu, Phys. Rev. Lett. **69**, 1741 (1992).

[68] G. Morigi, J. Eschner, and C. H. Keitel, Phys. Rev. Lett. **85**, 4458 (2000).

[69] A. Aspect, E. Arimondo, R. Kaiser, N. Vansteenkiste, and C. Cohen-Tannoudji, J. Opt. Soc. Am. B **6**, 2112 (1989).

[70] The ACME Collaboration, J. Baron, W. C. Campbell, D. DeMille, J. M. Doyle, G. Gabrielse, Y. V. Gurevich, P. W. Hess, N. R. Hutzler, E. Kirilov, I. Kozyryev, B. R. O'Leary, C. D. Panda, M. F. Parsons, E. S. Petrik, B. Spaun, A. C. Vutha, and A. D. West, Science **343**, 269 (2014).

[71] A. Dantan, J. Cviklinski, E. Giacobino, and M. Pinard, Phys. Rev. Lett. **97**, 023605 (2006).

[72] H. Schempp, G. Günter, C. S. Hofmann, C. Giese, S. D. Saliba, B. D. DePaola, T. Amthor, M. Weidemüller, S. Sevinçli, and T. Pohl, Phys. Rev. Lett. **104**, 173602 (2010).

[73] B. W. Shore, Adv. Opt. Photon. **9**, 563 (2017).

[74] D. Finkelstein-Shapiro, S. Felicetti, T. Hansen, T. Pullerits, and A. Keller, Phys. Rev. A **99**, 053829 (2019).

[75] R. Penrose, Math. Proc. Cambr. Philos. Soc. **51**, 406 (1955).

[76] G. H. Golub and C. F. van Loan, in *Matrix Computations*, 4th ed. (Johns Hopkins University Press, Baltimore, Maryland, 2013), p. 290.