Resummation for Nonequilibrium Perturbation Theory and Application to Open Quantum Lattices

Lattice models of fermions, bosons, and spins have long served to elucidate the essential physics of quantum phase transitions in a variety of systems. Generalizing such models to incorporate driving and dissipation has opened new vistas to investigate nonequilibrium phenomena and dissipative phase transitions in interacting many-body systems. We present a framework for the treatment of such open quantum lattices based on a resummation scheme for the Lindblad perturbation series. Employing a convenient diagrammatic representation, we utilize this method to obtain relevant observables for the open Jaynes-Cummings lattice, a model of special interest for open-system quantum simulation. We demonstrate that the resummation framework allows us to reliably predict observables for both finite and infinite Jaynes-Cummings lattices with different lattice geometries. The resummation of the Lindblad perturbation series can thus serve as a valuable tool in validating open quantum simulators, such as circuit-QED lattices, currently being investigated experimentally.


I. INTRODUCTION
Lattice models describe particles or spins residing on a set of sites, arranged in a regular fashion. Different types of interactions among these components are possible and can be included in the formulation of the model. For this reason, lattice models can cover a large arena of physical systems and phenomena. Prominent examples are the fermionic Hubbard model [1][2][3], the Bose-Hubbard model [4,5], and the Heisenberg model [6,7].
The study of open quantum lattices tends to be challenging. Analytical and numerical techniques for open lattices are currently less developed than their counterparts for closed lattices. While for a large class of open quantum lattices the Lindblad master equation provides an appropriate theoretical framework [39,40], numerical methods for solving this master equation directly, such as exact diagonalization [41,42], time evolution, or averaging of quantum trajectories [43], are computationally demanding and become quickly infeasible as lattice size increases. More sophisticated numerical techniques have been suggested and are further being developed, including matrix-product [44,45] and variational methods [46,47]. These methods can handle larger lattices to some degree, but come with their own specific drawbacks.
As a result, the development of quantum simulators based on photons is particularly intriguing. Photonic systems represent an interesting open-system complement to the wellestablished paradigm of ultracold-atom quantum simulators. Since photons do not possess a chemical potential (however, see Ref. [48] for a proposal to engineer a chemical potential), realistic photonic lattices will typically include coherent driving and photon loss [49]. Such systems will thus be a particular useful tool to better understand, gain intuition, and ultimately devise tractable effective models for open quantum lattices of interest. Experiments with photonic quantum simulators will shed definitive light on both dynamical and steadystate phenomena by employing well-defined artificial lattice structures and systematically controlling parameters including drive strength, photon frequency, and strength of the mediated photon-photon interaction. Very promising experimental progress in this direction has already been made in the circuit QED architecture [49][50][51].
An experimental quantum simulator requires careful initial steps of validation [52] to ensure that the given physical system is correctly implementing the intended model. The validation procedure naturally demands that, for specific parameter regimes, theoretical understanding and reliable quantitative predictions are available and enable a comparison between theory and the experimental data obtained from the quantum simulator. For the purpose of validation, we utilize the well-controlled approximation scheme of Lindblad perturbation theory [28,[53][54][55][56]. We are particularly interested in the steady-state behavior of open lattices which can be directly related to experimental observables such as microwave transmission in circuit-QED lattices [50,51]. Steady-state quantities are of paramount importance for the detection of dissipative phase transitions [11, 24-26, 28-30, 32, 34, 36, 57-61]. In the work presented here, we take a crucial step beyond finite-order perturbation theory by demonstrating a partial resummation of the perturbation series. We then employ this method to study an open Jaynes-Cummings lattice [Fig. 1] and establish that the resummation affords a significant im- provement of the approximation accuracy. We illustrate the method's versatility in handling both finite-size and infinite lattices as well as different geometries and dimensionalities in a natural way. The method is hence well suited for validating data from the first circuit-QED quantum simulators currently being investigated [62].
Our discussion is structured as follows. We set the stage with a general review of (Markovian) open quantum lattices in Section II, examining their theoretical description in terms of the Lindblad master equation (II A), and the use of thirdquantization methods [63][64][65] to obtain exact solutions for non-interacting open quantum lattices (II B). Section III forms the centerpiece of our theoretical framework: here, we introduce the resummation scheme by which we can include certain perturbative corrections up to infinite order. We then apply this technique to the open Jaynes-Cummings (JC) lattice model in Section IV. After a few preparatory steps (IV A), we perform the resummation and obtain results for steady-state observables (IV B). We finally compare our results to both finite-order perturbation theory and exact solutions where available, and discuss the role of finite-size effects and lattice structures (IV C). We conclude with a summary in Section Section V.

II. BACKGROUND
A. Lindblad formalism for driven, dissipative lattices Open quantum lattice models are widely used to study many-body physics under nonequilibrium conditions. There exists a large variety of such lattice models, including open photon lattices [12,13,66], Jaynes-Cummings lattices [26,29,30,60,67,68], Dicke-type models [24,25,[69][70][71] and lattices with different types of interactions [33,72,73]. Open lattices are not limited to bosons but may also involve fermions and spins, and may include both on-site and off-site interactions. We shall denote the Hamiltonians governing the unitary evolution from on-site and off-site terms by h r and V rr , re-spectively. The resulting generic system Hamiltonian is then given by In many cases, off-site interactions can be limited to nearestneighbor pairs r, r . For an open quantum lattice, we further account for the coupling to environmental degrees of freedom. Considering the effective dynamics of the lattice only, one finds that the environment generally induces non-unitary evolution which cannot be captured by an effective lattice Hamiltonian. The environment induces effects such as dissipation and decoherence, so that the time evolution of the reduced density matrix of the lattice ρ generally deviates from the unitary von-Neumann equationρ = −i [H, ρ].
In many relevant cases of weak coupling to the environment, the lattice system will undergo Markovian dynamics: the state of the system at the current time t fully determines the state at the slightly advanced time t + dt in the future. (Another way to express this is to say that there are no "memory effects.") Under a fairly general set of conditions [39], the time evolution of ρ is then described by the Lindblad master equation [39,40], The influence of the environment is thus encoded in the damp- Typically, each particular jump operator f j points to a particular non-unitary process. For example, photon decay commonly results from the photon annihilation operator a acting as the jump operator. (This statement is over-simplifying matters somewhat. In general, care must be taken to derive the appropriate jump operator for a particular system and environment coupling [39,74].) The prefactor γ j characterizes the rate of the damping process.
As a concrete example of an open quantum lattice and paradigm for interacting photon lattices, we consider the κ g Γ γ FIG. 2. Constituents of the Jaynes-Cummings lattice: two-level systems are driven coherently with strength , and exchange excitations with local harmonic oscillators at a rate set by g. Together, a pair of two-level system and oscillator form one site of the Jaynes-Cummings system. Nearest-neighbor sites are coupled by propagation of oscillator excitations with rate ∼ κ. Dissipation is included in the form of two-level relaxation (rate Γ) and oscillator relaxation (rate γ). In the circuit-QED realization, two-level systems are implemented by superconducting qubits, and oscillators by photonmodes of on-chip microwave resonators with typical frequencies in the range of a few GHz. driven, damped Jaynes-Cummings lattice [Fig. 2]. In this system, each site consists of a harmonic oscillator, such as the mode of an electromagnetic resonator, coupled coherently to a two-level system, referred to as "qubit" in the following. Each site can be driven by a coherent tone. For simplicity, we consider the situation of a global drive frequency ω d , identical on each site.
The single-site Hamiltonian for this lattice is the Jaynes-Cummings Hamiltonian plus drive term, h r = δω r a † r a r + δΩ r σ + r σ − r + g r a r σ + r + a † r σ − r (3) In the usual way, we have already eliminated the original time dependence of the drive by switching to a frame rotating with the drive frequency. Consequently, the photon and qubit terms involve frequency detunings relative to the drive. Specifically, we have δω r ≡ ω r − ω d for the photon mode with frequency ω r , and δΩ r ≡ Ω r − ω d for the qubit with frequency Ω r on site r. Photon and qubit excitations are created (annihilated) by the standard ladder operators a † r and σ + r (a r and σ − r ), respectively. We denote the strengths of the Jaynes-Cummings coupling and of the coherent drive tone by g r and r . Off-site terms of this lattice arise from the hopping of photons between sites r and r with rate κ rr : The full Hamiltonian consists of the appropriate sum of onsite and off-site contributions (as given above). Finally, we consider the damping induced by photon decay (rate γ r ) and qubit relaxation (rate Γ r ). If both processes occur due to coupling to separate zero-temperature baths, then the appropriate jump operators can be shown to be the photon annihilation operator a r and the pseudo-spin lowering opera-tor σ − r . Overall, we thus obtain the Lindblad master equatioṅ We will frequently find it convenient to write the Lindblad master equation in the short formρ = Lρ. Here, L is the so-called Liouville superoperator. The term "superoperators" designates an object mapping an ordinary Hilbert-space operators such as ρ to another Hilbert-space operator. (In our notation, we distinguish superoperators, operators, and real/complex numbers by using double-stroke, sans-serif, and regular lettering, respectively.) Using this shorthand of the master equation, we can easily characterize the steady state ρ s of the lattice: it is the stationary solution of this Lindblad master equation, i.e.ρ s = Lρ s = 0 with normalization tr ρ s = 1.
Mathematically, the equation Lρ s = 0 is a system of linear equations for the components of ρ s . It can also be interpreted as a special instance of the eigenvalue problem for the Liouville superoperator L, namely as the instance of eigenvalue λ µ = 0. (We will always take µ = 0 to denote this case in the following.) The steady state ρ s is thus the eigenstate u 0 = ρ s of L associated with the eigenvalue λ 0 = 0. While it is tempting to think of Eq. (7) as a superoperator-analogue of the stationary Schrödinger equation, it is important to note that the Liouville superoperator L is in general not hermitian, i.e. L † is not equal to L. Hence, its eigenvalues may be complex-valued and we have to distinguish right eigenstates u µ from left eigenstatesȗ µ , given byȗ † As long as we keep this in mind, however, it is useful to mimic bra-ket notation and allow ourselves the freedom to write operators and their adjoints also in the alternative form u µ ↔ |u µ ) andȗ † µ ↔ (ȗ µ |. We can then denote the Hilbert-Schmidt inner product between operators as (x|y) ≡ tr[x † y]. Linear algebra dictates that the right and left eigenstates of L are orthogonal and, by appropriate normalization, can be chosen orthonormal, Assuming that the eigenstates of L form a complete set [75], we can represent an arbitrary operator X as and an arbitrary superoperator A as Except for the matter of left vs. right eigenvectors, these expressions are familiar from the usual decomposition of states and operators in Hilbert space, and we will make use of them in Section III.

B. Non-interacting open lattices, third-quantization method
Due to competition of on-site interaction and off-site hopping, as well as non-conservation of excitation number, the open Jaynes-Cummings lattice is not exactly solvable in general. The model does become tractable once the coupling to qubits is removed, by which one obtains an open lattice of non-interacting photons subject to coherent driving and photon loss (and a disconnected set of qubits). Third-quantization methods can then be used to find the exact solutions to the stationary Lindblad equation (7) [63][64][65]. To sketch the procedure, we consider a uniform photon-hopping Hamiltonian without disorder, H P = r δωa † r a r + (a † r + a r ) + κ r,r (a † r a r +h.c.) and the Liouville superoperator We bring the hopping Hamiltonian into diagonal form in two steps. First, we introduce the appropriate Bloch states which are labeled by quasimomentum k, and describe the normal modes of the undriven photon lattice H P | =0 . Second, we note that uniform driving is synonymous with driving of the k = 0 mode. Hence, we can eliminate the drive term by a coherent displacement of this mode, a k=0 → a k=0 + α, where α depends on both drive strength and relaxation rate. These steps yield the simplified Hamiltonian H P = k δω k a † k a k and Liouville superoperator This Liouville superoperator can now be diagonalized analytically by using the third-quantization method [63][64][65] (or alternative techniques, e.g., [76]). We employ the superop- While b † † r and b r are not proper adjoints, the use of the unconventional " † †" symbol is motivated by the the fact that it leads to commutation relations of the ordinary form, and all other commutators vanish. Thanks to this commutator algebra, L P takes on the compact form [64], Here, b k (b † † k ) and β k (β † † k ) may be viewed as "normalmode" superoperators, and complex-valued "mode frequencies" t k = −i δω k − γ/2 and b * k . From this, it is straightforward to read off eigenvalues and eigenstates of L P , analogous to the way we find eigenvalues and eigenstates of a non-interacting boson Hamiltonian. For a given k-mode, the right and left "vacuum states" obey b k r k 00 = β k r k 00 = 0 andȓ k 00 b † † k =ȓ k 00 β † † k = 0. The right "vacuum state" is therefore the projector r k 00 = |0 k 0 k | onto the pure state without any photons in mode k. One can show that the left "vacuum states" must always coincide with the identity operator,ȓ k 00 = 1. Once the coherent displacement is reversed, we thus obtain ρ s = k =0 |0 k 0 k | ⊗ |α k=0 α k=0 | for the steady state of the original Liouville superoperator L P .
The "excited" eigenstates of L P are obtained by acting with the creation superoperators on the "vacuum states." For given k, this means The corresponding eigenvalues are λ k mn = mt k + nt * k . When forming the appropriate product states and summing eigenvalues over k, we thus obtain the entire spectrum of the Liouville superoperator.
While our brief review here has been limited to the simplest scenario, we emphasize that the third-quantization method applies to a much larger class of lattice models which may involve bosons, fermions as well as spins [77][78][79][80].

III. LINDBLAD PERTURBATION THEORY AND RESUMMATION
In this section, we present the formalism of Lindblad perturbation theory and its resummation. This section remains general and is applicable to Markovian open quantum systems of various types. The concrete application of the formalism to an open Jaynes-Cummings lattice follows in Section IV. Together, these two sections constitute the central result of our paper.
Consider the general case of an open quantum system with Hamiltonian H and Liouville superoperator L. We shall assume that L is amenable to a perturbative treatment and can be decomposed into a sum L = L 0 + L 1 , consisting of the unperturbed Liouville superoperator L 0 , and the perturbation L 1 . For L 0 to qualify as such, it is expected that we can obtain its spectrum exactly. We denote the resulting unperturbed eigenvalues by λ 0 µ , and the corresponding unperturbed right and left eigenstates by |u 0 µ ) and (ȗ 0 µ |, respectively. The spectra of L and L 0 differ, but corrections may be calculated by a perturbative series expansion if L 1 is "sufficiently small." The corrections to eigenvalues and eigenstates can then be determined recursively, order by order [56]. Our interest here primarily regards the steady state ρ s , and we apply Lindblad perturbation theory assuming the nondegenerate case in which the steady state is unique. Two remarks may be useful for clarification. First, we emphasize that non-degeneracy refers to the spectrum of L 0 , not to the Hamiltonian H; we make no assumptions about the spectrum of H. Second, we note that non-uniqueness of the steady state and resulting non-analyticities are crucial at the phase boundary of a dissipative phase transition. Perturbative series expansions will generally not hold directly at such a boundary, but may still be applicable in its vicinity. Turning now to the concrete series expansion |ρ s ) = j |ρ j ) of the steady state, we note that the j-th order contribution ρ j is obtained from the recursion relation Here, inversion of L 0 will always be understood as restricted to the space orthogonal [81] to the unperturbed steady state, i.e., L −1 With this, we obtain the formal expression where ρ 0 is the unperturbed steady state of L 0 , and we have introduced the shorthand U ≡ −L −1 0 L 1 . The matrix elements of the U-superoperator are for µ > 0 and (U) 0ν = 0. Using this shorthand, we write the j-th order contribution to the steady state in the form of the chain expression and represent it diagrammatically as shown in Fig. 3(a). As a diagrammatic rule, we choose dots to represent unperturbed eigenstates, |u 0 µj ), and interconnecting lines to represent factors of U. Reading from the left to right, the leftmost state is the unperturbed steady state |ρ 0 ) = |u 0 0 ), and the rightmost one the final state |ρ µj ) which appears explicitly in the expression (21). All intermediate states and the final state, µ 1 through µ j are subject to summation, but cannot coincide with the initial unperturbed steady state due to (U) 0ν = 0.
To facilitate our resummation scheme for the steady-state series we decompose the U-superoperator into two parts To make the definition parallel to expressions to come, we have explicitly recorded the exponent j = 1 on the left-hand side. We now specify the terms on the right-hand side in such a way that resummation takes a particular simple form. We define the first-order Σ-superoperator as the diagonal part of U 1 , i.e. (Σ 1 ) µν = δ µν (U 1 ) µν . Accordingly, Σ 1 and L 0 share the same set of right eigenstates, i.e., Σ 1 |u 0 µ ) = Σ 1;µ |u 0 µ ), and the eigenvalue is Σ 1;µ = (U 1 ) µµ . We will see that this simplicity of Σ 1 will be important for the resummation of the series, Eq. (22). The term T 1 in Eq. (23) is the off-diagonal remainder of the U-superoperator.
A simplification of the previous expressions occurs when making the natural assumption that the perturbation L 1 is itself off-diagonal with respect to the unperturbed eigenstates of L 0 . (Whenever L 1 does not satisfy this assumption, a simple re-definition of L 0 and L 1 can be used to turn L 1 off-diagonal.) Now, if L 1 is off-diagonal, so is U = −L −1 0 L 1 and we immediately obtain This may initially make the decomposition of U seem pointless, but we will see momentarily that this simplification does not carry over to higher orders j > 1, thus justifying our approach. We next consider the second-order term U 2 |ρ 0 ), which warrants a decomposition of U 2 = T 1 U into Analogous to our strategy above, we define the second-order Σ-superoperator, as the diagonal part of the left-hand side, Note that the off-diagonal character of U automatically leads to exclusion of the term with τ = ν. As before, T 2 represents the remaining off-diagonal part in Eq. (25). We represent Σ 2 by the loop diagram shown in Fig. 3(b). Since Σ 2 is diagonal, its initial and final state |u 0 µ ) must be identical. However, the intermediate state |u 0 τ ) involved in the expression (26) must differ from |u 0 µ ). For resummation of terms to infinite order, we need to formulate our decomposition strategy for arbitrary order j. It is natural to extend the definitions for diagonal and off-diagonal superoperators, Eqs. (25) and (26), by setting The recurrence relation is solved by where A denotes the off-diagonal part of A with respect to the unperturbed basis {|u 0 µ )}. We must note, however, that the definition (27) does not yet determine a unique separation scheme beyond second order. Consider, for instance, the case of the third-order term involving U 3 . While we know the decomposition of U 2 = T 1 U from Eq. (25), we still have the freedom to perform the substitution for either U 3 = U(U 2 ) or U 3 = (U 2 )U. Both forms are mathematically equivalent, but only the systematic usage of one replacement rule produces expressions for which resummation becomes simple. We will consistently employ the form where U j−1 signals that U j−1 is to be replaced by an expression composed of Σ, T, and U-superoperators. Multiple replacements, in some cases making use of the identity U = T 1 , may be necessary to reach the final decomposed form only involving Σ and T-superoperators.
For illustration, we consider the decompositions of U 3 , U 4 , and U 5 . For the third-order case, we first make use of Eq. (29) and then Eq. (27) to obtain The last term on the right-hand side cannot be simplified further (except for substituting U = T 1 ), the first term is further decomposed by using Eq. (27), leading to the final expression For the fourth order, we merely sketch the decomposition, We give the fifth-order result without showing substeps, Inspection of Eqs. (31)-(33) indicates a systematic structure underlying the expressions, namely Each term in this sum has one factor of T k of order 0 ≤ k ≤ j and a prefactor σ j−k consisting of all possible combinations of Σ-superoperators of total order j −k. A formal proof of this is given in Appendix A. (Recall from Eq. (24) that Σ 1 = 0, which reduces the number of terms significantly.) Using the decomposition (34) and regrouping terms according to each occurrence of T j , we can now rewrite the perturbation series for the steady state in the form Here, the superoperator f = f(Σ 1 , Σ 2 , · · · ) = ∞ n=0 σ n implements the resummation of terms that we have been aiming for. It is given by i.e., the sum of all possible products of Σ-superoperators (here explicitly shown up to fifth order). Diagrammatically, we represent f in the form shown in Fig. 3 The role of Σ j resembles that of an irreducible self-energy contribution of order j in closed-system perturbation theory. Indeed, if we define Σ = ∞ j=1 Σ j as the sum of all irreducible "self-energy" contributions, we can rewrite f = ∞ j=0 Σ j = (1 1 − Σ) −1 , and obtain for our resummed series expansion of the steady state. Due to the (1 1 − Σ) −1 prefactor, each term |ρ [j] ) = (1 1 − Σ) −1 T j |ρ 0 ) in this sum includes perturbative corrections up to infinite order. We therefore call the term |ρ [j] ) the rank-j term of the resummed series. We note that, formally, Eq. (37) is an exact expression for the steady state. Practical calculations will typically involve a truncation in both the maximum summation index j as well as the maximum order of irreducible selfenergy contributions taken into account. We represent individual terms |ρ [j] ) in the resummed series by the type of diagram shown in Fig. 3(c). The final-state dot is marked with a circle to indicate the inclusion of the self-energy correction. The diagrammatic rules are similar to the case without resummation, Fig. 3(a), except that the offdiagonal nature of T j in addition requires that all intermediate states be different from the final state. This is simple to infer when writing T j in the form of Eq. (28). Each diagram then translates to an expression with the following structure: × (U) µj νj−1 (U) νj−1νj−2 · · · (U) ν2ν1 (U) ν10 .

IV. APPLICATION TO THE OPEN JAYNES-CUMMINGS LATTICE
Lindblad perturbation theory and resummation as discussed in the previous section are applicable to a large class of open quantum systems. Here, we present its concrete use in studying the open Jaynes-Cummings lattice [Eqs. (1) and (3)-(5)] as a specific example of an open driven quantum lattice. The example is of particular interest due to its role as a minimal model for highly anticipated experiments with circuit-QED lattices [62]. The photonic backbone of the lattice has already been demonstrated in the experimental work by Underwood et al. [49].

A. Preparatory steps
We shall consider a uniform lattice, in which resonator frequencies, qubit frequencies, and related quantities have uniform values across the lattice. (Disorder levels, especially in qubit frequencies, may need to be considered carefully for detailed modeling of future circuit-QED experiments, but this consideration is beyond the scope of the present paper.) We find the modes of the photonic lattice structure by diagonalizing the N × N hopping matrix, N being the number of lattice sites [82]. For periodic lattices, diagonalization is achieved in the usual way by switching from real space to momentum space via the transformation a r = kã k e ik·r / √ N . Here, photons inside the mode with quasi-momentum k are annihilated byã k , and k runs over all reciprocal lattice vectors from the first Brillouin zone. Note that k = 0 corresponds to the uniform mode with identical amplitudes on all sites, which is the mode being excited by a global coherent drive.
Depending on the values of model parameters, it is beneficial to perform a displacement transformation which eliminates the coherent drive on the photon mode and converts it into an effective qubit drive instead. This is particularly helpful when the uniform photon mode is approximately in a coherent state with a large number of photons. The coherent displacement then serves as a tool incorporating this insight directly into the unperturbed Liouville superoperator and mitigates the need for large photon number cutoffs. (In a regime of low photon occupation, however, the displacement transformation can be skipped and the perturbative treatment carried out directly.) The displacement is applied to the k = 0 mode, i.e., b 0 =ã 0 − α, with α = − √ N /(δω k=0 − iγ/2). After this displacement, the resonator drive is converted into an effective qubit drive of strength q = gα/ √ N . The resulting Hamiltonian has the form where the three terms correspond to the photon, qubit, and photon-qubit coupling contributions. The resonator part now lacks the drive term, H k r = δω k b † k b k . (We define b k =ã k for k = 0 to unify notation.) The qubit Hamiltonian including the effective drive reads and the interaction Hamiltonian is given by Finally, the dissipator term simply transforms as In the absence of the interaction H kr rq , resonator modes and qubits decouple, and the associated master equation is exactly solvable. This presents us with an ideal starting point for a perturbative treatment of H kr rq , which physically is a very sensible treatment of the dispersive regime. The unperturbed Liouville superoperator is then L 0 = k L k r + r L r q with separate photon contribution, First, the photonic part L k r is directly amenable to the diagonalization procedure we discussed in Section II A. The resulting eigenvalues read λ k mn = mt k + nt * k with t k = −i δω k − γ/2 (m, n = 0, 1, . . .). The associated right and left eigenvectors |r k mn ) and (ȓ k mn | are those given in Eqs. (16) and (17). Second, the qubit Liouvillian r L r q can be diagonalized exactly since it decomposes into a direct product of 4 × 4 matrices. For each L r q , we denote the eigenvalues, and right/left eigenstates by r µ , |q r µ ) and (q r µ |, respectively (µ = 0, . . . , 3). Except for special parameters, analytical expressions for these quantities are too lengthy to provide much insight, and will hence not be recorded here.
Altogether, right and left eigenstates of the full-lattice Liouvillian L 0 thus take the form with corresponding eigenvalues The multi-indices m, n, and µ collect the sets of all "quantum numbers" m k , n k and µ r . With this, we are ready for the perturbative treatment of the Liouvillian L 1 capturing the Jaynes-Cummings interaction H kr rq .

B. Perturbative treatment and resummation
The perturbation superoperator decomposes into a sum L 1 = k,r L kr 1 , in which each term describes the interaction between an individual resonator mode k and a qubit at position r: It is therefore convenient to write the U-superoperator (see Section III) as an analogous sum, i.e., U = k,r U kr with U kr = −L −1 0 L kr 1 . Each U kr is off-diagonal with respect to the unperturbed basis, so that U kr = U kr = T kr 1 holds. The rank-1 term in the resummation [Eq. (37)] is given by Here, |ρ [1] kr s ) = |r k mn q r µ ) k =k |r k 00 ) r =r |q r 0 ) is an interacting cluster involving the photon mode k and qubit at position r in state s = (m, n, µ). In a similar manner, we see that the rank-2 term incorporates the resummation of a cluster |ρ [2] kr s,k r s ) composed of two photon modes and two qubits in states s and s . (The bracket notation in the corresponding summation signals that the allowed choices of these states is dictated by the off-diagonalism requirements in T 2 = U kr U k r .) In the definition of |ρ [2] kr s,k r s ), several cases must be distinguished according to whether k = k and/or r = r . For the case where both pairs are distinct, we have |ρ [2] kr s,k r s ) = |r k mn r k m n q r µ q r µ ) Analogous definitions hold in the other three cases. By inspection of Eqs. (45) and (46) we expect that, in general, the rank-j correction consists of a sum over all possible terms in which clusters of j photon modes and j qubits deviate from the unperturbed density matrix. Thanks to resummation, interaction within each cluster includes terms up to infinite order. We note that this cluster structure directly implies a hierarchy of correlations with increasing rank j. Specifically, every n-point correlation function with p photon and q qubit operators, b r1 · · · σ aq rq ss , does not trivially separate into a product of correlators if the rank j of the correction satisfies j ≥ max{p, q}. We also emphasize that clusters automatically include long-range correlations between distant qubits.
The essential tasks of determining the perturbative corrections and resummation consist of evaluating matrix elements of the form given in Eqs. (45) and (46), and computing the effect of the resummation superoperator Σ to a given order. We illustrate the procedure for the example of rank-1 corrections. Plugging in the definition of U kr and recalling that L 0 is diagonal with respect to the unperturbed basis states, we obtain ρ [1] kr s U kr r k 00 q r 0 = − ȓ k mnq r µ L kr 1 r k 00 q r 0 /λ kr s (47) where λ kr s = λ k mn + r µ . Once the commutator is opened, it is useful to note that the simple properties of the L reigenstates lead to vanishing overlaps tr(ȓ k † mn b k ) = 0, so that any application of U kr must switch to a different resonator-mode eigenstate. The same does not hold for traces of the qubit degrees of freedom, i.e., the overlaps tr(q k † µ σ + r q k µ ) etc. may indeed be non-zero. As a result, we obtain the two types of terms for the rank-1 correction which are diagrammatically depicted in Fig. 4(a). The evaluation of rank-2 corrections follows the same basic scheme. Unsurprisingly, it is more tedious and we only show two examples of corresponding diagrams in Fig. 4(b).
The effect of the resummation superoperator is to redistribute weights among cluster contributions. Since Σ is diagonal with respect to unperturbed Liouvillian eigenstates, = + (k, r) (k) + ... we can cast its contribution to a particular into the form Σ|u 0 s ) = |u 0 s )(ȗ 0 s |Σ|u 0 s ) and evaluate the occurring matrix element as follows. We choose an appropriate truncation for the series Σ = Σ 2 + Σ 3 + · · · of irreducible resummation operators (Σ j ) s s = δ s s (T j−1 U) s s . Matrix elements for Σ j are calculated in the same way as for |ρ [j] ) except that the final state of the chain must be identical to the initial state. Fig. 5 shows the resulting two diagrams for Σ 2 .
If we are merely interested in certain steady-state expectation values (rather than in the density matrix itself), then the calculation of perturbative results may be simplified. As an example, consider computing an expectation value of a local qubit operator σ a r up to corrections of rank j, σ a r ≈ j j =0 tr(σ a r ρ [j ] ). To effect the desired simplification, we recall that all eigenstates of L k r and L r q other than the steady state must be traceless, i.e., tr(q r µ ) = 0 for µ = 0 and tr(r k mn ) = 0 for non-zero m or n [83]. Therefore, any perturbative contribution ∼ tr(σ a r u 0 m n µ ) in which µ r = 0 for some r = r will immediately vanish since the partial trace over the qubit at position r is zero. Similarly, any term with (m k , n k ) = (0, 0) for some photon mode k will vanish. As a result, none of the rank-1 corrections [ Fig. 4(a)] contribute to local qubit expectation values. Only those diagrams that terminate in a state labeled (r) will yield a nonzero contribution to σ a r . Analogous diagrammatic rules apply for photon-mode operators.

C. Perturbative results for a near-resonant regime
We now illustrate the validity and improvement achieved by (partial) resummation of the Lindblad perturbation series. For this purpose, we compare perturbative results for finite-size Jaynes-Cummings chains [ Fig. 2(a)] to exact results computed via quantum trajectories methods. Perturbative resummation further allows us to treat periodic chains (or, open chains if desired) of sizes beyond the computational capabilities of exact quantum-trajectory solutions. Finally, we can carry out perturbative resummation even for an infinite system with chain or global-coupling geometry. This versatility enables us to predict finite-size effects, the approach to the thermodynamic limit, as well as differences according to distinct lattice geometries. (We discuss the quite moderate computational costs of the perturbative treatment in Appendix B.) In our treatment, we capture photon mediated qubit-qubit interactions by second-order Lindblad perturbation theory with resummation based on single-loop terms Σ 2 (i.e., corrections of rank 2 in the above terminology). The most natural regime for treating the Jaynes-Cummings coupling perturbatively in this way is the dispersive regime where the detuning ∆ = min k |Ω − ω k | between qubit and photon-mode frequencies is large compared to their mutual coupling strength g [84]. We have confirmed by exact numerics that the perturbation theory is indeed reliable in this regime and, over a wide range of parameters, we identify g 2 /Γ∆ as the relevant small parameter governing the series expansion.
In the following, we choose to present results from exploring a different parameter regime more directly based on the open-system nature of Jaynes-Cummings lattices. Specifically, we consider the case where photon hopping dominates over both photon decay and Jaynes-Cummings coupling, and where the latter two are permitted to be of the same order, i.e. κ g ∼ γ. The strong hopping κ shifts spectral weight of the photon modes away from the bare resonator frequency which is chosen degenerate with the qubit frequency, Ω = ω. This regime is not fully dispersive, so that nonlinearities are more pronounced, and the significance of resummation becomes easily visible.
We begin with the comparison between perturbative and exact results for the steady state of few-site Jaynes-Cummings chains with periodic boundary conditions. In our calculations, we have considered several qubit and resonator expectation values. Among those, we find that | σ − | = σ x 2 + σ y 2 /2 is a convenient choice for clearly resolving resonances. Representing the reduced steady-state density matrix for one of the qubits by means of the Bloch sphere picture, this quantity is directly proportional to the distance of the Bloch vector from the z axis, see Fig. 6(a). Computing exact steady-state solutions for Jaynes-Cummings chains even as small as 4 sites, is a non-trivial task which we accomplish by averaging of quantum trajectories. For instance, exact results presented in Fig. 6 for N = 4 were determined from stochastic time evolution of a quantum state vector of size 10 4 . Sufficient averaging of a single data point required a runtime of several days on one core.
The comparison between exact and rank-2 perturbative data (including resummation at the level of Σ 2 ) in Fig. 6(b) shows very good agreement, and indicates that the resummation procedure closely matches the exact solution. Plots in this figure show the steady-state value of | σ − | as a function of the detuning δΩ = Ω − ω d between the drive and the bare qubit frequency. Multiple resonances are visible over the chosen frequency range (the nature of which we will further discuss below). The only notable quantitative deviations occur for N = 4 in the vicinity of the bare qubit frequency where δΩ = 0. This deviation has a simple explanation: a look ahead at Fig. 7(b) shows that the 4-site chain has a photon mode with large spectral weight directly on resonance with the bare qubit frequency, so that there we must expect the perturbative treatment in g to break down. With the exception of this finding, we conclude that resummation of the perturbative series in the chosen parameter regime works very well. The improvement gained over the pure second-order approximation is illustrated in Fig. 6(c). In this panel, curves show the difference between approximate and exact results, for the case including resummation (solid lines) and lacking resummation (dashed lines). Excluding the pathological region for N = 4 around δΩ = 0, we observe that resummation consistently improves the results, reducing the deviation from the exact solution. The improvement is especially significant in the resonance region between δΩ/Γ ≈ −17 and −23. Here, the drive populates the uniform mode (centered at δΩ/Γ = −20) and renders photonmediated qubit-qubit interaction important, making resummation of corrections up to infinite order particularly fruitful.
We next turn to the discussion of the resonances visible in Fig. 6 and shown for additional site numbers N in Fig. 7(a). For small chain lengths up to N = 5 sites, we observe three resonances labeled A, B, and C in the range of drive frequencies spanning the photon-chain eigenmodes and qubit frequencies.
[The uniform photon mode has the lowest frequency, δΩ/Γ = −20, and the bare qubit frequency is at δΩΓ = 0, see vertical lines in Fig. 7(b)]. We find that the detailed positions and strengths of resonances depends on the number of sites, revealing systematic finite-size effects for chains of short lengths. Both the nature and N dependence of resonances can be explained, or at least motivated, by the following considerations.
The resonance marked A directly coincides with the frequency of the uniform photon mode. Equivalent interpretations of the resonance can be given based on the original Hamiltonian Eq. (3) with a coherent tone driving this particular mode with strength , or for the Hamiltonian following the displacement transformation Eqs. (39) and (40). Employing the language of the latter description, we note that the strength q = −g /(δω k=0 − iγ/2) of the effective qubit drive reaches its local maximum at the uniform-mode frequency (δΩ/Γ = −20). This peak in the off-resonant Rabi drive, modified by weak coupling to photon modes, is responsible for resonance A. Its dependence on the site number N is relatively weak and mainly affects the shoulders of the resonance. This is further confirmed by our results for the infinitesystem case with periodic-chain and global-coupling geometry [ Fig. 6(c)]. Both show the same resonance A, but differ in the resonance shoulders.
For resonance B, the N dependence of resonance position and strength is much more pronounced. The general region where B occurs is close to δΩ = 0, i.e., where bare qubit frequency and drive frequency match. Upon displacement, the drive transforms into an effective Rabi drive on each qubit [Eq. (40)]. Hence, the presence of a resonance is natural, and variations in its strength and precise position must be governed by the Jaynes-Cummings interaction playing the role of the perturbation in our treatment [Eq. (44)]. The importance of this interaction is influenced by the position of photon-mode resonances ω k relative to the bare qubit frequency. In Fig. 7(b) and (d), we depict width and position of photon modes in terms of the spectral function s(ω) = k γ/2 (ω−ω k ) 2 +(γ/2) 2 , which is the sum over individual Lorentzians for each photon mode, normalized such that dω s(ω) = 1. Inspection of resonance-B positions [ Fig. 7(a),(c)] and peaks in the spectral function [ Fig. 7(b),(d)] shows that peaks in s(ω) with significant weight in the region Ω ± g, shift B resonances towards the close-by photon mode (such as for N = 3, 5). Further, it is clear that strongly increased weight of the spectral function directly at the qubit frequency (such as for N = 4 and for the global-coupling geometry) endangers the validity of perturbation theory in the Jaynes-Cummings coupling. Above, we recognized this as the reason for the observed deviations between perturbation theory and exact solution close to δΩ = 0 in the N = 4 case. A look at the spectral function for the global-coupling geometry shows that the same issue occurs here. Accordingly, we show the perturbative result in Fig. 7(c) only with dashes in (c,d) show analogous plots for infinite lattices with periodic chain or global-coupling geometry. Resonances A and B are visible, but resonance C is (nearly) absent.
[All system parameters are identical to those in Fig. 6.] that region. We note that steady-state expectation values for infinite lattices are not always easily accessible by other methods. Thanks to the possibility of carrying out leading-order resummation analytically in the infinite-system case, our treatment gives direct access to the thermodynamic limit of different lattice geometries. Here, we have chosen two extreme cases: the infinite periodic chain with a minimum number of links between sites, and the global-coupling model with a maximum number of links. Fig. 7(c) depicts results for both lattice structures. We expect that the region close to δΩ = 0 is unproblematic for the infinite chain case, but potentially pathological for the global-coupling model which accumulates maximum spectral weight at the bare qubit frequency [ Fig. 7(d)]. Away from the δΩ = 0 range, the two geometries yield similar behavior of | σ − | versus drive detuning δΩ. As before, characteristic differences occur primarily in the shoulders of resonance A. Interestingly, the resonance C is absent for both infinite lattices, and we discuss the rather anomalous behavior of this finite-size resonance next.
The anomalous properties of resonance C are illustrated in Fig. 8. The position of this resonance close to δΩ/Γ ≈ 8 does not simply coincide with a resonance between photon modes and bare qubits. Panel (a) shows the decrease of the resonance strength for increasing chain length N . Moreover, both resonance position and strength depend sensitively on the drive power ∼ as shown in Fig. 8 for the dimer case, N = 2. We investigate this anomaly for the N = 2 case, where a semi-quantitative reduced model can shed light on the origin and nature of this resonance.
For N = 2 we can confirm analytically that the anomalous resonance C is closely related to an eigenstates |ψ of the displaced Hamiltonian H [Eq. (39)] without the effective drive. This eigenstate comprises two excitations distributed between the uniform mode and the qubits. Truncation to first order in g/κ (a small quantity for the chosen parameter set) yields with eigenenergy 2Ω + 2κ − g 2 /2κ. Here, |n 0 is the nphoton state of the uniform mode and |↑↑ etc. denote states of the qubits on the two sites. The effective drive Hamiltonian with strength q connects the ground state |g = |0 0 |↓↓ to the state |ψ via two intermediate states |r and |q , see Fig. 9. These intermediate states belong to the one-excitation manifold and primarily consist either of a photon in the uniform mode or or of a qubit excitation, respectively. Truncated again to first order in g/κ, |r and |q are given by Our description of the anomalous resonance in the following is based on the effective four-level model spanned by the states |g , |r , |q , and |ψ , see Fig. 9. Within this model, the effective drive Hamiltonian connects the ground state |g to the two-excitation state |ψ via |r and |q . Since the effective drive creates (or annihilates) qubit excitations only, there is stronger hybridization of |g with |q (strength ∼ q ) than with |r (strength ∼ q g/κ). An analogous argument applies to explain hybridization of |r with |ψ (again, strength ∼ q )). We thus obtain the picture of two pairs of hybridized states, |g ↔ |q and |r ↔ |ψ , with only small drive matrix elements connecting the pairs.
Due to the energy differences between states within each pair, hybridization is only partial. The two emerging hybridized states relevant to the anomalous resonance have significant overlap with the ground state |g in one case, and the two-excitation state |ψ in the other case. The resonance C can be approximately viewed as a resonance between these two hybridized states. Note that the degree of hybridization critically depends on the effective drive strength q , which is in turn proportional to the drive strength . As a consequence, the energy separation between the two relevant hybridized states depends on drive strength as observed in Fig. 8(b).
The generalization of the effective four-level model to periodic chains with larger number of sites N is difficult due to the proliferation of degeneracies among eigenstates of H in the absence of a drive. Based on our perturbative calculations, we find a clear trend of diminishing resonance strength with increasing number of sites [ Fig. 8(a)]. The anomalous resonance C hence provides an interesting example of an interaction-induced feature which is limited to a small N , which should be accessible in experiments with Jaynes-Cummings chains of only a few sites.

V. CONCLUSION
We have extended the general Lindblad perturbation framework to include resummation of an infinite subset of perturbative corrections. We have formulated the scheme at a general level, emphasizing that is not limited to a particular open quantum system, but benefits a variety of problems in Markovian quantum systems amenable to perturbative treatment. For the examples we have investigated, we find that the series resummation can significantly improve the accuracy of the perturbative treatment.
We have applied perturbation theory with resummation to a specific model of an open quantum lattice, the open Jaynes-Cummings lattice, and have introduced a diagrammatic representation systematically organizing the contributing terms. For small lattices, we find very good agreement with exact results which we obtained by extensive quantum trajectory simulations for an interesting parameter regime near resonance. Our perturbative treatment is capable of predicting steadystate observables for both finite and infinite Jaynes-Cummings lattices with different lattice geometries and dimensionalities -thus including settings that may not be easily accessible by other methods.
The capability of obtaining reliable results beyond exactly solvable limits of open quantum lattices is particularly promising as a method for validating experimental implementations of quantum simulators. Concrete realizations of open quantum lattices are currently being investigated in the circuit QED architecture [62].
Finally, questions that warrant future investigations regard the relation between mean-field approximations and suitable choices of resummation within the presented formalism. While Keldysh techniques tend to be more difficult for handling (pseudo-)spin degrees of freedom, it will be interesting to compare and relate results obtainable for simpler systems such as the open Bose-Hubbard lattice. Further investigations into spectra of Liouvillians, handling cases of degeneracies of the Liouvillian spectrum, and extending resummation to timedependent perturbation theory offer exciting perspectives for the study of dissipative phase transitions in open quantum systems.