{'fraction_non_alphanumeric': 0.05290129717610634, 'fraction_numerical': 0.017882880478300327, 'mean_word_length': 3.963963963963964, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Projection based adiabatic elimination of bipartite open quantum systems

Adiabatic elimination methods allow the reduction of the space dimension needed to describe systems dynamics which exhibits separation of time scale. For open quantum system, it consists in eliminating the fast part assuming it has almost instantaneously reached its steady-state and obtaining an approximation of the evolution of the slow part. These methods can be applied to eliminate a linear subspace within the system Hilbert space, or alternatively to eliminate a fast subsystems in a bipartite quantum system. In this work, we extend an adiabatic elimination method used for removing fast degrees of freedom within a open quantum system (Phys. Rev. A 2020, 101,042102) to eliminate a subsystem from an open bipartite quantum system. As an illustration, we apply our technique to a dispersively coupled two-qubit system and in the case of the open Rabi model.


I. INTRODUCTION
Adiabatic elimination is a method whereby the fast degrees of freedom of a system are removed while retaining an effective description of the slow degrees of freedom. This simplification can be very useful to obtain tractable and intuitive equations when only a coarse-grained or long times description is desired [1][2][3][4][5][6][7][8][9], depending on if the target system has a conservative [5,6,10] or a dissipative [11][12][13][14][15][16][17][18] evolution. There are two classes of manifolds on which adiabatic elimination has been applied, i) those that consist of levels within a subsystem, for example the excited states of an atom, and ii) those that consist of a separate subsystem, such as ancillary qubits or measuring devices. For a slow and fast parts described by Hilbert spaces H (A) and H (B) , the first case corresponds to the Hilbert space H = H (A) ⊕ H (B) (direct sum) while the second case corresponds to the Hilbert space H = H (A) ⊗ H (B) (tensor product). Adiabatic elimination is useful in developing protocols for dissipative state preparation in ion traps [19,20], reservoir engineering [21] and the description of measurement devices [22]. The simplicity of the resulting equations can also be computationally advantageous in the study of quantum phase transitions where the size of the system is cumbersomly large [23].
There are several approaches to obtain effective operators, ranging from perturbative expansions of the Liouville operator [17,18], the corresponding Kraus maps [11][12][13][14]16], the resolvent [24], or using stochastic methods [22]. Eliminating a fast subsystem (that forms a tensor product with the slow subsystem) is typically done with a partial trace over the fast subsystem. This can result in a set of hierarchical equations that allows error estimation and correcting the approximation as the slow * arne.keller@u-psud.fr and fast timescales get closer. Importantly, the expansion can be built to preserve the Lindblad structure and as a consequence the physicality of the map [12,25,26]. The procedure for eliminating sublevels within a subsystem (direct sum with the slow subsystem) is best carried out with Feshbach projections [17,24]. However, as the fast-slow separation breaks down, or when incoherent pumping channels exist, the population of the fast subsystem becomes non-negligible (i.e. there can be a finite fraction of population in the excited states). When this happens, the exact time evolution of the slow part becomes non-trace preserving. The loss of trace can be corrected using contour integral methods [24]. It would be however advantageous to have a method that can handle both classes of fast manifolds. This is more important considering that systems from atomic physics are inspiring a number of chemical versions that have much more complicated Hamiltonians and it would be ideal to transform them into effective operators for a direct comparison to the atomic physics counterparts [27][28][29].
In this work, we extend the methodology developed in Ref. [24] to bipartite open quantum systems whose dynamics are described by a Lindblad operator [30,31]. We use the projection operator method suggested by Knezevic and Berry [32] in order to derive equations for a slow subsystem A coupled to a fast subsystem B. The paper is organized as follows. We first recall the main results of Ref [24]. We then apply it to the general bipartite case to obtain a recipe for describing the slow subsystem. Finally, we illustrate the method in the case of a spin dispersively coupled to a second highly dissipative driven spin and to describe the dynamics of the open Rabi model.

A. Adiabatic elimination through projectors techniques
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 generated by a Lindblad operator L:ρ(t) = Lρ(t). We define the Hilbert space H of operator O on H, equipped with the scalar product tr O † 1 O 2 . We first recall the main results presented in Ref [24] related to projector techniques. Let P the projector such that ρ s (t) = Pρ(t) is the the slow part of the system and write Q = 1 1 − P, where 1 1 is the identity operator on H. Let G(z) = (z − L) −1 be the resolvent of the Lindblad operator L. Operators like P, Q, L or G are operators on H. They are sometimes called super-operators. They are here denoted with calligraphic letter, to distinguish them from operators on H (belonging to H), like the density matrix ρ.
We define the effective Lindblad operator L eff (z), such that PG(z)P = (z − L eff (z)) −1 . The effective Lindblad operator L eff (z) can be written as : where For any ρ(t = 0), such that Qρ(t = 0) = 0, the slow dynamics inside PH can be obtained with the inverse Laplace transform as: where ρ s (t = 0) = Pρ(t = 0) and the integral on the complex plane is performed on a straight line D = {z ∈ C; z = a > 0}. At this point no approximation has been made. Eq. (3) gives the exact dynamics inside PH, as long as the initial condition is also inside PH, that is Qρ(t = 0) = 0. Because L eff (z) captures the effect of the dynamics in QH, it is a nonlinear superoperator, in the sense that (z − L eff (z)) ν = 0 is a nonlinear eigenvalue problem. The approximation of a slow dynamics of Pρ(t), with respect to the fast dynamics of Qρ(t) is equivalent to considering the dynamics inside PH in the vicinity of the stationary state reached in the limit t → ∞. In this long time limit only the z → 0 limit will contribute to the inverse Laplace transform of Eq. (3). We thus approximate L eff (z) to the lowest relevant order: where L 0 = L eff (z = 0) and L n = 1 n! d dz L eff (z) z=0 . Using the expression of L eff (z) given by Eq. (1) allows to express L 0 and L n as: In this work, we consider the approximation given by Eq. (4) with n ≤ 1 only, which is a standard approximation for most of the effective operators that are calculated explicitly. The systematic study of higher order approximations (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 computed explicitly. We obtain: The stationary state ρ f of this dynamics, reached at t → +∞, is in the kernel of L 0 . We note that although the dynamics described by Eq. (6) is an approximation, the final reached stationary state ρ f is the exact one.
To conclude, after the adiabatic elimination of the fast part, the generator of the slow dynamics is approximatively given by where L 0 and L 1 can in principle be computed using Eq. (5). The hard part in these equations is the evaluation of the inverse (QLQ) −1 , which in the most general case, as we will see later, can only be achieved through a perturbative expansion.
Theses results are very general, and require only the definition of a projector P and that the initial state fulfills the condition Qρ(t = 0) = 0. We note that P doesn't have to be hermitian, that is the projection does not need to be orthogonal.
In Ref. [24], this formalism was applied to the case where the slow and fast degrees of freedom correspond to a partition of the underlying Hilbert space in two complementary sub-spaces, that is H = H (A) ⊕ H (B) . In the next section we will adapt this general formalism to the bipartite case where H = H (A) ⊗ H (B) .

B. Adiabatic elimination in a bipartite system
We suppose that the state of the bipartite system at time t is described by a density operator ρ (AB) (t) acting on the Hilbert space H = H (A) ⊗ H (B) . We consider that the dynamics of subsystem A is very slow compared to the dynamics of subsystem B. We suppose that the exact stationary state in H is unique and that it is a product state ρ 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 consider that at t = 0, ρ (AB) (t = 0) = ρ (A) 0 ⊗ρ b . In other word, we consider that B reaches its steady state instantaneously in the time scale of the subsystem A. We thus define P as where tr B [] denotes the partial trace over B. The reduced density operator ρ (A) (t) in H (A) can then be obtained as ρ (A) (t) = tr B Pρ (AB) (t) .
For the purpose of simplifying some expressions and calculations, it will be useful to use the operator-vector isomorphism [33] which maps each element of H to a vector in H ⊗ H as follows. An operators such as |a b| ∈ H is mapped to the vector |b ⊗ |a in the H ⊗ H Hilbert space, where |b is the complex conjugate of |b . Consequently, any n × n density matrix ρ ∈ H is mapped to a column vector ||ρ ∈ H ⊗ H, with n 2 elements, by stacking the columns of the ρ matrix. Under this isomorphism, super-operators on H are mapped to operators on H ⊗ H. In particular, the super-operator O performing In this way, the scalar product tr ρ † 1 ρ 2 between two operators ρ 1 and ρ 2 in H is equal to the usual scalar product ρ 1 |ρ 2 in H ⊗ H. Some useful remarks can be made. The identity operator 1 1 in H is mapped to the maximally entangled state We also note that the usual density matrix normalization tr [ρ] = 1 does not correspond to the normalization induced by the scalar product tr ρ 2 = 1 (except in the case of a pure state). Using the previous remark, we have that tr [ρ] = tr [1 1ρ] = 1 is mapped to 1 1|ρ = 1. For our bipartite case, H = H (A) ⊗ H (B) . Therefore, an operator in H as |a 1 a 2 | ⊗ |b Consequently, the operator Pρ (AB) ∈ H is mapped to the vector 1 1 (B) |ρ (AB) ⊗||ρ b ∈ H ⊗H. The projector P acting on H ⊗ H can thus be written as: where 1 1 (2A) is the identity operator on H (A) ⊗H (A) . The operator Q = 1 1 − P simply reads as: In general, the Lindblad operator of the system can be split in 3 terms as follows: 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: Where we have used the fact that 1 1 (B) ||L B = 0, as L B is a trace preserving operator. In some specific cases, it can be a good approximation to take as ρ B , a stationary state of L B . In that case L B P B = 0 and QLP in Eq. (14) can be simplified as For computing L 0 using Eq. (6), the main difficulty resides in the inversion of QLQ. In general this inversion can not be done explicitly, but a perturbative expansion can give a good approximation. In many cases QLQ can be divided into two terms, QLQ = D + V, where the computation of D −1 is easy, and V is small compared to D. In that case we can write Retaining the terms up to n = 0 or n = 1 may give a good approximation. In some cases, the term (15) is dominant over the two other terms. In that case, approximating the inverse of QLQ by can be sufficient as we will see in the examples in the next section. So, L 0 can be approximated by the following expression: This is the main result of this work.

III. EXAMPLES
We apply the formalism of the preceding section to two examples. We first address the case of a strongly dissipative driven qubit B dispersively coupled to a target qubit A. Then, as a second example, we consider the open Rabi model in the regime where the dynamics of the spin is very fast compared to the boson frequency.
A. A 2-qubit system This two qubits system has been considered previously by Azouit et in Ref. [12] to test another method of bipartite adiabatic elimination (note that in their work, it is the A spin which is the strongly dissipative spin). It consists in a strongly dissipative driven qubit B dispersively coupled to a target qubit A. This model is used in Ref. [15] to describe the continuous measurement of a harmonic oscillator excitation number (corresponding to system A) by a spin (corresponding system B). The Lindblad equation for the bipartite system can be written as: where σ + = |1 0|, σ − = σ † + and |0 , |1 are the eignvectors of the Pauli matrix σ z with eigenvalues −1, 1 respectively. Defining the new parameters τ = ut, χ = χ u , and γ = γ u and using the column-vector isomorphism, we write the superoperator form of the Liouvillian as: where we have used the relations From now on, we will drop the prime in all the parameters γ , χ . As the qubit A only comes into play through the Hamiltonian term −iχ[σ A z ⊗ σ B z , ρ] (see Eq. (20)), the kernel of the Lindblad operator is two dimensional, according to the two eigenvectors σ A z . Hence the kernel can be considered as the span of {ρ In order to avoid any unnecessary complications, we will assume that the qubit A has an extremely slow dissipation rate that we will omit from our calculations, but will ensure the uniqueness of the steady state. In other words, the steady state of the considered system will be ρ ss = ρ s0 . We thus assume that the initial state is of the form ρ 0 = ρ s0 , and we define the projector P (see Eq. (9)): where according to Eq. (8), ||1 1 (B) = |0 ⊗ |0 + |1 ⊗ |1 . This ensures that Q = 1 1 − P verifies Q|ρ 0 = 0 as it should be. In Appendix A, we calculate in detail all the quantities necessary to compute L 0 and L 1 as given by Eq. (5) : where X i,j , X ∈ {A, B, C, D}, and i, j ∈ 0, 1, are operators acting on H (B) and we have simplified the notation as |i, j = |i ⊗ |j ∈ H (A) ⊗ H (A) (i, j = 0, 1), see Appendix. A. From the block diagonal form of Eq. (28), it is relatively easy to invert QLQ exactly. In addition, we are only interested in quantities of the form PLQ(QLQ) −n QLP. Hence, from the form of PLQ in Eqs. (26), we only need to calculate B −1 0,1 and B −1 1,0 . Using Eqs. (25) to (28) in Eq. (5), we can write L 0 and L 1 as: Since P (B) is a projector of rank 1, we can write: where (see Appendix. A): and, see Appendix A: With these variables and truncating (4) to the zeroth order, the effective Liouville equation can be written in operator space as: Note that approximating the expression of ξ and ζ (given in Appendix. A) by their lowest order in χ gives for Eq. (34) the same result as the one obtained by Azouit et al. [12] using a completely different method. Finally, it is straightforward to calculate (1 1 − L 1 ) −1 exactly , given that L 1 (Eq. (30)) is block diagonal: Let us define the modified initial state of the qubit A as |ρ where ρ 0,0 , ρ 1,1 represent population in the state ρ (A) 0 while ρ 0,1 , ρ 1,0 represent initial coherences. The physical meaning of this redefined initial density matrix is the state of qubit A immediately after qubit B has reached its steady-state. In this case this corresponds to a rescaling of the coherences. Using Q B ρ (B) s0 = 0, we can rewrite Eq. (6) to describe the dynamics of the slow qubit as: whereL 0 = (1 − L 1 ) −1 L 0 = α 1+β |0, 1 0, 1| + α 1+β |1, 0 1, 0|. The evolution operator U (t) = eL 0t is simple to calculate: where we have defined The evolution of the approximated expectation value of the Pauli matrices for qubit A and B for γ = 1 and χ = 0.1 and with the initial state taken as are compared in Figure 1 to the ones obtained through an exact full numerical propagation. We see that the adiabatic elimination captures the exact dynamics faithfully and that indeed qubit B reaches its steady-state before any appreciable dynamics in A has taken place.
It is worth mentioning that when adiabatic elimination is valid, i.e. χ γ, the exact final state of the fast qubit B is very close to the steady state of L B : Defining the projector P B = ||ρ (B) ss 1 1 (2B) || leads to considerable simplifications in Eq. (14) where the first term becomes zero. Thus, taking the interaction QL AB Q to be small, we only need the term n = 0 in Eq. (17) and taking the zero order only in (4), one can check that it leads to the same Lindblad operator derived in [12] which is enough to obtain a very good approximation in the adiabatic limit.

B. Open Rabi model
The open Rabi model has been considered recently by Garbe et al. [34] in a quantum metrology context. It consist in a spin-1 2 (with frequency Ω) interacting with one bosonic mode (frequency ω 0 ) of a cavity described by the following Hamiltonian: where a (a † ) is the annihilation (creation) operator of the bosonic mode. The dynamics of the open Rabi model where the relaxation of the spin (at a rate Γ ) and the photon losses from the cavity (at a rate κ ) are taken into account is generated by the following Lindblad operator: where we have used the notation We also assume that Γ ∼ Ω and κ ∼ ω 0 . If we rescale the time by dividing the above equation by √ Ωω 0 , and define: we rewrite the Lindbladian: with : It has been shown that, in the limit where ω 0 Ω, this model exhibits a quantum phase transition when g increases [34][35][36][37]. The critical point corresponding to g = 1 separates a normal phase (g < 1) from a superadiant phase (g > 1). Here we show that our method can be used to obtain an effective Lindblad operator for the boson in the normal phase after the elimination of the fast spin. After rescaling, the adiabatic limit ω 0 Ω corresponds to η → 0.
In the normal phase (g < 1), the steady state of the system, which is the kernel of the Lindblad operator given by Eq. (42), is separable and unique [34]. It is straight forward to verify that the steady state of L B is: Following the same steps of the last example, see Ap-pendix. B, we define the projector P B = ||ρ (B) ss 1 1 (2B) || and we only calculate the term QL B Q, the inverse of which corresponds to QLQ −1 up to zeroth order in Eq. (17). Simple and straightforward calculations lead to the following Lindbladian evolution of the boson: which is exactly the formula derived in [34] using a completely different method, where one should take into consideration that the parameters in (42) are double those considered in [34].