Towards generic adiabatic elimination for bipartite open quantum systems

TL;DR

This paper develops a systematic adiabatic elimination method for bipartite open quantum systems, allowing derivation of reduced models that accurately describe slow subsystem dynamics with explicit second-order formulas.

Contribution

It introduces a general adiabatic elimination technique for bipartite open quantum systems using asymptotic expansion and geometric singular perturbation theory, ensuring physical consistency.

Findings

Derived explicit second-order formulas for reduced dynamics.

Ensured the reduced model is in Lindblad form.

Provided formulas for Hamiltonian and cascade couplings.

Abstract

We consider a composite open quantum system consisting of a fast subsystem coupled to a slow one. Using the time-scale separation, we develop an adiabatic elimination technique to derive at any order the reduced model describing the slow subsystem. The method, based on an asymptotic expansion and geometric singular perturbation theory, ensures the physical interpretation of the reduced second-order model by giving the reduced dynamics in a Lindblad form and the state reduction in Kraus map form. We give explicit second-order formulas for Hamiltonian or cascade coupling between the two subsystems. These formulas can be used to engineer, via a careful choice of the fast subsystem, the Hamiltonian and Lindbald operators governing the dissipative dynamics of the slow subsystem.

Equations

$$\frac{d\rho }{dt}=u[{\mathbf{\sigma }}_{\text{+}}^{\mathbf{A}}-{\mathbf{\sigma }}_{\text{-}}^{\mathbf{A}},\rho ]+\gamma ({\mathbf{\sigma }}_{\text{-}}^{\mathbf{A}}\rho {\mathbf{\sigma }}_{\text{+}}^{\mathbf{A}}-\frac{{\mathbf{\sigma }}_{\text{+}}^{\mathbf{A}}{\mathbf{\sigma }}_{\text{-}}^{\mathbf{A}}\rho +\rho {\mathbf{\sigma }}_{\text{+}}^{\mathbf{A}}{\mathbf{\sigma }}_{\text{-}}^{\mathbf{A}}}{2})-i\chi [{\mathbf{\sigma }}_{\mathbf{z}}^{\mathbf{A}}\otimes {\mathbf{\sigma }}_{\mathbf{z}}^{\mathbf{B}},\rho ]$$

$$\displaystyle\frac{d\rho}{dt}=i\chi\frac{\gamma^{2}}{\gamma^{2}+8u^{2}}[\boldsymbol{\sigma_{\!z}^{B}},\rho]+64\gamma\chi^{2}u^{2}\,\frac{(\gamma^{2}+2u^{2})}{(\gamma^{2}+8u^{2})^{3}}\;\big{(}\boldsymbol{\sigma_{\!z}^{B}}\rho\boldsymbol{\sigma_{\!z}^{B}}-\rho\big{)}\,.$$

$$\frac{d}{dt}\rho ={\mathcal{L}}_A(\rho )+ϵ{\mathcal{L}}_{\mathrm{int}}(\rho )+ϵ{\mathcal{L}}_B(\rho )$$

$${\mathcal{L}}_{\xi }(\rho )=-i[{\mathbf{H}}_{\xi },\rho ]+\sum _{\mu }{\mathcal{D}}_{{\mathbf{L}}_{\xi ,\mu }}(\rho ),\xi =A,B$$

$$\mathcal{D}_{\boldsymbol{X}}(\rho)=\boldsymbol{X}\rho\boldsymbol{X}^{\dagger}-\frac{1}{2}\big{(}\boldsymbol{X}^{\dagger}\boldsymbol{X}\rho+\rho\boldsymbol{X}^{\dagger}\boldsymbol{X}\big{)};$$

$$\rho$$

$$\frac{d}{dt}{\rho }_s$$

$${\mathcal{K}}_0({\rho }_s)={\overline{\rho }}_A\otimes {\rho }_s\text{ and }{\mathcal{L}}_0({\rho }_s)=0.$$

$${\mathcal{L}}_A(\mathcal{K}({\rho }_s))+ϵ{\mathcal{L}}_{\mathrm{int}}(\mathcal{K}({\rho }_s))+ϵ{\mathcal{L}}_B(\mathcal{K}({\rho }_s))=\frac{d}{dt}\rho =\mathcal{K}({\mathcal{L}}_s({\rho }_s))$$

$$\begin{gathered} {\mathcal{L}}_A\left(ϵ{\mathcal{K}}_1({\rho }_s)+{ϵ}^2{\mathcal{K}}_2({\rho }_s)+\dots \right)+ϵ{\mathcal{L}}_{\mathrm{int}}\left({\mathcal{K}}_0({\rho }_s)+ϵ{\mathcal{K}}_1({\rho }_s)+{ϵ}^2{\mathcal{K}}_2({\rho }_s)+\dots \right) \\ +ϵ{\mathcal{L}}_B\left({\mathcal{K}}_0({\rho }_s)+ϵ{\mathcal{K}}_1({\rho }_s)+{ϵ}^2{\mathcal{K}}_2({\rho }_s)+\dots \right)={\mathcal{K}}_0\left(ϵ{\mathcal{L}}_1({\rho }_s)+{ϵ}^2{\mathcal{L}}_2({\rho }_s)+\dots \right) \\ +ϵ{\mathcal{K}}_1\left(ϵ{\mathcal{L}}_1({\rho }_s)+{ϵ}^2{\mathcal{L}}_2({\rho }_s)+\dots \right)+{ϵ}^2{\mathcal{K}}_2\left(ϵ{\mathcal{L}}_1({\rho }_s)+{ϵ}^2{\mathcal{L}}_2({\rho }_s)+\dots \right)+\dots \end{gathered}$$

$${\mathcal{L}}_A({\mathcal{K}}_1({\rho }_s))+{\mathcal{L}}_{\mathrm{int}}({\mathcal{K}}_0({\rho }_s))+{\mathcal{L}}_B({\mathcal{K}}_0({\rho }_s))$$

$${\mathcal{L}}_A({\mathcal{K}}_2({\rho }_s))+{\mathcal{L}}_{\mathrm{int}}({\mathcal{K}}_1({\rho }_s))+{\mathcal{L}}_B({\mathcal{K}}_1({\rho }_s))$$

$${\mathcal{L}}_A({\mathcal{K}}_1({\rho }_s))+{\mathcal{L}}_{\mathrm{int}}({\overline{\rho }}_A\otimes {\rho }_s)+{\overline{\rho }}_A\otimes {\mathcal{L}}_B({\rho }_s)={\overline{\rho }}_A\otimes {\mathcal{L}}_1({\rho }_s).$$

$$\mathcal{L}_{1}(\rho_{s})={\rm{tr}}_{A}\Big{(}\mathcal{L}_{\!\rm{int}}(\overline{\rho}_{A}\otimes\rho_{s})\Big{)}+\mathcal{L}_{B}(\rho_{s})\;.$$

$$\mathbf{X}=\overline{\tau }{\overline{\mathcal{K}}}_A(\mathbf{W})+\lambda {\overline{\rho }}_A,\forall \lambda \in \mathbb{C}.$$

$${\mathcal{K}}_1({\rho }_s)$$

$$\begin{gathered} \mathcal{K}({\rho }_s)={\mathcal{K}}_0({\rho }_s)+ϵ{\mathcal{K}}_1({\rho }_s)+O({ϵ}^2) \\ ={\overline{\mathcal{K}}}_A\left({\overline{\rho }}_A\otimes {\rho }_s+ϵ\overline{\tau }{\mathcal{L}}_{\mathrm{int}}({\overline{\rho }}_A\otimes {\rho }_s)\right)+O({ϵ}^2) \\ ={\overline{\mathcal{K}}}_A\left(e^{ϵ\overline{\tau }{\mathcal{L}}_{\mathrm{int}}}({\mathcal{K}}_0({\rho }_s))\right)+O({ϵ}^2). \end{gathered}$$

$$\mathcal{K}(\rho_{s})=\overline{\mathcal{K}}_{A}\circ e^{\epsilon\overline{\tau}\mathcal{L}_{\!\rm{int}}}\circ\mathcal{K}_{0}\circ e^{-\epsilon\overline{\tau}{{\rm{tr}}_{A}\Big{(}\mathcal{L}_{\!\rm{int}}(\overline{\rho}_{A}\otimes\rho_{s})\Big{)}}}(\rho_{s})+O(\epsilon^{2}).$$

$$\frac{d}{dt}\rho=\mathcal{L}_{A}(\rho)+\epsilon\Big{(}-i[\boldsymbol{H}_{\!\rm{int}},\rho]+\mathcal{L}_{B}(\rho)\Big{)}$$

$${\mathcal{L}}_1({\rho }_s)$$

$$\mathcal{K}(\rho_{s})=\mathcal{K}_{0}(\rho_{s})+\epsilon\mathcal{K}_{1}(\rho_{s})+O(\epsilon^{2})=\\ \exp\big{(}\text{-}i\epsilon\boldsymbol{M}\big{)}\;(\overline{\rho}_{A}\otimes\rho_{s})\;\exp\big{(}i\epsilon\boldsymbol{M}^{\dagger}\big{)}\;\;+O(\epsilon^{2})\quad\text{ with }\boldsymbol{M}=\sum_{k=1}^{m}\boldsymbol{F}_{k}\otimes\boldsymbol{B}^{\dagger}_{k}\;,$$

$${\mathcal{L}}_2({\rho }_s)$$

$${\mathcal{L}}_2({\rho }_s)=\mathrm{tr}(\mathbf{F}{\overline{\rho }}_A\mathbf{A}+\mathbf{A}{\overline{\rho }}_A{\mathbf{F}}^{†})\mathbf{B}{\rho }_s\mathbf{B}-\mathrm{tr}(\mathbf{A}{\overline{\rho }}_A{\mathbf{F}}^{†}){\rho }_s\mathbf{B}\mathbf{B}-\mathrm{tr}(\mathbf{F}{\overline{\rho }}_A\mathbf{A})\mathbf{B}\mathbf{B}{\rho }_s.$$

$${\mathcal{L}}_2({\rho }_s)=-i\left[\frac{\mathrm{tr}(\mathbf{F}{\overline{\rho }}_A\mathbf{A}-\mathbf{A}{\overline{\rho }}_A{\mathbf{F}}^{†})}{2i}(\mathbf{B}{)}^2,{\rho }_s\right]+\mathrm{tr}(\mathbf{F}{\overline{\rho }}_A\mathbf{A}+\mathbf{A}{\overline{\rho }}_A{\mathbf{F}}^{†}){\mathcal{D}}_{\mathbf{B}}({\rho }_s).$$

$$\frac{d}{dt}{\rho }_s$$

$${\mathcal{L}}_2({\rho }_s)$$

$$X_{k,j}=\mathrm{tr}\left({\mathbf{F}}_j{\overline{\rho }}_A{\mathbf{A}}_k^{†}+{\mathbf{A}}_j{\overline{\rho }}_A{\mathbf{F}}_k^{†}\right)$$

$${\mathcal{L}}_2({\rho }_s)=-i\left[\sum _{k,j}{Y}_{k,j}{\mathbf{B}}_k{\mathbf{B}}_j^{†},{\rho }_s\right]+\sum _{k=1}^m{\mathcal{D}}_{{\mathbf{L}}_k}({\rho }_s)\text{ where }{\mathbf{L}}_k=\sum _{j=1}^m{\Lambda }_{j,k}^{\ast }{\mathbf{B}}_j^{†},X=\Lambda {\Lambda }^{†}.$$

$$\frac{d}{dt}{\rho }_s=-i\left[ϵ\sum _k\mathrm{tr}({\mathbf{A}}_k{\overline{\rho }}_A){\mathbf{B}}_k+{ϵ}^2\sum _{k,j}{Y}_{k,j}{\mathbf{B}}_k{\mathbf{B}}_j^{†},{\rho }_s\right]+ϵ{\mathcal{L}}_B({\rho }_s)+{ϵ}^2\sum _{k=1}^m{\mathcal{D}}_{{\mathbf{L}}_k}({\rho }_s).$$

$$\frac{d}{dt}{\rho }_s$$

$${\mathcal{L}}_A=-iu[{\mathbf{\sigma }}_{\mathbf{y}},\cdot ]+{\kappa }_{-}{\mathcal{D}}_{{\mathbf{\sigma }}_{\text{-}}}+{\kappa }_{+}{\mathcal{D}}_{{\mathbf{\sigma }}_{\text{+}}},$$

$$\frac{{\kappa }_{+}}{{\kappa }_{-}}\ll 1,\frac{|u|}{{\kappa }_{-}}\ll 1.$$

Full text

1 Introduction

As any quantum system interacts with an environment, the rigorous way to describe its evolution would be a Schrödinger equation including the environment. Studying this whole system is usually not possible due to the complexity of the environment. Therefore, using typical assumptions on the environment one can get rid of it by some Markov approximations and obtain a Lindblad master equation [6] governing only the system of interest. A similar situation arises naturally within a quantum system composed of several interacting subsystems: one may want to get the dynamics of a particular subsystem of interest by using some assumptions to get rid of the other ones. When the subsystem of interest has much slower dissipation rates (i.e. time scales) than the other ones and is weakly coupled to the other ones, such procedure is known as adiabatic elimination.

This kind of setting is related to reservoir engineering [27, 10, 22] where the goal is to design the dynamics of a target subsystem by properly choosing its interaction with the other subsystems. The interaction is engineered in such a way that the dynamics of the target subsystem is, after adiabatic elimination of the other subsystems, the desired one. Rigorously computing the reduced dynamics associated with the subsystem of interest is not straightforward. We here propose a systematic adiabatic elimination procedure that preserves the quantum structure. As illustrated by the following example, the resulting reduced dynamics may include some non-intuitive phenomena.

Consider a pair of qubits. The first qubit, labelled by $A$, is subject to strong energy loss (life time of the excited state $1/\gamma$ ) and is driven by a coherent resonant field of amplitude $u\geq 0$. The second qubit, labelled by $B$, is isolated from the environment but is weakly coupled to the first one by a dispersive Hamiltonian of strength $\chi$. This setting is an abstraction of a target qubit $B$ coupled to a spurious “TLS”, i.e. a strongly decoherent two-level system with $\gamma\gg|\chi|$, see e.g. [23]. The density operator $\rho$ of the joint system is governed by the following master equation:

[TABLE]

where $\boldsymbol{\sigma_{\!\text{\bf-}}^{A}}$ is the lowering operator for qubit $A$, $\boldsymbol{\sigma_{\!\text{\bf+}}^{A}}=(\boldsymbol{\sigma_{\!\text{\bf-}}^{A}})^{\dagger}$, and $\boldsymbol{\sigma_{\!z}^{A}}$, $\boldsymbol{\sigma_{\!z}^{B}}$ are the usual Pauli operators for qubits $A$ and $B$ respectively. Since $\gamma\gg|\chi|$, the effect of the “TLS” qubit $A$ on the target qubit $B$ can be obtained by adiabatically eliminating qubit $A$. A direct application of the formula (17) that we obtain below, shows that up to third order terms in $\epsilon=|\chi|/\gamma$, the reduced dynamics of qubit $B$ are:

[TABLE]

As one may intuitively expect, qubit $A$ induces on qubit $B$ a deterministic phase shift, corresponding to the mean value of $\boldsymbol{\sigma_{\!z}^{A}}$, and some phase loss related to the uncertainty in $\boldsymbol{\sigma_{\!z}^{A}}$. When $u$ is small, as expected, qubit $B$ suffers little phase loss as qubit $A$ is close to the ground state, which is an eigenstate of $\boldsymbol{\sigma_{\!z}^{A}}$. However, more surprisingly, in the presence of strong drive (i.e large $u$) the phase loss also vanishes. This is due to the fact that the system is then subject to a sort of dynamical decoupling [31].

The present paper thus provides a systematic method for deriving the Lindblad form of reduced dynamics like (1), up to third order terms, for general bipartite systems with two times-scales. For a Hamiltonian coupling between the two subsystems, this second order Lindblad form is given by the generic formula (19). For a cascade coupling, it is given by (38).

Throughout the literature, adiabatic elimination formulas for open quantum systems have been obtained for several particular examples separately. A specific atom-optics system is analyzed in [1]; for lambda systems see e.g. [7], [25], [28]; systems with Gaussian dynamics are investigated in [11]. In contrast, generic techniques for adiabatic elimination have attracted much less attention. Treating the Lindblad master as a standard linear system, Kessler has developed in [21] a generalization of the Schrieffer-Wolff formalism. In [18, 5], the authors derive a reduced dynamics for systems described by Quantum Stochastic Differential Equations in the limit where the speed of the fast system goes to infinity. Our aim is to go further and precisely characterise the order of the approximation. This is motivated by the fact that in some recent quantum experiments aiming via reservoir engineering at strong indirect stabilization of quantum systems [22], a better knowledge of the order of validity of such approximations is becoming necessary : the increase in the accuracy of experiments implies the need to compute higher order model in order to properly describe the slow dynamics. Our aim with the present paper is to provide formulas with improved direct applicability, control on the order of approximation, and explicit expression of how the reduced system is embedded in the full Hilbert space.

We focus on bipartite open quantum systems. Using a geometric approach, we treat the slow dynamics as a perturbation of the fast one. This leads to an asymptotic expansion in powers of the small parameter describing the time-scales separation and readily allows one to select the order of the approximation. Moreover, the specificity of our approach compared to standard abstract systems theory, is to preserve the structural properties that ensure the physical interpretation of the reduced model: first, the reduced density operator is governed by a Lindblad master equation; and second, the mapping from the reduced density operator to the density operator of the entire system is a Kraus map [12] i.e. completely positive trace preserving. These two features ensure the preservation of a density operator interpretation. The computation of the Kraus map allows to estimate for instance the residual entanglement, inside the protected subspace, between the two subsystems of a reservoir engineering setup, which may have important practical consequences. To our knowledge, such care about the mapping between the complete and reduced systems is new in adiabatic elimination techniques.

Our approach builds on the framework of center manifold theory [9] and geometric singular perturbations theory [15] to obtain recurrence relations between the approximations at different orders. At the first order, we readily retrieve the so-called Zeno Hamiltonian for any type of coupling between the fast and slow subsystems. Then, focusing on two standard forms of coupling namely Hamiltonian interaction and cascade interaction, we derive formulas for the second order approximations while ensuring a Lindblad form for the reduced dynamics and a mapping in Kraus map form. The formulas show that the Hamiltonian interaction with a decohering subsystem leads to decoherence at the second order. Before performing any detailed calculations, we obtain the general structural result that the number of decoherence channels is equal to the minimal number of tensor-product terms required to express the interaction Hamiltonian.

The paper is organized as follows. In section 2 we present the considered system types, our general adiabatic elimination approach, and the results on first order approximation. Section 3 is devoted to the computation of the second order approximation in a Lindblad form for a Hamiltonian coupling between the fast and slow subsystems. The typical cases of dispersive or resonant interaction are mentioned explicitly (for readers interested mainly in the final formulas: see equations (17) and (20)). In Section 3.2, the use of these formulas is illustrated by obtaining new explicit expressions for the residual imperfections on a reservoir engineering setup with two-photon pumping [13, 30]. Section 4 investigates a different kind of coupling, namely the cascade interaction. We derive the second order approximation in a Lindblad form provided some technical condition holds.

Acknowledgements:

The authors want to thank Zaki Leghtas and Mazyar Mirrahimi for useful discussions about this subject.

Notations:

The underlying Hilbert spaces are denoted by $\mathcal{H}_{\bullet}$ with different subscripts $\bullet$. Operators on Hilbert spaces are denoted with bold letters such as $\boldsymbol{a}$ or $\boldsymbol{X}$. Super-operators, i.e. operators acting e.g. on $\boldsymbol{X}$, are denoted via calligraphic capital letters such as $\mathcal{L}$ or $\mathcal{K}$.

2 Adiabatic elimination in bipartite systems

Consider on the Hilbert space $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ a Lindblad master equation with a two-time scale structure:

[TABLE]

where

• •

$\epsilon$ is a small positive parameter;

• •

$\mathcal{L}_{A}$ and $\mathcal{L}_{B}$ are Lindbladian super-operators acting only on $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$ respectively:

[TABLE]

with $\boldsymbol{H}_{\xi}$ a Hermitian operator on $\mathcal{H}_{\xi}$, $\boldsymbol{L}_{\xi,\mu}$ operators on $\mathcal{H}_{\xi}$ (not necessarily Hermitian) and with $\mathcal{D}_{\boldsymbol{X}}$ denoting the dissipation super-operator associated with operator $\boldsymbol{X}$,

[TABLE]

• •

$\mathcal{L}_{\!\rm{int}}$ is an interaction Lindbladian super-operator acting on both $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$.

We assume that for any initial density operator $\rho$ on $\mathcal{H}_{A}$, the solution $\frac{d}{dt}\rho=\mathcal{L}_{A}(\rho)$ converges exponentially towards a unique steady-state density operator $\overline{\rho}_{A}$.

Consider the solution of (2) starting from a density operator $\rho$ on $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$. For $\epsilon=0$, this solution converges towards the separable state $\overline{\rho}_{A}\otimes{\rm{tr}}_{A}(\rho)$, with ${\rm{tr}}_{A}(\rho)$ the partial trace over $\mathcal{H}_{A}$ of the initial density operator on $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$. For $\epsilon>0$ and small, this solution, after the relaxation time associated with the dynamics on $\mathcal{H}_{A}$ governed by $\mathcal{L}_{A}$, remains close to $\overline{\rho}_{A}\otimes\rho_{s}$ where $\rho_{s}$ is a slowly evolving density operator on $\mathcal{H}_{B}$. We explain here how to express quantitatively, to various orders in $\epsilon$, the fact that the solution of (2) remains close to $\overline{\rho}_{A}\otimes\rho_{s}$ and how $\rho_{s}$ evolves.

More precisely, we express the solution $\rho$ of (2) and the time derivative of $\rho_{s}$ as linear functions of $\rho_{s}$, by means the linear super-operators $\mathcal{K}$ and $\mathcal{L}_{s}$; afterwards, we develop these two super-operators in power of $\epsilon$, thus obtaining an asymptotic expansion of $\rho$ and the time derivative of $\rho_{s}$:

[TABLE]

with

[TABLE]

Plugging these series into (2):

[TABLE]

gives

[TABLE]

where the zero order term vanishes since $\mathcal{L}_{A}(\mathcal{K}_{0}(\rho_{s}))\equiv 0$. By identifying terms of same order versus $\epsilon$, it is possible to compute higher order terms from lower order terms as usual:

[TABLE]

Equation (6) corresponds to first order terms in $\epsilon$ and using the known expression of $\mathcal{K}_{0}$ it reads:

[TABLE]

Since ${\rm{tr}}_{A}(\mathcal{L}_{A})=0$, one takes the partial trace over $\mathcal{H}_{A}$ of this equation and gets

[TABLE]

To obtain $\mathcal{K}_{1}$, we must invert $\mathcal{L}_{A}$. Remark 1, in appendix, shows that the general solution $\boldsymbol{X}$ of $-\mathcal{L}_{A}(\boldsymbol{X})=\boldsymbol{W}-{\rm{tr}}(\boldsymbol{W})\overline{\rho}_{A}$ can be written as

[TABLE]

where $\overline{\mathcal{K}}_{A}$ is a linear trace-preserving and complete positive map on $\mathcal{H}_{A}$ and $\overline{\tau}$ a positive constant.

When $\ker\overline{\rho}_{A}\subset\ker\boldsymbol{W}$, i.e. $\boldsymbol{W}=\tilde{\boldsymbol{W}}\overline{\rho}_{A}$ for some operator $\tilde{\boldsymbol{W}}$, lemma 4 shows that $\ker\overline{\rho}_{A}\subset\ker\boldsymbol{X}$, i.e. that $\overline{\mathcal{K}}_{A}(\tilde{\boldsymbol{W}}\overline{\rho}_{A})=\tilde{\boldsymbol{X}}\overline{\rho}_{A}$ for some operator $\tilde{\boldsymbol{X}}$ (relevant when $\overline{\rho}_{A}$ is not of full rank). We will write as usual $\overline{\mathcal{K}}_{A}$ meaning $\overline{\mathcal{K}}_{A}\otimes\mathcal{I}$ when applied to operators on $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$. With this super-operator, using the just computed expression for $\mathcal{L}_{1}(\rho_{s})$, we readily get

[TABLE]

where $\boldsymbol{G}_{1,B}$ is any Hermitian operator on $\mathcal{H}_{B}$ (gauge degree of freedom). We will show in the following how the gauge choice enables to select between system propertie. We will show in the following how to fix the gauge in order to guarantee different properties (e.g. trace preservation, complete positivity, or a simple expression) of the expansion.

Lemma 5, in appendix, proves that $\rho_{s}\mapsto{\rm{tr}}_{A}(\mathcal{L}_{\!\rm{int}}(\overline{\rho}_{A}\otimes\rho_{s}))$ remains of Lindblad type when $\mathcal{L}_{\!\rm{int}}$ is of Lindblad type. Thus $\mathcal{L}_{1}$ in (8) is of Lindblad type and $\frac{d}{dt}\rho_{s}=\epsilon\mathcal{L}_{1}(\rho_{s})$ can be interpreted in this case as a Zeno dynamics. For $\mathcal{K}$, with the gauge $\boldsymbol{G}_{1,B}=0$ in (10), we have

[TABLE]

Since $\overline{\mathcal{K}}_{A}$, $e^{\epsilon\overline{\tau}\mathcal{L}_{\!\rm{int}}}$ and $\mathcal{K}_{0}$ are linear completely positive trace-preserving (CPTP) maps, the first order approximation $\mathcal{K}_{0}+\epsilon\mathcal{K}_{1}$ coincides, up to second order terms, with the linear CPTP map $\overline{\mathcal{K}}_{A}\circ e^{\epsilon\overline{\tau}\mathcal{L}_{\!\rm{int}}}\circ\mathcal{K}_{0}$. With the linear gauge \boldsymbol{G}_{1,B}(\rho_{s})=-\overline{\tau}{\rm{tr}}_{A}\Big{(}\mathcal{L}_{\!\rm{int}}(\overline{\rho}_{A}\otimes\rho_{s})\Big{)}, we get instead

[TABLE]

This expression is not always completely positive due to the backwards propagation with $\mathcal{L}_{1}$, but it is e.g. when $\mathcal{L}_{1}$ is Hamiltonian (thus in particular when $\mathcal{L}_{\!\rm{int}}$ is Hamiltonian).

The computations of the second order corrections $\mathcal{K}_{2}$ and $\mathcal{L}_{2}$ can be done via (2) along the same lines. The obtained expressions will depend on the gauge choice. We conjecture that, at any order $n$ versus $\epsilon$, we can choose $(\mathcal{K}_{j},\mathcal{L}_{j})_{1\leq j\leq n}$ such that all equations corresponding to orders less that $n$ are satisfied and additionally, such that $\sum_{j=0}^{n}\epsilon^{j}\mathcal{K}_{j}(\rho_{s})$ and $\sum_{j=1}^{n}\epsilon^{j}\mathcal{L}_{j}(\rho_{s})$ coincide, up to $n+1$ order terms, with a trace-preserving completely positive map and with a Lindbladian dynamics, respectively. We have seen on several examples that the related expressions can be quite involved and this is the topic of ongoing research. For instance, in general, the autonomous dynamics $\mathcal{L}_{B}$ undesirably appears in $\mathcal{L}_{2}$ if we take the gauge choice $\boldsymbol{G}_{1,B}=0$ which always guarantees the map $(\mathcal{K}_{0}+\epsilon\mathcal{K}_{1})$ to be CPTP.

We emphasize here that all gauge choices describe the same physical model i.e. for any gauge choice, the physical evolution of $\rho=\mathcal{K}(\rho_{s})$ (truncated at any order $\epsilon^{k}$) gives the same result (up to $\epsilon^{k+1}$ order terms). In this sense, the choice of the gauge degree of freedom has to be made in order to allow a better physical interpretation of $\rho_{s}$ and/or ensure physical properties (e.g. trace preservation, complete positivity $\dots$).

Therefore in this paper, we have carried out further analysis leading to more explicit expressions up to second order, for particular structures that are relevant in typical quantum systems. The associated conclusions, e.g. expressions proving that we retrieve trace-preserving completely positive maps and Lindblad type dynamics at this higher order as well, are presented in the next sections.

3 Hamiltonian interaction

One typical structure is when $\mathcal{L}_{\!\rm{int}}$ reduces to a Hamiltonian interaction:

[TABLE]

where $\boldsymbol{H}_{\!\rm{int}}=\sum_{k=1}^{m}\boldsymbol{A}_{k}\otimes\boldsymbol{B}_{k}^{\dagger}$ where $\boldsymbol{A}_{k}$ and $\boldsymbol{B}_{k}$ are (not necessarily Hermitian) operators on $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$, respectively. Such models are encountered e.g. in reservoir engineering, where $\mathcal{H}_{A}$ is the main dissipator allowing to stabilize the target system $\mathcal{H}_{B}$ (see example below). More generally, weak Hamiltonian coupling of a target system $\mathcal{H}_{B}$ to some environment is a standard model. The two typical cases leading to $\boldsymbol{H}_{\!\rm{int}}$ Hermitian are $\boldsymbol{A}_{k}$ and $\boldsymbol{B}_{k}$ Hermitian, as appears in “dispersive coupling” ($\boldsymbol{A}_{k}$ and $\boldsymbol{B}_{k}$ are diagonal in the energy bases of the respective subsystems); and $\boldsymbol{A}_{k+1}\otimes\boldsymbol{B}_{k+1}^{\dagger}=\boldsymbol{A}_{k}^{\dagger}\otimes\boldsymbol{B}_{k}$, as appears in “resonant coupling” ($\boldsymbol{A}_{k}$ and $\boldsymbol{B}_{k}$ are annihilation operators). We will treat these two cases separately for the second-order approximation.

For the first order approximation, we just plug this particular setting in the general formula and we get

[TABLE]

In particular, $\mathcal{L}_{1}$ is obviously of Lindblad type as the correction due to the coupling amounts to the Zeno Hamiltonian associated with $\boldsymbol{H}_{\!\rm{int}}$ (i.e. the Hamiltonian with $\mathcal{H}_{A}$ frozen to $\overline{\rho}_{A}$ by the fast dynamics).

For the remainder of this section, we take the gauge \boldsymbol{G}_{1,B}(\rho_{s})=-\overline{\tau}{\rm{tr}}_{A}\Big{(}\mathcal{L}_{\!\rm{int}}(\overline{\rho}_{A}\otimes\rho_{s})\Big{)} that will lead to simpler explicit expressions. First, for $\mathcal{K}_{0}+\epsilon\mathcal{K}_{1}$, by directly plugging the Hamiltonian interaction into (6) we obtain:

[TABLE]

where $\boldsymbol{F}_{k}\overline{\rho}_{A}=\overline{\tau}\overline{\mathcal{K}}_{A}\left(\,\boldsymbol{A}_{k}\;\overline{\rho}_{A}\,\right)-\overline{\tau}{\rm{tr}}(\boldsymbol{A}_{k}\;\overline{\rho}_{A})\overline{\rho}_{A}$ satisfies ${\rm{tr}}(\boldsymbol{F}_{k}\overline{\rho}_{A})=0$, for $k=1,2,...,m$. Thus the abstract inversion formula now has to be applied only on the given interaction operators $\boldsymbol{A}_{k}\overline{\rho}_{A}$ ($k=1,2,...,m$), instead of on any possible operator. Note that expression (14) differs from the general one (11) by second order terms and, while the expression without the $O(\epsilon^{2})$ terms is completely positive, $\boldsymbol{M}$ is not always Hermitian so corrections of $O(\epsilon^{2})$ might be necessary here to exactly preserve the trace.

We now turn to the second order.

3.1 Second order approximation

The gauge choice \boldsymbol{G}_{1,B}(\rho_{s})=-\overline{\tau}{\rm{tr}}_{A}\Big{(}\mathcal{L}_{\!\rm{int}}(\overline{\rho}_{A}\otimes\rho_{s})\Big{)} also facilitates the expression of $\mathcal{L}_{2}$. Taking the partial trace ${\rm{tr}}_{A}$ on (2), note that ${\rm{tr}}_{A}(\overline{\mathcal{K}}_{A}(\boldsymbol{X}))={\rm{tr}}_{A}(\boldsymbol{X})$ for all $\boldsymbol{X}$ on $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ thanks to trace preservation, and thus ${\rm{tr}}_{A}(\mathcal{K}_{1}(\boldsymbol{X}))=\overline{\tau}(\mathcal{L}_{1}(\boldsymbol{X})-\mathcal{L}_{B}(\boldsymbol{X}))+\boldsymbol{G}_{1,B}(\boldsymbol{X})$ for all $\boldsymbol{X}$ on $\mathcal{H}_{B}$. From there we get simply:

[TABLE]

3.1.1 Dispersive interaction

Let us first consider the case of a single interaction term, i.e. where $\boldsymbol{H}_{\!\rm{int}}=\boldsymbol{A}\otimes\boldsymbol{B}$ with $\boldsymbol{A},\boldsymbol{B}$ both Hermitian. Plugging this expression of $\mathcal{L}_{\!\rm{int}}$ and (14) for $\mathcal{K}_{1}$ into (15), we obtain

[TABLE]

By separating real and imaginary coefficients in the last two terms this can be rewritten as

[TABLE]

The first term contains a Hermitian Hamiltonian, second-order correction to the Zeno Hamiltonian of $\mathcal{L}_{1}$. The second term contains a dissipation in $\boldsymbol{B}$, at the rate ${\rm{tr}}(\boldsymbol{F}\overline{\rho}_{A}\boldsymbol{A}+\boldsymbol{A}\overline{\rho}_{A}\boldsymbol{F}^{\dagger})$. It expresses how the quantum uncertainty among eigenstates of $\boldsymbol{A}$ on $\mathcal{H}_{A}$, implies uncertainty in the Hamiltonian evolution with $\boldsymbol{B}$ on $\mathcal{H}_{B}$ when the two systems are dispersively coupled. There only remains to check that we will always get ${\rm{tr}}(\boldsymbol{F}\overline{\rho}_{A}\boldsymbol{A}+\boldsymbol{A}\overline{\rho}_{A}\boldsymbol{F}^{\dagger})\geq 0$, which is done by Lemma 7 in appendix.

Thus, the second order reduced dynamics for the dispersive interaction is given by:

[TABLE]

with $\boldsymbol{F}$ given by \boldsymbol{F}\overline{\rho}_{A}=-\mathcal{L}_{A}^{-1}\big{(}\boldsymbol{A}-{\rm{tr}}(\boldsymbol{A}\;\overline{\rho}_{A})\overline{\rho}_{A})\big{)}={+}\overline{\tau}\overline{\mathcal{K}}_{A}\big{(}\,\boldsymbol{A}\;\overline{\rho}_{A}-{\rm{tr}}(\boldsymbol{A}\;\overline{\rho}_{A})\overline{\rho}_{A}\big{)}.

3.1.2 General case and resonant interaction

For $\boldsymbol{H}_{\!\rm{int}}$ comprising several terms, we can still put $\mathcal{L}_{2}$ in a form with Hamiltonian and dissipation super-operator by following the same procedure. Instead of directly getting a dissipation rate ${\rm{tr}}(\boldsymbol{F}\overline{\rho}_{A}\boldsymbol{A}+\boldsymbol{A}\overline{\rho}_{A}\boldsymbol{F}^{\dagger})$ as for the dispersive interaction, we here get a positive semidefinite matrix in $\mathbb{C}^{m\times m}$ which determines dissipation operators consisting of linear combinations of the $\boldsymbol{B}_{k}$. The number of decoherence channels in $\mathcal{L}_{2}$ will be equal to the number $\leq m$ of linearly independent terms in $\boldsymbol{H}_{\!\rm{int}}$. Explicitly, plugging the expression (14) of $\mathcal{K}_{1}$ and the Hamiltonian interaction into (15), a few algebraic manipulations readily yield:

[TABLE]

Both matrices $X$ and $Y$ $\in\mathbb{C}^{m\times m}$ are Hermitian. Thus the first term contains a Hermitian Hamiltonian. Lemma 7 in appendix proves that $X$ is positive semidefinite. Then equation (3.1.2) is of standard Lindblad form. Indeed, writing $X=\Lambda\Lambda^{\dagger}$, where $\Lambda$ can be obtained for instance easily by Cholesky factorization, we get

[TABLE]

The second order reduced dynamics is then given by:

[TABLE]

The choice of $\Lambda$ is not unique, reflecting the non-uniqueness of Lindblad representations with several dissipation channels. Usually and as illustrated in the example of Section 3.2, a preferred decomposition, identifying physically meaningful components, can be imposed.

For a single resonant interaction $\boldsymbol{H}_{\!\rm{int}}=\boldsymbol{A}\otimes\boldsymbol{B}^{\dagger}+\boldsymbol{A}^{\dagger}\otimes\boldsymbol{B}$, we thus have $\boldsymbol{F}_{1}$ and $\boldsymbol{F}_{2}$ given essentially by $-\mathcal{L}_{A}^{-1}\equiv\overline{\tau}\overline{\mathcal{K}}_{A}$ acting on $\boldsymbol{A}\overline{\rho}_{A}-{\rm{tr}}{\boldsymbol{A}\overline{\rho}_{A}}\overline{\rho}_{A}$ and on $\boldsymbol{A}^{\dagger}\overline{\rho}_{A}-{\rm{tr}}{\boldsymbol{A}^{\dagger}\overline{\rho}_{A}}\overline{\rho}_{A}$ respectively; $X,Y,\Lambda$ $\in\mathbb{C}^{2\times 2}$ defined as above with no particular simplifications, unless $\overline{\rho}_{A}$ takes a particular form; and finally

[TABLE]

The following example considers a combination of dispersive and resonant interaction, where $m=3$. The result on this example is new to our knowledge, and with the number of decoherence channels that have to be combined, it is not easily obtained from intuitive reasoning.

3.2 Example: Two-photon pumping

We illustrate the result of Section 3.1.2 on the model of the experiment presented in [13, 30]. Such system, following the theoretical proposal [24], is a promising way towards dynamically protected quantum processors. The reservoir is based on a scheme that induces 2-photon loss at order $\epsilon^{2}$ on the target system, as a dominant effect. We show how our method can calculate a more precise reduced model taking into account all the effects at order $\epsilon^{2}$, and thereby quantify the effects of potential 2-photon excitation and cross-Kerr nonlinearity on this reservoir.

The system is composed of two interacting cavities: the first cavity (fast subsystem, $\mathcal{H}_{A}$) is driven by an electromagnetic field and exchanges energy with the environment, with an energy loss term dominant with respect to both the energy gain term and the electromagnetic field drive. This implies that only the lowest energy level is significantly populated, so we can model the cavity as a two level system, i.e. a qubit with energy levels $\left|g\right\rangle,\left|e\right\rangle$. The second cavity (slow subsystem, $\mathcal{H}_{B}$) weakly interacts with this qubit, with auxiliary pump fields matched such that the resonant interaction exchanges 1 energy quantum of $\mathcal{H}_{A}$ with 2 energy quanta of $\mathcal{H}_{B}$ [13, 30]; an additional residual dispersive interaction (“cross-Kerr”) is unavoidable in this setup.

The system model thus writes as follows. As usual, $\boldsymbol{\sigma_{\!x}},\boldsymbol{\sigma_{\!y}},\boldsymbol{\sigma_{\!z}}$ denote the standard Pauli matrices, $\boldsymbol{\sigma_{\!\text{\bf-}}}=|g\rangle\langle e|$ and $\boldsymbol{\sigma_{\!\text{\bf+}}}=|e\rangle\langle g|$ the energy loss operator and energy gain operators respectively, all on $\mathcal{H}_{A}$, and $\boldsymbol{b},\boldsymbol{b}^{\dagger}$ denote respectively the annihilation and creation operators, on $\mathcal{H}_{B}$. The dynamics of the fast subsystem are described by the Lindbladian operator

[TABLE]

where the coupling constants satisfy

[TABLE]

The interaction of the two systems is described by

[TABLE]

The fact that $\boldsymbol{H}_{\!\rm{int}}$ is much weaker than $\mathcal{L}_{A}$ is expressed by $|g|,|\chi|\ll\kappa_{+},u,\kappa_{-}$. Formally we can take units such that $\kappa_{+},|u|,\kappa_{-}$ are of order $1$ or larger, and $g,\chi$ are of order $\epsilon\ll 1$. In particular, if $g=\chi=0$, then the two systems are independent. The fast system converges to the steady state

[TABLE]

while the slow target system $B$ does not move.

When $g,\chi$ are nonzero, the dynamics of the slow dynamics, approximately corresponding to the second cavity $\mathcal{H}_{B}$, can be written in the form (5). By simple computations this yields the first-order dynamics of system $B$, given by the Zeno Hamiltonian:

[TABLE]

We now compute the second-order dynamics, choosing in equation (9) the gauge \boldsymbol{G}_{1,B}(\rho_{s})=-\overline{\tau}{\rm{tr}}_{A}\Big{(}\mathcal{L}_{\!\rm{int}}(\overline{\rho}_{A}\otimes\rho_{s})\Big{)}. Algebraic computations, by solving for $\mathcal{L}_{A}^{-1}$ with the Bloch equations for the qubit, yield $X$ and $Y$ in the form:

[TABLE]

[TABLE]

Towards interpreting these expressions, we take into account the relative strengths of the couplings (22). More precisely, with $\kappa_{-}=1+n_{\text{thermal}}$, $\kappa_{+}=n_{\text{thermal}}$ and $n_{\text{thermal}}\ll 1$, we define $\delta^{2}=\kappa_{+}/\kappa_{-}\approx n_{\text{thermal}}$ and $\eta=u/\kappa_{-}$, and we neglect the terms of order three or higher in $\delta,\eta$.

The Hamiltonian operator in (19) then reads

[TABLE]

up to terms of order at least three in $\delta,\eta$. For the dissipative part, up to the same terms,

[TABLE]

From this form, it is already clear that the dominant effect of the dissipation involves the two-photon annihilation operator. We next write $X=\Lambda\Lambda^{{\dagger}}$, choosing $\Lambda^{{\dagger}}$ as an upper triangular matrix:

[TABLE]

so that the three dissipation channels are given by the operators

[TABLE]

With the same approximation, the terms in $\mathcal{L}_{1}$ involve $(1+z_{\infty})/2=\delta^{2}+4\eta^{2}$ and $x_{\infty}/2=-2\eta$. Wrapping up, the effects on the slow subsystem are thus, with units such that $\kappa_{-}=1$:

• •

At order $(g\eta)$ and $(g^{2})$ respectively: a two-photon pumping Hamiltonian, in $\boldsymbol{b}^{{\dagger}2}+\boldsymbol{b}^{2}$, and a two-photon dissipation, with $\boldsymbol{L}_{1}\approx\boldsymbol{b}^{2}$, precisely as intended in [24];

• •

At order $(\chi\delta^{2},\chi\eta^{2})$: a Stark shift Hamiltonian $\boldsymbol{b}^{\dagger}\boldsymbol{b}$, which just shifts the cavity frequency and can be compensated for;

• •

At order $(g\chi\eta)$: a Hamiltonian in $(\boldsymbol{b}^{\dagger}\boldsymbol{b}^{3}-\boldsymbol{b}^{{\dagger}3}\boldsymbol{b})/i$, whose precise deformation effect would have to be investigated; and a modification of the two-photon dissipation channel $\boldsymbol{L}_{1}$ by some dephasing effect, leading to terms like $\boldsymbol{b}^{2}\rho\boldsymbol{b}^{\dagger}\boldsymbol{b}$.

• •

At order $(\chi^{2}\delta^{2},g^{2}\delta^{2},g^{2}\eta^{2})$: a Hamiltonian effect in $(\boldsymbol{b}^{4}-\boldsymbol{b}^{{\dagger}4})/i$, that is essentially a 4-photon drive; two new dissipation channels, namely in $\boldsymbol{b}^{{\dagger}2}$ and $\boldsymbol{b}^{\dagger}\boldsymbol{b}$; and a further modification of $\boldsymbol{L}_{1}$, now leading to terms like $\boldsymbol{b}^{2}\rho\boldsymbol{b}^{2}$.

These are the dominant ones provided $|g|,|\chi|\ll|\eta|,|\delta|\ll 1$; once $g,\chi$ become comparable to $\eta,\delta$, effects of order $\epsilon^{3}$ might become important as well. The next step of our asymptotic expansion method would address this point; carrying out the related computations goes beyond the present paper.

Finally, we must emphasize that the reduced dynamics describes the evolution of the slow subsystem $\rho_{S}$. This $\rho_{S}$ does not exactly coincide with the state of subsystem $B$. Indeed, the Kraus map $\mathcal{K}_{0}+\epsilon\mathcal{K}_{1}+...$ (not given here) expresses precisely how subsystem $B$ gets slightly hybridized with subsystem $A$ in order to obtain Markovian slow dynamics. The dynamics of just ${\rm{tr}}_{A}(\rho)$ would not be Markovian at this order of precision, and this is one reason that precludes a simple intuitive derivation of these results.

4 Cascade interaction

We consider in this section a different kind of interaction between the fast and the slow subsystem. Instead of the standard reciprocal Hamiltonian coupling of Section 3 we consider a unidirectional coupling which allows the output of the first system to feed the second system while forbidding the reverse process. The expression ”cascaded system” was first introduced in [8], see also e.g. [17] for more details on this type of interaction. The network extension of the cascade structure can be found in [19], see [20, 26] for its relation with adiabatic elimination.

The master equation of a fast subsystem $A$ on a Hilbert space $\mathcal{H}_{A}$ coupled via cascaded interactions to a slow subsystem $B$ on a Hilbert space $\mathcal{H}_{B}$ is given by:

[TABLE]

where $\mathcal{L}_{A}$ and $\mathcal{L}_{B}$ are Lindbladian super-operators acting only on $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$, respectively, and where $\boldsymbol{a}$ and $\boldsymbol{b}$ are operators acting only on $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$ respectively. We assume that there exists a unique steady state density operator $\overline{\rho}_{A}$ solution of $\mathcal{L}_{A}(\rho)+\mathcal{D}_{\boldsymbol{a}}(\rho)=0$. Note that unlike the Hamiltonian interaction case, the fast dynamics is given by the super-operator $\mathcal{L}_{A}(\rho)$ but also by $\mathcal{D}_{\boldsymbol{a}}(\rho)$ due to the unidirectional coupling. The completely positive map $\overline{\mathcal{K}}_{A}$ of Remark 1 is defined accordingly.

We detail now the computation leading to the adiabatic elimination of the fast variables associated with the subsystem $A$ in the case of a cascaded system whose dynamics is given by (24). The zero order approximation is readily given by $\mathcal{L}_{0}(\rho_{s})=0$ and $\mathcal{K}_{0}(\rho_{s})=\overline{\rho}_{A}\otimes\rho_{s}$.

The first order reduced dynamics can be directly computed using equation (8). A straightforward calculation yields:

[TABLE]

The super-operator $\mathcal{K}_{1}$ corresponding to the first order entanglement between the two subsystems is given by (10). From this expression, with the gauge choice \boldsymbol{G}_{1,B}(\rho_{s})=-\overline{\tau}{\rm{tr}}_{A}\bigl{(}\mathcal{L}_{int}(\overline{\rho}_{A}\otimes\rho_{s})\bigr{)} we get

[TABLE]

where $\overline{\boldsymbol{a}}=\boldsymbol{a}-{\rm{tr}}(\boldsymbol{a}\overline{\rho}_{A})\boldsymbol{I}_{A}$. We made this particular gauge choice in order to ensure ${\rm{tr}}_{A}(\mathcal{K}(\rho_{s}))=\rho_{s}$ and therefore $\rho_{s}$ corresponds to the density operator of the subsystem $B$. This gauge choice is different from the one of Section 2. However we show that $\mathcal{K}(\rho_{s})=\mathcal{K}_{0}(\rho_{s})+\epsilon\mathcal{K}_{1}(\rho_{s})$ is also a completely positive trace-preserving map up-to second order terms. To exhibit such property, we express this map in an explicit Kraus operator form.

Using Lemma 4 and Corollary 1 we get that $\overline{\mathcal{K}}_{A}\left(\overline{\boldsymbol{a}}\hskip 2.84526pt\overline{\rho}_{A}\right)=\overline{\mathcal{K}}_{A}\left(\overline{\boldsymbol{a}}\hskip 2.84526pt\overline{\rho}_{A}\right)\overline{\rho}_{A}^{-1}\overline{\rho}_{A}$ and $\overline{\mathcal{K}}_{A}\left(\overline{\rho}_{A}\overline{\boldsymbol{a}}^{\dagger}\right)=\overline{\mathcal{K}}_{A}\left(\overline{\rho}_{A}\overline{\boldsymbol{a}}^{\dagger}\right)\overline{\rho}_{A}^{-1}\overline{\rho}_{A}$ where $\overline{\rho}_{A}^{-1}$ stands for the Moore-Penrose pseudo-inverse of $\overline{\rho}_{A}$. Therefore we derive the following Kraus operator form for the first order entanglement:

[TABLE]

4.1 Second order approximation

In order to compute the second order dynamics we proceed in the same manner as the Hamiltonian interaction case: first, we choose the gauge \boldsymbol{G}_{1,B}(\rho_{s})=-\overline{\tau}{\rm{tr}}_{A}\bigl{(}\mathcal{L}_{int}(\overline{\rho}_{A}\otimes\rho_{s})\bigr{)} in order to facilitate the expression of $\mathcal{L}_{2}$. Then, take the partial trace ${\rm{tr}}_{A}$ on equation (2) leading to the cancellation of the unknown term $\mathcal{L}_{A}(\mathcal{K}_{2}(\rho_{s}))$. Using ${\rm{tr}}_{A}(\mathcal{L}_{A}(\bullet))\equiv 0$ with ${\rm{tr}}_{A}(\mathcal{K}_{0}(\bullet))\equiv\bullet$ and the fact that ${\rm{tr}}_{A}(\mathcal{K}_{1}(\bullet))\equiv 0$ due to our choice of the gauge degree of freedom we get:

[TABLE]

Note that the term $\mathcal{D}_{\boldsymbol{b}}$ is intrinsically related to the cascaded structure while $\mathcal{L}_{B}$ represents the slow dynamics on the subsystem $B$ and is independent of the structure (in particular, it could satisfy $\mathcal{L}_{B}\equiv 0$). Therefore the term $\mathcal{L}_{B}$ must be treated independently from the other terms. For this reason we impose $\mathcal{L}_{2}=\widetilde{\mathcal{L}_{2}}+\mathcal{L}_{B}$ with

[TABLE]

A direct expansion leads to the following expression:

[TABLE]

where $\{x_{1},x_{2},y_{1},y_{2}\}\in\mathbb{C}^{4}$ are solutions of the set of equations:

[TABLE]

with $\alpha=\overline{\tau}{\rm{tr}}\left(\boldsymbol{a}\overline{\mathcal{K}}_{A}\left(\overline{\rho}_{A}\overline{\boldsymbol{a}}^{\dagger}\right)\right)$ and $\beta=\overline{\tau}{\rm{tr}}(\boldsymbol{a}^{\dagger}\overline{\mathcal{K}}_{A}(\overline{\rho}_{A}\overline{\boldsymbol{a}}^{\dagger}))$.

This set of equation can be viewed as the norm and the scalar product between the two vectors $\{x_{1},x_{2}\}^{\intercal}\in\mathbb{C}^{2}$ and $\{y_{1},y_{2}\}^{\intercal}\in\mathbb{C}^{2}$. Therefore, there exists a (non-unique) solution to (35) if and only if holds:

[TABLE]

Using a slight adaptation of Lemma 7, we get that condition (36) holds for any system described by the master equation (24). We conjecture that (37) holds for any system. The non-uniqueness of the solutions is expected due to the fact that the Lindblad decomposition (31) is also not unique.

In the case of cascade interaction, the second order reduced dynamics is given by:

[TABLE]

4.2 Arbitrary system with a squeezed drive

To illustrate our results, we consider a driven linear cavity producing a squeezed output field, unidirectionally ”feeding” the slow subsystem. Such kind system was first studied in [16] considering the slow subsystem as a qubit and shows how one can engineer the coherence times of this qubit using a squeezed input field. It has recently been realized experimentally in [29] validating the theoretical results. Our aim is to show through this example, how to apply our method, emphasizing that it doesn’t need the specification of the slow subsystem and readily retrieve the known results of the adiabatic elimination of the linear cavity when the slow subsystem is a qubit.

Using [19] , we get the following dynamics for the system:

[TABLE]

where $\boldsymbol{a}$ is the annihilation operator on the cavity, $\boldsymbol{b}$ is an operator on the unspecified slow subsystem (usually the annihilation operator for a cavity and the energy loss operator for a qubit). The term $g[\boldsymbol{a}^{2}-{\boldsymbol{a}^{\dagger}}^{2},\rho]$ corresponds to the standard squeezing Hamiltonian by an appropriate phase choice for $\boldsymbol{a}$ ($g$ is real here). The parameter $\kappa$ is the dissipation rate of the cavity. The time-scale separation is given by $\kappa\gg\epsilon^{2}$. This model is valid under the assumption $\kappa>4g$, otherwise the cavity subsystem is unstable and its energy grows to infinity (and therefore additional phenomenon such as Kerr effect have to be taken into account). As any dynamics $\mathcal{L}_{B}$ on the slow subsystem plays no role for second order computations, we assume for simplicity $\mathcal{L}_{B}\equiv 0$.

To derive the reduced dynamics up to second order on this example, we use formulas (25) and (31) where the fast system $A$ is the cavity and the slow system $B$ is unspecified. Therefore we have only to compute the three coefficients ${\rm{tr}}(\sqrt{\kappa}\boldsymbol{a}\overline{\rho}_{A})$, $\alpha=\kappa\overline{\tau}{\rm{tr}}\left(\boldsymbol{a}\overline{\mathcal{K}}_{A}\left(\overline{\rho}_{A}\overline{\boldsymbol{a}}^{\dagger}\right)\right)$ and $\beta=\kappa\overline{\tau}{\rm{tr}}(\boldsymbol{a}^{\dagger}\overline{\mathcal{K}}_{A}(\overline{\rho}_{A}\overline{\boldsymbol{a}}^{\dagger})$. As the fast dynamics of (39) corresponds to the evolution of a linear quantum harmonic oscillator, these coefficients are easier to compute by taking their Heisenberg representation counterpart i.e:

[TABLE]

where $\rho_{0}$ is the initial state, ${\mathcal{L}_{A}}^{*}$ is the dual of $\mathcal{L}_{A}$. With a small abuse of notation, $e^{\infty{\mathcal{L}_{A}}^{*}}$ means $\lim_{t\to\infty}e^{t{\mathcal{L}_{A}}^{*}}$.

Then, some usual computations lead to (see Appendix C for the key elements of the computation):

[TABLE]

In this case one can check that $(\alpha+\alpha^{*}+1)(\alpha+\alpha^{*})=4|\beta|^{2}$ verifies condition (37). It is therefore possible to solve (35) and the simplest solution is given by $x_{2}=y_{2}=0,x_{1}=\sqrt{2\alpha+1},y_{1}=-\sqrt{2\alpha}$. Noting that $\alpha\in\mathbb{R}$, we get the second order reduced dynamics:

[TABLE]

Recovering the result of [16] when the slow subsystem is a qubit and $\boldsymbol{b}=\boldsymbol{\sigma_{\!\text{\bf-}}}$ the energy loss operator.

5 Conclusion

We have presented an adiabatic elimination technique for bipartite quantum systems characterized by one subsystem converging rapidly towards a unique steady state. By considering the slow dynamics as a perturbation of order $\epsilon$, we explicitly obtain, up to second order terms in $\epsilon$, the reduced slow dynamics as a Lindblad master equation. Our explicit formulas directly apply to reservoir engineering settings e.g. for computing perturbations on the engineered reservoir, or the residual dissipation added when first-order Hamiltonian Zeno dynamics is introduced towards performing gates on top of the reservoir.

The asymptotic expansion method, developed in this paper on bipartite systems, seems in fact to be applicable for any type of quantum system with two time-scales. In [2, 13] for instance it is exploited in the context of fast convergence towards a decoherence free space, with some preliminaries results that would have to be completed. In this case, the fast dissipative dynamics drives the system towards a particular subspace inside the Hilbert space. Hence the slow/fast decomposition involves a Cartesian product rather than a tensor product as presented here.

We conjecture that this geometric adiabatic elimination technique will lead to formulas conveying a physical interpretation of the reduced model at higher orders as well: by ensuring the Lindblad form of the reduced dynamics, and by providing a completely-positive trace-preserving mapping from reduced model to actual system states. The key element to do so would be to properly choose the gauge degree of freedom at each order. This is the topic of ongoing research.

Appendix A Inverse of Lindblad super-operators via Kraus maps

Consider the Lindblad master equation $\frac{d}{dt}\rho=\mathcal{L}(\rho)$ where $\rho$ is a density operator on a finite dimensional Hilbert space $\mathcal{H}$. Assume that for any operator $\boldsymbol{X}$ on $\mathcal{H}$, $e^{t\mathcal{L}}(\boldsymbol{X})$ converges exponentially towards a fixed point depending on $\boldsymbol{X}$. This means that there exists a complete-positive and trace-preserving map $\overline{\mathcal{R}}$ such that $\lim_{t\mapsto+\infty}e^{t\mathcal{L}}(\boldsymbol{X})=\overline{\mathcal{R}}(\boldsymbol{X})$. Thus we have $\mathcal{L}(\overline{\mathcal{R}}(\boldsymbol{X}))\equiv 0\equiv\overline{\mathcal{R}}(\mathcal{L}(\boldsymbol{X}))\equiv 0$ since $e^{t\mathcal{L}}(\overline{\mathcal{R}}(\boldsymbol{X}))\equiv\overline{\mathcal{R}}(\boldsymbol{X})\equiv\overline{\mathcal{R}}(e^{t\mathcal{L}}(\boldsymbol{X}))$.

Lemma 1**.**

There exists $\bar{\tau}>0$ such that the super-operator $\overline{\mathcal{K}}$ sending operator $\boldsymbol{X}$ to

[TABLE]

is a linear, trace-preserving and completely positive mapping with

[TABLE]

Remark 1**.**

When for any initial density operator the solution of $\frac{d}{dt}\rho=\mathcal{L}(\rho)$ converges toward a unique density operator $\overline{\rho}$, we have $\overline{\mathcal{R}}(\boldsymbol{X})={\rm{tr}}(\boldsymbol{X})\overline{\rho}$. In this case, for any given operator $\boldsymbol{W}$ with ${\rm{tr}}(\boldsymbol{W})=0$, the general solution of $-\mathcal{L}(\boldsymbol{X})=\boldsymbol{W}$ reads

[TABLE]

where $\lambda$ is an arbitrary complex number. Moreover $\boldsymbol{X}=\bar{\tau}\overline{\mathcal{K}}(\boldsymbol{W})$ is the unique solution with zero trace.

Proof.

Due to exponential convergence of $e^{t\mathcal{L}}(\boldsymbol{X})$ towards $\overline{\mathcal{R}}(\boldsymbol{X})$, the indefinite integral

[TABLE]

is absolutely convergent. Since \frac{d}{dt}e^{t\mathcal{L}}\Big{(}\boldsymbol{X}-\overline{\mathcal{R}}(\boldsymbol{X})\Big{)}=\mathcal{L}\left(e^{t\mathcal{L}}\Big{(}\boldsymbol{X}-\overline{\mathcal{R}}(\boldsymbol{X})\Big{)}\right), we have $\mathcal{L}(\overline{\mathcal{K}}(\boldsymbol{X}))=-\frac{\boldsymbol{X}-\overline{\mathcal{R}}(\boldsymbol{X})}{\bar{\tau}}$. Since, for each $t\geq 0$ the propagator $e^{t\mathcal{L}}$ is trace preserving and ${\rm{tr}}(\boldsymbol{X})={\rm{tr}}\left(\overline{\mathcal{R}}\boldsymbol{X}\right)$, simple computations yield ${\rm{tr}}(\overline{\mathcal{M}}(\boldsymbol{X}))=0$ and thus ${\rm{tr}}(\overline{\mathcal{K}}(\boldsymbol{X}))={\rm{tr}}(\boldsymbol{X})$.

To prove complete-positivity, consider the extension of $\mathcal{L}$, $\overline{\mathcal{K}}$, $\overline{\mathcal{R}}$ on the tensor product $\mathcal{H}\otimes\widetilde{\mathcal{H}}$ where $\widetilde{\mathcal{H}}$ is any Hilbert space of finite dimension. Let us prove that for $\bar{\tau}$ large enough, such extension of $\overline{\mathcal{K}}$ is non-negative, i. e., that for any $\left|\Phi\right\rangle,\left|\Psi\right\rangle\in\mathcal{H}\otimes\widetilde{\mathcal{H}}$, we have \left\langle\Psi\Big{|}\overline{\mathcal{K}}\left(\left|\Phi\right\rangle\!\left\langle\Phi\right|\right)\Big{|}\Psi\right\rangle\geq 0. Consider an Hilbert basis $\left(\left|n\right\rangle\right)_{1\leq n\leq d}$ of $\mathcal{H}$ whose dimension is denoted by $d$. Take

[TABLE]

where, for each $n$ and $\nu$ in $\{1,\ldots,d\}$, $\left|\phi_{n}\right\rangle,\left|\psi_{\nu}\right\rangle\in\widetilde{\mathcal{H}}$. Then standard computations give

[TABLE]

with $z_{n,\nu}=\left\langle\phi_{n}\right|\left.\psi_{\nu}\right\rangle$ and

[TABLE]

Similarly

[TABLE]

with R_{n^{\prime},\nu^{\prime},n,\nu}=\left\langle\nu^{\prime}\Big{|}\overline{\mathcal{R}}(\left|n^{\prime}\right\rangle\!\left\langle n\right|)\Big{|}\nu\right\rangle. This means that, defining the $d^{2}$-dimensional vector $z=(z_{n,\nu})_{n,\nu\in\{1,\ldots,d\}}$ and the $d^{2}\times d^{2}$ Hermitian matrices M=\Big{(}M_{n^{\prime},\nu^{\prime},~{}n,\nu}\Big{)} and R=\Big{(}R_{n^{\prime},\nu^{\prime},~{}n,\nu}\Big{)} , we have the following quadratic forms

[TABLE]

where $z$ depends on $\left|\Phi\right\rangle$ and $\left|\Psi\right\rangle$, where $M$ and $R$ depend only on $\overline{\mathcal{M}}$ and $\overline{\mathcal{R}}$. We have thus to prove that there exists $\bar{\tau}>0$ such that $M+\bar{\tau}R\geq 0$. Since $\overline{\mathcal{R}}$ is a complete-positive map, the $d^{2}\times d^{2}$ Hermitian matrix $R$ is non-negative.

Take $z$, such that $Rz=0$. Take $T>0$. We have

[TABLE]

since e^{t\mathcal{L}}\Big{(}\overline{\mathcal{R}}(\left|\Phi\right\rangle\!\left\langle\Phi\right|)\Big{)}=\overline{\mathcal{R}}(\left|\Phi\right\rangle\!\left\langle\Phi\right|) and  \left\langle\Psi\Big{|}\overline{\mathcal{R}}(\left|\Phi\right\rangle\!\left\langle\Phi\right|)\Big{|}\Psi\right\rangle=z^{\dagger}Rz=0. Since for each $t\geq 0$, $e^{t\mathcal{L}}$ is completely positive, then there exists a Kraus decomposition

[TABLE]

with operators $\boldsymbol{W}_{\mu,t}$ on $\mathcal{H}$ such that $\sum_{\mu}\boldsymbol{W}_{\mu,t}^{\dagger}\boldsymbol{W}_{\mu,t}=\boldsymbol{I}$. We have

[TABLE]

Consequently

[TABLE]

Since $z^{\dagger}Mz$ is the limit when $T$ tends to $+\infty$ of

[TABLE]

we have for any $T>0$

[TABLE]

that is, M is non-negative definite. In particular, if we assume that $z^{\dagger}Mz=0$. The above inequality implies that for any $t>0$

[TABLE]

Recall that, by assumption, \sum_{n,\nu}\left\langle\nu^{\prime}\Big{|}\overline{\mathcal{R}}(\left|n^{\prime}\right\rangle\!\left\langle n\right|)\Big{|}\nu\right\rangle~{}z_{n,\nu}=0 for any $n^{\prime},\nu^{\prime}\in\{1,\ldots,d\}$. Consequently

[TABLE]

Thus for any $z$ such that $Rz=0$ and $z^{\dagger}Mz=0$, we have necessarily $Mz=0$ by taking the limit for $T$ tending to $+\infty$.

To summarize we have shown that

$R\geq 0$ ; 2. 2.

if $z^{\dagger}Rz=0$ then $z^{\dagger}Mz\geq 0$; 3. 3.

if $z^{\dagger}Rz=z^{\dagger}Mz=0$ then $Mz=0$.

According to Lemma 2 proved below, there exists $\bar{\tau}>0$ such that $M+\bar{\tau}R\geq 0$.

∎

Lemma 2**.**

Consider two Hermitian matrices of same dimension $R$ and $M$ such that $R$ is non negative, such that $z^{\dagger}Rz=0$ implies that $z^{\dagger}Mz\geq 0$, and such $z^{\dagger}Rz=z^{\dagger}Mz=0$ implies $Mz=0$. Then for $\tau\geq 0$ large enough, $M+\tau R\geq 0$.

Proof.

Up to a unitary transformation, we have the block decomposition associated with $\ker R$ and $\ker R^{\perp}$:

[TABLE]

with $A$, $B$ and $D$ Hermitian matrices with $D>0$. For any $z$ such that $z^{\dagger}Rz=0$ we have $z^{\dagger}Mz\geq 0$, this means that $A\geq 0$. Up to some unitary transformation on $A$ only, we can always assume the following sub-block decomposition for $A$, $C$ and $x$

[TABLE]

with $\bar{A}>0$. According to these block decompositions, $z^{\dagger}Rz=z^{\dagger}Mz=0$ means that $y=0$ and $\bar{x}=0$ with $\tilde{x}$ arbitrary. But $Mz=0$ means that

[TABLE]

for all $\tilde{x}$. Thus $\tilde{C}=0$. To summarize, up to a unitary transformation, we have the following decomposition for $M$ and $R$

[TABLE]

with $\bar{A}>0$ and $D>0$. For $\tau$ large enough $\tau\bar{A}\geq\bar{C}^{\dagger}D^{-1}\bar{C}$ and thus $M+\tau R\geq 0$ (see, e.g, [4, Theorem 1.3.3, page 14]). ∎

We end this section by proving a factorization property with $\overline{\rho}_{A}$ for the result of $\mathcal{L}_{A}^{-1}$ used after (9) and then in Section 3. This property is a direct consequence of the following two short lemmas, first introduced in [3]. For the sake of self-completeness we also present their proof. Then, we use these lemmas to derive Corollary 1 in the particular case of composite systems with a cascaded structure, as presented in Section 4.

Lemma 3**.**

Let $\mathcal{L}_{A}(\overline{\rho}_{A})=0$, with $\mathcal{L}_{A}$ of the form (3) and $\overline{\rho}_{A}$ a density operator. Then for any $\left|\nu\right\rangle\in\ker(\overline{\rho}_{A})$ we have $\sqrt{\overline{\rho}_{A}}\boldsymbol{L}_{A,\mu}^{\dagger}\left|\nu\right\rangle=0$, for all $\mu$.

Proof.

For $\left|\nu\right\rangle\in\ker(\overline{\rho}_{A})$ we have $\left\langle\nu\right|\mathcal{L}_{A}(\overline{\rho}_{A})\left|\nu\right\rangle=\sum_{k}\left\langle\nu\right|\boldsymbol{L}_{A,\mu}\overline{\rho}_{A}\boldsymbol{L}_{A,\mu}^{\dagger}\left|\nu\right\rangle$. Since $\mathcal{L}_{A}(\overline{\rho}_{A})=0$ each term of this positive sum must vanish. ∎

Lemma 4**.**

Denote by $\rho=\overline{\rho}_{A}$ the unique density operator solution of $\mathcal{L}_{A}(\rho)=0$. For a traceless operator $\boldsymbol{Y}$ such that $\ker(\overline{\rho}_{A})\subseteq\ker(\boldsymbol{Y})$, the traceless solution to $\boldsymbol{X}=\mathcal{L}_{A}^{-1}(\boldsymbol{Y})$ also satisfies $\ker(\overline{\rho}_{A})\subseteq\ker(\boldsymbol{X})$.

Proof.

Note that the operators have such kernels if and only if they can be written $\boldsymbol{X}=\tilde{\boldsymbol{X}}\overline{\rho}_{A}$, $\boldsymbol{Y}=\tilde{\boldsymbol{Y}}\overline{\rho}_{A}$. Since $\mathcal{L}_{A}$ is a bijection on the space of traceless operators, the property is equivalent to show that $\boldsymbol{Y}\left|\nu\right\rangle=\mathcal{L}_{A}(\tilde{\boldsymbol{X}}\overline{\rho}_{A})\left|\nu\right\rangle=0$ for all $\left|\nu\right\rangle\in\ker(\overline{\rho}_{A})$. By using $\overline{\rho}_{A}\left|\nu\right\rangle=0$ and Lemma 3, we directly get

[TABLE]

Subtracting $0=\mathcal{L}_{A}(\overline{\rho}_{A})$ inside the bracket, applying $\overline{\rho}_{A}\left|\nu\right\rangle=0$ and Lemma 3 once again, we do get 0. ∎

Corollary 1**.**

Let the density operator $\overline{\rho}_{A}$ be solution of $\mathcal{L}_{A}(\overline{\rho}_{A})=\widetilde{\mathcal{L}}_{A}(\overline{\rho}_{A})+\mathcal{D}_{\boldsymbol{a}}(\overline{\rho}_{A})=0$ with $\widetilde{\mathcal{L}}_{A}$ of the form (3) . For the traceless operator $\boldsymbol{Y}=\overline{\rho}_{A}a^{\dagger}-{\rm{tr}}(\overline{\rho}_{A}a^{\dagger})\overline{\rho}_{A}$, the traceless solution to $\boldsymbol{X}=\mathcal{L}_{A}^{-1}(\boldsymbol{Y})$ satisfies $\ker(\overline{\rho}_{A})\subseteq\ker(\boldsymbol{X})$.

Proof.

From Lemma 3 we get that for any $\left|\nu\right\rangle\in\ker(\overline{\rho}_{A})$ we have $\sqrt{\overline{\rho}_{A}}\boldsymbol{a}^{\dagger}\left|\nu\right\rangle=0$. Therefore $\ker(\overline{\rho}_{A})\subseteq\ker(\boldsymbol{Y})$, then apply Lemma 4. ∎

Appendix B Lindblad form of super-operators

B.1 $\mathcal{L}_{1}$: partial trace of Lindblad super-operators

On the Hilbert space $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ consider the Lindblad super-operator

[TABLE]

with $\boldsymbol{H}$ an Hermitian operator on $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ and with $\boldsymbol{L}_{\mu}$ operators on $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$.

Lemma 5**.**

For any density operator $\rho_{A}$ on $\mathcal{H}_{A}$ the super-operator given by

[TABLE]

remains a Lindbaldian super-operator on $\mathcal{H}_{B}$.

Proof.

By linearity, it is enough to consider the following two cases:

[TABLE]

With $\boldsymbol{H}=\sum_{k}\boldsymbol{A}_{k}\otimes\boldsymbol{B}_{k}$, we have {\rm{tr}}_{A}\Big{(}[\boldsymbol{H},\rho_{A}\otimes\rho_{B}]\Big{)}=\left[\boldsymbol{H}_{B},\rho_{B}\right] where $\boldsymbol{H}_{B}=\sum_{k}{\rm{tr}}(\boldsymbol{A}_{k}\rho_{A})\boldsymbol{B}_{k}$ is an Hermitian operator on $\mathcal{H}_{B}$.

With $\boldsymbol{L}=\sum_{k}\boldsymbol{A}_{k}\otimes\boldsymbol{B}_{k}$, we have

[TABLE]

Since $\rho_{A}\geq 0$, the square matrix M=\Big{(}{\rm{tr}}(\boldsymbol{A}_{k}\rho_{A}\boldsymbol{A}_{k^{\prime}}^{\dagger})\Big{)} is Hermitian an non-negative (use the spectral decomposition of $\rho_{A}$ and the fact that, for any $\left|\psi\right\rangle\in\mathcal{H}_{A}$, the matrix \Big{(}{\rm{tr}}(\boldsymbol{A}_{k}\left|\psi\right\rangle\left\langle\psi\right|\boldsymbol{A}_{k^{\prime}}^{\dagger})\Big{)} is the Gram matrix associated with the vectors $\boldsymbol{A}_{k}\left|\psi\right\rangle$ of $\mathcal{H}_{A}$). Take the Cholesky decomposition $M=NN^{\dagger}$. We have

[TABLE]

with $\boldsymbol{X}_{k^{\prime\prime}}=\sum_{k}N_{k,k^{\prime\prime}}\boldsymbol{B}_{k}$. ∎

B.2 $\mathcal{L}_{2}$: positivity of quadratic form on $\{\boldsymbol{B}_{k},\boldsymbol{B}_{k}^{\dagger}\}$

The proof uses the following property.

Lemma 6**.**

: For any operators $\boldsymbol{X}$ and $\boldsymbol{Y}$ on $\mathcal{H}_{A}$, it holds that

[TABLE]

Proof.

To check this, just write down the expressions, cancel some terms and use $\mathcal{L}_{A}(\overline{\rho}_{A})=0$ for the remaining ones. ∎

Lemma 7**.**

The matrix $X$ in (3.1.2) is always positive semidefinite. In particular, for $\boldsymbol{H}_{\!\rm{int}}=\boldsymbol{A}\otimes\boldsymbol{B}$ with $\boldsymbol{A}$ and $\boldsymbol{B}$ Hermitian, the coefficient ${\rm{tr}}(\boldsymbol{F}\overline{\rho}_{A}\boldsymbol{A}+\boldsymbol{A}\overline{\rho}_{A}\boldsymbol{F}^{\dagger})$ is always nonnegative.

Proof.

The proof for the particular case is simple and introduces the main idea. Apply Lemma 6 with $\boldsymbol{X}=\boldsymbol{Y}=\boldsymbol{F}$; use that $-\mathcal{L}_{A}(\boldsymbol{F}\overline{\rho}_{A})=\boldsymbol{A}\overline{\rho}_{A}-{\rm{tr}}(\boldsymbol{A}\overline{\rho}_{A})\overline{\rho}_{A}$; take the trace over $\mathcal{H}_{A}$ and use that ${\rm{tr}}(\boldsymbol{F}\overline{\rho}_{A})=0={\rm{tr}}(\mathcal{L}_{A})$. This yields:

[TABLE]

where the right hand side is obviously nonnegative.

Applying the same idea with $(\boldsymbol{X},\boldsymbol{Y})=(\boldsymbol{F}_{j},\boldsymbol{F}_{k})$, shows that the components of $X$ can be re-expressed as

[TABLE]

For each $\mu$ and replacing $\overline{\rho}_{A}$ by any pure state $\left|\overline{\psi}\right\rangle\left\langle\overline{\psi}\right|$, we would thus have

[TABLE]

The corresponding matrix $X^{(\mu,\left|\overline{\psi}\right\rangle\left\langle\overline{\psi}\right|)}$ is a Gram matrix, which is always positive semi-definite, see e.g. [4, exercise 1.1.1, page 3]. The same then obviously holds for $X$, which is obtained by taking the sum over $\mu$ and a convex combination over different $\left|\overline{\psi}\right\rangle\left\langle\overline{\psi}\right|$. ∎

Appendix C Computations for the cascaded example

We detail here the computations of the coefficients involved in the second order reduced dynamics of a cavity with a squeezed drive coupled to an unspecified subsystem presented in Section 4.2. The fast dynamics is given by

[TABLE]

with $\kappa>4g$ to guarantee the convergence towards a steady state. The evolution of the operator $\boldsymbol{a}$ in the Heisenberg picture reads:

[TABLE]

where $\mathcal{D}_{\boldsymbol{a}}^{*}$ is the dual of $\mathcal{D}_{\boldsymbol{a}}$. We search $e^{t\mathcal{L}_{A}^{*}}(\boldsymbol{a})$ in the particular form $e^{t\mathcal{L}_{A}^{*}}(\boldsymbol{a})=f_{a}(t)\boldsymbol{a}+h_{a}(t)\boldsymbol{a}^{\dagger}$ for some scalar functions $f_{a}(t)$ and $h_{a}(t)$ verifying $f_{a}(0)=1$ and $h_{a}(0)=0$. Injecting this expression into (41) and solving the differential equations for $f_{a}(t)$ and $h_{a}(t)$ leads to $f_{a}(t)=\frac{e^{-t\frac{\kappa-4g}{2}}+e^{-t\frac{\kappa+4g}{2}}}{2}$, $h_{a}(t)=\frac{e^{-t\frac{\kappa+4g}{2}}-e^{-t\frac{\kappa-4g}{2}}}{2}$ and thus for any initial state $\rho_{0}$, {\rm{tr}}\bigl{(}e^{\infty\mathcal{L}_{A}^{*}}(\boldsymbol{a})\rho_{0}\bigr{)}=0. As the computation of the coefficients $\alpha$ and $\beta$ is similar, we detail here only the calculation of \alpha=\kappa{\rm{tr}}\biggl{(}e^{\infty\mathcal{L}_{A}^{*}}\bigl{(}\overline{\boldsymbol{a}}^{\dagger}\int_{0}^{\infty}e^{t\mathcal{L}_{A}^{*}}(\boldsymbol{a})dt\bigr{)}~{}\rho_{0}\biggr{)}. Using the previous computation and by linearity of the integral we get:

[TABLE]

Then pick $e^{t\mathcal{L}_{A}^{*}}(\boldsymbol{a}^{\dagger}\boldsymbol{a})=f_{N}(t){\boldsymbol{N}}+h_{N}(t)\boldsymbol{a}^{2}+l_{N}(t){\boldsymbol{a}^{\dagger}}^{2}+m_{N}(t)\boldsymbol{I}$, (similarly for $e^{t\mathcal{L}_{A}^{*}}(\boldsymbol{a}^{\dagger}\boldsymbol{a}^{\dagger})$) and standard computations lead to $\alpha=\frac{32\kappa^{2}g^{2}}{((\kappa+4g)(\kappa-4g))^{2}}$.