# Explicit formulas for adiabatic elimination with fast unitary dynamics

Angela Riva, Alain Sarlette, Pierre Rouchon

► To cite this version:

Angela Riva, Alain Sarlette, Pierre Rouchon. Explicit formulas for adiabatic elimination with fast unitary dynamics. CDC 2024 - IEEE Conference on Decision and Control, IEEE-CSS, Dec 2024, Milan, Italy. hal-04530303

HAL Id: hal-04530303

https://inria.hal.science/hal-04530303

Submitted on 3 Apr 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Explicit formulas for adiabatic elimination with fast unitary dynamics

Angela Riva, Alain Sarlette, and Pierre Rouchon

The so-called “adiabatic elimination” of fast decaying degrees of freedom in open quantum systems can be performed with a series expansion in the timescale separation. The associated computations are significantly more difficult when the remaining degrees of freedom (center manifold) follow fast unitary dynamics instead of just being slow. This paper highlights how a formulation with Sylvester’s equation and with adjoint dynamics leads to systematic, explicit expressions at high orders for settings of physical interest.

## I. INTRODUCTION

Using quantum effects in practical technological platforms relies on balancing two seemingly contradictory requirements: maintaining isolation from the environment to protect the fragile quantum coherences, and simultaneously allowing for fast manipulation towards performing operations. This is the fundamental problem that quantum control seeks to solve. The specific setups engineered in this context are often composed of several components, whose dimensions multiply and span several timescales. The “useful quantum information” represents only a small, ideally non-decaying part of this system. Fast decaying degrees of freedom (DOFs) are thus associated to auxiliary control signals, as e.g. when stabilizing a quantum system via dissipation engineering [15], [11], [12]. Model reduction techniques would thus seek to summarize the effect of these fast stabilizing DOFs on the dynamics of the target variables, preferably with operational principles and analytic expressions to guide system design.

The linearity of the master equation for Markovian open quantum systems allows, in principle, a spectral decomposition very much like Hamiltonian (block-)diagonalization. Specific approaches in the physics literature rather focus on interpreting dissipative systems as being subject to uncertain events [1, 4, 16]. While this may be natural in many situations, it fails to provide a general and systematic approach. Our research group has initiated a line of work with a more mathematical approach: identifying, via series expansion still, the variables and dynamics corresponding to an invariant subspace where dynamics is at the slowest timescale. The approach was applied to perform model reduction inside a single system (Cartesian product, [3]) and in most general composite systems (tensor product, [2]). It has recently been translated to the Heisenberg picture, with particular benefits for multi-partite systems [9, 10].

All previously mentioned studies address scenarios where a rapid relaxation, characterized by eigenvalues with negative real parts, surrounds a slower evolution on an invariant subspace associated with eigenvalues near zero. However, like in the center manifold theory, the same model reduction technique should also admit a rapid relaxation towards an invariant subspace featuring large, almost purely imaginary eigenvalues (unitary dynamics). This would fit practical situations where e.g. the stabilized state keeps turning at significant characteristic frequencies, prefiguring in/out of resonance control [13]. Such a case was first treated in [5, 6], establishing the abstract expressions up to second-order expansion, under technical conditions on the fast unitary dynamics. In the present work, we provide integral formulas removing these conditions (Section III). We then highlight how, for typical systems of interest — namely, eliminating a rapidly decaying harmonic oscillator environment — these integral formulas also lead more easily to interpretable explicit expressions for the long-term dynamics on the invariant manifold (see the physical examples of section IV). The two key technical ingredients are the use of stable Sylvester equation solutions (Proposition 1) and the treatment of the typical dynamics on the adjoint state (Proposition 2).

## II. MODEL AND METHODS

This section recalls the setting and adiabatic elimination approach as used e.g. in [2, 6, 9, 18].

We consider a density operator $\rho_t$ at time $t$, defined over a bipartite Hilbert space $\mathcal{H} = \mathcal{H}_A \otimes \mathcal{H}_B$. The time evolution of this operator is described by the linear dynamics:

$$ \dot{\rho}_t = \mathcal{L}(\rho_t) = -i(H_A \otimes I_B)^\times(\rho_t) + \mathcal{I}_A \otimes \mathcal{L}_B(\rho_t) - ig[H_I, \rho_t]. \quad (1) $$

Here $I_M$ is the identity operator and $\mathcal{I}_M$ the identity superoperator, associated to the Hilbert space $\mathcal{H}_M$, and we use the notation $H^\times(\bullet) := [H, \bullet]$ to represent the commutator with operator $H$, i.e. $[H, X] = HX - XH$. The Lindbladian term for the mode $B$ takes the general form:

$$ \mathcal{L}_B(\bullet) = -iH_B^\times(\bullet) + \kappa \sum_k \mathcal{D}[L_k](\bullet), \quad (2) $$
with $\kappa$ the damping rate and $\mathcal{D}[L](\bullet) := L \bullet L^\dagger - \frac{1}{2}(L^\dagger L \bullet + \bullet L^\dagger L)$ the dissipation channel associated to operator $L$. The coupling between $\mathcal{H}_A$ and $\mathcal{H}_B$ is weak and Hamiltonian:

$$gH_I = g \sum_k A_k \otimes B_k$$

with $A_k$ and $B_k$ operators on $\mathcal{H}_A$ and $\mathcal{H}_B$ respectively and $g$ a small positive parameter.

The adiabatic elimination procedure is based on the assumption that the system exhibits a clear separation of timescales, allowing for a perturbative expansion in the small parameter $\epsilon = g/\kappa \ll 1$. The dynamics for $g = \epsilon = 0$ is assumed “relatively easy to compute”, with $\mathcal{L}_B$ converging exponentially towards an equilibrium state $\bar{\rho}_B$ and subsystem A independently rotating with a Hamiltonian $H_A$; thus, $\mathcal{L}(\rho)$ then features an invariant subspace $S_0$ spanned by $\{\rho = \rho_s \otimes \bar{\rho}_B, \text{ for any } \rho_s \text{ on } \mathcal{H}_A\}$ and on which all eigenvalues are purely imaginary (center manifold). For $\epsilon \neq 0$, The interaction Hamiltonian $H_I$ acts as a perturbation on this situation. The goal is then to efficiently compute the “long-term dynamics”, taking place on $S_\epsilon$, an invariant subspace $\epsilon$-close to $S_0$ and characterized by the eigenoperators of $\mathcal{L}(\rho)$ associated to eigenvalues with real part close to zero.

When $\epsilon \neq 0$, the interaction Hamiltonian creates correlations between the two subsystems A and B. As a result, the invariant subspace $S_\epsilon$ does not take the same form as $S_0$ anymore. Still, it is convenient to parameterize states on $S_\epsilon$ by using $\rho_s \otimes \bar{\rho}_B \in S_0$ as coordinates. Indeed, the dynamics are then reduced to $\rho_s$ whose evolution can be interpreted in comparison to the case $\epsilon = 0$. In order to characterize the long-term behavior of the system, we then aim to find two linear, time independent maps: $\mathcal{L}_s$, describing the dynamics of the reduced system via its coordinates $\dot{\rho}_s = \mathcal{L}_s(\rho_s)$; and $\mathcal{K}$, mapping $\rho_s$ to the solution $\rho = \mathcal{K}(\rho_s) \in S_\epsilon$ of the complete dynamics (1). Expressing that the evolution of $\rho_s$ must mirror the one of $\mathcal{K}(\rho_s)$ leads to the condition:

$$\mathcal{K}(\mathcal{L}_s(\rho_s)) = \mathcal{L}(\mathcal{K}(\rho_s)). \qquad (3)$$

Following usual approximation theory, we solve (3) by expanding the maps $\mathcal{L}_s$ and $\mathcal{K}$ in powers of $\epsilon \ll 1$,

$$\mathcal{L}_s(\rho_s) = \sum_{j=0}^{\infty} \epsilon^j \mathcal{L}_{s,j}(\rho_s), \quad \mathcal{K}(\rho_s) = \sum_{j=0}^{\infty} \epsilon^j \mathcal{K}_j(\rho_s), \qquad (4)$$

and requiring to satisfy the condition (3) separately at each order $\epsilon^j$. Since the lowest-order contributions are the dominant ones, usually only the solutions $\mathcal{L}_{s,j}, \mathcal{K}_j$ for the first few orders are computed explicitly, e.g. $j \le 2$. The solution to (3) is certainly not unique, as it involves a coordinate choice on $S_\epsilon$, expressed via $\mathcal{K}$.

The most natural solution to (3) at order $\epsilon^0$, is the one corresponding to $\epsilon = 0$, thus

$$\mathcal{L}_{s,0}(\rho_s) = -iH_A^\times(\rho_s), \quad \mathcal{K}_0(\rho_s) = \rho_s \otimes \bar{\rho}_B. \qquad (5)$$

In this paper, we further specify the coordinate choice by imposing the *partial trace gauge*, namely $\rho_s = \text{Tr}_B(\rho) = \text{Tr}_B(\mathcal{K}(\rho_s))$ at all orders of approximation, where $\text{Tr}_B$ denotes the partial trace over $\mathcal{H}_B$. This choice implies that $\text{Tr}_B \mathcal{K}_1 = \text{Tr}_B \mathcal{K}_2 = \dots = 0$.

Note that when $H_A = 0$, we get $\mathcal{L}_{s,0} = 0$ and this simplifies all further orders (see e.g. [2]). Our contributions are meant to address the difficulties appearing in solving (3),(4) at orders $\epsilon^1, \epsilon^2$ for $H_A \neq 0$.

### III. RESULT

Obtaining according to (3),(4) the elements $\mathcal{L}_{s,0}, \mathcal{K}_0$ (see (5)) and

$$\epsilon \mathcal{L}_{s,1}(\rho_s) = -ig \sum_k \text{Tr}(B_k \bar{\rho}_B) [A_k, \rho_s] \qquad (6)$$

(see proof of Proposition 1 below) poses no particular problem. Our results concern the next order. In Proposition 1, we derive a new integral closed form for the first order map $\mathcal{K}_1$ by solving a Sylvester equation. In Proposition 2, the second-order reduced dynamics $\mathcal{L}_{s,2}$, involving the partial trace of the commutator between $H_I$ and $\mathcal{K}_1$ is re-expressed with this form, putting the integral part on a dual operator which is more systematically manageable in typical situations. Examples indeed illustrate the usefulness of these formulas in Section IV. For ease of notation, we start by defining some terms.

**Definition 1.** The operator $B_{0,k}$, acting on $\mathcal{H}_B$, is defined as

$$B_{0,k} = B_k - \text{Tr}_B(B_k \bar{\rho}_B) I_B. \qquad (7)$$
**Definition 2.** The operator $A_k^-(t)$, acting on $\mathcal{H}_A$ is defined as

$$A_k^-(t) = e^{-itH_A^\times}(A_k) = e^{-itH_A}A_k e^{itH_A}. \qquad (8)$$

(Note that this is the opposite propagation to the Heisenberg picture evolution under $H_A$, see Remark 1.)

**Proposition 1.** Consider model (1) and the adiabatic elimination expansion of Section II with partial trace coordinates. The first order map $K_1$ is given by

$$\epsilon K_1(\rho_s) = -ig \sum_k \int_0^{+\infty} \{A_k^-(t)\rho_s \otimes e^{t\mathcal{L}_B}(B_{0,k}\bar{\rho}_B) - \rho_s A_k^-(t) \otimes e^{t\mathcal{L}_B}(\bar{\rho}_B B_{0,k})\} dt. \qquad (9)$$

*Proof.* By inserting the expansions (4) in the general invariance equation (3) and grouping terms of order $\epsilon$, we obtain the condition

$$-i[H_A \otimes \mathcal{I}_B, \epsilon K_1(\rho_s)] + \mathcal{I}_A \otimes \mathcal{L}_B(\epsilon K_1(\rho_s)) - ig[H_I, K_0(\rho_s)] = K_0(\epsilon \mathcal{L}_{s,1}(\rho_s)) + \epsilon K_1(\mathcal{L}_{s,0}(\rho_s)). \qquad (10)$$

This equation has two unknowns: $\mathcal{L}_{s1}$ and $K_1$. Two terms involving $K_1$ here are due to $H_A \neq 0$ and therefore new with respect to [2]. Nevertheless, thanks to the coordinate choice with $\mathrm{Tr}_B(K_1) = 0$, we can much like in [2] apply $\mathrm{Tr}_B$ to the condition (10) to directly isolate the first order reduced dynamics in the form (6). To get this, we further observe that $\mathrm{Tr}_B(\mathcal{L}_B) = 0$ and $\mathrm{Tr}_B(H_A^\times(\bullet)) = H_A^\times(\mathrm{Tr}_B(\bullet))$.

Inserting (6) into (10) yields the expression:

$$-i[H_A \otimes \mathbb{I}_B, \epsilon K_1(\rho_s)] + \mathcal{I}_A \otimes \mathcal{L}_B(\epsilon K_1(\rho_s)) + i\epsilon K_1([H_A, \rho_s]) = ig \sum_k [A_k \otimes B_{0,k}, \rho_s \otimes \bar{\rho}_B], \qquad (11)$$

which should be solved for $K_1$. At the level of superoperators, (11) is in fact a Sylvester equation of the form $AX+XB = C$, where $\epsilon K_1$ plays the role of the unknown term $X$. To be explicit: $A \to -i(H_A \otimes \mathbb{I}_B)^\times(\bullet) + \mathcal{I}_A \otimes \mathcal{L}_B(\bullet)$, $B \to iH_A^\times(\bullet)$, $C \to ig \sum_k [A_k \otimes B_{0,k}, \bullet \otimes \bar{\rho}_B]$. The solution to this Sylvester equation can be written as

$$X = - \int_0^{+\infty} e^{tA} C e^{tB} dt, \qquad (12)$$

provided

$$e^{tA} C e^{tB} \to_{t \to \infty} 0. \qquad (13)$$

Using that (super)operators acting on different subsystems commute, as well as the explicit expression $e^{iH^\times}(Q) = e^{iH}Qe^{-iH}$, we obtain:

$$\begin{aligned} \frac{e^{tA} C e^{tB}(\rho_s)}{g} &= i \sum_k e^{t((-iH_A \otimes \mathbb{I}_B)^\times + \mathcal{I}_A \otimes \mathcal{L}_B)} [A_k \otimes B_{0,k}, e^{itH_A^\times}(\rho_s) \otimes \bar{\rho}_B] \\ &= i \sum_k \left( e^{-itH_A^\times}(A_k)\rho_s \otimes e^{t\mathcal{L}_B}(B_{0,k}\bar{\rho}_B) - \rho_s e^{-itH_A^\times}(A_k) \otimes e^{t\mathcal{L}_B}(\bar{\rho}_B B_{0,k}) \right). \end{aligned}$$

This last expression corresponds to the proposition statement. Since $e^{t\mathcal{L}_B(\bullet)} \to_{t \to \infty} \mathrm{Tr}(\bullet)\bar{\rho}_B$, we have $\lim_{t \to \infty} e^{t\mathcal{L}_B(B_{0,k}\bar{\rho}_B)} = \lim_{t \to \infty} e^{t\mathcal{L}_B(\bar{\rho}_B B_{0,k})} = 0$, such that condition (13) is satisfied. $\square$

**Remark 1.** Note that $A_k^-(t)$ in the statement of Proposition 1 follows the opposite dynamics to the Heisenberg evolution under $H_A$. We can understand this as follows. According to $\mathcal{L}_{s,0}$, $\rho_s$ already rotates with $H_A$. Inside the integral of (9), we must replace it by $(A_k \rho_s)$ rotating with $H_A$, much like the second tensor factor takes $(B_{0,k} \bar{\rho}_B)$ evolving under $\mathcal{L}_B$.

Compared to previous work [5, 6], The proposition 1 avoids any technical conditions and the need to solve several equations by treating the components of $H_A$ individually. In turn, it does not guarantee positivity — whose importance has anyways been re-evaluated since [18] — and it leaves a propagator equation to solve, for all $t$ and separately for each subsystem. Propagation on subsystem A, with a fixed Hamiltonian $H_A$, is usually no big deal. Often, $H_A$ can be diagonalized efficiently or commuted through typical operators to compute $A_k^-(t)$. Propagation on subsystem B, with a Lindbladian superoperator, is computationally less convenient. Indeed, even for simple steady-state cases like convergence towards a thermally broadened vacuum, it is not trivial to express the state at any time along the trajectory starting from rather arbitrary initial states like $B_{0,k}\bar{\rho}_B$. The following result mitigates this difficulty.
**Definition 3.** The operator $B_k(t)$, evolving under the action of $\mathcal{L}_B$ in the Heisenberg picture, and acting on $\mathcal{H}_B$, is defined as

$$B_k(t) = e^{t\mathcal{L}_B^*}(B_k), \qquad (14)$$

with the adjoint of the Lindbladian $\mathcal{L}^*$ defined in general as $\mathcal{L}^*(\bullet) = +i[H, \bullet] + \kappa \sum_k (L_k^\dagger \bullet L_k - \frac{1}{2}(L_k^\dagger L_k \bullet + \bullet L_k^\dagger L_k))$.

**Proposition 2.** Consider model (1) and the adiabatic elimination expansion of Section II with partial trace coordinates. The second order reduced dynamics is given by:

$$\epsilon^2 \mathcal{L}_{s,2} = -\frac{g^2}{\kappa} \sum_{k,l} \int_0^{+\infty} \left\{ \mathrm{Tr}\{B_l(t) B_{0,k} \bar{\rho}_B\} [A_l, A_k^-(t) \rho_s] - \mathrm{Tr}\{B_l(t) \bar{\rho}_B B_{0,k}\} [A_l, \rho_s A_k^-(t)] \right\} \kappa dt. \qquad (15)$$

*Proof.* In the second order invariance condition:

$$-i[H_A \otimes \mathcal{I}_B, \epsilon^2 \mathcal{K}_2(\rho_s)] + \mathcal{I}_A \otimes \mathcal{L}_B(\epsilon^2 \mathcal{K}_2(\rho_s)) - i[gH_I, \epsilon \mathcal{K}_1(\rho_s)] = \mathcal{K}_0(\epsilon^2 \mathcal{L}_{s,2}(\rho_s)) + \epsilon^2 \mathcal{K}_1(\mathcal{L}_{s,1}(\rho_s)) + \epsilon^2 \mathcal{K}_2(\mathcal{L}_{s,0}(\rho_s)), \qquad (16)$$

take the partial trace over subsystem B, like we did before Proposition 1 towards getting $\mathcal{L}_{s,1}$, in order to obtain:

$$\epsilon^2 \mathcal{L}_{s,2}(\rho_s) = -i \mathrm{Tr}_B\{[gH_I, \epsilon \mathcal{K}_1(\rho_s)]\}.$$

Substituting the expression of $H_I$ and of $\mathcal{K}_1$ from Proposition 1 yields:

$$\epsilon^2 \mathcal{L}_{s,2}(\rho_s) = (-ig)^2 \sum_{l,k} \left\{ \left[ A_l, \int_0^{+\infty} c_{k,l}(t) A_k^-(t) dt \rho_s \right] - \left[ A_l, \rho_s \int_0^{+\infty} \tilde{c}_{k,l}(t) A_k^-(t) dt \right] \right\}, \qquad (17)$$

where

$$\begin{align}
c_{k,l}(t) &= \mathrm{Tr}(B_l e^{t\mathcal{L}_B}(B_{0,k}\bar{\rho}_B)) = \mathrm{Tr}(e^{t\mathcal{L}_B^*}(B_l)B_{0,k}\bar{\rho}_B) \tag{18} \\
\tilde{c}_{k,l}(t) &= \mathrm{Tr}(B_l e^{t\mathcal{L}_B}(\bar{\rho}_B B_{0,k})) = \mathrm{Tr}(e^{t\mathcal{L}_B^*}(B_l)\bar{\rho}_B B_{0,k}). \nonumber
\end{align}$$

On the right-hand side of (18), we have transferred the Lindblad dynamics to the adjoint, in other words evolving the partner operator inside the trace in Heisenberg picture. Recalling Definition 3, this corresponds to the statement. □

**Remark 2.** The Proposition does not claim that the resulting reduced dynamics preserves positivity — thus taking the typical Lindblad form with positive dissipation rates. This property is proven in [2] for $H_A = 0$, and in [5] under more technical conditions.

The essential part of Proposition 2 is to replace, when computing $\mathcal{L}_{s,2}$, the Lindbladian trajectory starting at $(B_{0,k}\bar{\rho}_B)$ by an adjoint Lindbladian trajectory on $B_l$. As already noted in [7], the latter is often much easier to compute in typical situations. For instance, the operator propagation can be computed with little more effort than a classical system, when both $\mathcal{L}_B$ and $B_l$ on subsystem B correspond to a so-called linear quantum system [14].

The strategy followed by these two Propositions serves to address the same major issue in resolving the series expansion (3),(4), namely the inversion (in matrix inverse sense) of the fast dynamics. When pursuing the expansion at higher orders, the main computational issue remains the same, and we can reiterate the same procedure to obtain explicit expressions.

**Proposition 3.** Consider model (1) and the adiabatic elimination expansion of Section II with partial trace coordinates. Assume that the terms of the series expansion have been computed up to $n-1$. Then, for all $n \ge 1$,

$$\epsilon \mathcal{L}_{s,n}(\rho_s) = -ig \sum_k [A_k, \mathrm{Tr}_B((I_A \otimes B_k)\mathcal{K}_{n-1}(\rho_s))] \qquad (19)$$

and $\epsilon \mathcal{K}_n$ can be computed with an integral formula similar to (9).

*Proof.* The invariance condition at order $\epsilon^n$ writes:

$$\epsilon \mathcal{K}_0(\mathcal{L}_{s,n}(\rho_s)) = \epsilon \left(-i(H_A \otimes I_B)^\times + \mathcal{I}_A \otimes \mathcal{L}_B\right) (\mathcal{K}_n(\rho_s)) - \epsilon \mathcal{K}_n(\mathcal{L}_{s,0}(\rho_s)) - \epsilon \sum_{m=1}^{n-1} \mathcal{K}_m(\mathcal{L}_{s,n-m}(\rho_s)) - ig(H_I^\times)(\mathcal{K}_{n-1}(\rho_s)). \qquad (20)$$
The last line contains only known terms at stage n. Like in Proposition 1, by taking the partial trace over subsystem B, all the terms involving $K_n$ vanish with the coordinate choice $\rho_s = \text{Tr}_B(\rho)$, while $\text{Tr}_B K_0 = \mathcal{I}_s$. The terms in $K_m$ also vanish by partial trace. We thereby obtain the explicit expression (19) for $\mathcal{L}_{s,n}$, assuming all previous orders were known.

After this, the left-hand side of (20) is also known and the equation takes the Sylvester form $AX + XB = C$ like in the proof of Prop.1, with A and B unchanged, $\epsilon K_n$ playing the role of X, and C containing all known terms. There remains to show that the integral form solution (12) remains stable at this order. Like in the proof of Prop.1, this is ensured thanks to $e^{tA}(\bullet) \to_{t \to \infty} U_{t,A} \text{Tr}_B(\bullet)\bar{\rho}_B$ with $U_{t,A}$ a unitary evolution, and the annihilation of terms under $\text{Tr}_B$. $\square$

When one is only interested in the dynamics $\mathcal{L}_{s,n}$, it may be possible to also apply the trick of Proposition 2 iteratively at all orders. This would simplify computations significantly, compared to computing each $K_m$ explicitly. Investigating this possibility is left for future work.

## IV. EXAMPLES

In this section, we give some practical demonstrations of the applicability of the results presented in Section III. The starting point is a harmonic oscillator B subsystem undergoing so-called quantum linear dynamics [14]:

$$ \mathcal{L}_B(\bullet) = -i\omega_B[b^\dagger b, \bullet] + \kappa_\phi \mathcal{D}[b^\dagger b](\bullet) + \kappa(1 + n_{\text{th}})\mathcal{D}[b](\bullet) + \kappa n_{\text{th}}\mathcal{D}[b^\dagger](\bullet), \quad (21) $$

with $b(b^\dagger)$ the bosonic annihilation (creation) operator, $b^\dagger b$ the number operator, $\omega_B$ the undamped harmonic oscillator frequency, $\kappa$ the damping rate for the dissipator $\mathcal{D}$, and $n_{\text{th}}$ the residual thermal excitation. Note that we treat both the single photon loss relaxation channel $\mathcal{D}[b]$, and the noise induced by thermal excitation $\mathcal{D}[b^\dagger]$. We also note that the unique steady state such that $\mathcal{L}_B(\rho) = 0$ is the thermal state,

$$ \bar{\rho}_B = \sum_{n=0}^{+\infty} \frac{n_{\text{th}}^n}{(n_{\text{th}} + 1)^{n+1}} |n\rangle\langle n| \quad (22) $$

in the number-basis $|n\rangle$, i.e. where $b^\dagger b = \sum_{n=0}^{+\infty} n |n\rangle\langle n|$.

We couple this system to an arbitrary A subsystem, using a standard dipolar coupling which is also linear:

$$ \tilde{H}_I = \tilde{A} \otimes (b + b^\dagger). \quad (23) $$

In the vast majority of cases [8], the timescale separation between $\omega_B$ and $g$ is so large that significant coupling effects on A only happen when $H_A$ contains frequencies close to $\omega_B$. This is what lies behind resonance conditions invoked in more high-level descriptions of physical phenomena. More precisely, going to a rotating frame with $\omega_B$ and averaging the resulting time-dependent dynamics singles out the relevant dominant long-term behavior [17]. Note that this is mathematically justified only if $\omega_B \gg \kappa$ as well. In physics terms [8], the resonance has a width of order $\kappa$, thus frequency selection only works for $\kappa \ll \omega_B$.

With the framework developed here, we are able to compare the results of adiabatic elimination in both cases: [Section IV A] The physics traditional way, i.e. eliminating the large almost-resonant part of $H_A$ and of $\omega_B b^\dagger b$ by first averaging in a rotating frame;

[Section IV B] Maintaining all terms in so-called “lab frame”, with thus $\omega_B$ and $H_A$ of the same order and non-negligible compared to $\kappa$. The results should match when $\omega_B$ and $H_A$ are close and large compared to $\kappa$.

### A. Jaynes-Cummings interaction

The Jaynes-Cummings interaction, see e.g. [8, Chap.3.4], corresponds to the canonical theoretical and practical case of coupling subsystem B with a (quasi-) resonant qubit A. After first-order averaging in a frame where both subsystems rotate at the resonance frequency, the interaction boils down to an energy exchange term.

**Example 1.** Consider the bipartite system composed of a harmonic oscillator and a two-level system (qubit), whose dynamics is described by

$$ \begin{aligned} \mathcal{L}_B(\bullet) &= -i\Delta[b^\dagger b, \bullet] + \kappa(1 + n_{\text{th}})\mathcal{D}[b](\bullet) + \kappa n_{\text{th}}\mathcal{D}[b^\dagger](\bullet) + \kappa_\phi\mathcal{D}[b^\dagger b](\bullet), \\ H_A &= 0, \\ H_I &= \sigma_+ \otimes b + \sigma_- \otimes b^\dagger, \end{aligned} $$
with $\mathcal{H}_A = \text{span}\{|g\rangle, |e\rangle\}$ and $\sigma_- = (\sigma_+)^\dagger = |g\rangle\langle e|$ the qubit lowering operator (from excited to ground state). With Proposition 2, assuming $\kappa \gg g$, the second order reduced dynamics is given by

$$ \frac{d}{dt}\rho_s = \epsilon^2 \mathcal{L}_{s,2}(\rho_s) = -i(1+2n_{th})\frac{4\Delta g^2}{|\gamma|^2} \left[\frac{\sigma_z}{2}, \rho_s\right] + (1+n_{th})(\kappa + \kappa_\phi)\frac{4g^2}{|\gamma|^2}\mathcal{D}[\sigma_-](\rho_s) + n_{th}(\kappa + \kappa_\phi)\frac{4g^2}{|\gamma|^2}\mathcal{D}[\sigma_+](\rho_s), \quad (24) $$

where $\sigma_z = |g\rangle\langle g| - |e\rangle\langle e|$ and $\gamma = \kappa + \kappa_\phi + 2i\Delta$.

**Solution.** The coupling corresponds to our general setting with $B_1 = b$, $B_2 = b^\dagger$, $A_1 = \sigma_+$, $A_2 = \sigma_-$. The zeroth-order dynamics $\mathcal{L}_{s,0}$ vanishes since $H_A = 0$. Note that $\mathcal{L}_B$ in (24) still features the unique steady state (22), and since $\langle n | b | n \rangle = 0$ for all $n$ we also have $\mathcal{L}_{s,1} = 0$. To recover the second order reduced dynamics, we start by computing the Heisenberg representation of the coupling operators (Definition 3):

$$ \frac{d}{dt}b = \mathcal{L}_B^*(b) = -\frac{\gamma}{2}b, \quad \text{thus} \quad b(t) = e^{-\frac{\gamma}{2}t}b(0), \quad (25) $$

with $\gamma$ defined as in the statement. In this step, we highly benefit from assuming a linear quantum system [14] on B in order to obtain such easy closed-form solution. Indeed, for a general quantum system B, the computation of (14) in Definition 3 can be a significant bottleneck.

From there, with easy calculations involving the geometric series, we obtain the coefficients $c_{k,l}(t)$ and $\bar{c}_{k,l}(t)$ mentioned in (18). Since $H_A = 0$, the operators $A_k^-(t)$ are time independent, so the integral to be solved in Proposition 2 is simply $\int_0^{+\infty} e^{-s\gamma}ds = \frac{1}{\gamma}$. By rearranging the terms, $\mathcal{L}_{s,2}$ can be put in standard Lindblad form to obtain (24). $\square$

The reduced model (24) thus contains a unitary shift with B pulling on the frequency of A (Hamiltonian in $\sigma_z$); and it translates the thermal dissipations on B into corresponding thermal dissipations in $\sigma_-$ and $\sigma_+$ on the two-level system A, yet with decreasing effect as $\kappa$ increases. This result is well-known [2]. It readily generalizes to a subsystem A of higher dimension. One just replaces, in the model and in the resulting dissipator, $\sigma_-$ by A and $\sigma_+$ by $A^\dagger$; in the resulting slow Hamiltonian, one replaces $(1+2n_{th})\frac{\sigma_-}{2}$ by $(n_{th}AA^\dagger - (1+n_{th})A^\dagger A)$.

## B. Treatment with fast unitary dynamics in lab frame

We move on to consider now the same setting in lab frame, i.e. without going to a rotating frame and averaging, but thus with remaining Hamiltonian $\hat{H}_A$ (and $\omega_B b^\dagger b$) not small compared to $\kappa$. The model thus corresponds to (21), (23), with for a qubit $\hat{A} = (\sigma_- + \sigma_+) = \sigma_x$ and $H_A$ proportional to $\sigma_z$.

**Example 2.** Consider the bipartite system whose dynamics is described by

$$ \begin{aligned} \mathcal{L}_B(\bullet) &= -i\omega_B[b^\dagger b, \bullet] + \kappa(1+n_{th})\mathcal{D}[b](\bullet) + \kappa n_{th}\mathcal{D}[b^\dagger](\bullet) + \kappa_\phi\mathcal{D}[n](\bullet), \\ H_A &= -\omega_{eg}\frac{\sigma_-}{2}, \\ H_I &= \sigma_x \otimes (b+b^\dagger), \end{aligned} $$

with the energy gap $\omega_{eg} = \omega_e - \omega_g$, $\sigma_x = |e\rangle\langle e| + |g\rangle\langle g|$, $\sigma_z = |g\rangle\langle g| - |e\rangle\langle e|$. With proposition 2, assuming $\kappa \gg g$, the second order reduced dynamics is given by

$$ \frac{d}{dt}\rho_s = \mathcal{L}_{s,0}(\rho_s) + \epsilon^2 \mathcal{L}_{s,2}(\rho_s) = -i[-\omega_{eg}\frac{\sigma_-}{2}, \rho_s] - ig^2Y\left[\frac{\sigma_z}{2}, \rho_s\right] + g^2 \sum_{\ell, \ell' \in \{+, -\}} X_{\ell\ell'}\left(\sigma_{\ell'}\rho_s\sigma_\ell^\dagger - \frac{\rho_s\sigma_\ell^\dagger\sigma_{\ell'} + \sigma_\ell^\dagger\sigma_{\ell'}\rho_s}{2}\right), \quad (26) $$

with the hermitian matrix $X$ and coefficient $Y$ defined by:

$$ \begin{aligned} X_{\ell\ell'} &= r_{\ell'} + r_{\ell}^* + e_{\ell'} + e_{\ell}^* \\ Y &= \frac{1}{2i}(r_+ + e_+ - r_+^* - e_+^* - r_- - e_- + r_-^* + e_-^*), \end{aligned} \quad (27) $$

and the coefficients $r_\ell$ and $e_\ell$, $\ell \in \{+, -\}$:

$$ r_{\pm} = \frac{2(1 + n_{th})}{\gamma_{\pm}}, \quad e_{\pm} = \frac{2n_{th}}{\gamma_{\mp}^*}, \quad (28) $$

where $\gamma_{\pm} = \kappa + \kappa_\phi + 2i(\omega_B \pm \omega_{eg})$.
**Solution.** The coupling corresponds to the above setting with just $A_1 = \sigma_x$, $B_1 = b+b^\dagger$. The zero order, given by Eq 5, is $\mathcal{L}_{s,0} = -i\omega_{eg}\sigma_z^\times$. The first order $\mathcal{L}_{s,1} = 0$ for the same reason as in the previous example. To obtain the second order reduced dynamics via Definition 3, the time dependence for $b(t)$ is the same as in (25), with $\omega_B$ replacing $\Delta$ in the definition of $\gamma$. With it one can easily obtain according to (18) the sole coefficients $c_{11}(t) = (n_{th}+1)e^{-\gamma t/2}+n_{th}e^{-\gamma^* t/2}$ and $\tilde{c}_{11}(t) = c_{11}(t)^*$. Since we just have $A_1^-(t) = \sigma_+e^{-i\omega_{eg}t} + \sigma_-e^{+i\omega_{eg}t}$, we compute the time integrals

$$ \int_0^{+\infty} c_{11}(t)A_1^-(t)dt = \sum_{\ell \in \{+, -\}} (r_\ell + e_\ell)\sigma_\ell, \\ \int_0^{+\infty} \tilde{c}_{11}(t)A_1^-(t)dt = \sum_{\ell \in \{+, -\}} (r_\ell^* + e_\ell^*)\sigma_\ell^\dagger. \qquad (29) $$

The statement then follows after straightforward calculations and relabeling. $\square$

The result generalizes to higher-dimensional subsystems A without much difficulties. Replace $\sigma_x$ by $(A+A^\dagger)$ and $H_A = \sum_{s=1}^{d_A} \omega_s |s\rangle\langle s|$ in the problem statement. The computations of $c_{11}(t)$ and $\tilde{c}_{11}(t)$ involving exponentials remains unchanged. The time evolution of the operators $A_k^-(t)$ will be given by

$$ A_k^-(t) = \sum_{n,m=1}^{d_A} e^{-i(\omega_n - \omega_m)t} \langle n | A_k | m \rangle |n\rangle\langle m|. \qquad (30) $$

With these expressions, we can thus in principle easily compute the integrals

$$ \int_0^{+\infty} c_{k,l}(t)A_k^-(t)dt \quad \text{and} \quad \int_0^{+\infty} \tilde{c}_{k,l}(t)A_k^-(t)dt, \qquad (31) $$

and obtain the second order reduced dynamics.

**Interpretation.** The following observations are in order about the result in (26).

* One checks that for $\omega_{eg} = \omega_B$ large, averaging in a frame rotating with $U(t) = e^{i\omega_{eg}t\sigma_z/2}$, yields back the result of Example 1 as expected. Without averaging, terms in $X_{+-}$ and $X_{-+}$ remain. Note though that the averaging condition now relaxes to $\omega_B, \omega_{eg} \gg e^2\kappa$, thanks to confining ourselves to a manifold with slow dynamics.

* For $\omega_{eg} = 0$ instead, we have $r_+ = r_-$ and $e_+ = e_-$ such that $Y = 0$ and $X$ is proportional to the all-ones matrix. This singularity implies a single dissipation channel, in $\sigma_x$ i.e. proportional to the coupling operator, in agreement with the result of [2] when $H_A = 0$.

* The formula (26) thus allows us to capture all intermediate scaling cases.

* The matrix X has a positive trace, but its determinant is independent of $n_{th}$ and negative as soon as $\omega_{eg} \neq 0$. Then (26) cannot be put in completely positive Lindblad form, as also happened in [18].

* The corresponding Bloch equations for $\rho = \frac{I+x\sigma_x+y\sigma_y+z\sigma_z}{2}$ are:

$$ \begin{aligned} \dot{x} &= \left( \omega_{eg} - g^2 Y \frac{2n_{th}}{1+2n_{th}} \right) y \\ \dot{y} &= \left( -\omega_{eg} + g^2 Y \frac{2+2n_{th}}{1+2n_{th}} \right) x - g^2 (1+2n_{th}) r_z y \\ \dot{z} &= -g^2 (1+2n_{th}) r_z (z - \bar{z}) \end{aligned} \qquad (32) $$

where $r_z = 4(\kappa + \kappa_\phi)(\frac{1}{|\gamma_+|^2} + \frac{1}{|\gamma_-|^2})$ and $\bar{z} = \frac{|\gamma_+|^2 - |\gamma_-|^2}{(|\gamma_+|^2 + |\gamma_-|^2)(1+2n_{th})}$. For $\omega_{eg} \neq 0$ this system is exponentially stable and converges to $x = y = 0, z = \bar{z}$. For $\omega_{eg} = 0$, the $x$ coordinate remains invariant and the two others converge exponentially to $y = 0$ and $z = \bar{z}$.

* The Bloch equations (32) can thus be seen as coordinates, accurately describing the completely positive Lindblad form (26) when restricted to an invariant subspace inside $\mathcal{H}_A \otimes \mathcal{H}_B$. In [18], a necessary and sufficient inequality on the spectrum of Bloch equations is presented to decide whether another 3-dimensional coordinate choice (thus relaxing the choice $\rho_s = \text{Tr}_A(\rho)$) may yield a completely positive reduced model, identifiable with a qubit. Remarkably, for any parameter values, our particular result (32) appears to lie on the boundary of these inequalities. Higher orders of the expansion thus have to be examined before concluding on this possibility.
V. CONCLUSION

This paper leverages the integral solution of Sylvester equation to propose explicit formulas for quantum model reduction via adiabatic elimination, when the remaining subsytem undergoes fast unitary dynamics. From a formal viewpoint, this completes the picture of spectral block-decomposition, by assuming only a timescale separation on the real part of the eigenvalues, i.e. the degrees of freedom which vanish versus remain in the long term, similarly to allowing a central manifold with non-trivial motion in classical systems theory. From a practical viewpoint, we have illustrated how our new results avoid the need to move to the interaction picture and perform an averaging approximation before treating the model reduction.

Future work should allow us to generalize the computations to any linear quantum system [14], providing insights on how their various tunings can be integrated into quantum dissipation engineering. We would also use this tool to further explore corrections to averaging approximations, also called “Rotating Wave Approximation”, RWA, and ubiquitous in quantum engineering by making things on/off-resonant. Indeed, while averaging expansions can in principle be carried out at higher order, they are a priori not converging, unlike the block-spectral decomposition of the present paper.

Acknowledgment

The authors would like to thank Mazyar Mirrahimi, Jérémie Guillaud, Samuel Deléglise, Masaaki Tokieda, Lev-Arcady Sellem and other colleagues for useful discussions. This research has been funded by ANR grants HAMROQS and MECAFLUX (French Research Agency) and by Plan France 2030 through the project ANR-22-PETQ-0006. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 884762).

[1] T. Albash, S. Boixo, D. A. Lidar, and P. Zanardi. Quantum adiabatic markovian master equations. New Journal of Physics, 14(12):123016, 2012.

[2] R. Azouit, F. Chittaro, A. Sarlette, and P. Rouchon. Towards generic adiabatic elimination for bipartite open quantum systems. Quantum Science and Technology, 2(4):044011, 2017.

[3] R. Azouit, A. Sarlette, and P. Rouchon. Adiabatic elimination for open quantum systems with effective Lindblad master equations. In 2016 IEEE 55th Conference on Decision and Control (CDC), pages 4559-4565. IEEE, 2016.

[4] D. Finkelstein-Shapiro, D. Viennot, I. Saideh, T. Hansen, T. Pullerits, and A. Keller. Adiabatic elimination and subspace evolution of open quantum systems. Physical Review A, 101(4):042102, 2020.

[5] P. Forni, T. Launay, A. Sarlette, and P. Rouchon. A palette of approaches for adiabatic elimination in bipartite open quantum systems with hamiltonian dynamics on target. In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 1362-1368. IEEE, 2019.

[6] P. Forni, A. Sarlette, T. Capelle, E. Flurin, S. Deléglise, and P. Rouchon. Adiabatic elimination for multi-partite open quantum systems with non-trivial zero-order dynamics. In 2018 IEEE Conference on Decision and Control (CDC), pages 6614-6619. IEEE, 2018.

[7] J. Guillaud. Thermal adiabatic elimination with adjoint Lindbladian. inria Technical Note, 2020.

[8] S. Haroche and J.-M. Raimond. Exploring the quantum: atoms, cavities, and photons. Oxford University Press, 2006.

[9] F.-M. Le Régent and P. Rouchon. Heisenberg formulation of adiabatic elimination for open quantum systems with two timescales. In 2023 IEEE Conference on Decision and Control (CDC), pages -. IEEE, 2023.

[10] F.-M. Le Régent and P. Rouchon. Adiabatic elimination for composite open quantum systems: Reduced-model formulation and numerical simulations. Phys. Rev. A, 109:032603, Mar 2024.

[11] A. Metelmann and A. A. Clerk. Nonreciprocal photon transmission and amplification via reservoir engineering. Physical Review X, 5(2):021025, 2015.

[12] 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. New Journal of Physics, 16(4):045014, 2014.

[13] B.-L. Najera-Santos et al. High-sensitivity ac-charge detection with a mhz-frequency fluxonium qubit. Phys. Rev. X, 14:011007, 2024.

[14] H. I. Nurdin and N. Yamamoto. Linear dynamical quantum systems. Analysis, Synthesis, and Control, 2017.

[15] J. Poyatos, J. I. Cirac, and P. Zoller. Quantum reservoir engineering with laser cooled ions. Physical review letters, 77(23):4728, 1996.

[16] F. Reiter and A. S. Sørensen. Effective operator formalism for open quantum systems. Physical Review A, 85(3):032111, 2012.

[17] J. A. Sanders, F. Verhulst, and J. Murdock. Averaging methods in nonlinear dynamical systems, volume 59. Springer, 2007.
[18] M. Tokieda, C. Elouard, A. Sarlette, and P. Rouchon. Complete positivity violation in higher-order quantum adiabatic elimination. *arXiv preprint arXiv:2211.11008*, 2022.