Adiabatic elimination for multi-partite open quantum systems with non-trivial zero-order dynamics

# Adiabatic elimination for multi-partite open quantum systems with non-trivial zero-order dynamics

Paolo Forni, Alain Sarlette, Thibault Capelle, Emmanuel Flurin, Samuel Deleglise,
and Pierre Rouchon
Centre Automatique et Systèmes, Mines-ParisTech, PSL Research University. 60 Bd Saint-Michel, 75006 Paris, France; and INRIA Paris, 2 rue Simone Iff, 75012 Paris, France. paolo.forni@inria.fr, pierre.rouchon@mines-paristech.fr. INRIA Paris, 2 rue Simone Iff, 75012 Paris, France; and Ghent University / Data Science Lab, Technologiepark 914, 9052 Zwijnaarde, Belgium. alain.sarlette@inria.fr. Laboratoire Kastler Brossel, ENS-PSL, CNRS, Sorbonne Université et Collège de France, Paris, France. thibault.capelle@lkb.upmc.fr, emmanuel.flurin@lkb.upmc.fr, samuel.deleglise@lkb.upmc.fr.
###### Abstract

We provide model reduction formulas for open quantum systems consisting of a target component which weakly interacts with a strongly dissipative environment. The time-scale separation between the uncoupled dynamics and the interaction allows to employ tools from center manifold theory and geometric singular perturbation theory to eliminate the variables associated to the environment (adiabatic elimination) with high-order accuracy. An important specificity is to preserve the quantum structure: reduced dynamics in (positive) Lindblad form and coordinate mappings in Kraus form. We provide formulas of the reduced dynamics. The main contributions of this paper are (i) to show how the decomposition of the environment into components enables its efficient treatment, avoiding the quantum curse of dimension; and (ii) to extend the results to the case where the target component is subject to Hamiltonian evolution at the fast time-scale. We apply our theory to a microwave superconducting quantum resonator subject to material losses, and we show that our reduced-order model can explain the transmission spectrum observed in a recent pump probe experiment.

## 1 Introduction

The evolution of a quantum system interacting with an environment is rigorously described by a Schrödinger equation on the joint Hilbert space. However, the complexity of the environment hampers the study of the system as a whole and one often resorts to the Born-Markov approximation to obtain a Lindblad master equation [7] describing the target system alone, and the environment’s effect summarized by dissipation or “decoherence” operators. Similarly, when a quantum system consists of several interacting components, e.g. a main computing subsystem coupled to an ancillary subsystem expressing a measurement device, one often seeks to analyze a dynamical equation for the main subsystem alone, approximately including the effect of the ancillary subsystem. In this perspective, model reduction methods come to aid to the physicists interested in gaining better physical insights, in running simplified numerical simulations, and in designing the dynamics of a target subsystem by smartly engineering its interaction with other subsystems, as in the case of reservoir engineering [18].

A classical approach to model reduction for quantum systems makes use of the time-scale separation between a slow subsystem of interest and the fast auxiliary subsystems coupled to it, and eliminates the fast variables in a procedure denominated as adiabatic elimination. In closed quantum systems – where the evolution stays unitary under Hamiltonian dynamics – adiabatic elimination is performed by means of standard perturbation theory techniques [22]. in contrast, the treatment of open quantum systems – including decoherence under Lindbladian dynamics – is more involved. In the literature, adiabatic elimination in the latter case has been addressed for specific examples separately: lambda systems up to second-order [8], a specific atom-optics example [2], systems where excited states decay toward ground states [17, 21], systems with Gaussian dynamics and subject to continuous measurement [14].

However, general approaches to adiabatic elimination of Lindblad systems – and maintaining the positivity-preserving quantum structure, beyond a standard linear systems treatment via singular perturbation theory – have attracted much less attention. In [15], Kessler has developed a generalization of the Schrieffer-Wolff formalism; in [13, 6], the authors address quantum stochastic differential equations in the limit where the speed of the fast system goes to infinity. A geometric approach to adiabatic elimination has been introduced by [4, 3], where the authors explore an asymptotic expansion of the reduced dynamics by a careful application of center manifold techniques [9] and geometric singular perturbation theory [11]. In order to succesfully retain the physical interpretation, the reduced dynamics is expressed by Lindblad equations and is mapped to the original dynamics via a trace-preserving completely-positive (CPTP) map, also called Krauss map.

The present work builds upon the geometric approach of [4, 3] and brings forward two novel features. First, unlike in [4, 3] where the target system was assumed to be static in the ideal case, we here develop formulas for the case where the target system undergoes non-trivial fast Hamiltonian dynamics, when uncoupled from the environment. This appears in all practical situations where the target system is detuned from the reference frame, e.g. by a field to be measured in quantum metrology, or when it undergoes (in this paper constant) drives to implement quantum operations. Second, we consider environments that consist not of a single bulk system, but which can be decomposed into a not-necessarily-finite number of fast dissipative subsystems. Such situations often appear in practice when the target quantum system is corrupted by various imperfection sources [16]. We show how to take advantage of this decomposition towards more efficient model reduction computations. Indeed, the first-order approximation amounts to the sum of the contributions of each fast dissipative subsystem, and the same result holds for the second-order approximation under specific commutation properties of the operators involved in the computation. The proposed theory is applied to a model of a microwave superconducting resonator subject to dielectric losses due to a bath of many two-level-systems. We show how a reduced model resulting from our theory allows to explain the non-trivial transmission spectrum observed in a pump probe experiment.

The outline of the paper is as follows. Setting and main assumptions are introduced in Section 2. Section 3 provides our main results with the formulas of our adiabatic elimination for the case of many fast dissipative subsystem weakly coupled to the target one. Section 4 contains the application and comparison to experimental data. We conclude the paper with few final remarks. Proof and computation details are given in appendix.

## 2 Setting

### 2.1 K-partite systems with non-trivial zero-order dynamics

Open quantum systems are typically described by differential equations evolving on the manifold of density operators , namely the set of all linear Hermitian nonnegative operators from a Hilbert space to itself, whose trace equals one. The evolution of an open quantum system is then described by the Lindblad master equation [7]:

 dρdt=L(ρ)=−i[H,ρ]+∑μDLμ(ρ),

where each is a “decoherence” operator on , is a Hermitian “Hamiltonian” operator on , and is a superoperator defined by:

 DLμ(ρ):=LμρL†μ−12L†μLμρ−12ρL†μLμ.

In this paper, we consider the composite Hilbert space of a target quantum system on and its environment on . The dynamics on satisfies a two-time scale separation:

 dρdt=LA(ρ)+εLint(ρ)+εLB(ρ)+(−i)[~HB,ρ], (1)

where is a small positive parameter; and are Lindbladian super-operators acting exclusively on and respectively; is a Lindbladian superoperator which captures the interaction between and . Here we assume that this interaction is Hamiltonian and expressed as:

 Lint(ρ):=−i[A⊗B†+A†⊗B,ρ],

where and respectively are non-necessarily-Hermitian operators acting on and only. Finally, is a Hamiltonian operator on , thus expressing fast unitary dynamics on the target system; its presence is the first novelty in our paper. For a set of interesting situations, the dynamics of typical quantum systems can be expressed in a rotating frame where the term would vanish. However, several reasons can justify to keep this term. For instance, in many significant situations the vanishing of is not rigorous and involves an additional treatment of appearing fast time-varying parameters in the equation via averaging theory; or, can be a term of particular interest like a field to be measured with the quantum device or an actuation towards applying some operation on the target system.

As a second novelty, we consider a generalized setting where is composed of a non-necessarily-finite number of Hilbert spaces . Each subsystem on is strictly dissipative and interacts with only. Then, system (1) reads as:

 dρdt=∑k(L(k)A(ρ)+εL(k)int(ρ))+εLB(ρ)+(−i)[~HB,ρ] (2)

where acts on only and where

 L(k)int(ρ):= −i[A(k)⊗B†+A(k)†⊗B,ρ],

captures the Hamiltonian interaction between and , with non-necessarily-Hermitian operators acting on only.

For , the system is uncoupled and the solution trajectories stay separable for all times, namely for we have for all times, with each factor in the product following its independent dynamics. To apply adiabatic elimination, we assume that each part of the environment is highly dissipative and relaxes fast to a unique steady state, i.e.: for any initial state on , the solution of the uncoupled system converges to where, for each , is the unique solution of ; and satisfies with . For ease of presentation, we will also denote .

### 2.2 Asymptotic expansion

Both in the bi-partite and the -partite case, for the uncoupled system , there exists an asymptotically stable center manifold of same dimension as , on which the dynamics have imaginary eigenvalues. It thus follows from Fenichel’s Invariant Manifold Theorem [11] that, for small enough , there exists an invariant and attractive manifold which has the same dimension as and which is close to it. Furthermore, by virtue of Carr’s result [9], an approximation of can be computed up to arbitrary precision. The quantum particularity, as explained in [3], is that such approximation should retain a physical interpretation by preserving the quantum structure:

• the mapping from the reduced space to the complete space is a mapping between density operators, and can be parameterized by for some Hilbert space that has same dimension as , and where is a Kraus map111A Kraus map takes the form for some operators in order to express any completely positive superoperator [10], and with ensuring trace-preservation i.e. . ;

• the reduced dynamics on are Lindbladian, i.e. for some Lindbladian superoperator .

In other words, we aim to find a Kraus map and a Lindbladian such that the following invariance equation is satisfied for all small enough and for all :

 LA(K(ρs))+εLint(K(ρs))+εLB(K(ρs)) +(−i)[~HB,K(ρs)]=K(Ls(ρs)). (3)

We impose that both the Kraus map and the Lindbladian are parameterized as an infinite series:

 K(ρs):=+∞∑h=0εhKh(ρs),Ls(ρs):=+∞∑h=0εhLs,h(ρs).

Then, by identifying the terms of the same order of in the invariance equation (3), we obtain an invariance relation at all orders . At zero-order, we have:

 LA(K0(ρs))+(−i)[~HB,K0(ρs)]=K0(Ls,0(ρs)). (4)

Similarly, the first-order invariance condition reads as:

 LA(K1(ρs))+Lint(K0(ρs))+LB(K0(ρs)) −i[~HB,K1(ρs)]=K0(Ls,1(ρs))+K1(Ls,0(ρs)), (5)

whereas the second-order invariance condition reads as:

 LA(K2(ρs))+Lint(K1(ρs)) +LB(K0(ρs))−i[~HB,K2(ρs)] (6)

## 3 Reduced-model formulas

The aim of this Section is to provide explicit solutions to the zero-, first-, and second-order invariance equations (4)-(6) for the case of -partite systems as introduced in Section 2.1, i.e. for model (2). We immediately observe that the zero-order (4) is naturally solved by setting:

 Ls,0(ρs):=−i[~HB,ρs],K0(ρs):=(⨂k¯ρ(k)A)⊗ρs. (7)

At first order, let the Kraus map have the following structure inspired by [3]:

 K(ρs)=K0(ρs)+εK1(ρs):= (I−iεM)(¯ρA⊗ρs)(I+iεM†)+O(ε2), (8)

where , for any . This would immediately imply that:

 K1(ρs)=−iM(¯ρA⊗ρs)+i(¯ρA⊗ρs)M. (9)

The following assumption will be instrumental in establishing our main results.

###### Assumption 1

There exists such that:

 [~HB,B†]=cB†B†. (10)
###### Theorem 1

Consider model (2). Let Assumption 1 hold. Then, the first-order invariance equation (5) is satisfied by the Lindbladian and by a map of the form (9) where, for each , respectively are the unique solutions of:

 L(k)A(F(k)1¯ρ(k)A)+A(k)¯ρ(k)A−icB†F(k)1¯ρ(k)A= 0, (11a) L(k)A(F(k)2¯ρ(k)A)+A(k)†¯ρ(k)A+ic∗B†F(k)2¯ρ(k)A= 0. (11b)

Furthermore, is a CPTP map up to second-order terms.

###### Proof 1

see Appendix A.1.

###### Remark 1

The first-order (5) is also satisfied by the Lindbladian with , and by a map of the form (9) where respectively are the unique solutions of:

 L(k)A(F(k)1¯ρ(k)A)+S(k)(A(k)¯ρ(k)A)−icB†F(k)1¯ρ(k)A= 0, (12a) L(k)A(F(k)2¯ρ(k)A)+S(k)(A(k)†¯ρ(k)A)+ic∗B†F(k)2¯ρ(k)A= 0, (12b)

where, for an operator acting on , notation denotes . Furthermore, is a CPTP map up to second-order terms. The possibility of having alternative solutions to the first-order invariance equation hinges upon a gauge degree of freedom in the selection of the trace of terms and . It appears that gauge choices are instrumental for positivity-preservation in the solution of the second-order invariance equation, as we consider next.

###### Theorem 2

Consider model (2). Let Assumption 1 hold. Assume furthermore that , with selected according to Theorem 1. Then, the second-order invariance equation is satisfied by a Lindbladian:

 Ls,2(ρs)=∑k −iI(z(k)1)[BB†,ρs]−iI(z(k)2)[B†B,ρs] +2R(z(k)1)DB†(ρs)+2R(z(k)2)DB(ρs) +∑k>k′{−iδ(k,k′)[[B,B†],ρs]}, (13)

with:

 δ(k,k′) =−2R(z(k)0z(k′)∗0)cB†,z(k)0=Tr(A(k)¯ρ(k)A), (14a) z(k)1 =Tr(F(k)1¯ρ(k)AA(k)†),z(k)2=Tr(F(k)2¯ρ(k)AA(k)), (14b)

and by a map , obtained from formulas (23)-(31), such that is a CPTP map up to third-order terms.

###### Proof 2

see Appendix A.2.

###### Corollary 1

Under the same assumptions of Theorem 2, if for some , then the second-order invariance equation (6) is satisfied by the Lindbladian:

 Ls,2(ρs) = ∑k{−iI(z(k)1+z(k)2)[B†B,ρs] +2R(z(k)1)DB†(ρs)+2R(z(k)2)DB(ρs)}.

## 4 Application

Microwave superconducting resonators are an important component in various quantum devices, and in particular in the quantum electrodynamics circuits [5, 19] that are one of the most promising current technologies towards building a quantum computer [1]. Losses due to imperfections in amorphous materials constitute a dominant loss channel of such resonators [23, 12], and can be represented by a bath of two-level systems (TLSs). In many practical cases, strong microwave tones are applied with significant frequency detuning with respect to the resonance frequency [20] in order to activate a parametric interaction between the resonator mode and another circuit degree of freedom. Within this framework, the LKB team has performed a pump probe experiment on the microwave resonator in Figure 1: a strong “pump” drive, at a frequency far detuned from the resonator, is applied to essentially scramble the quantum behavior (“saturate”) of the TLS bath, whereas a weak probe tone, assumed not to disturb the bath behavior, is used to retrieve the transmission spectrum of the resonator. The latter allows to extract induced detuning and damping rate.

Let and respectively be the Hilbert space of the -th TLS=qubit and the resonator, and . Respectively denote with and the raising and lowering operator on the -th qubit, and with , , and the Pauli operators on the -th qubit. Let and be the annihilation and creation operators in the resonator mode. The experimental setup is modeled by the following system in Lindblad form:

 ddt~ρ=−i[H,~ρ]+Γ−∑kDσ(k)−(~ρ), (16)
 H = ωca†a+(veiωrt+v∗e−iωrt)(a†+a) +∑k(ω(k)q2σ(k)z+igσ(k)x(a†−a)).

Here , , and are the frequencies of the resonator, the pump drive, and the -th qubit respectively, is the strength of the pump, is the coupling strength between the resonator and each qubit, and is the dissipation rate associated to on each qubit. The goal would be to obtain a reduced order model for (16) which matches the transmission spectrum of this experiment.

For each , let and . Under the assumption that for any , we apply the standard rotating-wave approximation (i.e. first-order averaging) with the Hamiltonian corresponding to the rotating change of frame and the remaining Hamiltonian. The first-order RWA yields:

 ddtρrwa1=−i[Hrwa1,ρrwa1]+Γ−∑kDσ(k)−(ρrwa1),

where is the Jaynes-Cummings interaction Hamiltonian plus drive. We next apply a unitary coordinate change on the resonator state, to center it around its well-known steady state under off-resonant drive, namely by a complex field amplitude displacement . This yields:

 ddtρ=∑k{L(k)Q(ρ)+gL(k)int(ρ)}+(−i)[Δca†a,ρ], (17)

where

 L(k)Q(ρ) L(k)int(ρ) :=[σ(k)−a†−σ(k)+a,ρ].

The term with for now expresses an effective, indirect drive on the TLSs.

System (17) is in the form (2) with , , , . The hypothesis of Theorems 1 and 2 are satisfied since and Assumption 1 holds with . Let be a Hilbert space whose dimension matches the dimension of the resonator space , and the density operator on . By virtue of Theorem 1 and Corollary 1, the reduced model is given in Lindbladian form as follows:

 ddtρs= −i(Δc+g2∑kδ(k))[a†a,ρs] (18) δ(k) =I(z(k)1+z(k)2), Γ(k)a† =2R(z(k)1),Γ(k)a=2R(z(k)2), z(k)1 =Tr(−iF(k)1¯ρ(k)Qσ(k)+),z(k)2=Tr(iF(k)2¯ρ(k)Qσ(k)−)

and where, for each , matrices , satisfy equations (11). The solution of such equations can be computed directly since, on a qubit space , one can always parameterize operators in terms of Pauli matrices. We immediately find: and where

 W(k)1= −4g2v2(8ig2v2+ (Γ−+iΔc)Δc((Γ−+2iΔc)2+4Δ2q)), W(k)2= 32ig4v4+2i(Γ−−iΔc)Δ4c(Γ2+4Δ2q)⋅ ⋅(Γ−−2i(Δc+Δq))−4g2v2Δc(Γ3−−5iΓ2−Δc+ 4iΔc(Δ2c+2ΔcΔq−Δ2q)+4Γ−(−Δ2c+Δ2q)), Z(k)= (8g2v2+Δ2c(Γ2−+4Δ2q))(8g2v2(iΓ−+2Δc)+ Δ2c(iΓ−+Δc)((Γ−−2iΔc)2+4Δ2q)).

Coefficients , , and appearing in our reduced-order model (18) can be visualized for different values of pump detuning and intraresonator photon number . As depicted in Figure 1, we can compare the frequency shift of the resonator with experimental findings from the pump probe experiment. We find that, by properly calibrating the values of and by selecting a proper density function of the TLSs, we are able to match the resonance shift in the trasmission spectrum observed in the pump probe experiment. Quantitative agreement between data and our model should hence enable to extract characteristics about the TLS bath, pending other experimental features that will have to be calibrated.

## 5 Conclusions

We have studied adiabiatic elimination for open quantum systems in Lindblad form composed by a target subsystem weakly interacting with strongly dissipative subsystems. The key novel features of our approach are twofold: the decomposition of the environment into separately treated subsystems, and the presence of fast Hamiltonian dynamics on the target system. The time-scale separation between the uncoupled dynamics and the interaction allows model reduction via center manifold techniques and geometric singular perturbation theory. We have provided formulas for the first- and second-order expansion and shown that the asymptotic expansion of the center manifold retains a physical interpretation: the reduced model still evolves according to Lindbladian dynamics and can be mapped to the original model via Kraus map. Each strongly dissipative subsystem contributes linearly to the reduced model at first-order, and does the same at second-order if a specific commutation property about the interaction terms holds. We have successfully applied our proposed theory to the model of a microwave superconducting resonator subject to dielectric losses where our reduced-order model shows a trasmission spectrum whose shape matches experimental data. Future work will address the necessary conditions to satisfy the invariance equation at orders higher than two. The assumption about the commutator between the original Hamiltonian dynamics of the target system and the interaction terms might also be removed, yielding a full generalization of the proposed theory.

## Appendix A Proofs of Theorems

### a.1 Proof of Theorem 1

By plugging (27) and (29) into the first-order invariance condition (5) and by making use of Assumption 1, condition (5) reads as:

 ∑k¯ρ[k]A⊗{(−iL(k)A(F(k)1¯ρ(k)A)−iA(k)¯ρ(k)A −cB†F(k)1¯ρ(k)A) ⊗B†ρs+ (19a) (−iL(k)A(F(k)2¯ρ(k)A)−iA(k)†¯ρ(k)A +c∗B†F(k)2¯ρ(k)A) ⊗Bρs+ (19b) +\emph{herm.conj.}+(LB(ρs)−Ls,1(ρs)) ⊗¯ρA=0 (19c)

Case of Theorem 1. It can be immediately seen from (19) that one can select as long as each round parenthesis (19a),(19b) is set to zero. Taking the trace on (11a) and (11b) yields:

 cB†Tr(F(k)1¯ρ(k)A)=cB†Tr(F(k)2¯ρ(k)A)∗=−iTr(A(k)¯ρ(k)A), (20)

which solves the situation with the announced formulas.

Case of Remark 1. By taking the trace on equations (12), we observe that, for each :

 Tr(F(k)1¯ρ(k)A)=Tr(F(k)2¯ρ(k)A)∗=0. (21)

Then, by taking the partial trace over in (19), we immediately have with as in Theorem 1. Now, plugging in (19) yields:

 ∑k¯ρ[k]A⊗{(−iL(k)A(F(k)1¯ρ(k)A) −iS(k)(A(k)¯ρ(k)A)−cB†F(k)1¯ρ(k)A) ⊗B†ρs+ (22a) (−iL(k)A(F(k)2¯ρ(k)A) −iS(k)(A(k)†¯ρ(k)A)+c∗B†F(k)2¯ρ(k)A) ⊗Bρs+ (22b) \emph{herm.conj.}}=0

In order to solve (22) it is enough to set each round parenthesis of (22a) and (22b) to zero for each — see equations (12).

It is immediate to see from (8) that is a completely positive map, as long as one can neglect the terms of order . One concludes that it is also trace-preserving at order by checking that , thanks to (20) for the case of Theorem 1 and thanks to (21) for the case of Remark 1.

### a.2 Proof of Theorem 2

We use the notation for operators that, if their superscripts respectively are or , then they respectively apply to or only, possibly tracing out the remaining subsystems. Let denote . Define:

 M= ∑kM(k),M(k)=F(k)1⊗B†+F(k)2⊗B N:= ∑kN(k)+∑k≠k′N(k,k′) N(k):= U(k)BB⊗BB+U(k)BB†⊗BB†+ U(k)B†B⊗B†B+U(k)B†B†⊗B†B† N(k,k′):= U(k,k′)BB⊗BB+U(k,k′)BB†⊗BB†+ U(k,k′)B†B⊗B†B+U(k,k′)B†B†⊗B†B† P(k)W:= W(k)1⊗B†+W(k)2⊗B P(k)V:= V(k)1⊗B†+V(k)2⊗B P(k,k′)W:= W(k,k′)1⊗B†+W(k,k′)2⊗B P(k,k′)V:= V(k,k′)1⊗B†+V(k,k′)2⊗B, (23)

and define, for :

 ~W(k)ij:= F(k)i¯ρ(k)AF(k)†j+W(k)i¯ρ(k)AW(k)†j+V(k)i¯ρ(k)AV(k)†j ~W(k,k′)ij:= F(k,k′)i¯ρ(k,k′)AF(k)†j+W(k,k′)i¯ρ(k,k′)<