Projection based adiabatic elimination of bipartite open quantum systems
Ibrahim Saideh, Daniel Finkelstein-Shapiro, Camille Noûs, Tõnu Pullerits, Arne Keller
PProjection based adiabatic elimination of bipartite open quantum systems
Ibrahim Saideh, Daniel Finkelstein-Shapiro, Camille Noˆus, T˜onu Pullerits, and Arne Keller Laboratoire Mat´eriaux et Ph´enom`enes Quantiques,Universit´e Paris Diderot, CNRS UMR 7162, 75013, Paris,France. Universit´e Paris-Saclay, 91405 Orsay, France Division of Chemical Physics, Lund University, Box 124, 221 00 Lund, Sweden Cogitamus Laboratory ∗ Adiabatic elimination methods allow the reduction of the space dimension needed to describesystems dynamics which exhibits separation of time scale. For open quantum system, it consistsin eliminating the fast part assuming it has almost instantaneously reached its steady-state andobtaining an approximation of the evolution of the slow part. These methods can be appliedto eliminate a linear subspace within the system Hilbert space, or alternatively to eliminate a fastsubsystems in a bipartite quantum system. In this work, we extend an adiabatic elimination methodused for removing fast degrees of freedom within a open quantum system (Phys. Rev. A , ,042102) to eliminate a subsystem from an open bipartite quantum system. As an illustration, weapply our technique to a dispersively coupled two-qubit system and in the case of the open Rabimodel. I. INTRODUCTION
Adiabatic elimination is a method whereby the fast de-grees of freedom of a system are removed while retainingan effective description of the slow degrees of freedom.This simplification can be very useful to obtain tractableand intuitive equations when only a coarse-grained orlong times description is desired [1–9], depending on ifthe target system has a conservative [5, 6, 10] or a dissi-pative [11–18] evolution. There are two classes of man-ifolds on which adiabatic elimination has been applied,i) those that consist of levels within a subsystem, for ex-ample the excited states of an atom, and ii) those thatconsist of a separate subsystem, such as ancillary qubitsor measuring devices. For a slow and fast parts describedby Hilbert spaces H ( A ) and H ( B ) , the first case corre-sponds to the Hilbert space H = H ( A ) ⊕ H ( B ) (directsum) while the second case corresponds to the Hilbertspace H = H ( A ) ⊗ H ( B ) (tensor product). Adiabaticelimination is useful in developing protocols for dissipa-tive state preparation in ion traps [19, 20], reservoir engi-neering [21] and the description of measurement devices[22]. The simplicity of the resulting equations can also becomputationally advantageous in the study of quantumphase transitions where the size of the system is cumber-somly large [23].There are several approaches to obtain effective oper-ators, ranging from perturbative expansions of the Li-ouville operator [17, 18], the corresponding Kraus maps[11–14, 16], the resolvent [24], or using stochastic meth-ods [22]. Eliminating a fast subsystem (that forms a ten-sor product with the slow subsystem) is typically donewith a partial trace over the fast subsystem. This canresult in a set of hierarchical equations that allows errorestimation and correcting the approximation as the slow ∗ [email protected] and fast timescales get closer. Importantly, the expan-sion can be built to preserve the Lindblad structure andas a consequence the physicality of the map [12, 25, 26].The procedure for eliminating sublevels within a subsys-tem (direct sum with the slow subsystem) is best car-ried out with Feshbach projections [17, 24]. However, asthe fast-slow separation breaks down, or when incoher-ent pumping channels exist, the population of the fastsubsystem becomes non-negligible (i.e. there can be afinite fraction of population in the excited states). Whenthis happens, the exact time evolution of the slow partbecomes non-trace preserving. The loss of trace can becorrected using contour integral methods [24]. It wouldbe however advantageous to have a method that can han-dle both classes of fast manifolds. This is more importantconsidering that systems from atomic physics are inspir-ing a number of chemical versions that have much morecomplicated Hamiltonians and it would be ideal to trans-form them into effective operators for a direct comparisonto the atomic physics counterparts [27–29].In this work, we extend the methodology developed inRef. [24] to bipartite open quantum systems whose dy-namics are described by a Lindblad operator [30, 31]. Weuse the projection operator method suggested by Kneze-vic and Berry [32] in order to derive equations for a slowsubsystem A coupled to a fast subsystem B . The paperis organized as follows. We first recall the main resultsof Ref [24]. We then apply it to the general bipartitecase to obtain a recipe for describing the slow subsys-tem. Finally, we illustrate the method in the case of aspin dispersively coupled to a second highly dissipativedriven spin and to describe the dynamics of the openRabi model. a r X i v : . [ qu a n t - ph ] J un II. THEORYA. Adiabatic elimination through projectorstechniques
Let ρ ( t ) be the density operator on the Hilbert space H describing the quantum state of the system at time t . We suppose that the evolution of ρ ( t ) is generatedby a Lindblad operator L : ˙ ρ ( t ) = L ρ ( t ). We define theHilbert space H of operator O on H , equipped with thescalar product tr (cid:104) O † O (cid:105) . We first recall the main resultspresented in Ref [24] related to projector techniques. Let P the projector such that ρ s ( t ) = P ρ ( t ) is the the slowpart of the system and write Q = − P , where H . Let G ( z ) = ( z − L ) − bethe resolvent of the Lindblad operator L . Operators like P , Q , L or G are operators on H . They are sometimescalled super-operators. They are here denoted with cal-ligraphic letter, to distinguish them from operators on H (belonging to H ), like the density matrix ρ .We define the effective Lindblad operator L eff ( z ), suchthat PG ( z ) P = ( z − L eff ( z )) − . The effective Lindbladoperator L eff ( z ) can be written as : L eff ( z ) = PLP + PLQG ( z ) QLP , (1)where QG ( z ) Q = ( z − QLQ ) − . (2)For any ρ ( t = 0), such that Q ρ ( t = 0) = 0, the slowdynamics inside P H can be obtained with the inverseLaplace transform as: ρ s ( t ) = 12 πi (cid:90) D d ze zt P G ( z ) P ρ s ( t = 0) , (3)where ρ s ( t = 0) = P ρ ( t = 0) and the integral onthe complex plane is performed on a straight line D = { z ∈ C ; (cid:60) z = a > } . At this point no approximation hasbeen made. Eq. (3) gives the exact dynamics inside P H ,as long as the initial condition is also inside P H , thatis Q ρ ( t = 0) = 0. Because L eff ( z ) captures the effect ofthe dynamics in Q H , it is a nonlinear superoperator, inthe sense that ( z − L eff ( z )) (cid:126)ν = 0 is a nonlinear eigenvalueproblem.The approximation of a slow dynamics of P ρ ( t ), withrespect to the fast dynamics of Q ρ ( t ) is equivalent toconsidering the dynamics inside P H in the vicinity ofthe stationary state reached in the limit t → ∞ . Inthis long time limit only the z → L eff ( z ) to the lowest relevant order: L eff ( z ) (cid:39) L + z L + · · · + z n L n (4)where L = L eff ( z = 0) and L n = n ! dd z L eff ( z ) (cid:12)(cid:12) z =0 . Us-ing the expression of L eff ( z ) given by Eq. (1) allows to express L and L n as: L = PLP − PLQ ( QLQ ) − QLPL n = −PLQ ( QLQ ) − ( n +1) QLP . (5)In this work, we consider the approximation given byEq. (4) with n ≤ n >
1) will be considered in a future work.Within the approximation given by Eq. (4), with n = 1,the inverse Laplace transform of Eq. (3) can be computedexplicitly. We obtain: ρ s ( t ) = exp (cid:2) (1 − L ) − L t (cid:3) (1 − L ) − ρ s ( t = 0) (6)The stationary state ρ f of this dynamics, reached at t → + ∞ , is in the kernel of L . We note that although thedynamics described by Eq. (6) is an approximation, thefinal reached stationary state ρ f is the exact one.To conclude, after the adiabatic elimination of the fastpart, the generator of the slow dynamics is approxima-tively given by (1 − L ) − L , ˙ ρ ( t ) = (1 − L ) − L ρ ( t ),where L and L can in principle be computed usingEq. (5). The hard part in these equations is the evalua-tion of the inverse ( QLQ ) − , which in the most generalcase, as we will see later, can only be achieved through aperturbative expansion.Theses results are very general, and require only thedefinition of a projector P and that the initial state fulfillsthe condition Q ρ ( t = 0) = 0. We note that P doesn’thave to be hermitian, that is the projection does not needto be orthogonal.In Ref. [24], this formalism was applied to the casewhere the slow and fast degrees of freedom correspond toa partition of the underlying Hilbert space in two com-plementary sub-spaces, that is H = H ( A ) ⊕ H ( B ) . In thenext section we will adapt this general formalism to thebipartite case where H = H ( A ) ⊗ H ( B ) . B. Adiabatic elimination in a bipartite system
We suppose that the state of the bipartite system attime t is described by a density operator ρ ( AB ) ( t ) actingon the Hilbert space H = H ( A ) ⊗ H ( B ) . We consider thatthe dynamics of subsystem A is very slow compared tothe dynamics of subsystem B . We suppose that the exactstationary state in H is unique and that it is a productstate ρ f = ρ a ⊗ ρ b , where ρ a ∈ H ( A ) and ρ b ∈ H ( B ) with H ( i ) , i = A, B , the Hilbert space of operators on H ( i ) . As the the dynamics of B subsystem is very fast,we suppose that it is a good approximation to considerthat at t = 0, ρ ( AB ) ( t = 0) = ρ ( A )0 ⊗ ρ b . In other word, weconsider that B reaches its steady state instantaneouslyin the time scale of the subsystem A . We thus define P as P ρ ( AB ) ( t ) = tr B (cid:104) ρ ( AB ) ( t ) (cid:105) ⊗ ρ b , (7)where tr B [] denotes the partial trace over B . The re-duced density operator ρ ( A ) ( t ) in H ( A ) can then be ob-tained as ρ ( A ) ( t ) = tr B (cid:2) P ρ ( AB ) ( t ) (cid:3) .For the purpose of simplifying some expressions andcalculations, it will be useful to use the operator-vectorisomorphism [33] which maps each element of H to a vec-tor in H ⊗ H as follows. An operators such as | a (cid:105)(cid:104) b | ∈ H is mapped to the vector | b (cid:105) ⊗ | a (cid:105) in the H ⊗ H
Hilbertspace, where | b (cid:105) is the complex conjugate of | b (cid:105) . Con-sequently, any n × n density matrix ρ ∈ H is mappedto a column vector || ρ (cid:105)(cid:105) ∈ H ⊗ H , with n elements, bystacking the columns of the ρ matrix. Under this isomor-phism, super-operators on H are mapped to operators on H ⊗ H . In particular, the super-operator O performingthe operation ρ → O ρO † , with O and O operatorsin H , is mapped to || ρ (cid:105)(cid:105) → O ⊗ O || ρ (cid:105)(cid:105) , where O de-notes the complex conjugate of O ; that is O = (cid:0) O † (cid:1) T ,where O † is the adjoint and O T is the transpose of O . Inthis way, the scalar product tr (cid:104) ρ † ρ (cid:105) between two oper-ators ρ and ρ in H is equal to the usual scalar product (cid:104)(cid:104) ρ | ρ (cid:105)(cid:105) in H ⊗ H . Some useful remarks can be made.The identity operator H is mapped to the maximallyentangled state || (cid:105)(cid:105) = (cid:88) k | k (cid:105) ⊗ | k (cid:105) (8)in H ⊗ H , where {| k (cid:105)} is an orthonormal basis of H .We also note that the usual density matrix normaliza-tion tr [ ρ ] = 1 does not correspond to the normalizationinduced by the scalar product tr (cid:2) ρ (cid:3) = 1 (except in thecase of a pure state). Using the previous remark, wehave that tr [ ρ ] = tr [ ρ ] = 1 is mapped to (cid:104)(cid:104) | ρ (cid:105)(cid:105) = 1.For our bipartite case, H = H ( A ) ⊗ H ( B ) . Therefore,an operator in H as | a (cid:105)(cid:104) a | ⊗ | b (cid:105)(cid:104) b | , where | a i ( b i ) (cid:105) ∈H ( A ) ( H ( B ) )( i = 1 , | a (cid:105)⊗| a (cid:105)⊗| b (cid:105)⊗| b (cid:105) .The partial trace over H ( B ) , tr B [ ρ ] ∈ H ( A ) is mapped to (cid:104)(cid:104) ( B ) | ρ (cid:105)(cid:105) ∈ H ( A ) ⊗ H ( A ) , where | ( B ) (cid:105) = (cid:80) k | k (cid:105) ⊗ | k (cid:105) in H ( B ) ⊗ H ( B ) , where {| k (cid:105)} is now an orthonormal basis of H ( B ) .Consequently, the operator P ρ ( AB ) ∈ H is mapped tothe vector (cid:104)(cid:104) ( B ) | ρ ( AB ) (cid:105)(cid:105)⊗|| ρ b (cid:105)(cid:105) ∈ H ⊗H . The projector P acting on H ⊗ H can thus be written as: P = (2 A ) ⊗ P (2 B ) with P (2 B ) = || ρ b (cid:105)(cid:105)(cid:104)(cid:104) ( B ) || (9)where (2 A ) is the identity operator on H ( A ) ⊗H ( A ) . Theoperator Q = − P simply reads as: Q = (2 A ) ⊗ Q (2 B ) with Q (2 B ) = (2 B ) − P (2 B ) , (10)where (2 B ) is the identity operator on H ( B ) ⊗ H ( B ) .In general, the Lindblad operator of the system can besplit in 3 terms as follows: L = L A ⊗ (2 B ) + (2 A ) ⊗ L B + L AB (11)where L A ( B ) is a Lindblad operator acting on the A ( B )part only. The decomposition of L with the help of P and Q reads: PLP = L A ⊗ P B + PL AB P (12) PLQ = PL ( AB ) Q (13) QLP = (2 A ) ⊗ Q B L B P B + QL ( AB ) P (14) QLQ = L A ⊗ Q B + (2 A ) ⊗ Q B L B Q B + QL AB Q (15)Where we have used the fact that (cid:104)(cid:104) ( B ) ||L B = 0, as L B is a trace preserving operator. In some specific cases, itcan be a good approximation to take as ρ B , a stationarystate of L B . In that case L B P B = 0 and QLP in Eq. (14)can be simplified as
QLP = QL ( AB ) P (16)For computing L using Eq. (6), the main difficulty re-sides in the inversion of QLQ . In general this inversioncan not be done explicitly, but a perturbative expansioncan give a good approximation. In many cases
QLQ canbe divided into two terms,
QLQ = D + V , where thecomputation of D − is easy, and V is small compared to D . In that case we can write( QLQ ) − = D − ∞ (cid:88) n =0 (cid:0) VD − (cid:1) n (17)Retaining the terms up to n = 0 or n = 1 may give agood approximation.In some cases, the term (2 A ) ⊗ Q B L B Q B in Eq. (15)is dominant over the two other terms. In that case, ap-proximating the inverse of QLQ by(
QLQ ) − (cid:39) (2 A ) ⊗ ( Q B L B Q B ) − (18)can be sufficient as we will see in the examples in the nextsection. So, L can be approximated by the followingexpression: L = L A ⊗ P B + PL AB P + PL ( AB ) Q (cid:104) (2 A ) ⊗ ( Q B L B Q B ) − (cid:105) × (cid:104) (2 A ) ⊗ Q B L B P B + QL ( AB ) P (cid:105) (19)This is the main result of this work. III. EXAMPLES
We apply the formalism of the preceding section totwo examples. We first address the case of a stronglydissipative driven qubit B dispersively coupled to a targetqubit A . Then, as a second example, we consider theopen Rabi model in the regime where the dynamics ofthe spin is very fast compared to the boson frequency. A. A 2-qubit system
This two qubits system has been considered previouslyby Azouit et in Ref. [12] to test another method of bi-partite adiabatic elimination (note that in their work, itis the A spin which is the strongly dissipative spin). Itconsists in a strongly dissipative driven qubit B disper-sively coupled to a target qubit A . This model is usedin Ref. [15] to describe the continuous measurement ofa harmonic oscillator excitation number (correspondingto system A ) by a spin (corresponding system B ). TheLindblad equation for the bipartite system can be writtenas: dρdt = u [ σ B + − σ B − , ρ ] + γ (cid:18) σ B − ρσ B + − σ B + σ B − ρ + ρσ B + σ B − (cid:19) − iχ [ σ Az ⊗ σ Bz , ρ ] (20)where σ + = | (cid:105)(cid:104) | , σ − = σ † + and | (cid:105) , | (cid:105) are the eign-vectors of the Pauli matrix σ z with eigenvalues − , τ = ut , χ (cid:48) = χu ,and γ (cid:48) = γu and using the column-vector isomorphism,we write the superoperator form of the Liouvillian as: L = − i (cid:2) χ (cid:48) (cid:0) A ⊗ σ Az ⊗ B ⊗ σ Bz − σ Az ⊗ A ⊗ σ Bz ⊗ B (cid:1) − (cid:0) A ⊗ B ⊗ σ By + A ⊗ σ By ⊗ B (cid:1)(cid:3) + γ (cid:48) A ⊗ σ B − ⊗ σ B − − γ (cid:48) (cid:2) A ⊗ σ B + σ B − ⊗ B + A ⊗ B ⊗ σ B + σ B − (cid:3) where we have used the relations¯ σ y = − σ y , ¯ σ + = σ + , ¯ σ − = σ − (21)From now on, we will drop the prime in all the param-eters γ (cid:48) , χ (cid:48) . As the qubit A only comes into play throughthe Hamiltonian term − iχ [ σ Az ⊗ σ Bz , ρ ] (see Eq. (20)), thekernel of the Lindblad operator is two dimensional, ac-cording to the two eigenvectors σ Az . Hence the kernel canbe considered as the span of { ρ ( A ) s ⊗ ρ ( B ) s , ρ ( A ) s ⊗ ρ ( B ) s } ,where ρ ( A ) s = | (cid:105)(cid:104) | , ρ ( A ) s = | (cid:105)(cid:104) | and (see Appendix fordetails): ρ ( B ) s =
12 + 2 γ χ + γ + 8 σ ( B ) x + 8 χ χ + γ + 8 σ ( B ) y − χ + γ χ + 2 γ + 16 σ ( B ) z (22) ρ ( B ) s =
12 + 2 γ χ + γ + 8 σ ( B ) x − χ χ + γ + 8 σ ( B ) y − χ + γ χ + 2 γ + 16 σ ( B ) z (23)In order to avoid any unnecessary complications, wewill assume that the qubit A has an extremely slow dis-sipation rate that we will omit from our calculations, but will ensure the uniqueness of the steady state. In otherwords, the steady state of the considered system will be ρ ss = ρ ( A ) s ⊗ ρ ( B ) s . We thus assume that the initial state isof the form ρ = ρ ( A )0 ⊗ ρ ( B ) s , and we define the projector P (see Eq. (9)): P = A ⊗ || ρ ( B ) s (cid:105)(cid:105)(cid:104)(cid:104) ( B ) || , (24)where according to Eq. (8), || ( B ) (cid:105)(cid:105) = | (cid:105) ⊗ | (cid:105) + | (cid:105) ⊗ | (cid:105) .This ensures that Q = − P verifies Q| ρ (cid:105) = 0 as itshould be. In Appendix A, we calculate in detail all thequantities necessary to compute L and L as given byEq. (5) : PLP = | , (cid:105) (cid:104) , | ⊗ A , + | , (cid:105) (cid:104) , | ⊗ A , (25) PLQ = | , (cid:105) (cid:104) , | ⊗ C , + | , (cid:105) (cid:104) , | ⊗ C , (26) QLP = | , (cid:105) (cid:104) , | ⊗ D , + | , (cid:105) (cid:104) , | ⊗ D , + | , (cid:105) (cid:104) , | ⊗ D , + | , (cid:105) (cid:104) , | ⊗ D , (27) QLQ = | , (cid:105) (cid:104) , | ⊗ B , + | , (cid:105) (cid:104) , | ⊗ B , + | , (cid:105) (cid:104) , | ⊗ B , + | , (cid:105) (cid:104) , | ⊗ B , (28)where X i,j , X ∈ {A , B , C , D} , and i, j ∈ ,
1, are opera-tors acting on H ( B ) and we have simplified the notationas | i, j (cid:105) = | i (cid:105) ⊗ | j (cid:105) ∈ H ( A ) ⊗ H ( A ) ( i, j = 0 , QLQ exactly. In addi-tion, we are only interested in quantities of the form
PLQ ( QLQ ) − n QLP . Hence, from the form of
PLQ inEqs.(26), we only need to calculate B − , and B − , .Using Eqs. (25) to (28) in Eq. (5), we can write L and L as: L = PLP − PLQ ( QLQ ) − QLP = | , (cid:105) (cid:104) , | ⊗ (cid:2) A , − C , B − , D , (cid:3) + | , (cid:105) (cid:104) , | ⊗ (cid:2) A , − C , B − , D , (cid:3) (29) L = −PLQ ( QLQ ) − QLP = − | , (cid:105) (cid:104) , | ⊗ (cid:2) C , B − , D , (cid:3) − | , (cid:105) (cid:104) , | ⊗ (cid:2) C , B − , D , (cid:3) (30)Since P ( B ) is a projector of rank 1, we can write: A , − C , B − , D , = α , P ( B ) A , − C , B − , D , = α , P ( B ) C , B − , D , = β , P ( B ) C , B − , D , = β , P ( B ) (31)where (see Appendix. A): α , ≡ α = (cid:104)(cid:104) ( B ) || (cid:0) A , − C , B − , D , (cid:1) || ( B ) (cid:105)(cid:105) β , ≡ β = (cid:104)(cid:104) ( B ) || (cid:0) C , B − , D , (cid:1) || ( B ) (cid:105)(cid:105) α , = α , β , = β (32)and, see Appendix A: α = − ζ + iξ. (33)With these variables and truncating (4) to the zerothorder, the effective Liouville equation can be written inoperator space as: ddτ ρ A = i ξ (cid:2) σ Az , ρ A (cid:3) + ζ (cid:0) σ Az ρ A σ Az − ρ A (cid:1) (34)Note that approximating the expression of ξ and ζ (givenin Appendix. A) by their lowest order in χ gives forEq. (34) the same result as the one obtained by Azouit etal. [12] using a completely different method. Finally, itis straightforward to calculate ( − L ) − exactly , giventhat L (Eq. (30)) is block diagonal:( − L ) − = | , (cid:105) (cid:104) , | ⊗ (cid:20)
11 + β P B + Q B (cid:21) + | , (cid:105) (cid:104) , | ⊗ (cid:20)
11 + β P B + Q B (cid:21) + [ | , (cid:105) (cid:104) , | + | , (cid:105) (cid:104) , | ] ⊗ (2 B ) (35)Let us define the modified initial state of the qubit A as | ˜ ρ ( A )0 (cid:105) = (1 − L ) − ρ ( A )0 : | ˜ ρ ( A )0 (cid:105) = ρ , | , (cid:105) + ρ , | , (cid:105) + ρ , β | , (cid:105) + ρ , β | , (cid:105) , (36)where ρ , , ρ , represent population in the state ρ ( A )0 while ρ , , ρ , represent initial coherences. The physi-cal meaning of this redefined initial density matrix is thestate of qubit A immediately after qubit B has reachedits steady-state. In this case this corresponds to a rescal-ing of the coherences. Using Q B ρ ( B ) s = 0, we can rewriteEq. (6) to describe the dynamics of the slow qubit as: | ρ A ( t ) (cid:105) = e ˜ L t | ˜ ρ A (cid:105) (37)where ˜ L = (1 − L ) − L = α β | , (cid:105) (cid:104) , | + α β | , (cid:105) (cid:104) , | . The evolution operator U ( t ) = e ˜ L t issimple to calculate: U ( t ) = e − ζ (cid:48) t + iξ (cid:48) t | , (cid:105) (cid:104) , | + e − ζ (cid:48) t − iξ (cid:48) t | , (cid:105) (cid:104) , | + | , (cid:105) (cid:104) , | + | , (cid:105) (cid:104) , | (38)where we have defined ζ (cid:48) = −(cid:60) α β , ξ (cid:48) = (cid:61) α β (39)The evolution of the approximated expectation value ofthe Pauli matrices for qubit A and B for γ = 1 and χ = 0 . | φ + (cid:105) = 1 √ | (cid:105) + | (cid:105) ) are compared in Figure 1 to the ones obtained throughan exact full numerical propagation. We see that theadiabatic elimination captures the exact dynamics faith-fully and that indeed qubit B reaches its steady-statebefore any appreciable dynamics in A has taken place.It is worth mentioning that when adiabatic eliminationis valid, i.e. χ (cid:28) γ , the exact final state of the fast qubit B is very close to the steady state of L B : ρ ( B ) ss = lim χ → ρ ( B ) s . (40)Defining the projector P B = || ρ ( B ) ss (cid:105)(cid:105)(cid:104)(cid:104) (2 B ) || leads toconsiderable simplifications in Eq. (14) where the firstterm becomes zero. Thus, taking the interaction QL AB Q to be small, we only need the term n = 0 in Eq. (17) andtaking the zero order only in (4), one can check thatit leads to the same Lindblad operator derived in [12]which is enough to obtain a very good approximation inthe adiabatic limit. B. Open Rabi model
The open Rabi model has been considered recently byGarbe et al. [34] in a quantum metrology context. Itconsist in a spin- (with frequency Ω) interacting withone bosonic mode (frequency ω ) of a cavity describedby the following Hamiltonian: H R = Ω σ z + ω a † a + λ ( a + a † ) ⊗ σ x , (41)where a ( a † ) is the annihilation (creation) operator of thebosonic mode. The dynamics of the open Rabi modelwhere the relaxation of the spin (at a rate Γ (cid:48) ) and thephoton losses from the cavity (at a rate κ (cid:48) ) are taken intoaccount is generated by the following Lindblad operator: L ( ρ ) = − i [ H R , ρ ] + Γ (cid:48) D σ − ( ρ ) + κ (cid:48) D a ( ρ ) (42)where we have used the notation D X ( ρ ) = XρX † − X † Xρ − ρX † X. We also assume that Γ (cid:48) ∼ Ω and κ (cid:48) ∼ ω . If we rescalethe time by dividing the above equation by √ Ω ω , anddefine: g = λ √ Ω ω , η = (cid:114) ω Ω , κ = κ (cid:48) √ Ω ω , Γ = Γ (cid:48) √ Ω ω (43)we rewrite the Lindbladian: L ( ρ ) = L A ( ρ ) + L AB ( ρ ) + L B ( ρ ) (44)with : L A ( ρ ) = − iη [ a † a, ρ ] + κ D a ( ρ ) L B ( ρ ) = − i η [ σ z , ρ ] + Γ D σ − ( ρ ) L AB ( ρ ) = − ig [ (cid:0) a + a † (cid:1) ⊗ σ x , ρ ] . (45) FIG. 1. Top: Evolution of the expectations values of thePauli matrices (cid:104) σ Ax (cid:105) (initial value = 1), (cid:104) σ Ay (cid:105) (initial value=0)and (cid:104) σ Az (cid:105) (always zero) for the slow spin A . Dashed line:Adiabatic elimination, continuous line: exact. Bottom: Evo-lution of the fast spin B. The initial state has been taken as | φ + (cid:105) = √ ( | (cid:105) + | (cid:105) ). It has been shown that, in the limit where ω (cid:28) Ω, thismodel exhibits a quantum phase transition when g in-creases [34–37]. The critical point corresponding to g = 1separates a normal phase ( g <
1) from a superadiantphase ( g > ω (cid:28) Ω corre-sponds to η → g < L B is: || ρ ( B ) ss (cid:105)(cid:105) = | , (cid:105) (46)Following the same steps of the last example, see Ap- pendix. B, we define the projector P B = || ρ ( B ) ss (cid:105)(cid:105)(cid:104)(cid:104) (2 B ) || and we only calculate the term QL B Q , the inverse ofwhich corresponds to QLQ − up to zeroth order inEq. (17). Simple and straightforward calculations leadto the following Lindbladian evolution of the boson: L (cid:16) ρ ( A ) (cid:17) = − i [ ηa † a − g η Γ η + 16 (cid:0) a + a † (cid:1) , ρ ( A ) ]+ κ D a (cid:16) ρ ( A ) (cid:17) + 4 g η ΓΓ η + 16 D ( a + a † ) (cid:16) ρ ( A ) (cid:17) (47)which is exactly the formula derived in [34] using a com-pletely different method, where one should take into con-sideration that the parameters in (42) are double thoseconsidered in [34]. IV. CONCLUSION
We have derived a projection based adiabatic elimina-tion method that works for bipartite systems. This workprovides a direct connection to earlier work on adiabaticelimination of a subspace of the system Hilbert space [24]so that in principle now subsystems as well as sublevelscan be eliminated at the same time. We have illustratedthis with two simple examples of two dispersively cou-pled spins and the open Rabi model. In both cases, us-ing the lowest order approximations, we have obtain thesame expressions that have been previously obtained bycompletely different methods. We expect that this workwill find applications in the case of molecules in cavitieswhere the cavity and part of the molecular levels couldbe adiabatically eliminated.
Appendix A: Detailed Calculation for the 2-qubitsystem
Here we present in detail all the calculations involvedin the example presented in section III A: first we write || ρ ( B ) s (cid:105)(cid:105) in the standard basis: || ρ ( B ) s (cid:105)(cid:105) = 16 χ + γ + 416 χ + γ + 8 | , (cid:105) + 416 χ + γ + 8 | , (cid:105) + 2 γ + 8 iχ χ + γ + 8 | , (cid:105) + 2 γ − iχ χ + γ + 8 | , (cid:105) (A1)where in this appendix we use the notation | i, j (cid:105) = | i (cid:105)⊗| j (cid:105) to alleviate the complexity of mathematical expressions.Then, we define the necessary projectors of the partialtrace, namely: P ( B ) = || ρ ( B ) s (cid:105)(cid:105)(cid:104)(cid:104) ( B ) || (A2)and Q ( B ) = (2 B ) − P ( B ) χ + γ + 8 | Ψ − (cid:105) (cid:104) Ψ + | − | Ψ − (cid:105) (cid:104) , |− iχ + 2 γ χ + γ + 8 | , (cid:105)(cid:104) Ψ + | + | , (cid:105) (cid:104) , | + 8 iχ − γ χ + γ + 8 | , (cid:105)(cid:104) Ψ + | + | , (cid:105)(cid:104) , | (A3)where we have defined the following vectors: | Ψ + (cid:105) = | , (cid:105) + | , (cid:105)| Ψ − (cid:105) = | , (cid:105) − | , (cid:105) (A4) Later on, it will be useful to define | Φ + (cid:105) = | , (cid:105) + | , (cid:105)| Φ − (cid:105) = | , (cid:105) − | , (cid:105)| Θ + (cid:105) = 2 γ + 8 iχ χ + γ + 8 | , (cid:105) + 2 γ − iχ χ + γ + 8 | , (cid:105)| Θ − (cid:105) = 2 γ + 8 iχ χ + γ + 8 | , (cid:105) − γ − iχ χ + γ + 8 | , (cid:105)| Ω + (cid:105) = 16 χ + γ + 416 χ + γ + 8 | , (cid:105) + 416 χ + γ + 8 | , (cid:105)| Ω − (cid:105) = 16 χ + γ + 416 χ + γ + 8 | , (cid:105) − χ + γ + 8 | , (cid:105)| Ω ⊥ + (cid:105) = 16 χ + γ + 416 χ + γ + 8 | , (cid:105) − χ + γ + 8 | , (cid:105)| Ω ⊥− (cid:105) = 16 χ + γ + 416 χ + γ + 8 | , (cid:105) + 416 χ + γ + 8 | , (cid:105) (A5)and the unitary matrix: U SWAP = | , (cid:105) (cid:104) , | + | , (cid:105) (cid:104) , | + | , (cid:105) (cid:104) , | + | , (cid:105) (cid:104) , | (A6)as well. From Eqs. (A1), (21) and (A3), we find that: PLP = | , (cid:105) (cid:104) , | ⊗ A , + | , (cid:105) (cid:104) , | ⊗ A , PLQ = | , (cid:105) (cid:104) , | ⊗ C , + | , (cid:105) (cid:104) , | ⊗ C , QLP = | , (cid:105) (cid:104) , | ⊗ D , + | , (cid:105) (cid:104) , | ⊗ D , + | , (cid:105) (cid:104) , | ⊗ D , + | , (cid:105) (cid:104) , | ⊗ D , QLQ = | , (cid:105) (cid:104) , | ⊗ B , + | , (cid:105) (cid:104) , | ⊗ B , + | , (cid:105) (cid:104) , | ⊗ B , + | , (cid:105) (cid:104) , | ⊗ B , (A7)Where we have defined the following matrices on H (2 B ) A , = 2 iχ (cid:0) χ + γ (cid:1) (16 χ + γ + 8) P ( B ) , A , = U SWAP A , (A8) C , = − iχ | ρ ( B ) s (cid:105)(cid:104) Ω ⊥ + | , C , = U SWAP C , (A9) D , = 16 iχ (cid:0) χ + γ + 4 (cid:1) (16 χ + γ + 8) | Ψ − (cid:105) (cid:104) Ψ + | − iχ (cid:0) χ + γ (cid:1) (16 χ + γ + 8) | Θ + (cid:105) (cid:104) Ψ + | − iχ | Θ − (cid:105) (cid:104) Ψ + |D , = − iχ | Θ − (cid:105) (cid:104) Ψ + | , D , = U SWAP D , , D , = (A10) B , = 2 iχ (cid:0) χ + γ (cid:1) (16 χ + γ + 8) | Ψ − (cid:105) (cid:10) Ω ⊥ + (cid:12)(cid:12) + γ | Ψ − (cid:105) (cid:104) , | − γ | , (cid:105) (cid:104) , | − γ | , (cid:105) (cid:104) , | − | Ψ − (cid:105) (cid:104) Φ + | + | Φ + (cid:105) (cid:104) Ψ − | + 2 iχ (cid:0) χ + γ (cid:1) (16 χ + γ + 8) | Θ − (cid:105) (cid:104) Ψ − | + 64 iχ (4 iχ − γ )(16 χ + γ + 8) | , (cid:105) (cid:104) , | + 16 iχ (4 iχ + γ ) (cid:0) χ + γ + 4 (cid:1) (16 χ + γ + 8) | , (cid:105) (cid:104) , |B , = − | Ψ − (cid:105) (cid:104) Φ + | + | Φ + (cid:105) (cid:104) Ψ − | + γ | Ψ − (cid:105) (cid:104) , | + 4 iχ − γ | , (cid:105) (cid:104) , | − iχ + γ | , (cid:105) (cid:104) , |B , = B , + 4 iχ ( | Θ − (cid:105) (cid:104) Ψ + | + | , (cid:105) (cid:104) , | − | , (cid:105) (cid:104) , | ) , B , = U SWAP B , U SWAP (A11)With this diagonal form of
QLQ , it is straightfor-ward to compute (
QLQ ) − . It consists in computing B − i,j , i, j = 0 ,
1, which are 4 × PLQ ( QLQ ) − n QLP and
PLQ is of the form (26), weonly need to compute B , .To simplify this task, we define the unitary matrix U = 1 √ | , (cid:105) (cid:104) Ψ + |− √ | , (cid:105) (cid:104) Ψ − | + | , (cid:105) (cid:104) , | + | , (cid:105) (cid:104) , | , (A12) and we compute the pseudo-inverse of ˜ B , = UB , .Multiplying by U boils down to replacing the ket | Ψ − (cid:105) by | , (cid:105) in (A11), which implies that ˜ B , can be representedas a 3 × B − , . To find the pseudo inverse of ˜ B , ,we simply solve the set of equations corresponding to:˜ B , ˜ B − , = Π ran[ ˜ B , ] = (2 B ) − | , (cid:105) (cid:104) , | (A13)where Π ran[ ˜ B , ] is the hermitian projector to the rangeof ˜ B , . If we define the following quantities b = − γ (cid:0) χ + γ + 8 (cid:1) √ iχγ (16 χ + γ −
16) + ( γ + 8) (16 χ + γ + 8)) (16 χ + γ + 4) (A14) b = √ (cid:0) χ γ + 256 χ − iχγ (cid:0) χ + γ (cid:1) + γ + 16 γ + 64 (cid:1) (2 iχγ (16 χ + γ −
16) + ( γ + 8) (16 χ + γ + 8)) (16 χ + γ + 4) (A15) b = √ (cid:16) χ γ − iχγ (cid:0) χ + 3 γ + 16 (cid:1) + (cid:0) χ + γ + 8 (cid:1) (cid:17) (2 iχγ (16 χ + γ −
16) + ( γ + 8) (16 χ + γ + 8)) (16 χ + γ + 4) (A16) b = − χ γ + 64 χ γ + 576 χ γ + 1024 χ + 4 iχγ (cid:0) χ + γ (cid:1) + 2 γ + 36 γ + 192 γ + 256 γ (2 iχγ (16 χ + γ −
16) + ( γ + 8) (16 χ + γ + 8)) (16 χ + γ + 4) (A17) b = b − χ (cid:0) χ (cid:0) χ + γ + 4 (cid:1) − iγ (cid:0) χ + γ + 8 (cid:1)(cid:1) γ (2 iχγ (16 χ + γ −
16) + ( γ + 8) (16 χ + γ + 8)) (16 χ + γ + 4) (A18) b = b = 2 √ γ b , b = 2 √ γ b , b = 2 √ γ b (A19)and make the identification | , (cid:105) → | (cid:105) , | , (cid:105) → | (cid:105) and | , (cid:105) → | (cid:105) , then: ˜ B − , = (cid:88) i,j =1 b ij | i (cid:105) (cid:104) j | (A20)We can check that ˜ B − , verify all the Moore-Penrose con-ditions [38]:˜ B , ˜ B − , ˜ B , = ˜ B , , ˜ B − , ˜ B , ˜ B − , = ˜ B − , (cid:16) ˜ B , ˜ B − , (cid:17) † = ˜ B , ˜ B − , , (cid:16) ˜ B − , ˜ B , (cid:17) † = ˜ B − , ˜ B , (A21)Finally, we can easily see that the pseudo-inverse of B , is: B − , = (cid:16) U † ˜ B , (cid:17) − = ˜ B − , U (A22)with that, we have all the necessary ingredients to com-pute L and L . Using Eqs.(A7), we find that: L = PLP − PLQ ( QLQ ) − QLP = | , (cid:105) (cid:104) , | ⊗ (cid:2) A , − C , B − , D , (cid:3) + | , (cid:105) (cid:104) , | ⊗ (cid:2) A , − C , B − , D , (cid:3) (A23) L = −PLQ ( QLQ ) − QLP = − | , (cid:105) (cid:104) , | ⊗ (cid:2) C , B − , D , (cid:3) − | , (cid:105) (cid:104) , | ⊗ (cid:2) C , B − , D , (cid:3) (A24)Since P ( B ) is a projector of rank 1, we can write: A , − C , B − , D , = α , P ( B ) A , − C , B − , D , = α , P ( B ) C , B − , D , = β , P ( B ) C , B − , D , = β , P ( B ) (A25)Because P B is of the form | ρ ( B ) s (cid:105)(cid:104) Ψ + | , we find that: α , ≡ α = 12 (cid:104) Ψ + | (cid:0) A , − C , B − , D , (cid:1) | Ψ + (cid:105) β , ≡ β = 12 (cid:104) Ψ + | (cid:0) C , B − , D , (cid:1) | Ψ + (cid:105) (A26)where we have used the fact that (cid:104) Ψ + | ρ ( B ) s (cid:105) = 1 and (cid:104) Ψ + | Ψ + (cid:105) = 2. We also have: α , = 12 (cid:104) Ψ + | (cid:0) A , − C , B − , D , (cid:1) | Ψ + (cid:105) = 12 (cid:104) Ψ + | (cid:16) U SWAP A , − U SWAP C , U SWAP B − , D , (cid:17) | Ψ + (cid:105) = 12 (cid:104) Ψ + | (cid:16) A , − C , B − , D , (cid:17) | Ψ + (cid:105) = α (A27) where we have used the fact that U = (2 B ) , C , U SWAP = C , , and (cid:104) Ψ + |U SWAP = (cid:104) Ψ + | . In a sim-ilar way we can also show that: β , = β (A28)Let us define α = − ζ + iξ where ζ ≥
0. A tediouscalculation leads to: ζ = 128 χ γ (cid:0) γ + 8 (cid:1) (cid:0) χ + γ + 2 (cid:1) χ γ (16 χ + γ − + ( γ + 8) (16 χ + γ + 8) ξ = 2 χ (cid:0) χ + γ (cid:1) χ + γ + 8 + 256 χ γ (cid:0) χ + γ − (cid:1) (cid:0) χ + γ + 2 (cid:1)(cid:16) χ γ (16 χ + γ − + ( γ + 8) (16 χ + γ + 8) (cid:17) (16 χ + γ + 8) (A29)and β = x + iy x + iy where x = χ (cid:0) χ γ − χ + 6144 χ γ + 2048 χ γ − χ + 192 γ + 1152 γ − γ − (cid:1) y =256 χ γ (cid:0) χ + γ + 4 (cid:1) (cid:0) χ + γ + 8 (cid:1) x = (cid:0) χ + γ + 4 (cid:1) (cid:0) − χ γ + 16 χ γ + 128 χ − χγ + 32 χγ + γ + 16 γ + 64 (cid:1) × (cid:0) χ γ + 16 χ γ + 128 χ + 2 χγ − χγ + γ + 16 γ + 64 (cid:1) y =4 χγ (cid:0) γ + 8 (cid:1) (cid:0) χ + γ − (cid:1) (cid:0) χ + γ + 4 (cid:1) (cid:0) χ + γ + 8 (cid:1) (A30)Finally, it is straightforward to calculate ( − L ) − ex-actly , given that L (A24) is block diagonal:( − L ) − = | , (cid:105) (cid:104) , | ⊗ (cid:20)
11 + β P B + Q B (cid:21) + | , (cid:105) (cid:104) , | ⊗ (cid:20)
11 + β P B + Q B (cid:21) + [ | , (cid:105) (cid:104) , | + | , (cid:105) (cid:104) , | ] ⊗ (2 B ) (A31)Taking the fact Q B ρ ( B ) s = 0 into account and definingthe modified initial state of the qubit A as: | ˜ ρ ( A )0 (cid:105) = ρ , | , (cid:105) + ρ , | , (cid:105) + ρ , β | , (cid:105) + ρ , β | , (cid:105) , (A32)where ρ , , ρ , represent population in the state ρ ( A )0 while ρ , , ρ , represent initial coherences, we can sim-plify (6) to describe the dynamics of the slow qubit as: | ρ A ( t ) (cid:105) = e ˜ L t | ˜ ρ A (cid:105) , (A33) where ˜ L = α β | , (cid:105) (cid:104) , | + α β | , (cid:105) (cid:104) , | . The evolu-tion operator U ( t ) = e ˜ L t is simple to calculate: U ( t ) = e − ζ (cid:48) t + iξ (cid:48) t | , (cid:105) (cid:104) , | + e − ζ (cid:48) t − iξ (cid:48) t | , (cid:105) (cid:104) , | + | , (cid:105) (cid:104) , | + | , (cid:105) (cid:104) , | (A34)where we have defined ζ (cid:48) = −(cid:60) α β , ξ (cid:48) = (cid:61) α β (A35) Appendix B: open Rabi model
In this section, we carry out all the calculations neededto adiabatically eliminate a fast qubit interacting with aslow cavity mode according to the open Rabi model: L ( ρ ) = − i η [ σ z , ρ ] + Γ D σ − ( ρ ) − iη [ a † a, ]( ρ )+ κ D a ( ρ ) − ig [ (cid:0) a + a † (cid:1) ⊗ σ x , ρ ] (B1)The first step is to write L in the super-operator repre-sentation:0 L = − iη (cid:16) ( A ) ⊗ a † a ⊗ (2 B ) − a † a ⊗ ( A ) ⊗ (2 B ) (cid:17) − iη (cid:16) (2 A ) ⊗ ( B ) ⊗ σ z − (2 A ) ⊗ σ z ⊗ ( B ) (cid:17) − ig (cid:16) ( A ) ⊗ (cid:0) a + a † (cid:1) ⊗ ( B ) ⊗ σ x − (cid:0) a + a † (cid:1) ⊗ ( A ) ⊗ σ x ⊗ ( B ) (cid:17) + Γ (2 A ) ⊗ σ − ⊗ σ − − Γ2 (cid:104) (2 A ) ⊗ σ + σ − ⊗ ( B ) + (2 A ) ⊗ ( B ) ⊗ σ + σ − (cid:105) + κa ⊗ a ⊗ (2 B ) − κ (cid:104) a † a ⊗ ( A ) ⊗ (2 B ) + ( A ) ⊗ a † a ⊗ (2 B ) (cid:105) If we define: P B = || ρ Bf (cid:105)(cid:105)(cid:104)(cid:104) ( B ) || , Q B = (2 B ) − P B = | , (cid:105) (cid:104) , | − | , (cid:105) (cid:104) , | + | , (cid:105) (cid:104) , | + | , (cid:105) (cid:104) , | (B2)and P = (2 A ) ⊗ P B , Q = (2 A ) ⊗ Q B (B3)then we can compute the needed quantities for L PLP = (cid:104) − iη (cid:16) ( A ) ⊗ a † a − a † a ⊗ ( A ) (cid:17) + κa ⊗ a − κ (cid:104) a † a ⊗ ( A ) + ( A ) ⊗ a † a (cid:105)(cid:105) ⊗ P B PLQ = − ig (cid:104) ( A ) ⊗ (cid:0) a + a † (cid:1) − (cid:0) a + a † (cid:1) ⊗ ( A ) (cid:105) ⊗ | , (cid:105) (cid:104) Φ + |QLP = − ig (cid:16) ( A ) ⊗ (cid:0) a + a † (cid:1) ⊗ | , (cid:105) (cid:104) Ψ + | − (cid:0) a + a † (cid:1) ⊗ ( A ) ⊗ | , (cid:105) (cid:104) Ψ + | (cid:17) QL B Q = (2 A ) ⊗ (cid:20) − Γ | , (cid:105) (cid:104) , | + Γ | , (cid:105) (cid:104) , | − (cid:18) Γ2 + 2 iη (cid:19) | , (cid:105) (cid:104) , | − (cid:18) Γ2 − iη (cid:19) | , (cid:105) (cid:104) , | (cid:21) (B4) QL B Q represents the dominant term of QLQ which can be inverted quite easily: Q L B Q − = (2 A ) ⊗ − | , (cid:105) (cid:104) , | + 12Γ | , (cid:105) (cid:104) , | − (cid:18) Γ2 + 2 iη (cid:19) | , (cid:105) (cid:104) , | − (cid:18) Γ2 − iη (cid:19) | , (cid:105) (cid:104) , | (B5)Hence we can easily calculate L to be L = (cid:104) − iη (cid:16) ( A ) ⊗ a † a − a † a ⊗ ( A ) (cid:17) + κa ⊗ a − κ (cid:104) a † a ⊗ ( A ) + ( A ) ⊗ a † a (cid:105) − g Γ2 − iη Γ η (cid:104) ( A ) ⊗ (cid:0) a + a † (cid:1) − (cid:0) a + a † (cid:1) ⊗ (cid:0) a + a † (cid:1)(cid:105) − g Γ2 + 2 iη Γ η (cid:104)(cid:0) a + a † (cid:1) ⊗ ( A ) − (cid:0) a + a † (cid:1) ⊗ (cid:0) a + a † (cid:1)(cid:105) ⊗ P B (B6)From which we can deduce the reduced dynamics gov-erning the evolution of the slow system to be: L (cid:16) ρ ( A ) (cid:17) = − i [ H ( A ) , ρ ( A ) ]+ κ D a (cid:16) ρ ( A ) (cid:17) + 4 g η ΓΓ η + 16 D ( a + a † ) (cid:16) ρ ( A ) (cid:17) (B7) where we have defined: H ( A ) = ηa † a − g η Γ η + 16 (cid:0) a + a † (cid:1) (B8) [1] H. Haken, Z Physik B , 413 (1975). [2] Haken, Synergetics–An introduction (Springer Berlin, , 213 (1967).[4] C. Cohen-Tannoudji, Physics Reports , 153 (1992).[5] V. Paulisch, H. Rui, H. K. Ng, and B.-G. Englert, Eur.Phys. J. Plus , 12 (2014).[6] E. Brion, L. H. Pedersen, and K. Mølmer, J. Phys. A:Math. Theor. , 1033 (2007).[7] L. You, X. X. Yi, and X. H. Su, Phys. Rev. A , 032308(2003).[8] D. Nagy, G. K´onya, G. Szirmai, and P. Domokos, Phys.Rev. Lett. , 130401 (2010).[9] J. S. Douglas, H. Habibian, C.-L. Hung, A. V. Gorshkov,H. J. Kimble, and D. E. Chang, Nature Photonics , 326(2015).[10] A. Sinatra, F. Castelli, L. A. Lugiato, P. Grangier, andJ. P. Poizat, Quantum Semiclass. Opt. , 405 (1995).[11] R. Azouit, A. Sarlette, and P. Rouchon,arXiv:1603.04630 [quant-ph] (2016), arXiv: 1603.04630.[12] R. Azouit, F. Chittaro, A. Sarlette, and P. Rouchon,Quantum Sci. Technol. , 044011 (2017).[13] R. Azouit, Adiabatic elimination for open quantum sys-tems , Ph.D. thesis, PSL Research University (2017).[14] R. Azouit, F. Chittaro, A. Sarlette, and P. Rouchon,IFAC-PapersOnLine 20th IFAC World Congress, ,13026 (2017).[15] A. Sarlette, P. Rouchon, A. Essig, Q. Ficheux, andB. Huard, (2020), arXiv:2001.02550.[16] M. Mirrahimi and P. Rouchon, IEEE Transactions onAutomatic Control , 1325 (2009).[17] F. Reiter and A. S. Sørensen, Phys. Rev. A , 032111(2012).[18] E. M. Kessler, Phys. Rev. A , 012126 (2012).[19] Y. Lin, J. P. Gaebler, F. Reiter, T. R. Tan, R. Bowler,A. S. Sørensen, D. Leibfried, and D. J. Wineland, Nature , 415 (2013).[20] V. Albert, K. Noh, and F. Reiterr, arXiv:1809.07324(2019).[21] F. Pastawski, L. Clemente, and J. I. Cirac, Phys. Rev.A , 012304 (2011).[22] O. c. v. ˇCernot´ık, D. V. Vasilyev, and K. Hammerer,Phys. Rev. A , 012124 (2015).[23] F. Minganti, A. Biella, N. Bartolo, and C. Ciuti, PhysicalReview A , 042118 (2018).[24] D. Finkelstein-Shapiro, D. Viennot, I. Saideh, T. Hansen,T. o. Pullerits, and A. Keller, Phys. Rev. A , 042102(2020).[25] I. Lesanovsky and J. P. Garrahan, Phys. Rev. Lett. ,215305 (2013).[26] M. Marcuzzi, J. Schick, B. Olmos, and I. Lesanovsky,Journal of Physics A: Mathematical and Theoretical ,482001 (2014).[27] A. Castellini, H. R. Jauslin, B. Rousseaux, D. Dzsotjan,G. Colas des Francs, A. Messina, and S. Gu´erin, TheEuropean Physical Journal D , 223 (2018).[28] D. Finkelstein-Shapiro, P.-A. Mante, S. Sarizosen,L. Wittenbecher, I. Minda, S. Balci, T. Pullerits, andD. Zigmantas, arXiv:2002.05642.[29] R. F. Ribeiro, L. A. Mart´ınez-Mart´ınez, M. Du,J. Campos-Gonzalez-Angulo, and J. Yuen-Zhou, Chem.Sci. , 6325 (2018).[30] G. Lindblad, Communications in Mathematical Physics , 119 (1976).[31] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan,Journal of Mathematical Physics , 821 (1976). [32] I. Knezevic and D. K. Ferry, Phys. Rev. E , 016131(2002).[33] T. F. Havel, Journal of Mathematical Physics , 534(2003).[34] L. Garbe, M. Bina, A. Keller, M. G. A. Paris, and S. Fe-licetti, Phys. Rev. Lett. , 120504 (2020).[35] M.-J. Hwang, R. Puebla, and M. B. Plenio, PhysicalReview Letters , 180404 (2015).[36] R. Puebla, M.-J. Hwang, J. Casanova, and M. B. Plenio,Physical Review Letters , 073001 (2017), publisher:American Physical Society.[37] M.-J. Hwang, P. Rabl, and M. B. Plenio, Physical Re-view A , 013825 (2018), publisher: American PhysicalSociety.[38] R. Penrose, Mathematical Proceedings of the CambridgePhilosophical Society51