Perturbation Analysis of Quantum Reset Models

This paper is devoted to the analysis of Lindblad operators of Quantum Reset Models, describing the effective dynamics of tri-partite quantum systems subject to stochastic resets. We consider a chain of three independent subsystems, coupled by a Hamiltonian term. The two subsystems at each end of the chain are driven, independently from each other, by a reset Lindbladian, while the center system is driven by a Hamiltonian. Under generic assumptions on the coupling term, we prove the existence of a unique steady state for the perturbed reset Lindbladian, analytic in the coupling constant. We further analyze the large times dynamics of the corresponding CPTP Markov semigroup that describes the approach to the steady state. We illustrate these results with concrete examples corresponding to realistic open quantum systems.


Introduction
A major challenge when investigating small quantum systems is to assess their dynamics when coupled to several environments that put the system in an out-of-equilibrium situation. To do so, one often resorts to effective master equations governing the reduced density operator for the small system. Under the Born-Markov approximation (that involves weak system-bath coupling and short bath time-correlations), the evolution equation for the reduced density operator becomes linear, and is cast into the form of a Lindblad-type master equation [12,13] for the corresponding map to be CPTP (Completely Positive and Trace Preserving). A Hamiltonian approach using perturbation theory is probably the most standard way to derive such a (continuous in time) effective evolution equation for the reduced quantum system [5,28]. For an account of mathematical results, we refer the reader to the review [10] and Communicated by Eric A. Carlen.
B Géraldine Haack geraldine.haack@unige.ch 1 to the recent paper [22] which implements this procedure rigorously for a general class of systems. Alternatively, repeated-interaction schemes (discrete in time) have attracted lots of attention among both mathematicians [1,[7][8][9]14,18] and physicists [2,3,20,25,26,29,31], especially in the context of quantum thermodynamics. Exact solutions for the asymptotic steady states generated by both types of dynamics can in general be derived for quantum systems with low dimensional Hilbert space only.
Appealing master equations to investigate the dynamics of higher dimensional quantum systems are provided by a specific class of models, known as Quantum Reset Models (QRM hereafter). These models can be viewed as a natural extension of classical stochastic models, see [16] for a review and [11] for an example treating diffusion processes. Remarkably, QRM can be formulated in terms of Lindblad master equations so that they generate CPTP maps. This is achieved by making specific choices of dissipation channels (corresponding to a fully depolarized quantum channel), see [15,23,33] for examples in specific physical setups. These QRM, thanks to their structural simplicity, present the strong advantage to allow for analytical solutions for the reduced density operator of multipartite quantum systems and have been successfully exploited to assess the dynamics of specific quantum systems, namely small quantum thermal machines made of a few qubits, qutrits or higher-dimensional quantum systems [4,6,19,30,32,33].
In this work, we raise the question to which extent general properties of the dynamics generated by QRM can be analyzed mathematically. Our aim is to go beyond specific models to determine generic properties of the dynamics of QRM, i.e. induced by the mathematical structure itself of the QRM. A first step in that direction is performed in the recent work [27] where a single system driven by a Lindbladian subject to a reset process is considered. The spectral properties of the total Lindbladian perturbed by the reset processes are established, under the assumption that the unperturbed Lindbladian possesses a unique stationary state. Extensions to certain degenerate unperturbed Lindbladians are also discussed and examplified. In the present paper, we consider QRM describing the dynamics of more complex structures that are therefore intrinsically degenerate and not amenable to the cases dealt with above. We reach a two-fold objective. On the one hand, we show that those degenerate QRM nevertheless allow for a complete mathematical treatment revealing a rich structure. On the other hand, we demonstrate the relevance of our perturbative analysis to assess the dynamics of realistic multipartite quantum systems characterized by Hilbert spaces of dimension as high as 8.
More precisely, our generic model is made of a tripartite structure, A − C − B, where A and B are the two quantum systems subject to reset processes, and C is a central system with its own free evolution. The three subsystems are weakly interacting through a Hamiltonian. We first recall that QRM are always characterised by Lindblad generators, with explicit dissipators. Then we analyse the spectral properties of the resulting Lindbladians and the dynamics of the tri-partite system they generate, under generic hypotheses on the coupling term. We conduct this analysis first in absence of interaction between the A − C − B parts of the Hilbert space they are defined on, which gives rise to an uncoupled Lindbladian displaying large degeneracies, i.e. a large subspace of invariant states. Then, we introduce a generic interaction between these different parts and perform a perturbative analysis in the coupling constant. We prove uniqueness of an invariant steady state under the coupled dynamics, analytic in the coupling constant, and provide a description of the converging power series of this non-equilibrium steady state that develops in the small system. Building up on our spectral analysis, we elucidate the long time properties of the dynamics of the tri-partite system and its approach to the steady state. Finally, we focus on the case where the uncoupled system has no Hamiltonian drive and we describe in particular the emergence of a natural classical Markov process in the description of the large time behaviour of the coupled system. The paper closes with the study of two examples illustrating the key features of this analysis: the systems A and B are two qubits while the central system C is of arbitrary dimension N and the uncoupled dynamics has no Hamiltonian drive. For a rather general choice of QRM coupled dynamics, we compute the leading order of the steady state for N arbitrary and, for N = 2-when C is another qubit-we determine the steady state up to order three in the coupling constant as well as the associated classical Markov process.

Simple Hilbert Space Setup
As a warmup, we consider a single quantum system of finite dimension characterized by its Hamiltonian H defined on its Hilbert space H which is coupled to M reservoirs. QRMs assume the state of the quantum system to be reset to a given state τ l with probability γ l dt within each time interval dt. The QRM-type evolution equation is given by [4,15,19]: (2.1) The operator ρ is the reduced density operator of the system defined on H, and γ l characterizes the coupling rate to the reservoir l, l = 1, . . . , M.
For the sake of comparison with our main concern-tri-partite systems-and to set the notation, we discuss the dynamics of QRM defined in this simple setup, essentially along the lines of [27]. We provide a full description of its generic features, under the following assumptions.

Gen:
Let H be a Hilbert space, with dim H = N < ∞. The dissipative part of the generator is characterised by • {τ l } 1≤l≤M a collection of density matrices on H, i.e. τ l ∈ B(H), with τ l ≥ 0 and tr(τ l ) = 1, for all l ∈ 1, . . . , M, • γ l > 0, l ∈ 1, . . . , M, the collection of associated non-zero rates for the coupling to the M baths.
The generator of QRM is thus the (super-)operator L ∈ B(B(H)) defined by where ρ here is arbitrary in B(H), such that the dynamics of the QRM readṡ In case ρ ∈ DM(H), the set of density matrices DM(H) = {ρ ∈ B(H) | ρ ≥ 0, tr(ρ) = 1}, the trace factor in (2.2) disappears. Indeed, we will see below in a more general framework that the operator L enjoys further symmetries, being a Lindblad operator, see Proposition 3.2; in particular if ρ 0 ∈ DM(H), ρ(t) ∈ DM(H), for all t ∈ [0, ∞). However, we perform the full spectral analysis of L as an operator on B(H) and, accordingly, solve the equation (2.3) without resorting to these symmetries.
We first combine the density matrices τ l with corresponding rates γ l into a single density matrix T with corresponding rate . Setting In the sequel, we denote the matrix elements of any A ∈ B(H) in the basis {ϕ j } 1≤ j≤N by A jk = ϕ j |Aϕ k , and the operator |ϕ ψ| ∈ B(H), for ϕ, ψ ∈ H, is defined by |ϕ ψ| : η → ϕ ψ|η . is diagonalisable with spectrum given by All eigenvalues are simple, except − which has multiplicity N − 1.
Moreover, the solution to (2.3) reads (2.7) Expressed in the eigenbasis of H , this means that, with λ jk = i(e j − e k ) + , Remark 2.2 i) In the limit t → ∞ the steady state is independent of the initial condition and reads On the spectral side, the observation above immediately yields L(|ϕ j ϕ k |) = −λ jk |ϕ j ϕ k | for j = k, showing {−λ jk } j =k are simple eigenvalues by our genericity assumption. To compute the other nonzero eigenvalues of L, we note that if ρ is an eigenvector of L associated with an eigenvalue λ, then λtrρ = 0. Hence λ = 0 implies trρ = 0. Thus, considering the N − 1 dimensional subspace of diagonal traceless matrices in the eigenbasis of H , {ρ = 1≤ j≤N r j |ϕ j ϕ j | | 1≤ j≤N r j = 0}, and making use of the identity L(|ϕ j ϕ j |) = (T − |ϕ j ϕ j |), for any j, we see that it coincides with Ker (L + I).
Finally, the one-dimensional kernel of L is spanned by i[H , ·] + −1 (T ): the inverse is well defined thanks to (2.11), it has matrix elements T jk /λ jk , and trace one. Thus

Tri-partite Hilbert Spaces
We define here the tri-partite systems whose dynamical properties are studied in this paper. while H C is arbitrary at this point. In applications, the reset state τ # will typically be defined as a Gibbs state at some inverse temperature β # associated to H # ; i.e. τ # = e −β # H # /Z # which satisfies (2.13), where Z # is the corresponding partition function. In Sect.3, we perform the analysis of the uncoupled case (system A − C − B is non-interacting), and in Sect.4, we make use of analytic perturbation theory to treat the case where a weak interaction is added to the system A − C − B.

The Non-interacting Tripartite QRM
We define the uncoupled QRM by the generator where I # denotes the identity operator on H # and tr # denotes the operator on the tensor product of Hilbert spaces with indices different from #, obtained by taking the partial trace over H # . For later purposes, tr ## denotes the operator on the Hilbert space with index different from # and # obtained by taking the partial trace over H # ⊗ H # . For example, will be viewed as linear maps. We shall abuse notations and write H # for the Hamiltonian both on H # and H, the context making it clear what we mean. Also, we shall denote the non-Hamiltonian part of the generator by Let us start by a structural result saying that the QRM at time t, e tL (ρ 0 ), with ρ 0 a state, is a CPTP map, by recalling that its generator can be cast under the form of a Lindblad operator, see e.g. [4,15,19]. More precisely, the non-Hamiltonian part of their generator (3.1) takes the form of a dissipator, i.e. (3.5)

Remark 3.3 i)
This result applies to the non-Hamiltonian part of the generator of QRM defined on a simple Hilbert space as well, by considering H C = C, in which case tr A reduces to the scalar valued trace. ii) The operators A jk can be replaced by √ t j |ϕ j ψ k |⊗I C , where {ψ k } k is any orthonormal basis of H A without altering the result.

Spectrum of the Uncoupled QRM
We proceed by analysing the spectrum of the uncoupled QRM L (3.1) in the tri-partite case, making use of the fact that, by construction, the Hamiltonian part of the decoupled QRM commutes with the dissipator as we quickly check:  (3.8) since H B commutes with ϕ A j |⊗I and |ϕ A j ⊗I. Altogether, the dissipator and the Hamiltonian parts of L admit a basis of common eigenvectors that we now determine.
Let us start with the dissipator and its spectral properties.

Proposition 3.4 The dissipator, as an operator on B(H), admits the following spectral decomposition
where the spectral projectors Q # , # ∈ {0, A, B, AB} are given by Moreover, the different spectral subspaces in B(H) are there are only three distinct eigenvalues and the corresponding spectral projector is The dimensions referred to correspond to complex dimensions for B(H). iv) The result essentially follows from the observation that τ A ⊗ tr A (·) and tr B (·) ⊗ τ B are commuting projectors.

Proof
We start with point iv) of Remark 3.5. For any ρ in B(H), , and similarly for tr B (·) ⊗ τ B . Hence the dissipator is a linear combination of two commuting projectors to which we can apply the next Lemma.
The proof of the Lemma is immediate, and in case some eigenvalues coincide, the corresponding spectral projector is simply the sum of the individual projectors. The identifications P = I − τ A ⊗ tr A (·), Q = I − tr B (·) ⊗ τ B , α = −γ A , β = −γ B yield the announced spectral decomposition of the dissipator, together with the explicit spectral projectors. A direct verification then gives the corresponding spectral subspaces.
The eigenvectors of the Hamiltonian part of L are readily computed. For # ∈ {A, B, C}, let {ϕ # j } 1≤ j≤n # be an orthonormal basis of H # of eigenvectors of H # , with associated eigenvalues e # j , 1 ≤ j ≤ n # . The eigenvalues need not to be distinct at that point. We denote by P # j,k ∈ B(H # ), j, k ∈ {1, · · · , n # }, the operators P # j,k = |ϕ # j ϕ # k | that yield a basis of eigenvectors of the Hamiltonian part of (3.1) of B(H)): It remains to take into account the role of the trace in the spectral subspaces of the dissipator to get the sought for common basis of eigenvectors of (3.3). To do so, we introduce the n # − 1 dimensional basis of diagonal (w.r.t. to the eigenbasis of H # ) traceless matrices Together with τ # , the # j 's form a basis of diagonal matrices. Proposition 3.4 then provides the full spectral analysis of the uncoupled QRM . i) The uncoupled reset model Lindbladian L is thus diagonalisable, with eigenvalues located on the (generically) four vertical lines

Proposition 3.7 The vectors listed below form a basis of B(H) consisting in eigenvectors associated with the mentioned eigenvalue of the uncoupled QRM
in the complex plane, symmetrically with respect to the real axis. ii) In particular, the kernel of L is degenerate, since dim Ker L(·) ≥ n C . iii) It is straightforward to generalise this result to the case where the dissipator admits a reset part acting on H C as well, and to the case of a p-partite non interacting system, with p ∈ N arbitrary.
The spectral projectors of L can be constructed explicitly, making use of the next Lemma: (3.14) yield a complete set of rank one projectors onto the span of the corresponding basis vectors of (3.13) so that the composition of any two of them equals zero.

Remark 3.10
The spectral projectors of L corresponding to Proposition 3.7 are then given by the appropriate tensor products of projectors (3.14).

The Weakly-Interacting Tripartite QRM
We consider now the coupled QRM defined by the Lindblad generator on where H = H * ∈ B(H) is a Hamiltonian that effectively couples the different Hilbert spaces H # , while g ∈ R is a coupling constant. We focus on the determination of the kernel of L g , as g → 0, which describes the asymptotic state of the system driven by L g , under generic hypotheses. Then we turn to the consequences for the dynamics generated by L g . By generic hypotheses, we mean that all assumptions we make along the way ensure the coupling is effective enough to lift all degeneracies, so that all accidental degeneracies are eliminated order by order in g.

Leading Order Analytic Perturbation Theory
When g = 0, Proposition 3.7 shows that whatever the properties of the Hamiltonian H C . We shall consider below both cases H C = 0 and H C = 0, which give rise to different results. In case the Hamiltonian H C is trivial, has dimension n 2 C , and the corresponding spectral projector coincides with Q 0 , the spectral projector on Ker D, see Proposition 3.4. In order to avoid accidental degeneracies when H C = 0, we will assume H C satisfies the spectral hypothesis Under this assumption, we have which is of dimension n C . The corresponding spectral projector acts as follows where the projector Diag C : extracts the diagonal part of ρ C within the normalised eigenbasis of H C . Observe that , extracting the offdiagonal part of ρ C within the same basis, yields the complementary projector We also note, for later reference, that Q 0 on B(H) is trace preserving, so that Ran Analytic perturbation theory, see e.g. Chapter II §2 [17], allows us to compute the splitting of the degenerate eigenvalue zero of L 0 by the perturbation gL 1 . Recall here that L g being a Lindblad operator (Proposition 3.2), the following structural constraints hold: (4.8) Moreover, the eigenvalue 0 is semisimple, that is there is no eigennilpotent (Jordan block) corresponding to that eigenvalue in the spectral decomposition of L g . The same is actually true for all eigenvalues sitting on the imaginary axis. Let {λ j (g)} 1≤ j≤m be the set of eigenvalues of L g that stem from the eigenvalue 0 of L 0 , with m = n 2 They form the so-called λ−group for λ = 0, and for g ∈ C \ {0} with |g| small enough, {λ j (g)} 1≤ j≤m are analytic functions of a (fractional) power of g that tend to zero as g → 0. These eigenvalues may be permanently degenerate. For the structural reasons recalled above, one of these eigenvalues, denoted by λ 0 (g), is identically equal to zero, λ 0 (g) ≡ 0, ∀g ∈ C \ 0, and in case λ 0 (g) is degenerate, it is semisimple.
Let us denote by Q 0 (g) the analytic spectral projector of L g corresponding to the set of eigenvalues in the 0-group . It writes for |g| small where 0 is a circle of small radius centered at the origin. Also, since 0 is a semisimple eigenvalue of L 0 , where S 0 is the reduced resolvent of L 0 at 0, satisfying The analytic reduced operator in the corresponding subspace which describes the splitting reads

Lemma 4.1 Under assumption Spec(H
As a consequence, when H C = 0 the splitting is generically described by the order g 2 correction, while in case H C = 0, the non-zero first order correction imposes that the elements of the kernel of Q 0 (g) commute with H τ which, generically, decreases the degeneracy from n 2 C to n C . In both cases, the eigenvalue zero of Q 0 L 1 Q 0 is semisimple. Proof We first compute for any ρ C ∈ B(H C ), using (4.13), One gets the explicit expression for H τ by expressing the partial trace within the eigenbases of τ # . Therefore (4.16) The fact that H C = 0 implies Q 0 L 1 Q 0 = 0 then follows from (4.17) and the identity Let us investigate the next order correction in order to analyse the splitting from the eigenvalue zero. Following [17] we consider the analytic matrix where we observe with (4.10) that Let Q 0 be the eigenprojector onto Ker L 0 . Then the spectrum of Q 0 L 1 Q 0 describes the splitting to order g 2 , see [17], Thm 5.11: forλ In order to proceed, we shall also assume in the sequel that the operator H τ appearing in Lemma 4.1 has generic spectral properties.

Spec(H τ ):
The spectrum of H τ ∈ B(H C ) is simple and the corresponding Bohr frequencies are distinct.
We denote the normalised eigenvectors and eigenvalues of H τ by ϕ τ j and e τ j , 1 ≤ j ≤ n C . Under Spec(H τ ), we get from (4.5) and Lemma 4.1 where Diag τ is the projector that extracts the diagonal part of the matrices expressed in the orthonormal eigenbasis {ϕ τ j }. Therefore where Diag stands here for Diag C (resp. Diag τ ) if H C = 0 (resp. H C = 0). Equivalently, Q 0 L 1 Q 0 is fully characterised by the following linear map. Set Note that is well defined and takes the form Then, the restriction of to Diag B(H C ), which has dimension n C , satisfies We shall abuse notations in the sequel and simply write identifying operators defined on Q 0 B(H) and Diag B(H C ). Hence , a subspace of dimension n C − 1, in keeping with the fact that Ker L g is never trivial. Hence, for the zero eigenvalue of L g to be non-degenerate at second order perturbation in g, we assume the coupling satisfies the assumption.

Coup:
The linear map

Remark 4.2 Assumption
Coup is equivalent to the statement As a consequence,

Remark 4.4
Under assumption Spec(H τ ), the non-zero eigenvalues of L 0 are all simple, of The next order correction, given by the eigenvalue of the operator − Q λ jk (4.32)

Dynamics
We push here the spectral analysis a bit further in order to get sufficient information to analyse the behaviour of the dynamics of the coupled QRM L g (·), as g → 0. We first discuss the richer case H C = 0 and then describe the modifications required for the case H C = 0. Let Q 0 (g) be the spectral projector of L g given by (4.9), and Q 0 (g) = I − Q 0 (g). We have accordingly Since the spectrum of L g is a positive distance away from the imaginary axis, uniformly in g small enough, functional calculus yields the existence of > 0, independent of g, such that where, for H C = 0 under assumption Spec(H τ ), with simple non zero eigenvalues, see Remark 4.4. In case H C = 0 under hypothesis Spec(H C ), L 0 = 0 by Lemma 4.1, so that (4.36) holds with Q 0 = Q 0 and Q λ jk = 0.
in case H C = 0), the long time behaviour of e tL g is controlled by the first term in (4.34) when g is small. This requires addressing the behaviour of the non self-adjoint spectral projectors associated to eigenvalues of L g that vanish as g goes to zero.
Assuming H C = 0, Spec(H C ) and Coup, the same statement holds with Q 0 (g) Moreover, assuming Coup and Spec(H τ ), (respectively Spec(H C )), if H C = 0, (respec- and the corresponding spectral projector Q S 0 (g) is analytic for |g| < g 0 , and satisfies Here and

Remark 4.6
The spectral constraints on Lindblad operators imply, We give conditions ensuring λ (1) jk < 0 in case the model has no leading order Hamiltonian drive, L 0 = D, that we analyse in more details in Sect. 6.

Proof
We consider H C = 0 only, the other case being similar. Thanks to (4.35) and (4.36), perturbation theory applies to L g and yields the analytic projectors Q 0 (g) and Q λ jk (g) converging to Q 0 and Q λ jk respectively, and the analytic simple eigenvalues λ jk (g), such that (4.37) holds. Expanding the first term using Q 0 L 0 = L 0 Q 0 = 0, one gets thanks to (4.26) (4.41) Assumption Coup implies that Q 0 L 1 Q 0 has one dimensional kernel, with associated spectral projector we write Q S 0 . Hence, perturbation theory again ensures the existence of an analytic one dimensional spectral projector Q S 0 (g) of Q 0 (g) L g Q 0 (g) corresponding to the simple zero eigenvalue of Q 0 L 1 Q 0 at g = 0. Necessarily, Q S 0 (g) coincides with the spectral projector onto the nontrivial kernel of L g for all g small enough, which proves (4.38).
Let us turn to the dynamical implications.

Corollary 4.7
Under the hypotheses for H C = 0 above, the following holds for all t ≥ 0 and g real small enough: (4.42) Further assuming max j =k { λ (1) jk } < 0, there exists δ > 0 such that for all t ≥ 0, where the constant in the O is uniform in t ≥ 0 and g small.
where the constants in all O are uniform in t ≥ 0, g small. Under the hypotheses for H C = 0 above, for all t ≥ 0 and g real small enough,

45)
and there exists δ > 0 such that for all t ≥ 0, where the constant in the O is uniform in t ≥ 0 and g small. Moreover, where the constants in all O are uniform in t ≥ 0, g small.

Remark 4.8 0)
The identical statements (4.43) and (4.46) show that 1/g 2 is the time scale of the approach to the asymptotic state, as expected. i) The full evolution can be approximated by the restriction of e tg 2 τ A ⊗ D ⊗τ B to Ran Q 0 , (provided η is larger than the absolute value of the real part of the eigenvalues of τ A ⊗ D ⊗ τ B in case H C = 0). ii) In case L 0 = D, we provide in Sect. 6 an interpretation of the approximate evolution e tg 2 τ A ⊗ D ⊗τ B as a classical continuous time Markov process. iii) Set F = max{| λ| λ ∈ σ ( D )}. When H C = 0, the explicit term in (4.44) is the leading term if F < η, and for times which satisfy 0 ≤ t < 1 +F | ln(g)|/g 2 , as g → 0, for any > 0. When H C = 0, the same is true for the explicit term in (4.47), without constraint on F. iv) This corollary is relevant for an analysis along the lines of [21].
Proof Again we prove the statements for H C = 0 only, the other case being similar. The first two statements follow from functional calcul, and Proposition 4.5, taking into account the analyticity of the spectral data involved. To get the last statement, we observe that since the CPTP map e tL g has a norm which is uniformly bounded in t ≥ 0 and g small enough, the norm of is bounded above by a constant C > 0 which uniform in t ≥ 0 and g small enough. Thus, by Duhamel formula (4.50) Moreover, η = min j =k {| λ (1) jk |} > 0 immediately implies upon expanding Q 0 (g), where the constants in all O are uniform in t ≥ 0 and g small. Finally, Q 0 L 1 Q 0 = τ A ⊗ D ⊗ τ B allows us to express the exponential in terms of that of D .

Construction of the Asymptotic State
We now turn to the determination of the state ρ 0 (g) ∈ Ker L g where L g = L 0 + gL 1 ∈ B(B(H)) given by a power series in g where tr(ρ 0 ) = 1 and tr(ρ j ) = 0, ∀ j > 0. Expanding L 0 (ρ 0 (g)) + gL 1 (ρ 0 (g)) ≡ 0, and equating like powers of g we get The way to solve this set of equations, in principle, is as follows. Note that the spectral decomposition of L 0 yields 3) The first equation is solved by picking a trace one element R 0 in Ker L 0 = Q 0 (B(H)), described in Proposition 3.7. The addition of any traceless vector r 0 ∈ Ker L 0 yields an equally good solution for ρ 0 := R 0 + r 0 at that order. The next equation amounts to solve which determines r 0 = Q 0 r 0 up to the addition of an element of Ker Q 0 L 1 Q 0 (Q 0 L 1 Q 0 viewed as an operator on Q 0 (B(H))). Let us assume for the discussion here that Again, the addition of any traceless vector r 1 = Q 0 r 1 ∈ Ker L 0 to that R 1 yields an equally good solution ρ 1 := R 1 + r 1 to that equation. The next order requires L 1 (r 1 . This equation will then determine r 0 completely, under generic hypotheses, as we shall see. Then we proceed by induction.
The case H C = 0 is slightly different, see Lemma 4.1, but is approached in the same spirit. We start by working out the first few steps and then give the general statements about this construction in Theorem 5.2 for H C = 0 and Theorem 5.4 for H C = 0.
Again, the inverse of L 0 on its range is the reduced resolvent S 0 = L −1 To express S 0 , it is enough to consider the spectral decomposition L 0 = k>0 λ k Q k , where λ k = 0 and Q k are the spectral projectors corresponding to Proposition 3.7, while λ 0 = 0 corresponds to the projector Q 0 .

H C = 0
We consider here that H C = 0 and work under the spectral assumption Spec(H τ ) on the self-adjoint operator defined by (4.13). We first work out the orders g 0 and g 1 terms, i.e. ρ 0 and ρ 1 , and then state an abstract result on the full perturbation series in Theorem 5.2. The first equation C can be added to that choice so that Then we compute Q 0 L 1 (R 0 + r 0 ): The condition to solve the equation for R 1 requires r where Diag τ (·) extracts the diagonal part of r (0) C in the normalised eigenbasis of H τ . Thanks to our assumption, we set This is equivalent to the equation on B(H C ) where we note that Diag τ (r (1) C ) is arbitrary. Our hypotheses on H τ imply that which fixes Offdiag τ r (1) C and leaves Diag τ r (1) C open for now. At this point, the formula which defines R 2 makes sense, where R 2 depends parametrically on Diag τ r C . At order two, the contribution is R 2 + r 2 , where r 2 = Q 0 (r 2 ) = τ A ⊗ r (2) C ⊗ τ B is arbitrary. The term Diag τ r (1) C is determined by the requirement that Q 0 (L 1 (R 2 + r 2 )) = 0 necessary to solve for R 3 , i.e.

Remark 5.1 The fact that ρ (0)
C ∈ Ker D implies trρ (0) C = 0, so the assumption that ρ C is a state in the initial step amounts to set a normalisation.
Let us formulate a general result that summarises the foregoing and guarantees the process can be pursued:

20)
for all g ∈ C with |g| < g 0 . We have, see (4.29) and (5.4), and for all j ≥ 1. Moreover, there exists a linear map R :

Remark 5.3 0) Replacing R j and Offdiag τ r
where μ k , M k , N k and m k are the eigenvalues, eigenprojectors, eigennilpotents and algebraic multiplicities appearing in the spectral decomposition of R = N k=1 μ k M k + N k . Hence the radius of convergence is g 0 = 1/ max 1≤k≤N (|μ k |).
ii) In case σ (R) ∩ R * ± = ∅, the steady state ρ 0 (g) is well defined for all g ∈ R * ± . iii) The iteration terminates if and only if R has a zero eigenvalue and ρ 0 belongs to the corresponding eigenspace; see Sect. 8 for examples. iv) The restriction of the invariant state to H C is given by tr AB (ρ 0 (g)) = ρ (0) We provide necessary and sufficient conditions in Proposition 6.1 for Coup to be satisfied in case L 0 = D and H C = 0.
Proof Recall that dim Ker L g = 1 is proven in Theorem 4.3.
We solve the higher orders equations for ρ j = R j + r j of (5.2) with for all j by induction. Let j ≥ 2 and assume R k , r k are given traceless matrices satisfying (5.28) for 1 ≤ k ≤ j − 1 as well as This is the situation we arrived at for j = 2. Consider Q 0 (L 1 (R j+1 +r j+1 )) = 0, a necessary condition to compute R j+2 , which yields Splitting the equation into its diagonal and offdiagonal parts gives The first equation determines C is fully determined and therefore the second equation yields where Diag τ r ( j+1) C remains free, while r j is determined. This finishes the proof of the induction.

H C = 0
We consider here H C = 0 and the necessary modifications to compute the series (5.1) due to the identities Q 0 (·) = τ A ⊗ Diag C (tr AB (·)) ⊗ τ B and Q 0 L 1 Q 0 ≡ 0. (5.36) The first equation in (5.2) yields C ∈ Diag C B(H C ) is free. The condition to solve the second equation is Q 0 L 1 (ρ 0 ) = Q 0 L 1 Q 0 (ρ 0 ) = 0 which is trivially satisfied. Thus, writing ρ 1 = R 1 + r 1 with R 1 = (I − Q 0 )ρ 1 and r 1 = Q 0 ρ 1 , we can solve partially the equation setting The next equation L 0 (ρ 2 ) = −L 1 (ρ 1 ) requires Q 0 L 1 (R 1 ) + Q 0 L 1 (r 1 ) = 0. Thanks to r 1 = Q 0 r 1 and the identity (5.36), this equation reduces to where we used the expression for R 1 and ρ 0 = Q 0 ρ 0 . Thanks to assumption Coup for Thus R 1 is now determined, while the traceless part r 1 = τ A ⊗ Diag C r (1) C ⊗ τ B is not. With the familiar decomposition ρ 2 = R 2 + r 2 with respect to the projector Q 0 , we set and turn to the equation for where we used (5.36) and r 2 = Q 0 r 2 . With (5.40), this is equivalent to Thanks to Coup, we can thus determine In turn R 2 is fully determined while r 2 = τ A ⊗ Diag C r (2) C ⊗ τ B remains to be computed, and for all g ∈ C with |g| < g 0 . We have, see (4.29) and (5.39), and i) The map R can be expressed as Recall that {ϕ τ j } 1≤ j≤n C denotes the normalized eigenbasis of H τ with respect to which the projectors Diag τ and Offdiag τ are defined, and set P τ j = |ϕ τ j ϕ τ j |. Given the definition (4.29) of D , we need to compute for all j, k ∈ {1, . . . , n C } Thanks to Proposition 3.4, we can express L −1 0 = D −1 in a compact way. Letρ 0 ∈ B(H) such that tr AB (ρ 0 ) = 0, so that Q 0 (ρ 0 ) = 0. Thus Therefore, introducing and making use of Then we note using the cyclicity of the trace that where the operator in the first trace reads while the second trace yields the j j element of its partial tr AB . Hence, Similar considerations can be made for the traces of the other two operators in (6.5): (6.9) and we eventually obtain where D is viewed as a matrix on C n C , and any diagonal matrix r = n C k=1 r k P τ k ∈ Diag τ B(H C ) is viewed as a vector r 1 r 2 · · · r n C t of C n C .
We provide a necessary and sufficient condition on the coupling Hamiltonian H in terms of the diagonal matrix elements of h(k), 1 ≤ k ≤ n C for assumption Coup to hold, i.e. that D restricted to diagonal traceless matrices is invertible.
and consider the non negative operators {h(k)} 1≤k≤n C defined by (6.11). Assumption Coup holds if and only if there exists j ∈ {1, . . . , n C } such that h j j (k) > 0 for all 1 ≤ k = j ≤ n C .

Remark 6.2 i) Since h(k)
is a sum of non negative operators, it is sufficient to check the condition on any of its constituants. ii) Explicit computations show that for dim H C = 2, assumption Coup holds as soon as s = k and (r , s) = (k, j).
Proof Within the framework introduced above we identify D with its matrix ( D ) jk . We need to show it admits zero as a simple eigenvalue, which amounts to showing that Rank D = n C − 1.
We use the short hand notations h j (k) = h(k) j j ≥ 0 for j = k and h j ( j) = k = j h k ( j) ≥ 0 to express the matrix elements of − D . The proof follows once we establish the following Lemma

Remark 6.4
It is possible that Rank h = n − 1 and one diagonal element h j ( j) = 0, in which case he j = 0, where e j is the j th canonical basis vector of C n . We can associate to h a stochastic matrix p the elements of which are

in its kernel
Proof We know 0 ∈ σ (h) and by Jacobi's formula, where com(A) is the comatrix of A andÂ jk is obtained by deleting the j th row and k th column of A. In our casê is real valued so that σ (ĥ j j ) = σ (ĥ j j ). Moreover, by definition, for all k = j so that by Gershgorin Theorem where the circle G k centered at h k (k) of radius l =k l = j h l (k) intersects the imaginary axis if and only if h j (k) = 0, in which case the intersection reduces to the origin. Since the determinant ofĥ j j is the product of its complex conjugate eigenvalues, (6.18) yields This ends the proof of the Proposition.

Emergence of a Classical Markov Process
Coming back to Corollary 4.7, we know that for times s.t. 0 ≤ t ≤ 1 F+ | ln(g)|/g 2 , the evolution semigroup e t(D(·)−ig[H ,·]) can be approximated by In the case at hand, D is expressed in the orthonormal basis {|ϕ τ j ϕ τ j |} 1≤ j≤n C as the matrix (6.12) denoted by h in Lemma 6.3. The negative of the transpose h T of h is thus a transition rate matrix or Q-matrix, associated to a classical continuous time Markov chain with finitely many states, see [24]. Therefore we can associate to our quantum probleṁ ρ = D(ρ) − ig[H , ρ] a classical continuous time Markov chain (X t ) t≥0 on the state space {|ϕ τ j ϕ τ j |} 1≤ j≤n C identified with {1, 2, . . . , n} with n = n C , as follows. Let us recall the general framework. The Markov process (X t ) t≥0 is characterised by the probability to find the process in state j at time t ≥ 0, given the process at time 0 is in state i, is denoted by These transition probabilities P(t) = ( p i j (t)) 1≤i, j≤n are solutions to the matrix form forward and backward equations where Q = (q i j ) 1≤i, j≤n is a transition rate matrix such that q ii ≤ 0, q i j ≥ 0 and n j=1 q i j = 0. Hence, with the identification Q = −h T we get the following interpretation Then, the operator e tg 2 D arising in the approximation of e tL g provided in (4.44), describes a (rescaled) continuous time Markov process (X t ) t≥0 on the state space {|ϕ τ j ϕ τ j |} 1≤ j≤n C ≡ {1, . . . , n} such that for all S ≥ 0, Remark 6.6 Therefore, for any S ≥ 0, the transpose of e s D is a stochastic matrix.
Let us note that appearance of a classical Markov process on the eigenstates of the leading order driving Hamiltonian within the derivation of Lindblad generators for open quantum systems is well known. By contrast, in absence of leading order driving Hamiltonian, the state space of the Markov process into play is determined by the eigenstates of the averaged first order Hamiltonian H τ , which takes into account the effects of the reset matrices.
Finally, let us address the computation of the order g 2 corrections (4.32) of the simple eigenvalues λ jk (g) of L g (·) = D(·) − ig[H , ·] given bỹ λ (1) We prove in Appendix that Then, the eigenvalues λ jk (g) of L g , see Proposition 4.5, satisfy where 0 < t A , t B < 1. We consider again a case without leading order Hamiltonian drive, that is H A = H B = H C = 0, while the order g Hamiltonian reads In other words, j |ϕ j ϕ j |, # ∈ {↓, ↑}, and β = N k=1 β k ϕ k . On the one hand, this example shows our hypotheses can be checked for arbitrary N and, on the other hand, it can lead to physically relevant models under additional assumptions, see for instance Sect. 8 where we deal with qubits (N = 2) subject to inter-qubit Coulomb interaction and flip-flop type interaction Hamiltonian.
With these definitions we compute which yields We can choose the real parameters a j so that the generic assumption Spec H τ holds for any choice of 0 < t A , t B < 1.

Leading Order Term
The next step consists in determining the diagonal elements of the nonnegative operators h(k) defined in (6.11), 1 ≤ k ≤ N ; more precisely h j (k) := ϕ j |h(k)ϕ j , for j = k. We first compute Since we do not need the elements ϕ k |h(k)ϕ k , we do not make explicit their contribution, that we generically denote below by c i (k)P k , where c i (k) ≥ 0, i = 1, 2, 3, 4. With this convention, we get for the different elements h(k) is made of Eventually, The offdiagonal elements h j (k), j = k, of the matrix − D immediately follow: let Therefore the matrix form (8.17) of the operator D reads Note that α j = 0 ⇔ S j = 0 and U j = 0 , while β j = 0 ⇔ T j = 0 and V j = 0. Hence, looking at the first row of (7.12), one sees that Coup holds for this model when α 2 α 3 . . . α N −1 = 0 and |β 1 | 2 + |α N | 2 = 0, (7.13) or, looking at the last row, when β 2 β 3 . . . β N −1 = 0 and |β 1 | 2 + |α N | 2 = 0. (7.14) In either cases, this validates the conclusions of Theorem 5.2 on the invariant state and the way to compute it. From now on, we assume that either (7.13) or (7.14) holds. The leading term of the invariant state is determined by the one dimensional kernel of D which turns out to be computable explicitly. We have, noting that S j + T j > 0 for 2 ≤ j ≤ N − 1, The corresponding faithful leading order ρ 0 , i.e. ρ 0 > 0, of the invariant state of the QRM thus reads x k . (7.16) Actually, the following more explicit expressions are true. With we can write Note in particular the generic nontrivial dependence on j of the populations of (the reduced) leading order ρ 0 of the invariant state. Further remarks are in order: • For non zero coefficients α j and β j , • In case we consider thermal states for τ # on C 2 , # ∈ {A, B}, such that t # = 1 1+e −β # E # , with excitation energy E # > 0. We get that t # → 1 when β # → ∞, while t # → 1/2 when β # → 0, which shows that at high temperature, the populations tend to be constant.

Example on
With the previous example considering H C ∈ C N , we could derive the exact expressions of the map D and of the leading order solution. However, going to first order correction and beyond requires considerable effort and would not be enlightening for the reader. This motivates this second example, where we restrict H C to be in C 2 and consider an interaction Hamiltonian H that is appropriate to describe effective physical systems. The goal of this section is twofold. First, we derive explicitly higher order corrections illustrating the theorems of Sect. 5, showing that we can capture the main features of the dynamics with relatively little effort as compared to the complexity of the system. Second, we make a clear connection between a tri-partite quantum reset model and models suitable to describe realistic physical systems.

Model
Explicitly, we consider here a chain of three qubits characterized by their bare energies e A , e B , e C entering H 0 . They are interacting through H . The two Hamiltonians are given by Without loss of generality, we assume the interaction strengths U , J α , J β to be real. This model could be effective for instance for three qubits subject to nearest-neighbour interactions: a Coulomb interaction (set by U ) whenever two adjacent qubits are occupied and to a flip-flop interaction term of the form |01 10| + h.c. that conserves the number of excitations (set by J α , J β with J α = J β ). In the ordered computational basis of the three qubits This model corresponds exactly to the previous example with N = 2 and setting: The ground state for the three qubits is now simply given by |000 and corresponds to |g ⊗ ϕ 2 ⊗ ↑ in the previous example with N = 2. For clarity, we provide the expression of H tot in the form introduced in (7.2) The two ends (A and B) of the 3-qubit chain are weakly coupled to their own thermal baths at inverse temperatures β A and β B with coupling strengths γ A and γ B respectively. Dissipation takes place following QRM . The reset states are assumed to be thermal states defined by the Maxwell-Boltzmann distribution with their respective inverse temperature β # = 1/T # (k B = 1 in the following) in the basis {|0 , |1 }: Note that since the ground state |0 in the C part of the Hilbert space corresponds to | ↑ , the substitution t B → (1 − t B ) is in order to use the results of Sect. 7. Let us remark that this model for a tri-partite open quantum system differs from previous works on reset models in the context of quantum thermodynamics, studying in particular quantum absorption refrigerators and entanglement engines, Refs. [4,30,32]. These models consist of a chain of 2, 3 or N qubits, each of them being coupled to its own thermal bath. Dissipation due to the presence of environments is captured through QRM . In Ref. [30], the steady-state solution for 3 qubits with three different environments is derived analytically, whereas the case of two qubits is fully solved in Ref. [4]. In contrast, in this work, we derive the steady-state solution considering an arbitrary system C only coupled to the two ends A and B of the chain, as long as H C satisfies generic assumptions.

Generic Assumptions
We first check the assumptions for H C and H . The condition Spec(H C ), is trivially satisfied in this case as the spectrum σ (H C ) = {0, e C } is simple with e C = 0. We can then verify is not sufficient to ensure the required non-degeneracy conditions in the 0-subspace of L 0 . We easily verify that the kernel of L 0 has dimension n 2 C = 4 if H C = 0 and n C = 2 if H C = 0.
In the following, we will restrict the derivation of the steady-state solution up to the second order correction assuming no drive, i.e. H A = H B = H C = 0. Let us note that in two dimensions, there is no loss of generality to consider the reset states τ A and τ B defined as thermal states with respect to H A and H B .

Leading Order Solution, No Drive
Under Spec(H τ ) and Lemma 4.1, the first-order-correction projectorQ 0L1Q0 in the 0eigenvalue subspace is fully characterized by the map acting onto H C , see Eq. (4.25) and Theorem 4.3 (8.10) In contrast to the previous example, we can compute explicitly here the map and not only D . To this end, we consider ρ C to be initially in an arbitrary diagonal state (with respect to the eigenbasis of H τ ) Defining the linear form on C 2 (8.12) we find the matrix (ρ C ) ∈ B(H C ) to be given by (with respect to the eigenbasis of H τ ) . (8.13) Note that (ρ C ) is diagonal, so that for this example we have (·) = D (·). In particular is one dimensional, so that Assumption Coup is satisfied. Then Ker D provides the leading order steady-state solution ρ 0 = τ A ⊗ ρ (0) Interestingly, the zeroth order solution is the exact solution in the equilibrium situation, i.e. when τ A = τ B = τ , the state ρ 0 = τ ⊗ τ ⊗ τ satisfies for any g ∈ R (or C) L g (ρ 0 ) = 0, an instance of Remark ii) 5.3.

Remark 8.1
In this example, the matrix D can also be derived directly from the previous example with N = 2, starting from the positive operator h(k): In the basis {|0 0|, |1 1|}, given (8.5), the substitution t B → 1−t B , and according to (7.12), the matrix D reads (8.17) whose kernel in this same basis is generated by the two-dimensional vector Let us note that D , when written as a superoperator acting onto diagonal matrices, takes a diagonal form, see Eq. (8.13).

Underlying Markov Process
We have enough information here to determine the natural two-state classical continuous Markov process associated to the model, according to Theorem 6.5. The state space is denoted by {0, 1} ≡ {|0 0|, |1 1|}, and by (6.24) we need to compute e s D to determine the transition probabilities of the process The spectral decomposition of D in the matrix form (8.17) is easily obtained. Introducing we have with eigenvector associated to the non zero eigenvalue proportional to 1 −1 T . Hence, with spectral projectors Therefore, withs = 2s(ϕ + +ϕ − ) In turn this eventually yields the sought for transition probabilities We stress that in absence of leading order driving Hamiltonian, the state space of the Markov process into play is determined by the eigenstates of H τ , that takes into account the effects of the reset matrices.

Higer-Order Corrections, No Drive
We now illustrate Theorem 5.2 by deriving the converging expansion of the unique invariant state of L g ρ 0 (g) = ρ 0 + g ρ 1 + g 2 ρ 2 + . . . with, ρ 0 = τ A ⊗ ρ (0) and We recall the definitions for convenience For the first-order correction, we start computing R 1 = iL −1 0 ([H , τ A ⊗ρ (0) C ⊗τ B ]) which can be expressed with F = i(|01 10| − |10 01|) (acting on H A ⊗ H C or H C ⊗ H B depending on the context) as We first note that since = D , the expression for Offdiag τ r (1) C reduces to zero: Then, it remains to determine Diag τ r (1) C to get the first order correction in g. Thanks to (8.29) and using (8.28) for R 1 , we compute Hence, the first order correction is simply given by R 1 , ρ 1 = R 1 and we obtain We proceed with the second-order correction and compute R 2 = iL −1 0 [H , R 1 ] . The matrix R 2 is rather complex and we provide the expressions for its diagonal and off-diagonal elements separately. Its 8 diagonal elements in the ordered basis (8.3) are proportional to by For its off-diagonal elements, we introduce F 2 = |01 10| + |10 01| and the coefficient matrices . (8.33) The matrix R 2 can then be written in a compact form For Offdiagr (2) C , we find that it is equal to 0. This leads us to: Diag r The solution up to the second-order correction is then given by We note that coulomb-interaction term like in U starts playing a role when considering the second-order correction.
Acknowledgements GH acknowledges support from the Swiss National Science Foundation through the starting grant PRIMA PR00P2_179748 and the National Center of Competence in Research SwissMap for a stimulating research environment. AJ is partially supported by the Agence Nationale de la Recherche through the grant NONSTOPS (ANR-17-CE40-0006-01), and he wishes to thank the Université de Genève for hospitality during the first stages of this work. Both authors acknowledge support from the Banff International Research Station which hosted the 2019 meeting "Charge and Energy Transfer Processes: Open Problems in Open Quantum Systems" where this project started.
Funding Open Access funding provided by University of Geneva.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
Here A jk denotes the jk of the matrix A ∈ B(H C ) with respect to the basis {ϕ τ j }. Note that the operators in the full traces are non negative, whereas the last term is a priori complex valued.
Similarly, 4) and the analogous formula holds for the term involving H τ B . These expressions allow us to bound below their real part by a non negative quantity, as the next lemma shows.
Lemma 9.1 Under the hypotheses above, we compute ≥ tr (I A ⊗ (I C − P τ j j ))H τ B (τ A ⊗ P τ j j )H τ B (I A ⊗ (I C − P τ j j )) + same with k ↔ j. (9.7)

Remark 9.2 Since
is a non negative operator, we get from (9.1), (9.2) and the Lemma that which proves Proposition 6.7.
Proof We prove the first inequality, the others are similar. Let G = H (τ