---
abstract: 'We consider an open quantum system described by a Lindblad-type master equation with two times-scales. The fast time-scale is strongly dissipative and drives the system towards a low-dimensional decoherence-free space. To perform the adiabatic elimination of this fast relaxation, we propose a geometric asymptotic expansion based on the small positive parameter describing the time-scale separation. This expansion exploits geometric singular perturbation theory and center-manifold techniques. We conjecture that, at any order, it provides an effective slow Lindblad master equation and a completely positive parameterization of the slow invariant sub-manifold associated to the low-dimensional decoherence-free space. By preserving complete positivity and trace, two important structural properties attached to open quantum dynamics, we obtain a reduced-order model that directly conveys a physical interpretation since it relies on effective Lindbladian descriptions of the slow evolution. At the first order, we derive simple formulae for the effective Lindblad master equation. For a specific type of fast dissipation, we show how any Hamiltonian perturbation yields Lindbladian second-order corrections to the first-order slow evolution governed by the Zeno-Hamiltonian. These results are illustrated on a composite system made of a strongly dissipative harmonic oscillator, the ancilla, weakly coupled to another quantum system.'
author:
- 'R. Azouit[^1], A. Sarlette[^2], P. Rouchon'
title: '**Adiabatic elimination for open quantum systems with effective Lindblad master equations [^3]** '
---

Introduction
============

Solving the equation of evolution for a open quantum system - the Lindblad master equation [@BreuerPetruccioneBook] - is generally tedious. To gain better physical insight and/or for numerical simulations, it is of wide interest to compute rigorous reduced models of quantum dynamical systems. In a typical case, a system of interest is coupled to an ancillary system expressing a measurement device or a perturbing environment [@haroche-raimondBook06]. The quantum dynamics describes the joint evolution of both systems and in order to focus only on the system of interest we want to determine a dynamical equation for the system of interest only, from which we have “eliminated” the ancillary system.

A standard tool for model reduction is to use the different timescales of the complete system to separate the quantum dynamics into fast and slow variables and then eliminate the fast ones. This technique is known as adiabatic elimination. In quantum Hamiltonian systems, regular perturbation theory can be easily applied as the propagator remains unitary, and the construction of the reduced model to various orders of approximation is standard [@sakurai2011modern]. In contrast, for open quantum systems, described by a Lindblad master equation [@BreuerPetruccioneBook], the case is much more complicated and involves singular perturbation theory. Several particular examples have been treated separately. In [@BrionJPA07] different methods are proposed to perform an adiabatic elimination up to second-order on a lambda system. In [@mirrahimi-rouchonieee09] and [@ReiteS2012PRA] the problem of excited states decaying towards $n$ grounds state is treated. A specific atom-optics dynamics is investigated in [@AtkinsPRA03]. In the presence of continuous measurement, [@CernotikPRA15] presents a method of adiabatic elimination for systems with Gaussian dynamics.

However, no general and systematic method has been developed yet for adiabatic elimination in systems with Lindblad dynamics. More precisely, treating the Lindblad master equation as a usual linear system, or applying the Schrieffer-Wolff formalism which is generalized in [@Kessl2012PRA] to Lindblad dynamics, requires the inversion of super-operators which can be troublesome both numerically and towards physical interpretation. In [@AzouitCDC15] we made a first attempt to circumvent this inversion, using invariants of the dynamics. This provides a first-order expansion only, and in general linear form — i.e. not necessarily with the structure of a Lindblad equation.

In the present paper, we propose a geometric method to perform an adiabatic elimination for open quantum systems, whose key feature is that *the resulting reduced model is explicitly described by an effective Lindblad equation. The reduced system is parameterized by a reduced density operator and the mapping from the reduced model to the initial system state space is expressed in terms of Kraus operators, ensuring a trace-preserving completely positive map.* By preserving these structural properties of open quantum dynamics, we obtain a reduced model that directly conveys a physical interpretation. As far as we know, combining asymptotic expansion with completely positive map and Lindbladian formulation has never been addressed before. This work is a first attempt to investigate the interest of such combination with lemmas \[lem:Ls1\], \[lem:K1\] and \[lem:Fs2\] underlying the conjecture illustrated on figure \[fig:GeomAdiabElimLindblad\].

Our method applies to general open quantum systems with two timescales, described by two general Lindbladian super-operators , and where the fast Lindbladian makes the system converge to a decoherence-free subspace of the overall Hilbert space. We then use a geometric approach based on center manifold techniques [@carr-book] and geometric singular perturbation theory [@Fenichel79] to obtain an expansion of the effect of the perturbation introduced by the slow Lindbladian on this decoherence-free subspace. For general Lindbladians satisfying this setting, we get explicit formulas for the Lindblad operators describing the first-order expansion. In the particular case of a Hamiltonian perturbation, we retrieve the well known Zeno effect. Furthermore, for a fast Lindbladian described by a single decoherence operator and subject to a Hamiltonian perturbation, we derive explicit formulas for the first-order effect on the location of the center manifold and for Lindblad operators describing the second-order expansion of the dynamics. This allows to highlight how a first-order Zeno effect is associated to second-order decoherence.

We apply our method to a quantum system coupled to a highly dissipative quantum harmonic oscillator (ancilla). Our general formulas directly provide an effective Lindblad master equation of the reduced model where this ancilla is eliminated. The result for this example is well known [@CarmichaelBook07], which allows us to emphasize how the correct results are obtained also on infinite-dimensional systems, and to appreciate the computational simplicity of our method in comparison with previous ones.

The paper is organized as follows. Section \[sec:system\] presents the structure of two-timescales master equation for open quantum systems, as well as the assumptions and properties of the unperturbed system used for deriving our results. In section \[sec:firstorder\] we present a geometric approach for performing the adiabatic elimination and derive a first-order reduced model for arbitrary perturbations. In section \[sec:secondorder\] we develop the second-order expansion for a class of systems. In section \[sec:example\] we illustrate the method where the ancilla is a highly dissipative harmonic oscillator.

A class of perturbed master equations {#sec:system}
=====================================

Denote by $\mathcal{H}$ a Hilbert space of finite dimension, by $\mathcal{D}$ the compact convex set of density operators $\rho$ on $\mathcal{H}$ ($\rho$ is Hermitian, nonnegative and trace one). We consider a two-time scale dynamics on $\mathcal{D}$ described by the following master differential equation $$\label{eq:DynEpsilon}
 {\frac{d}{dt}}\rho = {\mathfrak{L}}_0(\rho) + \epsilon {\mathfrak{L}}_1(\rho)$$ where $\epsilon$ is a small positive parameter and the linear super-operators ${\mathfrak{L}}_0$ and ${\mathfrak{L}}_1$ are of Lindbladian forms [@BreuerPetruccioneBook]. That is, there exist two finite families of operators on $\mathcal{H}$, denoted by $(L_{0,\nu})$ and $(L_{1,\nu})$, and two Hermitian operators $H_0$ and $H_1$ (called Hamiltonians) such that, for $r=0,1$, we have $$\label{eq:GenLindblad}
{\mathfrak{L}}_r(\rho)= -i[H_r,\rho] + \sum_{\nu} L_{r,\nu}\rho L_{r,\nu}^\dag - \tfrac{1}{2} L_{r,\nu}^\dag   L_{r,\nu}\rho -\tfrac{1}{2} \rho L_{r,\nu}^\dag   L_{s,\nu}
.$$ It is well known that if the initial condition $\rho(0)$ belongs to $\mathcal{D}$, then the solution $\rho(t)$ of  remains in $\mathcal{D}$ and is defined for all $t\geq 0$. It is also a known fact that the flow of  is a contraction for many distances, such as the one derived from the nuclear norm $\|\cdot\|_1$: for any trajectories $\rho_1$ and $\rho_2$ of , we have $\| \rho_1(t) -\rho_2(t) \|_1 \leq \| \rho_1(t') -\rho_2(t') \|_1$ for all $t\geq t'$; see e.g. [@nielsen-chang-book Th.9.2].

We assume that, for $\epsilon=0$, the unperturbed master equation ${\frac{d}{dt}}\rho ={\mathfrak{L}}_0(\rho)$ converges to a stationary regime. More precisely, we assume that the unperturbed master equation admits a sub-manifold of stationary operators coinciding with the $\Omega$ limit set of its trajectories. Denote by $$\mathcal{D}_0=\left\{\rho\in\mathcal{D}~\big|~{\mathfrak{L}}_0(\rho)=0\right\}$$ this stationary manifold: it is compact and convex. We thus assume that for all $\rho_0\in\mathcal{D}$, the solution of ${\frac{d}{dt}}\rho={\mathfrak{L}}_0(\rho)$ with $\rho(0)=\rho_0$ converges for $t$ tending to $+\infty$ towards an element of $\mathcal{D}_0$ denoted by $R(\rho_0)= \lim_{t\rightarrow +\infty} \rho(t)$. Since for any $t\geq 0$, the propagator $ e^{t{\mathfrak{L}}_0}$ is a completely positive linear map [@nielsen-chang-book Chap.8], $R$ is also a completely positive map. By Choi’s theorem [@ChoiLAA1975] there exists a finite set of operators on $\mathcal{H}$ denoted by $(M_\mu)$ such that $$\label{eq:Rkraus}
   R(\rho_0)= \sum_{\mu} M_\mu \rho_0 M^\dag_\mu$$ with $\sum_{\mu} M_\mu^\dag M_\mu = I$, the identity operator on $\mathcal{H}$. The form is called a Kraus map. We thus assume that $$\mathcal{D}_0= \left\{ R(\rho)~\big|~\rho \in\mathcal{D} \right\}\text{ and }
\forall \rho\in\mathcal{D}_0, ~R(\rho)=\rho.$$

An invariant operator attached to the dynamics ${\frac{d}{dt}}\rho = {\mathfrak{L}}_0(\rho)$ is an Hermitian operator $J$ such that for any time $t \geq 0$ and any initial state $\rho_0=\rho(0)$, we have ${\text{Tr}\left(J \rho(t)\right)} = {\text{Tr}\left(J \rho_0\right)}$. Such invariant operators $J$ are characterized by the fact that ${\mathfrak{L}}_0^*(J)=0$ where the adjoint map to ${\mathfrak{L}}_0$ is given by $${\mathfrak{L}}^*_0(A)= i[H_0,A] +\sum_{\nu} L_{0,\nu}^\dag A L_{0,\nu} - \tfrac{1}{2} L_{0,\nu}^\dag   L_{0,\nu}A -\tfrac{1}{2} A L_{0,\nu}^\dag   L_{0,\nu}$$ for any Hermitian operator $A$.

Thus by taking the limit for $t$ tending to $+\infty$ in ${\text{Tr}\left(J \rho(t)\right)} = {\text{Tr}\left(J \rho_0\right)}$, we have, for all Hermitian operators $\rho_0$, ${\text{Tr}\left(J R(\rho_0)\right)} = {\text{Tr}\left(J \rho_0\right)}$. Denote by $R^*$ the adjoint map associated to $R$: $$\label{eq:Radjoint}
  R^*(A)= \sum_{\mu} M_{\mu}^\dag A M_\mu$$ for any Hermitian operator $A$ on $\mathcal{H}$. Then, ${\text{Tr}\left(R^*(J)\rho_0 \right)} = {\text{Tr}\left(J \rho_0\right)}$ for any Hermitian operator $\rho_0$ is equivalent to $R^*(J)= J $. I.e. invariant operators $J$ are characterized by ${\mathfrak{L}}_0^*(J)=0$ and satisfy $R^*(J)=J$.

We assume additionally that $\mathcal{D}_0$ coincides with the set of density operators with support in $\mathcal{H}_0$, a subspace of $\mathcal{H}$. In other words the unperturbed master equation features a decoherence-free space $\mathcal{H}_0$. Denote by $P_0$ the operator on $\mathcal{H}$ corresponding to orthogonal projection onto $\mathcal{H}_0$. Consequently, for any Hermitian operator $\rho$, we have $P_0 R(\rho)=R(\rho) P_0=R(\rho)$. Thus ${\text{Tr}\left(R^*(P_0) \rho\right)}={\text{Tr}\left(R(\rho)\right)}={\text{Tr}\left(\rho\right)}$ for all $\rho$ which implies: $$\label{eq:RP0}
R^{*}(P_0)=I
.$$ Moreover, for any vector ${\ensuremath{|c\rangle}\xspace}$ in $\mathcal{H}_0$, $R({\ensuremath{|c\rangle}\xspace}{\ensuremath{\langle c|}\xspace})= {\ensuremath{|c\rangle}\xspace}{\ensuremath{\langle c|}\xspace}$. This implies that, for the Kraus map , there exists a family of complex numbers $\lambda_\mu$ such that $\sum_{\mu} |\lambda_\mu|^2=1$ and $$\label{eq:MmuC}
  \forall {\ensuremath{|c\rangle}\xspace}\in\mathcal{H}_0,~M_\mu {\ensuremath{|c\rangle}\xspace} = \lambda_\mu {\ensuremath{|c\rangle}\xspace}
  .$$

First order expansion for arbitrary perturbations {#sec:firstorder}
=================================================

We consider here the perturbed master equation  whose unperturbed part ${\frac{d}{dt}}\rho = {\mathfrak{L}}_0(\rho)$ satisfies the assumptions of Section \[sec:system\]: any trajectory converges to a steady-state; the set of steady-states $\mathcal{D}_0$ coincides with the set of density operators with support on a subspace $\mathcal{H}_0$ of $\mathcal{H}$. This section develops a first-order expansion versus $\epsilon$ of  around $\mathcal{D}_0$.

Denote by $\mathcal{H}_s$ (subscript $s$ for slow), an abstract Hilbert space with the same dimension as $\mathcal{H}_0$. Denote by $\mathcal{D}_s$ the set of density operators on $\mathcal{H}_s$. Denote by $\{{\ensuremath{|\nu\rangle}\xspace}\}$ (resp. $\{{\ensuremath{|c_\nu\rangle}\xspace}\}$) a Hilbert basis of $\mathcal{H}_s$ (resp. $\mathcal{H}_0$). Consider the Kraus map $K_0$ defined by $$\label{eq:K0}
  \forall \rho_s\in\mathcal{D}_s,~K_0(\rho_s) = S_0 \rho_s S_0^\dag \in\mathcal{D}$$ where $S_0=\sum_\nu {\ensuremath{|c_\nu\rangle}\xspace} {\ensuremath{\langle \nu|}\xspace}$. We have $S_0 S_0^\dag =P_0$, the orthogonal projector onto $\mathcal{H}_0$ and $S_0^\dag S_0=I_s$, the identity operator on $\mathcal{H}_s$.

![Adiabatic elimination, based on geometric singular perturbation theory, of the fast relaxation dynamics ${\frac{d}{dt}}\rho = {\mathfrak{L}}_0(\rho)$ for an open-quantum system governed by the Lindbladian master equation ${\frac{d}{dt}}\rho = {\mathfrak{L}}(\rho)= {\mathfrak{L}}_0(\rho) + \epsilon {\mathfrak{L}}_1(\rho)$ where $\epsilon$ is a small positive parameter. It provides two asymptotic expansions of the slow dynamics. The parameterization via density operators $\rho_s$ of the slow invariant attractive sub-manifold (related to ${\mathfrak{L}}_0(\rho)=0$) is based on the map $\rho=K(\rho_s)=K_0(\rho_s)+\epsilon K_1(\rho_s) + \ldots$. The slow dynamics corresponds to ${\frac{d}{dt}}\rho_s = {\mathfrak{L}}_s(\rho_s)= \epsilon {\mathfrak{L}}_{s,1}(\rho_s)+ \epsilon^2 {\mathfrak{L}}_{s,2}(\rho_s)+\ldots$. The super-operators $K_0$, $K_1$, …and ${\mathfrak{L}}_{s,1}$, ${\mathfrak{L}}_{s,2}$, …are obtained by identifying terms of identical order versus $\epsilon$ in the geometric invariance condition\
$
    {\mathfrak{L}}_0\big(K_0+\epsilon K_1\big)+ \epsilon {\mathfrak{L}}_1\big(K_0+\epsilon K_1+\ldots\big)$\
$\hspace*{2.em} = K_0\big(\epsilon {\mathfrak{L}}_{s,1} + \epsilon^2 {\mathfrak{L}}_{s,2}+\ldots\big)
   + \epsilon K_1 \big(\epsilon {\mathfrak{L}}_{s,1} + \epsilon^2 {\mathfrak{L}}_{s,2}+\ldots\big)  + \ldots~.
  $\
We conjecture that, at any order versus $\epsilon$, the super-operator $K$ is a Kraus map (up to higher-order corrections) and the slow evolution $ {\frac{d}{dt}}\rho_s={\mathfrak{L}}_s(\rho_s)$ is of Lindbladian type. []{data-label="fig:GeomAdiabElimLindblad"}](GeomAdiabElimLindblad "fig:"){width="45.00000%"}\

As illustrated on figure \[fig:GeomAdiabElimLindblad\], we are looking for an expansion based on linear super-operators $\{K_m\}_{m\geq 0}$ between $\mathcal{D}_s$ and $\mathcal{D}$ and on Lindblad dynamics $\left\{{\mathfrak{L}}_{s,m}\right\}_{m\geq 0}$ on $\mathcal{D}_{s}$ such that any solution $t \mapsto \rho_s(t)\in\mathcal{D}_s$ of the reduced Lindblad master equation $$\label{eq:DynRed}
    {\frac{d}{dt}}\rho_s = {\mathfrak{L}}_s(\rho_s)=\sum_{m\geq 0}   \epsilon^m {\mathfrak{L}}_{s,m}(\rho_s)$$ yields, via the map $$\label{eq:KrausRed}
   \rho(t)  = K(\rho_s(t))=\sum_{m\geq 0} \epsilon^m K_m(\rho_s(t)) \; ,$$ a trajectory of the perturbed system . We combine here geometric singular perturbation theory [@Fenichel79] with center manifold techniques based on Carr asymptotic expansion lemma [@carr-book] to derive recurrence relationships for $K_m$ and ${\mathfrak{L}}_{s,m}$. These recurrences are obtained by identifying the terms of the same order in the formal invariance condition $${\mathfrak{L}}_{0}\big(K(\rho_s)\big)+ \epsilon {\mathfrak{L}}_{1}\big(K(\rho_s)\big) =  {\frac{d}{dt}}\rho = K\left({\frac{d}{dt}}\rho_s\right) = K \big( {\mathfrak{L}}_s(\rho_s)\big)
.$$ This means that, for any $\rho_s\in\mathcal{D}_s$, we have $$\begin{gathered}
\label{eq:invariance}
  {\mathfrak{L}}_{0}\left( \sum_{m\geq 0} \epsilon^m K_m(\rho_s) \right) + \epsilon {\mathfrak{L}}_{1}\left( \sum_{m\geq 0} \epsilon^m K_m(\rho_s) \right)
  \\
  = \sum_{m} \epsilon^m K_{m}\left(\sum_{m'} \epsilon^{m'} {\mathfrak{L}}_{s,m'}(\rho_s) \right)
  .\end{gathered}$$ For $m=0$ we have $$\label{eq:order0}
  {\mathfrak{L}}_0\big(K_0(\rho_s)\big) = K_0\big({\mathfrak{L}}_{s,0}(\rho_s)\big)
  .$$ With $K_0$ defined in , we have ${\mathfrak{L}}_0\big(K_0(\rho_s)\big) \equiv 0$ and thus ${\mathfrak{L}}_{s,0}(\rho_s)=0$. Consequently, for $m\geq 1$, we have $$\begin{gathered}
\label{eq:orderm}
{\mathfrak{L}}_{0}\left(  K_m(\rho_s) \right) + {\mathfrak{L}}_{1}\left(  K_{m-1}(\rho_s) \right)
\\
=
\sum_{m'=1}^m K_{m-m'}\left( {\mathfrak{L}}_{s,m'}(\rho_s) \right)
.
 \end{gathered}$$ When $m=1$, $K_1$ and ${\mathfrak{L}}_{s,1}$ are defined by $$\label{eq:order1}
  {\mathfrak{L}}_0\big(K_1(\rho_s)\big) + {\mathfrak{L}}_1\big(K_0(\rho_s)\big) = K_0\big({\mathfrak{L}}_{s,1}(\rho_s)\big)
  .$$ The following lemma proves that the super-operator ${\mathfrak{L}}_{s,1}(\rho_s)$ defined by this equation is always of Lindblad form.

\[lem:Ls1\]

Assume that ${\mathfrak{L}}_1(\rho)= - i[H_1,\rho]$ for some Hermitian operator $H_1$ on $\mathcal{H}$. Then, if ${\mathfrak{L}}_{s,1}$ satisfies , we have ${\mathfrak{L}}_{s,1}(\rho_s) = - i [H_{s,1},\rho_s]$ where $H_{s,1}= S_0^\dag H_1 S_0$ is a Hermitian operator on $\mathcal{H}_s$.

Assume that ${\mathfrak{L}}_1(\rho)= L_1 \rho L_1^\dag - \tfrac{1}{2} \big(L_1^\dag L_1 \rho +\rho L_1^\dag L_1  \big)$ for some operator $L_1$ on $\mathcal{H}$. Then, if ${\mathfrak{L}}_{s,1}$ satisfies , we have $$\label{eq:Ls1}
  {\mathfrak{L}}_{s,1}(\rho_s)= \sum_\mu A_\mu \rho_s  A_\mu^\dag - \tfrac{1}{2} \big(A_\mu^\dag A_\mu \rho_s +\rho_s A_\mu^\dag A_\mu  \big)$$ with $A_\mu= S_0^\dag M_\mu L_1 S_0$ and the Kraus operators $M_\mu$ defined by .

The result for a general Lindbladian dynamics for $r=1$ follows by linearity.

Since $R\circ {\mathfrak{L}}_{0} =  0$ and $R\circ K_0=K_0$ , we have $R\left({\mathfrak{L}}_1\big(K_0(\rho_s)\big) \right) =K_0\big({\mathfrak{L}}_{s,1}(\rho_s)\big)$. Left multiplication by $S_0^\dag$ and right multiplication by $S_0$ yields $${\mathfrak{L}}_{s,1}(\rho_s) = S^\dag_0 R\left({\mathfrak{L}}_1\Big(S_0\rho_sS^\dag_0 \Big) \right)  S_0$$ since $S^\dag_0 S_0=I_s$ is the identity operator on $\mathcal{H}_s$.

For ${\mathfrak{L}}_1(\rho)= - i[H_1,\rho]$ we have, exploiting the fact that $M_\mu S_0= \lambda_\mu S_0$ and $S_0^\dag S_0=I_s$: $$\begin{gathered}
\label{eq:Trick}
  S^\dag_0  R\left({\mathfrak{L}}_1\Big(S_0\rho_sS^\dag_0 \Big)\right)  S_0
  \\
 = - i \sum _\mu S_0^\dag M_\mu \big( H_1 S_0\rho_sS^\dag_0- S_0\rho_sS^\dag_0 H_1\big) M_\mu^\dag S_0
 \\
 = - i \sum _\mu \left(\lambda_\mu^* S_0^\dag \right)M_\mu  H_1 S_0\rho_s + i \sum_\mu  \rho_s S^\dag_0 H_1M_\mu^\dag  \big(\lambda_\mu S^0\big)
 \\
  = - i \sum _\mu  S_0^\dag  M_\mu^\dag M_\mu  H_1 S_0\rho_s + i \sum_\mu  \rho_s S^\dag_0 H_1M_\mu^\dag  M_\mu  S^0
  \\
  = - i \left[  S_0^\dag  H_1  S_0~,~\rho_s \right]\end{gathered}$$ since $\sum_\mu M_\mu^\dag M_\mu = I$. We get the Zeno Hamiltonian $H_{s,1}= S_0^\dag  H_1  S_0$.

For ${\mathfrak{L}}_1(\rho)= L_1 \rho L_1^\dag - \tfrac{1}{2} \big(L_1^\dag L_1 \rho +\rho L_1^\dag L_1  \big)$, similar computations yield $$\begin{gathered}
  S^\dag_0  R\left({\mathfrak{L}}_1\Big(S_0\rho_sS^\dag_0 \Big)\right)  S_0
  \\ = \sum _\mu S_0^\dag M_\mu L_1 S_0\rho_sS^\dag_0 L_1^\dag  M_\mu^\dag S_0
 \\ - \tfrac{1}{2}\sum _\mu S_0^\dag M_\mu  \big(L_1^\dag L_1 S_0\rho_sS^\dag_0 -  S_0\rho_sS^\dag_0 L_1^\dag L_1 \big) M_\mu^\dag S_0
 \\= \left(\sum _\mu A_\mu \rho_s A_\mu^\dag \right) - \tfrac{1}{2}S_0^\dag L_1^\dag L_1 S_0\rho_s-  \rho_sS^\dag_0 L_1^\dag L_1 S_0\end{gathered}$$ with $A_\mu=S_0^\dag M_\mu L_1 S_0$. It remains to prove that $\sum_\mu A_\mu^\dag A_\mu= S_0^\dag L_1^\dag L_1 S_0$ for showing that we indeed have a Lindblad formulation. This results from the following computations: $$\begin{gathered}
  \sum_\mu A_\mu^\dag A_\mu= \sum_\mu S_0^\dag  L_1^\dag M_\mu^\dag S_0S_0^\dagger M_\mu L_1 S_0
  \\= S_0^\dag L_1^\dag R^*(S_0 S_0^\dag) L_1 S_0 = S_0^\dag L_1^\dag L_1 S_0,\end{gathered}$$ where we use that $S_0 S_0^\dag =P_0$ and $R^*(P_0)=I$.

Second order expansion for Hamiltonian perturbations {#sec:secondorder}
====================================================

We assume here that ${\mathfrak{L}}_0$ is defined by a single operator $L_0$, $
{\mathfrak{L}}_0(\rho)= L_0 \rho L_0^\dag - \tfrac{1}{2}\big( L_0^\dag L_0 \rho  + \rho L_0^\dag L_0 \big)
,
$ and that the perturbation ${\mathfrak{L}}_1$ is Hamiltonian, $
{\mathfrak{L}}_1(\rho)= - i |H_1,\rho],
$ where $H_1$ is a Hermitian operator. The following lemma gives a simple expression for $K_1(\rho_s)$ solution of .

\[lem:K1\] Assume that ${\mathfrak{L}}_0(\rho)= L_0 \rho L_0^\dag - \tfrac{1}{2}\big( L_0^\dag L_0 \rho  + \rho L_0^\dag L_0 \big) $ and ${\mathfrak{L}}_1(\rho)=- i[H_1,\rho]$. Then ${\mathfrak{L}}_{s,1}(\rho_s)= -i \big[S_0^\dag H_1 S_0~,~\rho_s\big]$ and $K_1(\rho_s)=- i\big[C_1, S_0 \rho_s S_0^\dag \big] $ satisfy  where $C_1$ is the Hermitian operator $$C_1= 2 (L_0^\dag L_0)^{-1} H_1 P_0 +  2  P_0 H_1 (L_0^\dag L_0)^{-1}$$ with $P_0$ the orthogonal projector onto $\mathcal{H_0}$ and $(L_0^\dag L_0)^{-1}$ standing for the Moore-Penrose pseudo-inverse of the Hermitian operator $L_0^\dag L_0$.

The associated first order $\rho_s$-parametrization of the slow invariant attractive manifold, $$\begin{gathered}
  K_0(\rho_s)+ \epsilon K_1(\rho_s)
=
\\
\left( I - i\epsilon (L_0^\dag L_0)^{-1} H_1\right) S_0 \rho_s S_0^\dag  \left( I + i\epsilon (L_0^\dag L_0)^{-1} H_1\right)
+ 0(\epsilon^2)
,\end{gathered}$$ corresponds, up to second-order terms, to a trace-preserving completely positive map.

With $S_0 {\mathfrak{L}}_{s,1}(\rho_s) S_0^\dag = -i\big[P_0 H_1 P_0,S_0\rho_s S_0^\dag\big]$, reads $$\begin{gathered}
  {\mathfrak{L}}_0(K_1(\rho_s))=-i\big[P_0 H_1 P_0,S_0\rho_s S_0^\dag\big] + i\big[H_1,S_0\rho_s S_0^\dag\big]
\\
= -i \left[ P_0 H_1 P_0 -H_1~, ~S_0\rho_s S_0^\dag\right]
.\end{gathered}$$ With $K_1(\rho_s)= -i\big[C_1, S_0 \rho_s S_0^\dag \big] $ we have also $$\begin{gathered}
 {\mathfrak{L}}_0(K_1(\rho_s)) = - i L_0 \left[ C_1, S_0\rho_s S_0^\dag\right] L_0^\dag
 \\
 + \tfrac{i}{2}\left( L_0^\dag L_0 \left[ C_1, S_0\rho_s S_0^\dag\right]  +  \left[ C_1, S_0\rho_s S_0^\dag\right]  L_0^\dag L_0 \right)
 .\end{gathered}$$ Since $L_0 S_0=0$ and $S_0^\dag L_0^\dag=0$ we have $$L_0 \left[ C_1, S_0\rho_s S_0^\dag\right] L_0^\dag =0.$$ Since additionally, $P_0 S_0=S_0$, $L_0^\dag L_0 P_0=0$ and $L_0^\dag L_0 (L_0^\dag L_0)^{-1} = I-P_0$, we have $$L_0^\dag L_0 \left[ C_1, S_0\rho_s S_0^\dag\right] = 2 (I-P_0) H_1 P_0 S_0 \rho_s S_0^\dag
.$$ Thus $$\begin{gathered}
 {\mathfrak{L}}_0(K_1(\rho_s)) =   i (I-P_0) H_1 P_0 S_0 \rho_s S_0^\dag -  i S_0 \rho_s S_0^\dag   P_0 H_1 (I-P_0)
  \\
   =
- i \left[ P_0 H_1 P_0 -H_1~, ~S_0\rho_s S_0^\dag\right]
 .\end{gathered}$$

The second order term ${\mathfrak{L}}_{s,2}(\rho_s)$ is solution of  for $m=2$: $${\mathfrak{L}}_{0}\left(  K_2(\rho_s) \right) + {\mathfrak{L}}_{1}\left(  K_{1}(\rho_s) \right)
=
K_{0}\left( {\mathfrak{L}}_{s,2}(\rho_s) \right) + K_{1}\left( {\mathfrak{L}}_{s,1}(\rho_s) \right)
.$$ Using, once again, $R\circ {\mathfrak{L}}_0\equiv 0$ and $R \circ K_0 = K_0$, we get $$\label{eq:Ls2}
{\mathfrak{L}}_{s,2}(\rho_s)
=
S_0^\dag  R\Big({\mathfrak{L}}_{1}\left(  K_{1}(\rho_s) \right) -  K_{1}\left( {\mathfrak{L}}_{s,1}(\rho_s) \right)\Big) S_0
.$$ The following lemma shows that ${\mathfrak{L}}_{s,2}(\rho_s)$ admits a Lindbladian form.

\[lem:Fs2\] The super-operator ${\mathfrak{L}}_{s,2}$ defined by  admits the following Lindbladian formulation $${\mathfrak{L}}_{s,2}(\rho_s)= \sum_\mu B_\mu \rho_s  B_\mu^\dag - \tfrac{1}{2} \big(B_\mu^\dag B_\mu \rho_s +\rho_s B_\mu^\dag B_\mu  \big)$$ with $B_\mu= 2 S_0^\dag M_\mu  L_0 (L_0^\dag L_0)^{-1} H_1 S_0$, $M_\mu$ defined by  and $(L_0^\dag L_0)^{-1}$ standing for the Moore-Penrose pseudo-inverse of $L_0^\dag L_0$.

We have $ R\Big(  K_{1}\left( {\mathfrak{L}}_{s,1}(\rho_s) \right)\Big) =0$. This results from (${K_0^s}$ stands for $S_0 \rho_s S_0^\dag = K_0(\rho_s)$) $$\begin{gathered}
  K_{1}\left( {\mathfrak{L}}_{s,1}(\rho_s) \right)
  = -i[C_1,\, -i\,S_0[S_0^\dag H_1 S_0,\,\rho_s]S_0^\dag] \\
  = -[C_1, P_0H_1 {K_0^s}- {K_0^s}H_1 P_0]
   \\
   \label{eq:Trick2}
  = -2 (L_0^\dag L_0)^{-1} H_1 (P_0H_1{K_0^s}- {K_0^s}H_1P_0) \\
    +2 (P_0H_1{K_0^s}- {K_0^s}H_1P_0)  H_1 (L_0^\dag L_0)^{-1}\end{gathered}$$ where we have used Lemma \[lem:K1\] and $P_0 K_0 = K_0$.

Repeating computations similar to , we see that for any operator $A$ on $\mathcal{H}$, $R(AP_0)=R(P_0A)=P_0 A P_0$. Since $P_0 {K_0^s}= {K_0^s}P_0={K_0^s}$ we moreover have $R(A{K_0^s})=P_0 A {K_0^s}$ and $R({K_0^s}A)={K_0^s}A P_0$. This gives the result of applying $R$ on all the terms in , and since $P_0(L_0^\dag L_0)^{-1} = (L_0^\dag L_0)^{-1}P_0=0$, we conclude that $R(K_1({\mathfrak{L}}_{s,1}))= 0$.

Thus ${\mathfrak{L}}_{s,2}(\rho_s) = S_0^\dag  R\Big({\mathfrak{L}}_{1}\left(  K_{1}(\rho_s) \right) \Big) S_0$. Exploiting similar simplifications, we have $$\begin{gathered}
  {\mathfrak{L}}_{1}\left(  K_{1}(\rho_s) \right) = -H_1(C_1{K_0^s}-{K_0^s}C_1) +(C_1{K_0^s}-{K_0^s}C_1) H_1
  \\
  = H_1 {K_0^s}C_1 + C_1 {K_0^s}H_1 - (H_1 C_1 {K_0^s}+{K_0^s}C_1 H_1)
  \\
  = 2H_1 {K_0^s}H_1(L_0^\dag L_0)^{-1} + 2 (L_0^\dag L_0)^{-1} H_1 {K_0^s}H_1
  \\ - 2 H_1 (L_0^\dag L_0)^{-1} H_1 {K_0^s}-2  {K_0^s}H_1 (L_0^\dag L_0)^{-1} H_1\end{gathered}$$ and, using $S_0^\dag R(A{K_0^s})=S_0^\dag P_0 A {K_0^s}= S_0^\dag A {K_0^s}$ and the definition ${K_0^s}= S_0 \rho_s s_0^\dag$, we get $$\begin{gathered}
 {\mathfrak{L}}_{s,2}(\rho_s)
  =
  \\2 S_0^\dag R \Big( H_1 {K_0^s}H_1(L_0^\dag L_0)^{-1} + (L_0^\dag L_0)^{-1} H_1 {K_0^s}H_1\Big) S_0
  \\ - 2 S^\dag_0H_1 (L_0^\dag L_0)^{-1} H_1 S_0 \rho_s -2 \rho_s S^\dag_0H_1 (L_0^\dag L_0)^{-1} H_1 S_0
  .\end{gathered}$$ Since for all $A$, $R({\mathfrak{L}}_0(A))=0$, we have the identity $$R(L_0 A L_0^\dag) = R \big( \tfrac{1}{2} \big(L_0^\dag L_0 A + A L_0^\dag L_0\big) \big)
.$$ With $A=(L_0^\dag L_0)^{-1}  H_1 {K_0^s}H_1 (L_0^\dag L_0)^{-1}$ we get $$\begin{gathered}
2 R\Big(L_0 (L_0^\dag L_0)^{-1}  H_1 {K_0^s}H_1 (L_0^\dag L_0)^{-1} L_0^\dag\Big)  =
  \\
 R\Big( (I-P_0) H_1 {K_0^s}H_1 (L_0^\dag L_0)^{-1} +  (L_0^\dag L_0)^{-1}  H_1 {K_0^s}H_1 (I-P_0)\Big)
 \\
  = R\Big(  H_1 {K_0^s}H_1 (L_0^\dag L_0)^{-1} +  (L_0^\dag L_0)^{-1}  H_1 {K_0^s}H_1 \Big)\end{gathered}$$ since $
R\Big( P_0 H_1 {K_0^s}H_1 (L_0^\dag L_0)^{-1} \Big)=P_0 H_1 {K_0^s}H_1 (L_0^\dag L_0)^{-1}P_0$ and $(L_0^\dag L_0)^{-1}P_0=0$. Thus $$\begin{gathered}
 {\mathfrak{L}}_{s,2}(\rho_s)
  =
  \\4 S_0^\dag  R\Big(L_0 (L_0^\dag L_0)^{-1}  H_1 {K_0^s}H_1 (L_0^\dag L_0)^{-1} L_0^\dag\Big) S_0
  \\ - 2 S^\dag_0 H_1 (L_0^\dag L_0)^{-1} H_1 S_0 \rho_s -2 \rho_s S^\dag_0H_1 (L_0^\dag L_0)^{-1} H_1 S_0
  .\end{gathered}$$ Using the decomposition  of $R$ we have $$4 S_0^\dag  R\Big(L_0 (L_0^\dag L_0)^{-1}  H_1 {K_0^s}H_1 (L_0^\dag L_0)^{-1} L_0^\dag\Big) S_0 \\
= \sum_\mu B_\mu \rho_s B_\mu^\dag \, .$$ We conclude by the following computations: $$\begin{gathered}
 \tfrac{1}{2}\sum_\mu B_\mu^\dag B_\mu=
 \\ 2 \sum_\mu  S_0^\dag H_1(L_0^\dag L_0)^{-1}L_0^\dag  M_\mu^\dag S_0 S_0^\dag M_\mu  L_0 (L_0^\dag L_0)^{-1} H_1 S_0
 \\
 = 2   S_0^\dag H_1(L_0^\dag L_0)^{-1}L_0^\dag  R^*(P_0)  L_0 (L_0^\dag L_0)^{-1} H_1 S_0
 \\
 = 2   S_0^\dag H_1(L_0^\dag L_0)^{-1}L_0^\dag   L_0 (L_0^\dag L_0)^{-1} H_1 S_0
 \\
 = 2   S_0^\dag H_1 (L_0^\dag L_0)^{-1} H_1 S_0
 .\end{gathered}$$

Illustrative example: low-Q cavity coupled to another quantum system {#sec:example}
====================================================================

The developments above are rigorous in finite dimension, but they can be formally applied also on infinite- dimensional systems, as illsutrated in the following example.

We consider a strongly dissipative driven harmonic oscillator (low-Q cavity) coupled to another, undamped quantum system with the same transition frequency (“target” system). Denote $\mathcal{H}_A$ (resp. $\mathcal{H}_B$) the infinite-dimensional Hilbert space of the strongly dissipative harmonic oscillator (resp. the target system), spanned by the Fock states $\{{\ensuremath{|n_A\rangle}\xspace}\}_{n\in \mathbb{N}}$ (resp. a possibly infinite basis $\{{\ensuremath{|n_B\rangle}\xspace}\}_{n_B}$); $\rho$ is the density operator of the composite system, on $\mathcal{H}=\mathcal{H}_A \otimes \mathcal{H}_B$.

In the frame rotating at the common frequency of the two systems, their coupled evolution is described by the standard master differential equation: $$\begin{aligned}
\label{eq:ExDyn1}
\begin{aligned}
\frac{d}{dt}\rho = [u\tilde{\boldsymbol{a}}^\dagger - u^*\tilde{\boldsymbol{a}}, \rho] + \kappa \left(\tilde{\boldsymbol{a}}\rho \tilde{\boldsymbol{a}}^\dagger-\frac{1}{2}\left(\tilde{\boldsymbol{a}}^\dagger \tilde{\boldsymbol{a}}\rho + \rho \tilde{\boldsymbol{a}}^\dagger \tilde{\boldsymbol{a}}\right) \right) \\ - ig\left[\tilde{\boldsymbol{a}}^\dagger \tilde{\boldsymbol{b}}+ \tilde{\boldsymbol{a}}\tilde{\boldsymbol{b}}^\dagger , \rho \right] \, .\\
\end{aligned}\end{aligned}$$ Here $\tilde{\boldsymbol{a}}={\boldsymbol{a}}\otimes I_B $ and $\tilde{\boldsymbol{b}}= I_A \otimes {\boldsymbol{b}}$ are the annihilation operators respectively for the harmonic oscillator $A$ and for the quantum system $B$ (possibly generalized if $B$ is not a harmonic oscillator; e.g. if $B$ is a qubit, we have ${\boldsymbol{b}}= {\ensuremath{|g\rangle}\xspace}{\ensuremath{\langle e|}\xspace}$ the transition operator from excited to ground state). The first line describes the driven and damped evolution of harmonic oscillator A, while the second line describes the exchange of energy quanta between the two quantum systems. The constants $(\kappa, g) \in \mathbb{R}^2$ govern the speed of these dynamics. We here consider $\kappa \gg g$, with the goal to adiabatically eliminate the fast dynamics of the low-Q cavity and compute its effect on the other quantum system. The dynamics is then equivalent to $$\begin{aligned}
\label{eq:mast_eq_ex}
\frac{d}{dt}\rho = \mathfrak{L}_0(\rho) + \epsilon\mathfrak{L}_1(\rho)\end{aligned}$$ with $L_0 = \sqrt{\kappa}(\tilde{\boldsymbol{a}}- \alpha)$, $\alpha = 2u/\kappa$ and $\epsilon\mathfrak{L}_1(\rho)=- ig\left[\tilde{\boldsymbol{a}}^\dagger \tilde{\boldsymbol{b}}+ \tilde{\boldsymbol{a}}\tilde{\boldsymbol{b}}^\dagger , \rho \right]$. For this typical example, the results of $\mathfrak{L}_{s,1}$ and $\mathfrak{L}_{s,2}$ are well known (see e.g. [@CarmichaelBook07 chap.12]). Our results allow to readily retrieve their expression and thus completely circumvent the trouble of the usual calculation.

In the absence of coupling between the two subsystems ($\epsilon=0$), the overall system trivially converges towards $R(\rho_0) = {\ensuremath{|\alpha\rangle}\xspace}{\ensuremath{\langle \alpha|}\xspace}_A \otimes \text{Tr}_A(\rho(0))$. Here $\text{Tr}_A$ is the partial trace over $\mathcal{H}_A$ and ${\ensuremath{|\alpha\rangle}\xspace}$ denotes the coherent state of amplitude $\alpha \in \mathbb{C}$, towards which a classically driven and damped harmonic oscillator is known to converge. Therefore we have $\mathcal{H}_0 = {\ensuremath{|\alpha\rangle}\xspace}{\ensuremath{\langle \alpha|}\xspace} \otimes \mathcal{H}_B$, $P_0 = {\ensuremath{|\alpha\rangle}\xspace}{\ensuremath{\langle \alpha|}\xspace} \otimes I_B$, and $M_\mu = {\ensuremath{|\alpha\rangle}\xspace}{\ensuremath{\langle \mu_A|}\xspace} \otimes I_B$ with $\mu$ spanning $\mathbb{N}$. We will naturally describe $\rho_S$ on the Hilbert space $\mathcal{H}_s \equiv \mathcal{H}_B$ and with basis $\{{\ensuremath{|n_s\rangle}\xspace}\}_{n_s}$, so $S_0 = \sum_{n}{\ensuremath{|\alpha\rangle}\xspace}{\ensuremath{|n_B\rangle}\xspace}{\ensuremath{\langle n_s|}\xspace}$.

For the first-order perturbation, using the property $\tilde{\boldsymbol{a}}{\ensuremath{|\alpha\rangle}\xspace} = \alpha {\ensuremath{|\alpha\rangle}\xspace}$, Lemma \[lem:Ls1\] readily yields $$H_{s,1} = \alpha {\boldsymbol{b}}_s^\dagger + \alpha^* {\boldsymbol{b}}_s \; ,$$ denoting by $\boldsymbol{q}_s$ the operator on $\mathcal{H}_s$ equivalent to $\boldsymbol{q}$ on $\mathcal{H}_B$. This standard result shows that the oscillator $A$ can be approximated as a classical field of amplitude $\alpha$. Indeed, $H_{s,1}$ describes e.g. Rabi oscillations for a qubit driven by a classical field ($\mathcal{H}_B=\text{span}\{{\ensuremath{|g\rangle}\xspace},{\ensuremath{|e\rangle}\xspace}\}$); or, when $\mathcal{H}_B$ describes another harmonic oscillator, $H_{s,1}$ is the same Hamiltonian in fact as in the first line of , with classical drive amplitude $i u$ replaced by $\alpha$.

Next, using ${\mathfrak{D}_{\alpha}}$ the unitary displacement operator on $\mathcal{H}_A$, which satisfies ${\mathfrak{D}_{\alpha}}{\boldsymbol{a}}{\mathfrak{D}_{\alpha}^\dagger}= {\boldsymbol{a}}-\alpha I$, we compute $(L_0^\dag L_0)^{-1} = {\mathfrak{D}_{\alpha}}{{\boldsymbol{N}}}_A^{-1} {\mathfrak{D}_{\alpha}^\dagger}/\kappa$, where ${{\boldsymbol{N}}}_A = {\boldsymbol{a}}^\dag {\boldsymbol{a}}= \sum_{n \in \mathbb{N}} \; n \, {\ensuremath{|n_A\rangle}\xspace}{\ensuremath{\langle n_A|}\xspace}$ and the Moore-Penrose pseudo-inverse of ${{\boldsymbol{N}}}_A$ is just ${{\boldsymbol{N}}}_A^{-1}=\sum_{n \geq 1} \; \frac{1}{n} \, {\ensuremath{|n_A\rangle}\xspace}{\ensuremath{\langle n_A|}\xspace}$. We then compute $$\begin{gathered}
C_1 = \frac{2}{\kappa}\, {\mathfrak{D}_{\alpha}}{{\boldsymbol{N}}}_A^{-1} {\mathfrak{D}_{\alpha}^\dagger}(\tilde{\boldsymbol{a}}^\dag \tilde{\boldsymbol{b}}+ \tilde{\boldsymbol{a}}\tilde{\boldsymbol{b}}^\dag) {\ensuremath{|\alpha\rangle}\xspace}{\ensuremath{\langle \alpha|}\xspace}\otimes I_B \;\; + h.c. \\
= \frac{2}{\kappa}\, {\mathfrak{D}_{\alpha}}{{\boldsymbol{N}}}_A^{-1} ((\tilde{\boldsymbol{a}}^\dag+\alpha^* I) \tilde{\boldsymbol{b}}+ (\tilde{\boldsymbol{a}}+\alpha I) \tilde{\boldsymbol{b}}^\dag) {\mathfrak{D}_{\alpha}^\dagger}{\ensuremath{|\alpha\rangle}\xspace}{\ensuremath{\langle \alpha|}\xspace}\otimes I_B \;\; + h.c. \\
= \frac{2}{\kappa}\, {\mathfrak{D}_{\alpha}}{{\boldsymbol{N}}}_A^{-1} ((\tilde{\boldsymbol{a}}^\dag+\alpha^* I) \tilde{\boldsymbol{b}}+ (\tilde{\boldsymbol{a}}+\alpha I) \tilde{\boldsymbol{b}}^\dag) {\ensuremath{|0\rangle}\xspace}{\ensuremath{\langle \alpha|}\xspace}\otimes I_B \;\; + h.c. \\
= \frac{2}{\kappa}\, {\mathfrak{D}_{\alpha}}{{\boldsymbol{N}}}_A^{-1} ((\alpha^*\tilde{\boldsymbol{b}}+ \alpha \tilde{\boldsymbol{b}}^\dag){\ensuremath{|0\rangle}\xspace} + \tilde{\boldsymbol{b}}{\ensuremath{|1\rangle}\xspace}){\ensuremath{\langle \alpha|}\xspace} \;\; + h.c. \\
= \frac{2}{\kappa}\, {\mathfrak{D}_{\alpha}}{\ensuremath{|1\rangle}\xspace}{\ensuremath{\langle \alpha|}\xspace}\otimes {\boldsymbol{b}}\;\; + h.c. \, .\end{gathered}$$ From Lemma \[lem:K1\], we see that a pure state ${\ensuremath{|\psi_S\rangle}\xspace} \in \mathcal{H}_s$ gets mapped at order zero to ${\ensuremath{|\alpha\rangle}\xspace} \otimes {\ensuremath{|\psi_B\rangle}\xspace}$ with ${\ensuremath{|\psi_B\rangle}\xspace} \equiv {\ensuremath{|\psi_S\rangle}\xspace}$, but at order one to a slightly rotated state ${\ensuremath{|\alpha\rangle}\xspace} \otimes {\ensuremath{|\psi_B\rangle}\xspace} - \frac{i g}{\kappa} ({\mathfrak{D}_{\alpha}}{\ensuremath{|1\rangle}\xspace}) \otimes ({\boldsymbol{b}}{\ensuremath{|\psi_B\rangle}\xspace})$. This expresses that the coupled low-Q cavity $A$ contains slightly more energy than a coherent state, to the detriment of system B.

For the second order perturbation, from Lemma \[lem:Fs2\] we must compute $B_\mu= 2 S_0^\dag M_\mu  L_0 (L_0^\dag L_0)^{-1} H_1 S_0$. The computations made for $C_1$ above can be used, writing: $$\begin{gathered}
B_\mu = S_0^\dag M_\mu  L_0 \; \big( \frac{2}{\kappa}\, {\mathfrak{D}_{\alpha}}{\ensuremath{|1\rangle}\xspace}{\ensuremath{\langle \alpha|}\xspace}\otimes {\boldsymbol{b}}\big) S_0 \\
= \frac{2}{\sqrt{\kappa}} S_0^\dag M_\mu \big( {\mathfrak{D}_{\alpha}}{\boldsymbol{a}}{\ensuremath{|1\rangle}\xspace}{\ensuremath{\langle \alpha|}\xspace}\otimes {\boldsymbol{b}}\big) \, S_0 \\
= \frac{2}{\sqrt{\kappa}} \sum_{n,m} {\ensuremath{|n_s\rangle}\xspace}{\ensuremath{\langle n_B|}\xspace} {\ensuremath{\langle \mu_A|}\xspace} \; {\ensuremath{|\alpha\rangle}\xspace}{\ensuremath{\langle \alpha|}\xspace}\otimes {\boldsymbol{b}}\; {\ensuremath{|m_B\rangle}\xspace}{\ensuremath{\langle m_S|}\xspace}\\
= \frac{2}{\sqrt{\kappa}}\; {\langle \mu_A|\alpha \rangle}{\boldsymbol{b}}_s \, .\end{gathered}$$ All the obtained $B_\mu$ are in fact identical up to a scalar factor, so they may be combined into a single operator: $$\begin{gathered}
\epsilon^2  {\mathfrak{L}}_{s,2}(\rho_s)=g^2 \sum_\mu B_\mu \rho_s  B_\mu^\dag - \tfrac{1}{2} \big(B_\mu^\dag B_\mu \rho_s +\rho_s B_\mu^\dag B_\mu  \big) \\
= \frac{4g^2}{\kappa} \sum_\mu |{\langle \mu_A|\alpha \rangle}|^2 \left( {\boldsymbol{b}}_s\rho_s {\boldsymbol{b}}_s^\dagger - \frac{1}{2}\left({\boldsymbol{b}}_s^\dag {\boldsymbol{b}}_s \rho_s + \rho_s {\boldsymbol{b}}_s^\dagger {\boldsymbol{b}}_s \right) \right) \\
= \frac{4g^2}{\kappa}\left( b\rho_s b^\dagger - \frac{1}{2}\left(b^\dagger b \rho_s + \rho_s b^\dagger b \right)\right) \, .\end{gathered}$$ (Note that $\{ |{\langle \mu_A|\alpha \rangle}|^2 \}_{\mu \in \mathbb{N}}$ just corresponds to the expansion of the coherent state ${\ensuremath{|\alpha\rangle}\xspace}$, of unit norm, in the Fock basis.) We thus get the expected reduced dynamics: $$\begin{gathered}
{\frac{d}{dt}}\rho_s = -ig\left[\alpha {\boldsymbol{b}}_s^\dag + \alpha^*{\boldsymbol{b}}_s, \rho_s \right] \\+ \frac{4g^2}{\kappa}\left( {\boldsymbol{b}}_s\rho_s {\boldsymbol{b}}_s^\dagger - \frac{1}{2}\left({\boldsymbol{b}}_s^\dagger {\boldsymbol{b}}_s \rho_s + \rho_s {\boldsymbol{b}}_s^\dagger {\boldsymbol{b}}_s \right)\right) \; ,\end{gathered}$$ which expresses that the B system is subject to slow damping due to the presence of the low-Q cavity.

#### Remark {#remark .unnumbered}

Note that if the slow dynamics includes a Hamiltonian that acts only on the $B$ system, i.e. of the form $\tilde{\boldsymbol{H}}_B=I_A \otimes {\boldsymbol{H}}_B$ (acting only on B), then $C_1$ features an additional term $$\begin{gathered}
\frac{2}{\kappa}\, {\mathfrak{D}_{\alpha}}{{\boldsymbol{N}}}_A^{-1} {\mathfrak{D}_{\alpha}^\dagger}\big(I_A \otimes {\boldsymbol{H}}_B\big) {\ensuremath{|\alpha\rangle}\xspace}{\ensuremath{\langle \alpha|}\xspace}\otimes I_B \;\; + h.c.\\
= \frac{2}{\kappa}\, \big({\mathfrak{D}_{\alpha}}{{\boldsymbol{N}}}_A^{-1} {\mathfrak{D}_{\alpha}^\dagger}{\ensuremath{|\alpha\rangle}\xspace}{\ensuremath{\langle \alpha|}\xspace} \big) \otimes {\boldsymbol{H}}_B \\
= \frac{2}{\kappa}\, \big({\mathfrak{D}_{\alpha}}{{\boldsymbol{N}}}_A^{-1} {\ensuremath{|0\rangle}\xspace}{\ensuremath{\langle \alpha|}\xspace} \big) \otimes {\boldsymbol{H}}_B = 0 \, .\end{gathered}$$ Thus the second-order correction vanishes and the Zeno dynamics is the only addition up to second order: $$\begin{gathered}
{\frac{d}{dt}}\rho_s = -ig\left[\alpha {\boldsymbol{b}}_s^\dag + \alpha^*{\boldsymbol{b}}_s + {\boldsymbol{H}}_B, \rho_s \right] \\+ \frac{4g^2}{\kappa}\left( {\boldsymbol{b}}_s\rho_s {\boldsymbol{b}}_s^\dagger - \frac{1}{2}\left({\boldsymbol{b}}_s^\dagger {\boldsymbol{b}}_s \rho_s + \rho_s {\boldsymbol{b}}_s^\dagger {\boldsymbol{b}}_s \right)\right) \; .\end{gathered}$$

Conclusion
==========

We have shown how to eliminate the fast dynamics in an open quantum system (Lindblad equation) with two timescales, and obtain the resulting reduced dynamics explicitly in Lindblad form. This is important to guarantee that the approximate model preserves the structure of quantum states (positivity and trace). The slow system is hence parameterized explicitly with a quantum state on a lower-dimensional Hilbert space, and mapped to the complete Hilbert space by a completely positive trace preserving map (Kraus map). We have illustrated on a benchmark system (highly dissipating quantum oscillator resonantly coupled to another quantum system) how our explicit formulae directly retrieve the results previously obtained with lengthy ad hoc computations.

We have obtained explicit formulae for the second-order corrections only in the particular case of a fast Lindbladian with single-channel damping $L_0$, and a slow “perturbation” in Hamiltonian form. Conceptually there should be no obstacle to extending this theory to any Lindbladians, the key point being an appropriate way to generalize the pseudo-inversion $(L_0^\dag L_0)^{-1}$. However, the special case completed here will already allow to answer currently open questions about the influence of small Hamiltonian perturbations on stable open quantum systems built e.g. with engineered reservoirs [@MirrahimiCatComp2014; @LeghtTPKVPSNSHRFSMD2015S].

Acknowledgement {#acknowledgement .unnumbered}
===============

The authors thank Benjamin Huard, Zaki Leghtas and Mazyar Mirrahimi for many useful discussions.

[10]{}

D. J. Atkins, H. M. Wiseman, and P. Warszawski. Approximate master equations for atom optics. , 67:023802, Feb 2003.

R. Azouit, A. Sarlette, and P. Rouchon. Convergence and adiabatic elimination for a driven dissipative quantum harmonic oscillator. In [*2015 54th IEEE Conference on Decision and Control (CDC)*]{}, pages 6447–6453, Dec 2015.

H.-P. Breuer and F. Petruccione. . , 2006.

E Brion, L H Pedersen, and K M[ø]{}lmer. Adiabatic elimination in a lambda system. , 40(5):1033, 2007.

H. Carmichael. . Spinger, 2007.

J. Carr. . , [1981]{}.

O. Cernotik, D.V. Vasilyev, and K. Hammerer. Adiabatic elimination of gaussian subsystems from quantum dynamics under continuous measurement. , 92(1):012124–, July 2015.

M.D. Choi. Completely positive linear maps on complex matrices. , 10:285–290, 1975.

N. Fenichel. Geometric singular perturbation theory for ordinary differential equations. , 31:53–98, 1979.

S. Haroche and J.M. Raimond. Oxford University Press, 2006.

E. M. Kessler. Generalized [Schrieffer-Wolff]{} formalism for dissipative systems. , 86(1):012126–, July 2012.

Z. Leghtas, S. Touzard, I. M. Pop, A. Kou, B. Vlastakis, A. Petrenko, K. M. Sliwa, A. Narla, S. Shankar, M. J. Hatridge, M. Reagor, L. Frunzio, R. J. Schoelkopf, M. Mirrahimi, and M. H. Devoret. Confining the state of light to a quantum manifold by engineered two-photon loss. , 347(6224):853–857, February 2015.

M. Mirrahimi, Z. Leghtas, V.V. Albert, S. Touzard, R.J. Schoelkopf, L. Jiang, and M.H. Devoret. Dynamically protected cat-qubits: a new paradigm for universal quantum computation. , 16:045014, 2014.

M. Mirrahimi and P. Rouchon. Singular perturbations and [Lindblad-Kossakowski]{} differential equations. , 54(6):1325–1329, 2009.

M.A. Nielsen and I.L. Chuang. . Cambridge University Press, 2000.

F. Reiter and A.S. S[ø]{}rensen. Effective operator formalism for open quantum systems. , 85(3):032111–, March 2012.

J.J. Sakurai and J. Napolitano. . Addison-Wesley, 2011.

[^1]: Centre Automatique et Systèmes, Mines-ParisTech, PSL Research University. 60 Bd Saint-Michel, 75006 Paris, France.

[^2]: INRIA Paris, 2 rue Simone Iff, 75012 Paris, France; and Ghent University / Data Science Lab, Technologiepark 914, 9052 Zwijnaarde, Belgium.

[^3]: This work was partially supported by the Projet Blanc ANR-2011-BS01-017-01 EMAQS.
