Dispersive regime of the Jaynes-Cummings and Rabi lattice

Photon-based strongly-correlated lattice models like the Jaynes-Cummings and Rabi lattices differ from their more conventional relatives like the Bose-Hubbard model by the presence of an additional tunable parameter: the frequency detuning between the pseudo-spin degree of freedom and the harmonic mode frequency on each site. Whenever this detuning is large compared to relevant coupling strengths, the system is said to be in the dispersive regime. The physics of this regime is well-understood at the level of a single Jaynes-Cummings or Rabi site. Here, we extend the theoretical description of the dispersive regime to lattices with many sites, for both strong and ultra-strong coupling. We discuss the nature and spatial range of the resulting qubit-qubit and photon-photon coupling, demonstrate the emergence of photon- pairing and squeezing, and illustrate our results by exact diagonalization of the Rabi dimer.


Introduction
Bosonic and fermionic lattice systems, realized experimentally with ultracold atoms or trapped ions, have long served as a paradigm for quantum simulators of strongly-correlated many-body systems [1,2,3]. Recently, Jaynes-Cummings and Rabi lattices have been suggested as an exciting variant of this idea [4,5,6]. The basic physics of these lattices may enable new types of photon-based quantum simulators suitable for exploring many-body physics in and out of equilibrium, and has stirred lively interest [7,8,9,10].
In the Jaynes-Cummings and Rabi lattice, each lattice site consists of a photon mode interacting locally with a two-level system, and is described by the ordinary Jaynes-Cummings [11] or Rabi model [12]. In addition, photons are allowed to hop between nearestneighbor lattice sites. A key difference between the traditional boson and fermion lattices and the less conventional Jaynes-Cummings and Rabi lattice is that the latter provide an interesting additional tunable parameter: the detuning ∆. Its origin lies in the two-component nature of the Jaynes-Cummings and Rabi model, comprising an electromagnetic field component and a matter component. Each component is associated with its own characteristic energy scale. The possible energy mismatch between the two is quantified by the detuning ∆.
By changing the detuning, three qualitatively different regimes of the Jaynes-Cummings lattice and Rabi lattice can be accessed. In the resonant regime, the detuning is small and photons and matter excitations readily hybridize to form polaritons. In the dispersive regime, the magnitude |∆| of the detuning between the two-level splitting and the photon frequency ω is large. According to the sign of ∆ = − ω, we distinguish between the dispersive regime with negative detuning, where low-energy physics predominantly involves matter excitations, and the dispersive regime with positive detuning, in which photons govern the behavior at low energies. In this paper, we focus primarily on the two dispersive regimes where all interaction strengths are small compared to |∆|.
The circuit QED architecture [13,14,15] constitutes one of the most promising experimental platforms for the realization of dispersive Jaynes-Cummings and Rabi lattice systems [9,10]. In circuit QED lattices, photons can hop between transmission line resonators (which are coupled capacitively) and locally interact with a superconducting qubit. The qubit energy (and hence the detuning parameter ∆) can be tuned in-situ by an externally applied magnetic flux. Recently, the coherent photon exchange via hopping between several coupled resonators has also been realized experimentally [16].
Theoretically, the dispersive regime for a single Jaynes-Cummings [13] or Rabi system [17] is well understood. Its usual description is based on employing a perturbative Schrieffer-Wolff transformation, which eliminates the Jaynes-Cummings or Rabi coupling by switching to an appropriate dressed-state basis.
Here, we extend this procedure to an entire lattice of sites and systematically discuss all contributions in second-order perturbation theory (section 2). We find that effective qubitqubit interactions ‡, qubit-state dependent photon hopping terms and photon pairing terms (squeezing terms) emerge. All inter-site interactions are short-range but not limited to nearestneighbor sites. We explore the implications of the derived effective Hamiltonian for the lowenergy physics of the Jaynes-Cummings and Rabi lattices in section 3. For negative detuning, we show that the system reduces to an effective spin model with XY-type interaction for the Jaynes-Cummings lattice, and to a transverse Ising model for the Rabi lattice. For positive detuning, we discuss in detail the one-mode and two-mode vacuum squeezing relevant in the ultra-strong coupling regime as described by the dispersive Rabi model. Due to the ultrastrong coupling, approached indeed in a recent experiment for a single site in circuit-QED architecture [18], the non-trivial nature of the ground state makes the Rabi lattice particularly interesting. In section 4, we study the application of the dispersive regime to the Rabi dimer and confirm the validity of the derived effective Hamiltonians by comparison with results from exact diagonalization. Finally, we summarize and give an outlook on questions of future interest in section 5.

Derivation of the effective Hamiltonian
The Rabi lattice is described by the model Hamiltonian Here, each Rabi site is labeled by an index j, and consists of a harmonic mode coupled to a pseudo spin-1/2. As usual, excitations of the harmonic mode and spin on site j are created by a † j and σ + j , and annihilated by their Hermitean-conjugate counterparts. The energy necessary for an excitation of either type is fixed by ω and , respectively. The coupling strength is set by the parameter g. From the beginning on, we include the counter-rotating terms characteristic of the Rabi model and the ultra-strong coupling regime, and only drop those terms in our discussion of the Jaynes-Cummings limit where the rotating-wave approximation (RWA) is applicable (see text below). For the entirety of this paper, we will focus on the disorder-less case, assuming that all lattice sites have identical parameters. In the following and for the sake of brevity, we simply refer to the pseudo-spin degree of freedom as "qubit".
In addition to the onsite Rabi Hamiltonians, the last term in equation (1) introduces coupling between resonators. We consider weak resonator coupling, |t| ω, , and treat it within the RWA. As written in equation (1), the coupling term then allows for photons to hop from one resonator to its nearest neighbors, as indicated by the angular brackets in the summation over resonator pairs. The sign of the hopping amplitude t depends on the specific system realization. In the circuit QED architecture, for example, it depends on whether fullwavelength modes (t > 0) or half-wavelength modes (t < 0) are used [19,10].
The Rabi lattice (1) enters the dispersive regime when the detuning |∆| = | − ω| between qubit and resonator frequency is large compared to the coupling strength g. While the sum frequency Σ = + ω ≥ |∆| obeys Σ |∆| whenever the RWA holds, a hallmark of the dispersive regime of the Rabi model beyond RWA is that detuning and frequency sum may be of the same order, Σ ∼ |∆| [17]. In this case, two sub-regimes are of particular interest: the dispersive Rabi regime with negative detuning, g ω, and the dispersive Rabi regime with positive detuning, ω g . Both sub-regimes will be investigated in section 3. In any of the mentioned cases, the idea underlying the dispersive regime remains the suppression of interconversion between qubit and photon excitations due to a large energy mismatch which is not overcome by the coupling g. Under these conditions, the Rabi interaction terms ∼ g can be treated perturbatively.
To obtain a simple effective Hamiltonian describing the dispersive regime, it is convenient to carry out the perturbation theory in the form of a unitary Schrieffer-Wolff transformation,H = e iS He −iS [20,21]. For the cases of a single Rabi site and multiple qubits strongly coupled to a single resonator, this procedure was previously discussed by Zueco et al. [17]. Here, we extend it to a lattice of coupled Rabi sites. Using the appropriate Hermitian generator S, the unitary transformation switches to a dressed-state basis in which the original Rabi interaction is eliminated. Being unitary, the transformation clearly preserves the spectrum of the original Hamiltonian.
In our case, the unperturbed Hamiltonian H 0 consists of the terms in equation (1) with the exception of the Rabi coupling term: This corresponds to a photonic tight-binding model and uncoupled qubit. To avoid distractions, we limit our discussion to one-dimensional Rabi lattices. The generalization to lattices with more complex structure is straightforward. Thus, considering a one-dimensional photon tight-binding model with Born-von Karmann periodic boundary conditions over a number N of lattice sites, we define the itinerant photon operators by Here, for a large number of sites (N → ∞), the quasi-momentum k spans the entire first Brillouin zone, −π < k ≤ π. Employing this basis, the unperturbed Hamiltonian can be rewritten in the diagonal form as Here, the dispersion of the itinerant photons for a one-dimensional lattice is simply ω k = ω + 2t cos k. The perturbation V is the onsite Rabi interaction on each lattice site and can be expressed in terms of the itinerant photon operators as For itinerant photons, the energy gap towards the qubit energy now depends on the quasimomentum k. We thus define the k−dependent detuning parameter ∆ k = − ω k and the frequency sum Σ k = + ω k . Ensuring that the qubits remain detuned from all itinerant photons, we recognize that constitutes a necessary condition for the dispersive regime of the Rabi lattice model. Figure 1 illustrates this condition for both positive and negative detuning.
In the perturbative approach, the Schrieffer-Wolff transformation is carried out to a certain order in the interaction strength g. In this paper, we present results for the Rabi lattice model to second order in g. The generator necessary to achieve this must be of first order in g [21], and is given by We then construct the effective Hamiltonian by expanding the exponentials in e iS1 He −iS1 and collecting terms up to second order. § Due to the form of the interaction V , the generator S 1 splits into 4 terms, resulting in 4 2 = 16 second-order terms which contribute to the effective Hamiltonian. (Note that the first-order term vanishes since [H, S 1 ] = 0.) Figure  2 visualizes all 16 contributions in the form of "ladder-type" diagrams , similar to those previously introduced in Ref. [22]. Each diagram can be directly translated into a contribution to the effective Hamiltonian according the following rules: § Recall that S 2 does not contribute to the second-order effective Hamiltonian since [H, S 2 ] = 0, see e.g., Ref. [21].
Note, that the term "ladder-type" solely refers to the appearance of these contributions as visualized in Fig. 2 and should not be confused with other types of Feynman "ladder" diagrams. . The diagrams XIII-XVI vanish when summed over k, k , and hence do not contribute to the effective Hamiltonian. Figure 3. Paths from Fig. 2 for the special case of j = j . As an example, the paths I and II are shown here. In this case, initial and final states differ, no interference occurs and photon operators are not cancelled. Such paths contribute to all those second-order terms in the effective Hamiltonian (11) that involve photon operators.

I II
(i) Each ladder step corresponds to one unperturbed eigenstate of H 0 , with well-defined number of photons and qubit excitations. Virtual intermediate states are marked by dashing.
(ii) Diagrams are read from right to left. Arrows show virtual transitions produced by a specific operator term from S 1 (label on the arrow), e.g., a † k σ − j . The contribution to the effective Hamiltonian contains the same operator combination (including the operator ordering) as shown by the arrow label. For each occurrence of photon operators, traveling-wave factors should be included according to a k → a k e ikj and a † k → a † k e −ikj . (iii) Each contributing diagram involves summation over all intermediate labels j, j , k, k with a 1/N prefactor.
(iv) The energy coefficient for each contribution is given by g 2 As an example, the analytical expressions obtained for the first two diagrams read: We now turn to the systematic discussion of the contributions obtained from the diagrams shown in figure 2. All contributions from the paths XIII-XVI vanish. For each of these paths, summation over k and k leads to complete cancellation due to the opposite signs of the energy denominators involved. Several paths in Fig. 2, including paths I and II, are shown as pairs. Each member in such a pair has the same initial and final states but undergoes its two virtual transitions in opposite order, thus leading to interference and partial cancellation. This cancellation originates from the opposite signs of the prefactors, see, e.g., the example in equations (8) and (9).
We first discuss the case of distinct site indices j = j . Since Pauli operators on different sites commute and itinerant photon operators obey the canonical commutation relation [a k , a † k ] = δ kk , we find that all terms inH I +H II with k = k and j = j exactly cancel. In the remaining (k = k ) term, all photon operators are eliminated: Equation (10) represents photon-mediated flip-flop (XY) interaction between qubits on different sites. Further contributions of this type are produced by paths III and IV. By contrast, the paths IX-XII produce an effective qubit interaction of the type σ + j σ + j + σ − j σ − j . Finally, the j = j contributions of the paths V-VIII exactly vanish.
Next, we consider the onsite case, i.e., the situation of j = j in each path. By comparing with the expressions in equations (8) and (9), we note that the initial and final states must differ in this special case. For the example of the paths I and II, the resulting proper diagrams are shown in Fig. 3. They involve a single qubit flip accompanied by the creation and annihilation of one itinerant photon each. The j = j contributions of the paths V-VIII produce photon pair creation and annihilation terms. Finally, the j = j contributions of the paths IX-XII identically vanish due to the occurrence of two Pauli raising or lowering operators on the same site, where (σ + j ) 2 = (σ − j ) 2 = 0. By transforming all itinerant photon operators back into real space, we obtain the effective second-order Hamiltonian for the dispersive Rabi regime: where a global constant has been dropped. The coupling constantsC ± m depend on the distance m between lattice sites and are defined as Since the condition (6) for the dispersive regime of the Rabi lattice also implies that the inequalities |2t/∆| < 1 and |2t/Σ| < 1 must hold, we can ascertain that g/∆, g/Σ, t/∆ and t/Σ are all small parameters. The effective Hamiltonian (11) obtained from the perturbative Schrieffer-Wolff transformation is a series expansion in both g/∆ and g/Σ. Due to exact diagonalization of the photon tight-binding model, the coupling constantsC ± m contain terms to all orders in t/∆ and t/Σ. We can further elucidate this fact by re-expressing the denominators in equation (12) in terms of convergent geometrical series, namelỹ Applying the binomial theorem to the last factor, using the relation N −1 k e ikl = ∞ z=−∞ δ l,zN , and keeping only the leading order terms in t/∆ and t/Σ, we find the approximationC for the coupling constants. Equation (14) is valid for 0 ≤ m ≤ N/2. Whenever m is outside this range, the exponents in equation (14) should be replaced by |m mod N |. Equations (11), (12) and (14) constitute the main results of this first part of our paper. In the following we will explore the implications of this effective Hamiltonian and discuss the physics of the dispersive regime of the Rabi lattice and the Jaynes-Cummings lattice in the subsequent section.
As expected, the dispersive approximation to the Rabi lattice Hamiltonian involves energy shifts and interaction terms which do not inter-convert between qubit and photon excitations. Based on equation (11), we now discuss shift and interaction terms one by one. The second term on the right-hand side of equation (11) captures the Lamb shift for the on-site qubit-resonator system. The shift is identical to the one obtained for a single qubitresonator system in the dispersive regime [13] when using the weak-coupling approximatioñ The third term in equation (11) produces photon-mediated qubit-qubit interactions of the transverse-Ising type. ¶ The contributions responsible for going beyond the bare flip-flop (XY) interaction are the additional counter-rotating terms σ + j σ + j and σ − j σ − j . The strength of this interaction is set by the coupling constantC − m which, according to equation (14), decays exponentially with increasing distance m between lattice sites.
The terms 4 -7 on the right-hand side of equation (11) all directly involve photons. Term 4 with the form a † j a j σ z j produces the well-known AC Stark shift on each lattice site [13]. Term 5 with the structure a † j a † j σ z j is an onsite term as well but goes beyond a mere energy shift: here, photon pairs are created or annihilated on a single site. At the same time, the amplitude sign for this process depends on the state of the local qubit. Counter-rotating terms like this one are characteristic of ultra-strong coupling and reflect, naturally, that the total excitation number is not conserved in the case of the Rabi lattice.
The terms 6 and 7 involve two sites and describe conditional photon hopping and twomode photon pair creation or annihilation. Remarkably, in both cases the amplitude for these processes depends on the two-qubit operator (σ z j + σ z j ) including the z-projections of the qubits on the two sites involved in the hopping or pair creation. Assuming qubit configurations composed of σ z eigenstates, hopping and pair creation can be enhanced or suppressed by choosing qubits on the corresponding sites to be aligned or anti-aligned. Again, overall coupling strengths are fixed byC + m which favors hopping and pair creation across small distances m.

Reduction to the Jaynes-Cummings limit
When we reduce the strength g of the Rabi coupling sufficiently to reach the limit g ω, typical of the Jaynes-Cummings model, we can apply the RWA and drop counter-rotating terms. It is instructive to consider how, in this case, the effective Hamiltonian (11) reduces to the dispersive regime of the Jaynes-Cummings lattice.
Neglecting counter-rotating terms, the onsite interaction simplifies to the Jaynes-Cummings interaction Neglecting all counter-rotating terms in an analogous derivation of the dispersive Hamiltonian, we only obtain contributions from paths I and II. Note that, generally, smaller energy differences between the intermediate and initial/final levels in figure 2 lead to larger effective coupling. Further, Λ and V -shaped paths have larger effective coupling due to a constructive sum of the two inverse energy differences. Finally, the energy difference between the initial and final states for each path, when compared to the magnitude of the effective coupling, determines whether or not a path contributes within the RWA.
The resulting effective Hamiltonian describing the Jaynes-Cummings lattice in the dispersive regime reads where we have again dropped a global constant. We define and approximate the involved coupling constants C m by where we assume 0 ≤ m ≤ N/2. (Outside this range, the same substitution m → |m mod N | applies.) Note that the difference betweenC + m andC − m disappears once the counter-rotating term ∼ 1/Σ k is dropped.
The photon-mediated qubit-qubit interaction captured by the second term on the righthand side of equation (16) now has the typical flip-flop (XY) form reminiscent of the wellknown "quantum bus" interaction in the context of multiple qubits in a single resonator [23]. According to the form of the coupling constants C m , this interaction is again short-range and decreases exponentially with the distance between lattice sites. Mediation of this interaction requires a virtual photon to hop across m lattice sites, thus explaining the factor (t/∆) m in C m responsible for the short-range nature. The third term describes the same conditional photon hopping discussed for the Rabi lattice above. The fourth and final term combines the AC Stark and Lamb shifts of the on-site qubit-photon system, in agreement with the results in Ref. [13] when using the weak-coupling approximation C 0 ≈ 1/∆.

Physics of the Jaynes-Cummings and Rabi lattice model in the dispersive regime
Several previous studies of the single-site Rabi and Rabi lattice models have shown that interesting ground-state and steady-state properties emerge in the ultra-strong coupling regime [24,25,26,27,28]. As a particular example, the ground state of the Rabi lattice is believed to undergo a quantum phase transition involving symmetry breaking of the Z 2 parity in the ultra-strong coupling regime when qubit and photons are on resonance [27]. Here, we report that also the off-resonant, dispersive Rabi regime shows interesting ground-state properties when either the qubit or the photon frequency is comparable to the interaction strength g. Depending on the sign of the detuning ∆ = − ω, we distinguish between the dispersive Rabi regime with negative and positive detuning, respectively.
In the following three subsections, we discuss these two regimes along with the dispersive regime of the Jaynes-Cummings lattice. The most interesting aspect of the dispersive Jaynes-Cummings regime is the effective qubit-qubit interaction of XY-type and the possibility of next-nearest-neighbor frustration. In the negative-detuning dispersive Rabi regime, this interaction turns into an effective transverse-Ising model [27,28], which predicts the same kind of phase transition as in Ref. [27]. The effective Hamiltonian we obtain includes additional non-nearest-neighbor Ising-type interactions, which were not considered in Ref. [28]. In the positive-detuning dispersive Rabi regime, we obtain an effective photonic Hamiltonian which shows interesting one-mode and two-mode squeezing in its ground state.

Qubit-qubit interaction in the dispersive Jaynes-Cummings regime
In the Jaynes-Cummings regime, the interaction described by equation (15) only swaps qubit and resonator excitations. As a consequence, the Hamiltonian has a U (1) symmetry and the total number of excitations N tot = j (a † j a j + σ + j σ − j ) is conserved. In the dispersive Jaynes-Cummings regime, the inter-conversion of qubit and photon excitations is suppressed. By switching to a dressed-state basis, the Schrieffer-Wolff transformation eliminates this interaction term. In that dressed-state basis, the total numbers of photon excitations, N ph = j a † j a j , and of qubit excitations, N qu = j σ + j σ − j , are conserved separately and the effective Hamiltonian, equation (16), hence possesses a U (1) × U (1) symmetry. In other words, the effective Hamiltonian separates into blocks with fixed N ph and N qu , thus greatly simplifying the numerical diagonalization.
The interaction terms emerging in the second-order treatment of the Jaynes-Cummings lattice are ordinary AC Stark shifts, conditional photon hopping terms, and qubit-qubit interaction of flip-flop (XY) type, see equation (16). We now focus on the qubit-qubit interaction. As usual, the effective flip-flop interaction can be rewritten as an XY interaction between (pseudo) spins, with an interaction strength given by J m . Keeping the leading order term in t/∆, we can approximate the interaction strength as An interesting fact to note is that J m can be tuned with the detuning ∆ and the photon hopping strength t. It is thus conceivable to engineer both "ferromagnetic" (FM) or "antiferromagnetic" (AF) qubit-qubit interactions, including the possibility of terms leading to non-nearest-neighbor frustration. The signs for J 1 (nearest-neighbor), J 2 (next-nearestneighbor), and the presence or absence of frustration is summarized in Table 1 for the configurations that can occur.
It is well known that the 1D antiferromagnetic J 1 -J 2 XY and Heisenberg models show a phase transition related to frustration and spontaneous dimerization [29,30,31,32,33]. According to Refs. [32,33], the critical point for the J 1 -J 2 XY model is given by J 2 /J 1 = 0.32, which could indeed be accessible with the dispersive Jaynes-Cummings lattice where J 2 /J 1 = t/∆. [Recall: the necessary condition |t/∆| < 1/2 mentioned subsequent to equation (12) is indeed weaker and compatible with this value.] The frustration physics for the Jaynes-Cummings lattice, however, is more intricate: next-nearest neighbor frustration only occurs in the positive-detuning regime where the participation of photons is unavoidable. In order to determine the fate of the phase transition, one thus needs to investigate the relevance of the photon terms in equation (11) for the infinite chain near criticality, which is beyond the scope of this paper.

Dispersive Rabi regime for negative detuning (g ∼ ω)
In the case of negative detuning, the qubit frequency is small compared to the photon frequency ω and we may perform a series expansion in the small parameter /ω 1. Referring to equation (14), we find that the resulting coupling constants scale as  (16) and (18). Analogous statements hold for the dispersive Rabi lattice, equations (11) and (22), where Jm should be replaced byJm.
detuning hopping nearest neighbor next-nearest neighbor frustration?
from which we infer thatC + m C − m as long as m is sufficiently small. By inspection, we note that there are only two terms in the effective Hamiltonian (11) which do not conserve photon number. These correspond to photon pair creation contributions with individual strengths set by g 2C + m with m ≥ 1. For large photon energies ω , t g, it is clear that the inequality is satisfied automatically. Consequently, photon pair creation is strongly suppressed. Therefore, the resulting ground state is expected to be essentially free of photons in the dressed-state basis. These arguments lead us to the following dispersive Rabi model at negative detuning, valid in the N ph = 0 manifold: Here, the Ising coupling and Lamb shift parameter are given bỹ As before, the interaction strengthJ m decreases exponentially with the lattice site distance m. The relevant signs of the coupling parameters can be inferred from Table 1 (substituting J m →J m ). The renormalized qubit frequency = +K 0 plays the role of the effective "magnetic field" (up to a factor of 2), and tends to align the pseudo spin in negative z direction. However, the counter-rotating terms of the σ x j σ x j interaction compete with this tendency by tilting the pseudo spin towards the xy-plane.
For nearest-neighbor coupling only (i.e.,J m ≈ 0 for m > 1), we obtain a simple transverse-Ising model, which is exactly solvable via Jordan-Wigner transformation [34]. In the thermodynamic limit (infinite chain), one recovers the usual quantum phase transition between a paramagnetic and a ferromagnetic phase. The transition occurs atJ 1 = /2, which here implies a critical coupling of g = ω 2 t . (Corrections from weak non-nearestneighbor interaction are expected to slightly shift this transition point.) Below g , the system is in the "paramagnetic" (PM) phase -with the ground state approximately given by the trivial vacuum state |g PM ≈ |0 ph ⊗ | ↓↓↓ · · · ↓ , where |0 ph denotes the photon vacuum and | ↓ the low-energy σ z eigenstate of each qubit. Above g , the system is in the "ferromagnetic" (FM) phase (assuming t > 0). The two symmetry-broken ground-state wavefunctions of the system are approximately given by |g R FM ≈ |0 ph ⊗ |→→→ · · · → and |g L FM ≈ |0 ph ⊗ |←←← · · · ← , where |→ and |← are the two σ x eigenstates. The latter states, marked with "˜", denote eigenstates of the effective Hamiltonian. To obtain the eigenstates of the original Hamiltonian, we perform the inverse Schrieffer-Wolff transformation |g = e −iS1 |g with the previously used generator from equation (7). We illustrate the inverse transformation for the weak-hopping limit, t → 0. In this case, the generator reduces to and the inverse transformation decouples into mere onsite terms. Further approximating ∆ ≈ −ω and Σ ≈ ω, we obtain For the first of the two "FM" ground states, we evaluate Since each qubit is in a σ x eigenstate, we recognize the remaining operator acting on the photon vacuum as a displacement operator D(∓ g ω ) [35], producing a coherent photon state on each site: In a similar way, we find for the paramagnetic ground state in the t → 0 limit. We note that these results are consistent with those obtained by a different method in Ref. [28]. In the general case, and particularly for quantitative comparison, the inverse transformation must be carried out for arbitrary hopping strength t, leading to further corrections to equations (27) and (28). We will properly account for this in our discussion in section 4. Finally, we note that for the N ph > 0 manifolds, the situation is slightly more complicated. According to equation (20), when keeping terms with coefficientC − 1 , we can ignore all terms with coefficientsC + m except for those withC + 0 , which give rise to the Lamb shifts and the AC Stark shifts. Thus, our effective Hamiltonian turns into a transverse-Ising model and a photon tight-binding model, with the only coupling between the two being the AC Stark shifts term, namelỹ

Dispersive Rabi regime for positive detuning (g ∼ ω )
In the dispersive Rabi regime with positive detuning, the photon frequency ω is small compared to the qubit frequency . Thus, we may Taylor-expand in ω/ and find that the coupling constantsC ± m now depend on whether m is even or odd in the following way: ForC − m exactly the same expressions apply, except for an interchange of the roles of 'even' vs. 'odd'. Thus, in addition to the overall decrease in coupling with increasing lattice site distance m, we see thatC + m is further suppressed for odd m, whereasC − m is further suppressed for even m.
Following a line of argument analogous to the one we employed for negative detuning, we note that now the effective Ising coupling is small compared to the qubit frequency, J m . As a result, we may neglect the counter-rotating terms σ + j σ + j + σ − j σ − j of the Ising interaction. The remaining qubit-qubit interaction is then of XY-type, and conserves the total number of qubit excitations N qu . Since we are here interested in the ground state and low-lying states only, we can restrict our discussion to the N qu = 0 subspace. In this subspace, the effective dispersive Hamiltonian consequently takes the form This corresponds to a photon tight-binding model with additional on-site and off-site photon pair creation/annihilation as well as additional, second-order photon hopping terms. Comparison with the expressions forC + m in equations (30) and (31) shows that onsite terms dominate the second-order contributions. Nearest neighbor and next-nearest neighbor terms are suppressed by factors of ωt/ 2 and t 2 / 2 , respectively.
By rewriting the effective Hamiltonian (32) in k space, we obtaiñ The photon dispersion ω k and photon pairing amplitude δ k can be expressed as a Fourier series with coefficients determined byC + m . We approximate them by truncating the series at C + 0 and neglecting higher-order hopping and pairing. That way, we find Here, the pairing amplitude has become k-independent since only on-site pairing is taken into account. We solve this Hamiltonian by performing the Bogoliubov transformation where the coefficients u k and v k can be chosen real-valued and must satisfy u k = u −k , v k = v −k and u 2 k − v 2 k = 1, such that canonical bosonic commutation relations hold for the new operators. As usual, we ensure these conditions by expressing the coefficients in the form u k = cosh r k , v k = sinh r k . Defining the remaining r k parameter via tanh (2r k ) = δ k ω k , the effective Hamiltonian is rendered diagonal, i.e.,H eff Nqu=0 = k E k b † k b k , and the spectrum is given by Next, we show that the ground state (i.e., the 'vacuum' of Bogoliubov excitations) corresponds to a squeezed vacuum state of photons. To see this, we recall the two-mode squeezing operators where ξ k = r k e iϕ k is the squeezing parameter [36,35]. Note that the Bogoliubov transformation is then equivalent to a two-mode squeezing transformation according to Comparing to the results from the Bogoliubov transformation above, we find that the squeezing parameters are given by ϕ k = 0 and For the two special cases of k = 0 or k = π (center and edge of the first Brillouin zone), one obtains one-mode squeezing instead of two-mode squeezing. The ground state of the Rabi lattice in this regime can hence be expressed as a squeezed vacuum state, involving entangled dressed photon pairs with opposite quasi-momenta. + It is useful to point out that equations (36) and (39) also reveal the necessary breakdown of perturbation theory when g exceeds a critical value g c . Specifically, when |δ k | > ω k , the squeezing parameter is ill-defined and E k becomes imaginary. The resulting critical value thus marks the maximum possible value for the convergence radius of the perturbative expansion. In the context of a single Rabi site, the critical coupling strength g c = 1 2 √ ω has previously been derived in Ref. [37]. We note that for positive detuning, g g c is a more restrictive condition than the condition g ∆ min . The inequality g g c thus replaces equation (6) as the necessary condition for the dispersive regime at positive detuning. As g is increased beyond g c , perturbation theory is no longer valid. Entering this quasi-resonant regime of the Rabi lattice, it is plausible that the system will undergo the same type of phase transition with Z 2 symmetry breaking that was studied by Schiro et al. [27] for the case of exact resonance, = ω.
We illustrate the dressed-photon ground state by calculating several observables related to photon numbers and pairing amplitudes. The ground state expectation value for the photon number in mode k is given by Similarly, we find that the pair amplitude for dressed photons of opposite quasi-momenta is + Caveat: the photon pairs mentioned here are indeed dressed photon pairs. To assess the situation in the basis of the original Hamiltonian, the ground state |g = e −iS 1 |g should be calculated, and we will do so in section 4.
No such pairing occurs for photons with identical quasi-momenta, as long as the quasi-momenta k and −k are distinct (i.e., k = 0 and k = π are excluded). All these are the usual properties of two-mode squeezed states [36].

Application to the Rabi dimer and comparison with results from exact diagonalization
In order to illustrate the utility of the effective dispersive Hamiltonian, and to demonstrate its validity by comparison with results from exact numerical diagonalization, we now turn to the specific example of the Rabi dimer, i.e., two coupled Rabi sites. We choose this example to make exact diagonalization of the full Rabi lattice Hamiltonian (1) as tractable as possible, and emphasize that results obtained from the effective Hamiltonian easily carry over to larger lattices, as evident from our equation (11). For verification of our approximations, we select representative observables, and calculate their expectation values by exact diagonalization of the original Rabi lattice Hamiltonian (1). We then compare these results with those obtained from our effective Hamiltonians (29) and (32), which are simplified approximations to Hamiltonian (11) in the negative-and positive-detuning regimes respectively using the full expressions for the coupling constantsC ± m from equation (12). * Figure 4 shows the comparison for the dispersive regime with negative detuning, where example parameters have been chosen as ω = 5 and t = . With this choice of t, the inequality |t/∆| < 1 holds and the dispersive condition (6) can be satisfied. For the Rabi dimer with positive t, equation (6) takes the simple form ∆ min = ω − |t| − g , where ∆ min corresponds to the detuning of the antisymmetric photon mode. Perturbation theory and the effective Hamiltonian are expected to work reasonably well as long as g/∆ min is sufficiently small. This is indeed confirmed by figure 4. The individual results from the four panels are as follows.

The dispersive regime with negative detuning
Panel (a) shows the lowest seven excitation gaps E j − E 0 (E 0 being the ground-state energy and E j the j-th excited-state energy) as a function of the Rabi coupling strength g. The lower three gaps correspond to states in the N ph = 0 manifold, the remaining four to the N ph = 1 manifold. The different curves correspond to exact numerical diagonalization on one hand, and results using the effective Hamiltonian (29) on the other hand. We find very good agreement between approximate and exact results, up to values as high as g/∆ min 0.8. As expected, lower gaps in manifolds with lower photon number match comparatively better. Perturbation theory must break down for g/∆ min ≥ 1, as the qubit reaches quasi-resonance with the anti-symmetric photon mode at this point. Note that for small g/∆ min the 5th to 7th excitation gap differs from the 1st to 3rd excitation gap by a frequency 4 , which corresponds to the photon frequency of the antisymmetric mode. This recurrence illustrates the approximate decoupling of the transverse Ising and tight-binding model in the effective Hamiltonian (29).
Panel (b) shows the fidelity of the approximate ground-state wavefunctions, as obtained from the effective Hamiltonian (29) [or, for the ground state, equivalently to the transverse-Ising model, equation (22)] and subsequent inverse Schrieffer-Wolff (SW) transformation * An additional replacement t → t/2 is performed to account for the fact that a two-site Rabi ring is equivalent to a dimer except for a factor of 2 in the hopping amplitude. |g = e −iS1 |g . The generator S 1 we choose here is the exact expression from Eq. (7) and we do not take the weak hopping (t → 0) limit. For careful comparison, we perform the inverse transformation with two different schemes: (1) we apply the transformation directly in its full exponential form e −iS1 (full inverse SW transformation), and (2) we apply the transformation but keep only terms up to second order in g, namely e −iS1 = 1 − iS 1 + 1 2 (−iS 1 ) 2 + O(g 3 ) (truncated inverse SW transformation). (This truncation would be used for consistently keeping only terms up to second order.) For comparison, we also show the fidelity of the trivial vacuum state |0 ph ⊗ | ↓↓↓ · · · ↓ . We observe that for the g/∆ min → 0 limit, all three fidelities are similar and are very close to a 100% value. For large g, the vacuum fidelity decreases significantly, as expected in the ultra-strong coupling regime. Using the full inverse SW transform, the fidelity of the approximate state is very good, exceeding a 99% value over the full range shown and gives better results as obtained with the truncated transform.
In panels (c) and (d) we show plots for the ground-state expectation values of the qubit excitation σ + L σ − L and the photon number a † L a L (both measured on one of the two Rabi sites). The agreement between exact and approximate results is excellent over the entire range of g in the plot. As before, results are slightly better when using the full inverse SW transform instead of the truncated version. It is interesting to note that the number of qubit excitations on each site rises to about 0.25 around g/ = 2.5, indicating that the qubits are tilting up towards the xy-plane. The tilting is mainly induced by the transverse Ising interaction between qubits, namely the σ x L σ x R term. We have confirmed numerically that this tilting is negligible for a single-site Rabi system, in which no transverse Ising interaction is present. . Although the effective Hamiltonian (22) produces a ground state with zero dressed photons in the negative detuning regime, the undressing from the inverse Schrieffer-Wolff transformation induces small corrections. This effect can be motivated from equations (25)- (27), which indicate that the σ x qubit eigenstate is always accompanied by a photon coherent state on the same site, namely |α j = − g ω j ⊗ | → j . For the finite-size Rabi dimer, the transverse-Ising interaction cannot give rise to an actual phase transition to an ordered spin state. Such a phase transition may, however, manifest in the thermodynamic limit (infinite lattice size) but is beyond the scope of our current paper.

The dispersive regime with positive detuning
We next turn to the dispersive regime of the Rabi dimer with positive detuning. Figure 5 shows a comparison analogous to that presented in figure 4, now with model parameters fixed to = 10ω and t = 0.3ω. Here, we compare the original Rabi-lattice Hamiltonian (1) to the effective Hamiltonian in the N qu = 0 manifold, namely equation (32). We investigate the same set of observables as in the previous subsection and plot them as a function of g, here in units of the critical coupling strength g c = 1 2 (ω − |t|) [here we have already carried out the t → t/2 replacement relative to equation (41)] which is the relevant quantity marking the breakdown of perturbation theory for positive detuning. The critical interaction strength for our choice of parameters is given by g c = 1.32ω.
The lowest five excitation gaps plotted in panel 5(a) show very good agreement between exact and approximate results using the effective Hamiltonian over a wide coupling range, with best agreement for the lowest gap. Here, solid lines and circles represent results from exact diagonalization and diagonalizing the low-lying effective Hamiltonian (32) respectively, while the dot-dashed lines represent results using the analytical expression equation (36) for the quasiparticle energies (the quasimomentum k = 0 and k = π correspond to the symmetric and anti-symmetric photon mode respectively). In the g = 0 limit, original and Bogoliubov operators coincide, b k = a k , and excitations correspond to photon Fock states in the two modes. Specifically, for our parameters the 1st, 3rd and 5th excitation gaps correspond to creation of 1, 2 and 3 photons in the anti-symmetric mode with frequency ω − t = 0.7 . The 2nd excitation gap corresponds to creation of 1 photon in the symmetric mode, which here has frequency ω + t = 1.3 etc.. For g > 0, the Bogoliubov operators b k involve both photon annihilation (a k ) and creation (a † k ) operators. The eigenstates are no longer pure Fock states and the excitation energies decrease as a function of g, as expected from equation (36).
Panel 5(b) shows the ground-state fidelities for our perturbative approximations as well as the trivial vacuum state. For g near 0, the ground state remains quite close to the trivial vacuum state. However, as g is further increased, the expected squeezing of the ground state becomes more significant and the fidelity of the vacuum state drops significantly below the fidelities of our perturbative approximations, which remain very close to 1 until g approaches the critical value g c . The expected breakdown of the perturbative treatment close to g = g c is shown in the inset of panel (b).
Panels 5(c) and (d) show the expected photon number and qubit excitation on one of the two Rabi sites for the Rabi dimer ground state. The exact-diagonalization results based on equation (1) and our perturbative results based on the effective Hamiltonian (32) show good agreement in the relevant range of coupling strengths. Note that, by contrast to the situation of negative detuning, the photon number here exceeds the expected value σ + L σ − L of qubit excitation. Since the ground state is in the N qu = 0 manifold, the creation of qubit excitations is merely due to the qubit-photon dressing and is recovered as a correction from the inverse Schrieffer-Wolff transformation. We note that due to the qubit-photon dressing, an additional increase in the photon number arises from the qubit excitation, as dictated by the inverse Schrieffer-Wolff transformation.
4.2.1. Visualization of squeezing and photon pairing in the Rabi dimer In order to elucidate the effect of the photon pairing or squeezing terms in (32), we briefly discuss results for the Fock-state probability distribution and the Wigner function describing the Rabi dimer. Figure 6 shows the Fock-state probability distribution of the Rabi dimer ground state, where n L , n R denote the photon numbers on each Rabi site and we have traced out the qubit degrees of freedom. The distribution confirms that the pure vacuum state, even though dominant in the distribution, is not the true ground state, as expected for ultra-strong coupling. Expressed in terms of dressed photon states, equation (40) predicts that the ground state should only involve even-number Fock states, n L + n R = 2N , while the probability for odd-number Fock states should vanish. Due to the undressing by the inverse Schrieffer-Wolff transform, corrections to this simple picture emerge. As seen in the comparison of panels (a) and (b), these corrections become more significant as the detuning is decreased. In both cases, however, we clearly observe the fingerprint of photon pairing: even-number Fock states with n L + n R = 2N , have higher probability than their neighboring odd-number Fock states with n L + n R = 2N ± 1. By comparing P 2,0 and P 1,1 , we also observe that onsite pair creation Note that the distribution Pn L ,n R for even-number Fock states, n L + n R = 2N , is much larger than their neighboring odd-number Fock states, n L + n R = 2N ± 1, which evidently serves as a signature of photon pairing. This pairing signature for parameter set (a) is more significant than the one for set (b), due to the larger detuning of set (a) and hence smaller dressing effect which breaks photon parity.
is more significant than offsite pair creation. This agrees well with the fact thatC + 0 is larger thanC + 1 , as we demonstrated in equation (30). An alternative way of visualizing the squeezing predicted by equations (32) and (33) is to calculate and plot the Wigner function for the reduced density matrix representing the photon state on one of the two Rabi sites. Given the Hilbert space structure of the dimer, H dimer = H ph L ⊗ H qu L ⊗ H ph R ⊗ H qu R , we can obtain the desired reduced density matrix directly from the calculated ground state, ρ = tr n R tr σ L tr σ R |g g|.
Therefore, we perform a partial trace of the ground state density matrix over the Hilbert space except for the subspace H ph L , namely the photon Fock space on the left site. The Wigner function is then obtained via W (x, p) = 2π −1 tr [D(−α)ρ D(α) P], where D(α) is the usual displacement operator, α = x + ip and P = exp[iπa † a] is the photon number parity operator [38].
The resulting Wigner function is plotted in figure 7 for several values of g/ω, with the first row showing the result obtained from the full Rabi lattice Hamiltonian (1) and the second row the corresponding result calculated from the dispersive Hamiltonian (11). For g = ω (g/g c = 0.76), the exact and approximate Wigner functions are in good agreement and show  the expected squeezing. Increasing the coupling further to g = 1.3ω (g/g c = 0.98), squeezing becomes more significant and deviations between approximate and exact result are visible as the breakdown point g/g c = 1 is approached. For g = 1.4ω (g/g c = 1.06), we exceed the critical coupling and the Wigner functions differ significantly, signaling the breakdown of perturbation theory.
Overall, our comparison between exact and approximate results for the Rabi dimer in the last two subsections nicely confirms the validity of the perturbative approach in the expected parameter regimes.

Conclusion and Outlook
In summary, we have derived the effective Hamiltonian for the dispersive regime of the Rabi lattice by employing a Schrieffer-Wolff transformation to second order in the Rabi interaction ∼ g. Our results generalize the well-established treatment of the dispersive limit for a single Rabi site to the case of the Rabi lattice, which has enjoyed substantial recent interest in the context of photon-based quantum simulation. The effective interaction terms emerging from our treatment include transverse-Ising interaction between qubits (XY interaction in the Jaynes-Cummings limit), photon pairing terms and conditional photon hopping terms. We have presented analytical expressions for the resulting coupling constants and demonstrated that they are short-range but not restricted to nearest-neighbor sites.
For negative detuning, we found that the effective spin physics is given by a transverse-Ising model which includes interaction terms beyond nearest-neighbor spins. We showed how to recover the dressing of these spin states by photon coherent states via the inverse Schrieffer-Wolff transformation. For positive detuning, we studied the manifold of states without qubit excitations and discussed the effects of one-mode and two mode squeezing. We confirmed the validity of our effective model numerically and discussed in detail the Rabi dimer as the simplest non-trivial example of a Rabi lattice model.
Interesting future work should include extending the perturbative treatment to fourth order, where additional photon-photon interaction is expected due to self-Kerr and cross-Kerr terms. An interesting open question is whether the phase transition discussed in Ref. [27] can be accessed from within the dispersive Rabi regime with positive detuning when including such higher-order terms. Another interesting question concerns the fate of frustration induced phase transitions in the presence of spin-photon dressing. Finally, consideration of the open-system aspect including dissipation and external driving forms an important theoretical challenge in the near future.