Dispersive regime of the Jaynes–Cummings and Rabi lattice

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

$$\begin{equation} V=g \sum_j \left( a_j \sigma^+_j + \mathrm{h.c.} \right). \end{equation} \tag{ 15 }$$

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

$$\begin{eqnarray} &&\fl \tilde{H}_{\mathrm{eff}} =H_0+ g^2 \sum_{j\not=j'} C_{j-j'} \sigma_{j'}^+ \sigma_j^- + \frac{g^2}{2} \sum_{j\not=j'}C_{j-j'} a^{\dagger}_j a_{j'}( \sigma^z_j + \sigma^z_{j'})+ \frac{g^2}{2 } C_0 \sum_{j}( 2 a^{\dagger}_j a_{j}+1) \sigma^z_j \end{eqnarray} \tag{ 16 }$$

where we have again dropped a global constant. We define and approximate the involved coupling constants C_m by

$$\begin{equation} C_m =\frac{1}{N}\sum_k \frac{1 }{\Delta_k}{\rm e}^{{\rm{i}}mk}=\frac{1}{N \Delta}\sum_k \frac{ {\rm e}^{{\rm{i}}mk}}{1-\frac{2t}{\Delta} \cos{k}} \approx \frac{1}{\Delta} \left(\frac{t}{\Delta}\right)^m, \end{equation} \tag{ 17 }$$

where we assume 0 ⩽ m ⩽ N/2. (Outside this range, the same substitution m → |m mod N| applies.) Note that the difference between $\tilde {C}^+_m$ and $\tilde {C}^-_m$ disappears once the counter-rotating term ∼1/Σ_k is dropped.

The photon-mediated qubit–qubit interaction captured by the second term on the right-hand side of equation (16) now has the typical flip–flop (XY) form reminiscent of the well-known '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 [13] when using the weak-coupling approximation C₀ ≈ 1/Δ.

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_{\mathrm {tot}}=\sum _j(a^{\dagger }_j a_j + \sigma ^+_j \sigma ^-_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 SW transformation eliminates this interaction term. In that dressed-state basis, the total numbers of photon excitations, $N_{\mathrm {ph}}=\sum _j a^{\dagger }_j a_j$, and of qubit excitations, $N_{\mathrm {qu}}=\sum _j \sigma ^+_j \sigma ^-_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,

$$\begin{equation} \tilde{H}_{\mathrm{qubit\hbox{-}qubit}}=g^2\sum_{j\not=j'} C_{j-j'} \sigma_{j'}^+ \sigma_j^- = {\sum_{j>j'}} J_{j-j'} (\sigma_{j'}^x \sigma_j^x + \sigma_{j'}^y\sigma_j^y), \end{equation} \tag{ 18 }$$

with an interaction strength given by J_m. Keeping the leading order term in t/Δ, we can approximate the interaction strength as

$$\begin{equation} J_m= \frac{g^2}{2 } C_m\approx \frac{g^2 }{2\Delta}\left(\frac{t}{\Delta}\right)^m. \end{equation} \tag{ 19 }$$

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₁ (nearest-neighbor), J₂ (next-nearest-neighbor), 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 AF J₁–J₂ XY and Heisenberg models show a phase transition related to frustration and spontaneous dimerization [29–33]. According to [32, 33], the critical point for the J₁–J₂ XY model is given by J₂/J₁ = 0.32, which could indeed be accessible with the dispersive Jaynes–Cummings lattice where J₂/J₁ = 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.

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

$$\begin{equation} \fl \tilde{C}^+_m = -\frac{2(m+1)}{\omega}\frac{\epsilon}{\omega}\left(\frac{-t}{\omega}\right)^m + \mathcal{O}(\epsilon^3/\omega^3), \quad \tilde{C}^-_m = -\frac{2}{\omega}\left(\frac{-t}{\omega}\right)^m + \mathcal{O}(\epsilon^2/\omega^2), \end{equation} \tag{ 20 }$$

from which we infer that $\tilde {C}^+_m\ll \tilde {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^2 \tilde {C}^+_m$ with m ⩾ 1. For large photon energies ω ≫ , t g, it is clear that the inequality

$$\begin{equation} \omega \gg g^2 \tilde{C}^+_m \sim \epsilon \bigg(\frac{g}{\omega}\bigg)^2\left(\frac{t}{\omega}\right)^m \end{equation} \tag{ 21 }$$

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:

$$\begin{equation} \tilde{H}_{\mathrm{eff}}\bigg|_{N_{\mathrm{ph}}=0} \approx \frac{\epsilon+\tilde{K}_0}{2}\sum_j \sigma^z_j +\sum_{j> j'} \tilde{J}_{j - j'} \sigma^x_{j} \sigma^x_{j'}. \end{equation} \tag{ 22 }$$

Here, the Ising coupling and Lamb shift parameter are given by

$$\begin{equation} \tilde{J}_m=g^2 \tilde{C}^-_m \approx -\frac{2g^2}{\omega} \left(\frac{-t}{\omega}\right)^m \quad{\rm and}\quad \tilde{K}_0=g^2 \tilde{C}^+_0 \approx -2\frac{g^2\epsilon}{\omega^2}. \end{equation} \tag{ 23 }$$

As before, the interaction strength $\skew3\tilde{J}_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\to \skew3\tilde{J}_m$). The renormalized qubit frequency $\epsilon '=\epsilon +\tilde {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. $\skew3\tilde{J}_m\approx 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 (PM) and a FM phase. The transition occurs at $\skew3\tilde{J}_1=\epsilon /2$, which here implies a critical coupling of $g^\star = \frac {\omega }{2}\sqrt {\frac {\epsilon }{t}}$. (Corrections from weak non-nearest-neighbor interaction are expected to slightly shift this transition point.) Below g^⋆, the system is in the PM phase – with the ground state approximately given by the trivial vacuum state $|\tilde {g}\rangle _{\mathrm {PM}} \approx |0\rangle _{\mathrm {ph}}\otimes |\!\!\downarrow \downarrow \downarrow \cdots \downarrow \rangle$, where |0〉_ph denotes the photon vacuum and |↓〉 the low-energy σ^z eigenstate of each qubit. Above g^⋆, the system is in the FM phase (assuming t > 0). The two symmetry-broken ground-state wavefunctions of the system are approximately given by $|\tilde {g}_{\rm{R}}\rangle _{\mathrm {FM}} \approx |0\rangle _{\mathrm {ph}}\otimes |\!\!\rightarrow \rightarrow \rightarrow \cdots \rightarrow \rangle$ and $|\tilde {g}_{\rm{L}}\rangle _{\mathrm {FM}} \approx |0\rangle _{\mathrm {ph}}\otimes |\!\!\leftarrow \leftarrow \leftarrow \cdots \leftarrow \rangle$, where | → 〉 and |←〉 are the two σ^x eigenstates.

The latter states, marked with '$\tilde{\phantom{g}}$', denote eigenstates of the effective Hamiltonian. To obtain the eigenstates of the original Hamiltonian, we perform the inverse SW transformation $|g\rangle = {\mathrm { e}}^{-{\mathrm { i}}S_1} |\tilde {g}\rangle$ 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

$$\begin{equation} S_1 \approx \mathrm{i} \sum_j \left[ \frac{g}{\Delta}(a_j^{\dagger} \sigma^-_j -\mathrm{h.c.}) + \frac{g}{\Sigma}(a_j \sigma^-_j -\mathrm{h.c.}) \right], \end{equation} \tag{ 24 }$$

and the inverse transformation decouples into mere onsite terms. Further approximating Δ ≈ −ω and Σ ≈ ω, we obtain

$$\begin{equation} |g\rangle = {\rm e}^{-{\rm i}S_1} |\tilde{g}\rangle \approx \prod_j \exp\bigg[ -\frac{g}{\omega}(a_j^{\dagger}-a_j)\sigma^x_j\bigg]|\tilde{g}\rangle. \end{equation} \tag{ 25 }$$

For the first of the two 'FM' ground states, we evaluate

$$\begin{equation} |g_{\rm{R}}\rangle_{\mathrm{FM}} \approx \prod_j \exp\bigg[ -\frac{g}{\omega}(a_j^{\dagger}-a_j)\sigma^x_j\bigg] |0\rangle_{\mathrm{ph}}\otimes|\!\!\rightarrow \rightarrow \rightarrow \cdots \rightarrow\rangle. \end{equation} \tag{ 26 }$$

Since each qubit is in a σ^x eigenstate, we recognize the remaining operator acting on the photon vacuum as a displacement operator $D(\mp \frac {g}{\omega })$ [35], producing a coherent photon state on each site

$$\begin{equation} |g_{\rm{R}}\rangle_{\mathrm{FM}} \approx \prod_j |\alpha_j=-{\textstyle\frac{g}{\omega}}\rangle_j \otimes |\!\!\rightarrow\rangle_j \quad {\rm and} \quad |g_{\rm{L}}\rangle_{\mathrm{FM}} \approx \prod_j |\alpha_j={\textstyle\frac{g}{\omega}}\rangle_j \otimes |\!\!\leftarrow\rangle_j. \end{equation} \tag{ 27 }$$

In a similar way, we find

$$\begin{equation} |g\rangle_{\mathrm{PM}} \approx \prod_j \bigg[|\alpha_j=-{\textstyle\frac{g}{\omega}}\rangle_j \otimes |\!\!\rightarrow\rangle_j\,\,-\,\, |\alpha_j={\textstyle\frac{g}{\omega}}\rangle_j \otimes |\!\!\leftarrow\rangle_j \bigg]. \end{equation} \tag{ 28 }$$

for the PM ground state in the t → 0 limit. We note that these results are consistent with those obtained by a different method in [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 coefficient $\tilde {C}^-_1$, we can ignore all terms with coefficients $\tilde {C}^+_m$ except for those with $\tilde {C}^+_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, namely

$$\begin{eqnarray} &&\fl \tilde{H}_{\mathrm{eff}} \approx \frac{\tilde{K}_0 + \epsilon}{2}\sum_j \sigma^z_j +\sum_{j} \tilde{J}_{1} \sigma^x_{j} \sigma^x_{j+1} + \omega \sum_j a^{\dagger}_j a_j +t\sum_{\langle j,j' \rangle}\left(a_j^{\dagger} a_{j'} + \mathrm{h.c.} \right) \nonumber\\ &&+ g^2 \tilde{C}^+_0 \sum_{j} a^{\dagger}_j a_{j} \sigma^z_j. \end{eqnarray} \tag{ 29 }$$

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 constants $\tilde {C}^{\pm }_m$ now depend on whether m is even or odd in the following way:

$$\begin{equation} {\rm odd}\ m:\quad \tilde{C}^+_m = \frac{2(m+1)}{\epsilon}\frac{\omega}{\epsilon}\left(\frac{t}{\epsilon}\right)^m + \mathcal{O}(\omega^{3}/\epsilon^3), \end{equation} \tag{ 30 }$$

$$\begin{equation} {\rm even}\ m:\quad \tilde{C}^+_m = \frac{2}{\epsilon}\left(\frac{t}{\epsilon} \right)^m + \mathcal{O}(\omega^2/\epsilon^2). \ \end{equation} \tag{ 31 }$$

For $\tilde {C}^-_m$ exactly the same expressions apply, except for an interchange of the roles of 'even' versus 'odd'. Thus, in addition to the overall decrease in coupling with increasing lattice site distance m, we see that $\tilde {C}^+_m$ is further suppressed for odd m, whereas $\tilde {C}^-_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, $\skew3\tilde{J}_m \ll \epsilon$. 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

$$\begin{equation} \fl \tilde{H}_{\mathrm{eff}}\big|_{N_{\mathrm{qu}}=0} \approx \omega \sum_j a^{\dagger}_j a_j + t\sum_ j \left(a_j^{\dagger} a_{j+1} + \mathrm{h.c.} \right) - \frac{g^2}{2} \sum_{j, j'} \tilde{C}^+_{j-j'} (a^{\dagger}_{j}a_{j'}+a^{\dagger}_{j}a^{\dagger}_{j'} + \mathrm{h.c.}). \end{equation} \tag{ 32 }$$

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 for $\tilde {C}^+_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/² and t²/², respectively.

By rewriting the effective Hamiltonian (32) in k space, we obtain

$$\begin{equation} \tilde{H}_{\mathrm{eff}}\big|_{N_{\mathrm{qu}}=0} = \sum_k \left[\omega_k \mathsf{a}^{\dagger}_k \mathsf{a}_k + {\textstyle \frac{1}{2}}\delta_k ( \mathsf{a}_k \mathsf{a}_{-k} + \mathsf{a}^{\dagger}_k \mathsf{a}^{\dagger}_{-k}) \right]. \end{equation} \tag{ 33 }$$

The photon dispersion ω_k and photon pairing amplitude δ_k can be expressed as a Fourier series with coefficients determined by $\tilde {C}^+_m$. We approximate them by truncating the series at $\tilde {C}^+_0$ and neglecting higher-order hopping and pairing. That way, we find

$$\begin{equation} \omega_k \approx \omega+ 2t \cos{k} - 2g^2/\epsilon, \quad \delta_k \approx -2g^2/\epsilon. \end{equation} \tag{ 34 }$$

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

$$\begin{equation} b_k =u_k \mathsf{a}_k + v_k \mathsf{a}^{\dagger}_{-k}, \quad b^{\dagger}_{-k} = v_k \mathsf{a}_k + u_k \mathsf{a}^{\dagger}_{-k}, \end{equation} \tag{ 35 }$$

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²_k − v²_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)}=\frac {\delta _k}{\omega _k}$, the effective Hamiltonian is rendered diagonal, i.e. $\tilde{H}_{\mathrm {eff}}\big |_{N_{\mathrm {qu}}=0} = \sum _k E_k b^{\dagger }_k b_k$, and the spectrum is given by

$$\begin{equation} E_k = \sqrt{\omega_k^2 - \delta_k^2} \approx \sqrt{\left(\omega+2t\cos{k} - 2g^2/\epsilon\right)^2 - 4g^4/\epsilon^2}. \end{equation} \tag{ 36 }$$

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 $\mathcal {S}_2(\xi _k)=\exp [\xi ^*_k \mathsf {a}_k \mathsf {a}_{-k} - \xi _k \mathsf {a}^{\dagger }_k \mathsf {a}^{\dagger }_{-k}]$ where ξ_k = r_ke^iφ_k is the squeezing parameter [35, 36]. Note that the Bogoliubov transformation is then equivalent to a two-mode squeezing transformation according to

$$\begin{equation} b_k = \mathcal{S}_2(\xi_k) \mathsf{a}_k \mathcal{S}^{\dagger}_2(\xi_k) = \mathsf{a}_k \cosh{r_k} +{\rm e}^{{\rm{i}} \varphi_k} \mathsf{a}^{\dagger}_{-k} \sinh{r_k}, \end{equation} \tag{ 37 }$$

$$\begin{equation} b_{-k} = \mathcal{S}_2(\xi_k) \mathsf{a}_{-k} \mathcal{S}^{\dagger}_2(\xi_k) = \mathsf{a}_{-k} \cosh{r_k} +{\rm e}^{{\rm{i}} \varphi_k} \mathsf{a}^{\dagger}_k \sinh{r_k}. \end{equation} \tag{ 38 }$$

Comparing to the results from the Bogoliubov transformation above, we find that the squeezing parameters are given by φ_k = 0 and

$$\begin{equation} r_k=\frac{1}{2}\tanh^{-1}(\delta_k/\omega_k) \approx \frac{1}{2}\tanh^{-1}\left(\frac{-2g^2/\epsilon}{\omega+ 2t \cos{k} -2g^2/\epsilon}\right). \end{equation} \tag{ 39 }$$

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,

$$\begin{equation} |\tilde{g}\rangle= \prod_{k \geqslant 0} \mathcal{S}_2(r_k) |0\rangle = \prod_{k \geqslant 0} \exp\bigg[r_k \mathsf{a}_k \mathsf{a}_{-k}-r_k \mathsf{a}^{\dagger}_k \mathsf{a}^{\dagger}_{-k} \bigg] \,|0\rangle. \end{equation} \tag{ 40 }$$

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 $\left |\delta _k\right |>\omega _k$, the squeezing parameter is ill-defined and E_k becomes imaginary. The resulting critical value

$$\begin{equation} g_{\rm{c}} = \frac{1}{2}\min_k \sqrt{\epsilon(\omega-2t\cos{k})} = \frac{1}{2}\sqrt{\epsilon(\omega-2|t|)} \end{equation} \tag{ 41 }$$

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_{\rm{c}}=\frac {1}{2}\sqrt {\epsilon \omega }$ has previously been derived in [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₂ 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

$$\begin{equation} \langle\, \tilde{g} \,|\, \mathsf{a}^{\dagger}_k \mathsf{a}_k \,|\, \tilde{g} \,\rangle=\langle\, 0 \,|\, \mathcal{S}^{\dagger}_2(r_k)\mathsf{a}^{\dagger}_k \mathcal{S}_2(r_k)\mathcal{S}^{\dagger}_2(r_k)\mathsf{a}_k \mathcal{S}_2(r_k) \,|\, 0 \,\rangle = \sinh^2(r_k). \end{equation} \tag{ 42 }$$

Similarly, we find that the pair amplitude for dressed photons of opposite quasi-momenta is

$$\begin{equation} \langle\, \tilde{g} \,|\, \mathsf{a}_k \mathsf{a}_{-k} \,|\, \tilde{g} \,\rangle=\langle\, \tilde{g} \,|\, \mathsf{a}^{\dagger}_k \mathsf{a}^{\dagger}_{-k} \,|\, \tilde{g} \,\rangle = -\sinh{r_k} \cosh{r_k}. \end{equation} \tag{ 43 }$$

No such pairing occurs for photons with identical quasi-momenta,

$$\begin{equation} \langle\, \tilde{g} \,|\, \mathsf{a}_k \mathsf{a}_{k} \,|\, \tilde{g} \,\rangle=\langle\, \tilde{g} \,|\, \mathsf{a}^{\dagger}_k \mathsf{a}^{\dagger}_k \,|\, \tilde{g}\,\rangle = 0, \end{equation} \tag{ 44 }$$

as long as the quasi-momenta k and −k are distinct (i.e. k = 0 and π are excluded). All these are the usual properties of two-mode squeezed states [36].

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.

Zoom In Zoom Out Reset image size

Panel (a) shows the lowest seven excitation gaps E_j − E₀ (E₀ being the ground-state energy and E_j the jth 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 fifth to seventh excitation gap differs from the first to third 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 SW transformation $|g\rangle = {\mathrm { e}}^{-{\mathrm { i}}S_1} |\tilde {g}\rangle$. The generator S₁ we choose here is the exact expression from equation (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^−iS₁ (full inverse SW transformation), and (2) we apply the transformation but keep only terms up to second order in g, namely ${\mathrm { e}}^{-{\mathrm { i}}S_1}=1-\mathrm {i}S_1 +\frac {1}{2} (-\mathrm {i}S_1)^2 + \mathcal {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^†_La_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 SW 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 $|\alpha _j={-\textstyle \frac {g}{\omega }}\rangle _j \otimes |\!\!\rightarrow \rangle _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.

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_{\rm{c}}=\frac {1}{2}\sqrt {\epsilon (\omega -|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 ω.

Zoom In Zoom Out Reset image size

The lowest five excitation gaps plotted in figure 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 quasi-momentum k = 0 and π 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 first, third and fifth excitation gaps correspond to creation of 1, 2 and 3 photons in the anti-symmetric mode with frequency ω − t = 0.7. The second 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).

Figure 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 SW 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 SW 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,

$$\begin{equation} P_{n_{\rm{L}},n_{\rm{R}}}=|\mathrm{tr}_{\sigma_{\rm{L}},\sigma_{\rm{R}}}\langle n_{\rm{L}},\sigma_{\rm{L}};n_{\rm{R}},\sigma_{\rm{R}}|g\rangle|^2, \end{equation} \tag{ 45 }$$

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 SW 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 is more significant than offsite pair creation. This agrees well with the fact that $\tilde {C}^+_0$ is larger than $\tilde {C}^+_1$, as we demonstrated in equation (30).

Zoom In Zoom Out Reset image size

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, $\mathcal {H}_{\mathrm {dimer}}= \mathcal {H}^{\mathrm {ph}}_{\rm{L}} \otimes \mathcal {H}^{\mathrm {qu}}_{\rm{L}} \otimes \mathcal {H}^{\mathrm {ph}}_{\rm{R}} \otimes \mathcal {H}^{\mathrm {qu}}_{\rm{R}}$, we can obtain the desired reduced density matrix directly from the calculated ground state,

$$\begin{equation} \rho = \mathrm{tr}_{n_{\rm{R}}}\,\mathrm{tr}_{\sigma_{\rm{L}}}\,\mathrm{tr}_{\sigma_{\rm{R}}} \;|g\rangle\langle g|. \end{equation} \tag{ 46 }$$

Therefore, we perform a partial trace of the ground state density matrix over the Hilbert space except for the subspace $\mathcal {H}^{\mathrm {ph}}_{\rm{L}}$, namely the photon Fock space on the left site. The Wigner function is then obtained via $W(x,p)=2\pi ^{-1}\mathrm {tr}\,[D(-\alpha)\rho \, D(\alpha)\,\mathcal {P}]$, where D(α) is the usual displacement operator, α = x + ip and $\mathcal {P}=\exp [{\mathrm { i}}\pi a^{\dagger } 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.

Zoom In Zoom Out Reset image size

Overall, our comparison between exact and approximate results for the Rabi dimer in the last two sections nicely confirms the validity of the perturbative approach in the expected parameter regimes.