Two-photon Rabi-Hubbard and Jaynes-Cummings-Hubbard models: photon pair superradiance, Mott insulator and normal phases

TL;DR

This paper investigates the phase diagrams of two-photon Rabi-Hubbard, Jaynes-Cummings-Hubbard, and Dicke models, revealing unique superradiant, Mott insulator, and normal phases through various computational methods.

Contribution

It provides a comprehensive comparison of mean field predictions with exact QMC results and uncovers novel phase behaviors in two-photon models.

Findings

Excellent agreement between mean field and QMC for Dicke model.

JCH model shows a single Mott lobe with photon pair superfluid phase.

No superradiant transition in the Rabi-Hubbard model with QMC.

Abstract

We study the ground state phase diagrams of two-photon Dicke, the one-dimensional Jaynes-Cummings-Hubbard (JCH), and Rabi-Hubbard (RH) models using mean field, perturbation, quantum Monte Carlo (QMC), and density matrix renormalization group (DMRG) methods. We first compare mean field predictions for the phase diagram of the Dicke model with exact QMC results and find excellent agreement. The phase diagram of the JCH model is then shown to exhibit a single Mott insulator lobe with two excitons per site, a superfluid (SF, superradiant) phase and a large region of instability where the Hamiltonian becomes unbounded. Unlike the one-photon model, there are no higher Mott lobes. Also unlike the one-photon case, the SF phases above and below the Mott are surprisingly different: Below the Mott, the SF is that of photon {\it pairs} as opposed to above the Mott where it is SF of simple photons. The mean field phase diagram of the RH model predicts a transition from a normal to a superradiant phase but none is found with QMC.

Equations

$$H_D$$

$$H_{RH}$$

$$\mathcal{R}(\theta )=\exp \left(i\theta {\hat{a}}_j^{†}{\hat{a}}_j^{}+i2\theta {\sigma }_j^{+}{\sigma }_j^{-}\right),$$

$$G_{\alpha ,\beta }(r)\equiv \frac{1}{2N}\sum _i⟨{\alpha }_i{\beta }_{i+r}+h.c.⟩,$$

$$G_{a^{†},a^{}}(r)=\frac{1}{2N}\sum _i⟨{\hat{a}}_i^{†}{\hat{a}}_{i+r}^{}+{\hat{a}}_{i+r}^{†}{\hat{a}}_i^{}⟩.$$

$$G_{{\sigma }^{-},{\sigma }^{+}}(r)=\frac{1}{2N}\sum _i⟨{\sigma }_i^{-}{\sigma }_{i+r}^{+}+{\sigma }_{i+r}^{-}{\sigma }_i^{+}⟩,$$

$$G_{{\sigma }^{+},a}(r)=\frac{1}{2N}\sum _i⟨{\hat{a}}_i^{}{\sigma }_{i+r}^{+}+{\sigma }_{i+r}^{-}{\hat{a}}_i^{†}⟩,$$

$$G_{{\sigma }^{+},a^2}(r)=\frac{1}{2N}\sum _i⟨{\hat{a}}_i^2{\sigma }_{i+r}^{+}+{\sigma }_{i+r}^{-}{\hat{a}}_i^{†2}⟩.$$

$$N_{exc}=\sum _i⟨{\hat{a}}_i^{†}{\hat{a}}_i^{}+2{\sigma }_i^{+}{\sigma }_i^{-}⟩,$$

$${\rho }_s=\frac{⟨W^2⟩}{2t\beta L^{d-2}},$$

$${\hat{J}}^{+}={\hat{b}}^{†}\sqrt{N-{\hat{b}}^{†}{\hat{b}}^{}},{\hat{J}}^{-}=\sqrt{N-{\hat{b}}^{†}{\hat{b}}^{}}{\hat{b}}^{},{\hat{J}}_z={\hat{b}}^{†}{\hat{b}}^{}-\frac{N}{2},$$

$$H_D^{MF}=\omega {\hat{a}}^{†}\hat{a}+g_{\bar{b}}({\hat{a}}^2+{\hat{a}}^{†2})+{\omega }_q|\bar{b}{|}^2-\frac{{\omega }_{q}N}{2},$$

$$g_{\bar{b}}=\frac{g}{N}(\bar{b}+{\bar{b}}^{\ast })\sqrt{N-|\bar{b}{|}^2}.$$

$$\bar{b}$$

$$\lambda$$

$$\mu$$

$$P(x)=\frac{a+bx+c{x}^2+d{x}^3}{1+ex+f{x}^2+h{x}^3};$$

$$H_{JC}$$

$$|n,+⟩=\sin {\theta }_n|n,g⟩+\cos {\theta }_n|n-2,e⟩$$

$$|n,-⟩=\cos {\theta }_n|n,g⟩-\sin {\theta }_n|n-2,e⟩$$

$$E_{n\pm }=n\omega +\Delta /2\pm R_n(\Delta )/2.$$

$$\begin{split}&2\omega+(R_{2}-R_{3})-4Jt^{(3)2}_{--}+8J^{2}\Biggl{[}t_{--}^{(2)2}\\ &\left(-\sum_{\sigma=\pm}\frac{2t_{\sigma-}^{(3)2}}{\Delta-2R_{2}-\sigma R_{3}}+\frac{t_{--}^{(3)2}}{\Delta-2R_{2}+R_{3}}\right)\\ &+\sum_{\sigma=\pm}\frac{t_{\sigma-}^{(4)2}}{\Delta-R_{2}-R_{3}-\sigma R_{4}}\\ &+\frac{t_{+-}^{(3)2}t_{-+}^{(3)2}}{-2R_{2}-2R_{3}}-\frac{t_{+-}^{(3)2}t_{--}^{(3)2}}{R_{3}}-\frac{t_{--}^{(3)2}t_{-+}^{(3)2}}{-2R_{2}}\Biggr{]}=0,\end{split}$$

$$t^{(n)}_{\sigma\nu}=\mbox{$\langle n,\sigma|$}{\hat{a}}^{\dagger}\mbox{$|n-1,\nu\rangle$},$$

$$\begin{split}4\omega+\Delta-R_{2}+16J^{2}t_{--}^{(2)}\biggl{[}\frac{t_{--}^{(3)2}}{\Delta+R_{3}-2R_{2}}\\ +\frac{t_{+-}^{(3)2}}{\Delta-R_{3}-2R_{2}}-\frac{1}{\Delta-R_{2}}\biggr{]}=0.\end{split}$$

$$4w+\Delta -R_2+\frac{16J^2{t}_{--}^{(2)2}}{\Delta -R_2}=0.$$

$$E_G(\bar{b}\ne 0,g,J,\omega ,{\omega }_q)=E_G(\bar{b}=0,g,J,\omega ,{\omega }_q),$$

$$E_G$$

$$\omega -2J\ge 2g.$$

$$H_{RH}$$

$${\sigma }_i^{+}$$

$${\sigma }_i^{-}$$

$${\bar{b}}^{\ast }$$

Full text

I Introduction

Continued progress in controlling and tuning interactions between photons and individual atoms has made possible the realization of elementary quantum electrodynamics building blocksschoelkopf08 using two-level atoms in cavitiesharoche06 or Josephson junctions on solid state chipswallraff04 ; chen07 ; lambert09 ; nataf10 ; viehmann11 (circuit QED). Generally speaking, such couplings of photons to two-level systems are well described by the Rabi modelrabi36 ; rabi37 . In what is referred to as the strong coupling limit, $g/\omega\lesssim 0.1$, (see below) one can apply the random wave approximation (RWA) and ignore “counter rotating” (CR) terms which do not conserve the number of excitons (the number of photons plus the number of excited atoms). When these terms are ignored, we obtain the Jaynes-Cummings modeljaynes63 which, due its $U(1)$ symmetry, conserves the number of excitons and can be solved exactly for a single cavity. An ensemble of such cavities can be connected (here we consider a one-dimensional chain) by tuning the tunneling rate of the photon modes between near neighbor cavities resulting in a lattice model, the Jaynes-Cummings-Hubbard model, which consists of itinerant bosons (photons) hopping between near neighbor sites and interacting with localized two-level atoms (“spins” or qubits). The properties of this model were shown to be very similar to those of the one-dimensional Bose-Hubbard model (BHM)batrouni90 exhibiting a superfluid phase (of excitons) and incompressible Mott insulator (MI) lobesgreentree06 ; hartmann06 ; rossini07 ; koch08 ; hohenadler08 ; hartmann08 ; schmidt09 ; zhao08 . The MI is essentially a product single site states composed of a superposition of photons and excited atomsschmidt09 and exhibits behavior similar to photon blockadebirnbaum05 where there is a finite energy gap opposing the addition of a photon. Long range frustrated hoppinghohenadler12 and dynamic properties of this model have also elicited much interestpippan09 ; tomadin10 ; tomadin10b ; hohenadler11 ; carusotto13 . When the coupling becomes of the order of the photon frequency, $g/\omega$, the “ultra-strong coupling” regime which has been achieved experimentallyniemczyk10 ; forndiaz10 ; chen17 ; forndiaz17 ; yoshihara17 , the RWA is no longer valid and the counter rotating terms come into play reducing the symmetry from $U(1)$ to $Z_{2}$. This has several consequences: (a) CR terms cause the exciton number not to be conserved thus excluding the possibility of any Mott phases, (b) the system is now in the universality class of the Ising model with two phases, a disordered and an ordered (coherent) phase, (c) due to the discrete nature of the symmetry, it can break spontaneously in one dimension. Consequently, the ordered coherent phase is a photon Bose-Einstein condensate (BEC) and the transition between this phase and the disordered phase exhibits Ising critical exponentsschiro13 ; zheng11 ; schiro13 ; kumar13 ; flottat16 . This transition resembles the incoherent/coherent (normal/superradiant) phase transition in the Dicke modelhepp73 ; rotondo15 .

Multi-photon processes have also attracted scrutiny as more experimental systems are being realized where such physics enters into play. For example, the two-photon Rabi model has been used as an effective model to describe second order processes in Rydberg atoms in cavitiesbertet02 and quantum dotsstufler06 ; delvalle10 , and mechanisms have been proposed to realize it in circuit QEDfelicetti18 . Unlike the one-photon case, the two-photon Rabi model undergoes spectral collapse, where the Hamiltonian is no longer bounded from below, when the coupling exceeds a certain valueemary02 ; dolya09 ; travenec12 ; maciejewski15 ; travenec15 ; peng17 ; chen12 ; felicetti15 . As in the one-photon case, in the strong coupling limit where the CR terms can be ignored, the Hamiltonian has $U(1)$ symmetry. In the ultra strong coupling limit, where the CR terms must be restored, the symmetry is reduced to $Z_{2n}$ for the $n$-photon case, i.e. $Z_{4}$ in the two-photon case. The ground state of the many-body two-photon model was studied in the context of the Dicke modeldicke54 using mean field methodsgarbe17 . A quantum phase transition between a normal (disordered) and superradiant phase was predicted. In this model, however, the term superradiant is used to indicate a macroscopic change in the average number of photons but which remains relatively small. This is to be contrasted with the one-photon case where there is a very large number of photons in the superradiant phase and they form a BECschiro13 ; zheng11 ; schiro13 ; kumar13 ; flottat16 . The spectral collapse region was also found.

In this paper we will first use exact QMC simulations to determine the phase diagram of the two-photon Dicke model and compare with mean field resultsgarbe17 . This also serves to verify our QMC method. We then use our QMC algorithm, and also DMRG calculations, to study the phase diagram of another many-body system, namely the one-dimensional two-photon Jaynes-Cummings-Hubbard (JCH) model. The numerical results are compared with perturbation and mean field calculations. Finally, we use our numerical methods to study the full Rabi-Hubbard model.

Our results below show that the mean field phase diagram of the Dicke model is in agreement with our exact QMC simulations. In addition, we find that, surprisingly, the Jaynes-Cummings-Hubbard model exhibits only one Mott insulating lobe before the systems becomes unstable. Furthermore, doping above the MI yields a photon superfluid phase, while doping below the MI yields an unexpected photon pair superfluid phase. QMC simulations suggest that the Rabi-Hubbard model does not exhibit a phase transition from a normal to a superradiant phase.

The paper is organized as follows. In section II we describe the models and the physical quantities we will study. In section III we review the mean field calculation for the Dicke model and discuss our exact QMC results. In section IV we show our perturbation and exact QMC results for the JCH model followed by section V where we discuss the RH model. We present some conclusions in section VI followed by appendices A, B, and C where we show details of some of our calculations.

II Models

We will first study the two-photon Dicke model where $N$ two-level systems (qubits) couple to a single photon mode and which is governed by the Hamiltonian

[TABLE]

This model was examined with mean field in Ref. garbe17, . Here, $\omega$ is the photon frequency, $\omega_{q}$ the qubit energy spacing, $g$ the coupling constant, $\hat{a}$ (${\hat{a}}^{\dagger}$) is the photon destruction (creation) operator, $\sigma^{z}_{i}$ and $\sigma^{x}_{j}=\sigma^{+}_{j}+\sigma^{-}_{j}$ are Pauli matrices acting on the $j$th qubit, $\sigma^{+}_{j}$ ($\sigma^{-}_{j}$) is the corresponding raising (lowering) operator. A related model is the Rabi-Hubbard (RH) model,

[TABLE]

where now $N$ is the number of sites (or cavities) and ${\hat{a}}_{i}$ is the photon mode in the $i$th cavity. Note that photon modes can tunnel between sites. Ignoring in Eq. (2) the CR terms, $(\sigma^{+}_{i}{\hat{a}}^{\dagger 2}_{i}+\sigma^{-}_{i}{\hat{a}}^{2}_{i}$), yields the Jaynes-Cummings-Hubbard (JCH) model in which the number of excitons is conserved.

In both models, Eq. (1) and Eq. (2), when the CR terms are dropped, the system is invariant under the generalized rotation operator,

[TABLE]

with ${\cal R}(\theta)^{\dagger}{\hat{a}}_{j}{\cal R}(\theta)^{\phantom{\dagger}}={\rm e}^{i\theta}{\hat{a}}_{j}$, and ${\cal R}(\theta)^{\dagger}\sigma^{-}_{j}{\cal R}(\theta)^{\phantom{\dagger}}={\rm e}^{i2\theta}\sigma^{-}_{j}$ for any $\theta$, thus exhibiting $U(1)$ symmetry. However, the action of ${\cal R}(\theta)$ on the CR terms introduces a phase ${\rm exp}(i4\theta)$ thus reducing the symmetry to $Z_{4}$, in other words, the system is left invariant by the rotation only for $\theta=n2\pi/4$, $n=0,1,2,3$. Therefore, when the CR terms are ignored, the $U(1)$ symmetry results in the conservation of the number of excitons, $N_{exc}=N_{photon}+2N_{+}$ where $N_{+}$ is the number of qubits in the excited state. De-exciting a qubit generates two photons and vice versa. On the other hand, in the presence of the CR terms, $N_{exc}$ is not conserved, and quantum phase transitions would be expected to reflect the discrete $Z_{4}$ symmetry. Note that this discrete symmetry can break spontaneously for a one-dimensional quantum system in its ground state.

To characterize the various possible phases, we calculate several Green functions,

[TABLE]

where $\alpha$ and $\beta$ denote creation and annihilation operators of the photons (${\hat{a}}^{\dagger}_{i}$ and ${\hat{a}}^{\phantom{\dagger}}_{i}$) or the qubits ($\sigma^{+}_{j}$ and $\sigma^{-}_{j}$). For example, the photon Green function at equal time is given by,

[TABLE]

The qubit Green function is,

[TABLE]

and the following two functions will be particularly useful:

[TABLE]

and

[TABLE]

Power law decay of one of these Green functions would indicate quasi-long range order for the corresponding quantity. We also measure the average number of excitons,

[TABLE]

and the superfluid density,

[TABLE]

where $W$ is the winding number of an exciton and $d$ the dimensionality. To study the phase diagrams of these models, we use several methods: Mean field, perturbation expansion, QMC and DMRG.

III The Dicke Model

First, we will review briefly the mean field results of Ref. garbe17, and then discuss our QMC results. First the total angular momentum operators are defined, ${\hat{J}}_{z}\equiv\sum_{i}\sigma^{z}_{i}/2$ and ${\hat{J}}^{\pm}\equiv\sum_{i}\sigma^{\pm}_{i}$. Then using the Holstein-Primakoff (HP) transformation, they define the bosonic operators:

[TABLE]

where $[{\hat{b}},{\hat{b}}^{\dagger}]=1$. Then the bosonic operators are replaced by the average value in the ground state, ${\hat{b}}\to{\bar{b}}=\langle GS|{\hat{b}}|GS\rangle$. Making the substitutions in the Hamiltonian, Eq. (1), yields,

[TABLE]

with

[TABLE]

The Hamiltonian is now a quadratic in the photon operators and can be diagonalized exactly giving the ground state energy, $E_{G}$, which is to be minimized with respect to ${\bar{b}}$. This determines the value of the order parameter, ${\bar{b}}$ as a function of the other parameters $\omega$, $\omega_{q}$ and $g$. It was found that for $g<g_{c}=\sqrt{\omega\omega_{q}N/4}$, $E_{G}$ is minimum for ${\bar{b}}={\bar{b}}^{*}=0$. For $g>g_{c}$, the ground state is twofold degenerate and the order parameter, acquires a nonzero value, ${\bar{b}}={\bar{b}}^{*}$:

[TABLE]

Therefore, the system exhibits two phases, a normal phase (${\bar{b}}=0$) and a superradiant phase (${\bar{b}}\neq 0$) separated by the transition line $g_{c}=\sqrt{\omega\omega_{q}N/4}$. Furthermore, for $g>\omega/2$, the argument of the square root in Eq. (14) becomes negative indicating that the Hamiltonian is unbounded. These results map out the phase diagram as shown in Ref. garbe17, .

It is interesting to examine the accuracy of mean field calculations in various situations. We therefore performed QMC simulations using the stochastic Green function (SGF) methodrousseau08 ; rousseau08b , both to verify our QMC simulations and obtain the numerically exact ground state phase diagram. To map out the phase diagram, we choose $\omega=1$ to set the energy scale, and we calculate the order parameter for many values of $g$ at fixed $N\omega_{q}$. Making several such cuts for different $N\omega_{q}$ yields the phase diagram in the $(N\omega_{q},g)$ plane. In the mean field calculation, the order parameter was ${\bar{b}}=\langle GS|b|GS\rangle$; in the QMC simulations, we take the order parameter to be $\langle N_{q}\rangle$, the average number of excited qubits which corresponds to $|{\bar{b}}|^{2}$ in the mean field case. Figure 1 shows such a cut in $g$ for a $16$-particle system at $N\omega_{q}=0.1$ and several values of the inverse temperature, $\beta$. It is seen that the finite temperature effects are very pronounced for small values of the coupling $g$, and that to detect the transition properly, the temperature must be very low, and gets lower with increasing $N$.

Next, we study the finite size effects on the behavior of the order parameter. This is presented in Fig. 2 where the top panel shows the dependence of $\langle N_{q}\rangle/N$ on $N$ and $g$. The values of $\beta$ are chosen large enough so that the system is in its ground state. It is seen that there is a change in curvature as $g$ increases, and that the point of maximum slope shifts to smaller values of $g$ as $N$ increases. On the other hand, the lower panel shows that the average photon density does not exhibit any sudden changes; it increases mildly with $g$. This was remarked with the mean field results in Ref. garbe17, .

In order to determine the transition point for each system size, we fit Padé approximants to the order parameter as a function of $g$ for each system size. The approximants we use are of the form,

[TABLE]

the maximum of the derivative yields the transition point for that system size. We then extrapolate the value of the critical $g$ to the thermodynamic limit. To determine the boundary between the superradiant phase and the unstable region, we calculate $\langle N_{ph}\rangle/N$ as a function of $g$ for a fixed $N\omega_{q}$. We find that as $g\to 1/2$, the photon density increases very rapidly and becomes unmanageable at $g=1/2$ indicating the instability. This way we map the phase diagram shown in Fig. 3 and which agrees very well with the MF solution. This excellent agreement can be understood as being due to the fact that the photon mode is coupled to all particles and introduces an effective long range interaction.

IV The Jaynes-Cummings-Hubbard Model

We now consider the JCH model given by Eq. (2) but ignoring the CR terms,

[TABLE]

In what follows, we set the energy scale by fixing $g=1$.

IV.1 Perturbation

We start in the $J\ll 1$ limit where perturbation in $J$ can be expected to give accurate results.

In the $J=0$ limit, the model can be solved exactly since the eigenstates are dressed excitons labeled by the exciton number, $n$, and upper or lower branch, $\pm$, which can be written as a superposition of a Fock state with $n$ photons plus atomic ground state, $|n,g\rangle$, and a Fock state with $n-2$ photons plus atomic excited state, $|n-2,e\rangle$,

[TABLE]

with the angle $\tan\theta_{n}=2g\sqrt{n(n-1)}/(\Delta+R_{n})$, $R_{n}(\Delta)=\sqrt{4g^{2}n(n-1)+\Delta^{2}}$, and the detuning parameter $\Delta=\omega_{q}-2\omega$.

The corresponding eigenvalues are

[TABLE]

The zero-exciton state $|0-\rangle$ is a special case with $E_{0-}=0$. The energy difference $R_{n}(\Delta)$ defines the Rabi splitting of the upper and the lower branches. By comparing the energy of neighboring exciton number $n=2$ , we obtain the border $\omega\omega_{q}<1$ between vacuum and the $n=2$ Mott state. Leaving the $n=2$ Mott state or the vacuum at $\omega_{q}>1$ by increasing exciton number, renders the Hamiltonian unbounded from below and makes the system unstable at $\omega=1$. The $J=0$ phase diagram is shown in Fig. 4.

At small nonzero tunneling, $J\ll 1$, we ignore the upper exciton branch and perform a perturbation calculation to second order in $J$ to obtain the phase diagram of the JCH model. The ground state energy of the Mott phase, $E_{M}$, that of the state with one additional exciton, $E_{d}$, and that of the state with two additional excitons $E_{2d}$ are calculated. We find that the energy at the line determined by $E_{M}=E_{d}$ is smaller than $E_{2d}$ at the same $\omega$. This tells us that the upper Mott boundary is determined by doping with a single exciton (adding one photon). The boundary is given by,

[TABLE]

where the matrix elements

[TABLE]

and $\sigma,\nu=\pm$, can be expressed in terms of the angle $\theta_{n}$ and exciton number $n$. See Appendix C for details.

To determine the lower Mott boundary, we remove one exciton (single holon doping) and calculate $E_{h}$, and also remove two excitons (double holon doping) and calculate $E_{2h}$. The energy $E_{h}$ is always higher than $E_{2h}$ at $\omega$ determined by $E_{M}=E_{2h}$ indicating that the lower boundary of the Mott phase is given by $E_{M}=E_{2h}$,

[TABLE]

It is very interesting that, whereas single exciton doping determines the upper boundary, two excitons should be removed to determine the lower Mott boundary. This is confirmed by numerical calculations and has other consequences which will be discussed below. The natures of the phases above and below the Mott lobe will be discussed below.

The boundary between the vacuum and a state with at least two excitons is determined by the equation $E_{2e}=E_{0-}=0$, where $E_{2e}$ is the energy of the 2-exciton state:

[TABLE]

The resulting phase diagram, with the choice $\omega_{q}=0.4$, is shown in Fig. 5. Details of our perturbation calculation are shown in Appendix C. In addition, Fig. 5 shows the mean field boundary (see Appendix B for details), Eq. (56), between stable thermodynamic phases and the unstable region where the Hamiltonian is unbounded from below. For comparison, we also show in Fig. 5 the phase boundaries obtained with DMRG. As expected, we see very good agreement between perturbation and DMRG for small values $J$.

IV.2 Numerical results

To obtain the phase diagram with exact numerical methods, we employ the SGFrousseau08 ; rousseau08b quantum Monte Carlo method and DMRGwhite92 ; white93 using the ALPSalps library. The SGF and DMRG methods offer complementary advantages and allow us to map out the phase diagram more precisely. Some typical running times for the SGF QMC simulations are: $216$ hours for $N=24$ at $\beta=768$ and $192$ hours for $N=20$ at $\beta=640$. These long running times and large values of $\beta$ are necessary to ensure convergence to the ground state. As for DMRG, we typically took the maximum number of photons/site, $N_{max}=5$, $500$ states, and $500$ sweeps. These values are typical but the details depend on the couplings and the phase of the system. In all cases, we verified that increasing these values did not change the results.

As in the previous section, here we set the energy scale by fixing $g=1$ and take $\omega_{q}=0.4$. We start with grand canonical QMC simulations using the SGF algorithm. Figure 6 shows the density of excitons as a function of the photon frequency, $\omega$, for several values of the inverse temperature $\beta$. This allows us to determine how large $\beta$ needs to be to probe the ground state. Note that in the grand canonical ensemble, $\omega$ is minus the photon chemical potential. We see a clear incompressible Mott insulator (MI) plateau surrounded by compressible regions. We note that the MI becomes fully formed only for very large $\beta$. This is due to the small value of $J$ we chose for the figure; larger values of $J$ do not require such high values of $\beta$, but, of course, at larger $J$, the MI is not as wide (see Fig. 5).

To delineate the nature of the compressible regions surrounding the MI, we show in Figs. 7 and 8 the exciton density as a function of $\omega$ in the top panel, and the superfluid density, $\rho_{s}$, in the bottom panel. We see in both figures that, as soon as the system leaves the MI phase, it becomes superfluid. Hole doping, Fig. 7, results in a very dilute SF before the system eventually becomes empty as $\omega\approx 2.8$. Particle doping, Fig. 8, leads to a SF with $\rho_{s}$ approaching $2$ before the system starts to become unstable at $\omega\approx 2.23$. By making several cuts of this type, we map out the phase diagram shown in Fig. 9 obtained numerically with very good agreement between SGF and DMRG. The figure shows a single Mott lobe with $2$ excitons/site; below this lobe there is a rather narrow region of SF before the system becomes empty. Above the Mott lobe, there is an even narrower strip of SF before the system becomes unstable. It is hard to pinpoint precisely with SGF and DMRG the boundary where the system becomes unstable because the number of photons/site increases very rapidly near the instability. Neither DMRG nor QMC performs well under such conditions and, for that reason, we show in the figure the stable/unstable boundary given by MF calculations. Unlike the case of the 1-photon JCMrossini07 , there are no higher Mott lobes in this 2-photon case due to the instability triggered by the Hamiltonian becoming unbounded from below. The Mott lobe terminates in a cusp because the critical point at the tip is in the BKT universality class.

The question arises as to the nature of the SF phase. In the 1-photon systemrossini07 , the photon Green function, $\langle{\hat{a}}^{\dagger}_{0}{\hat{a}}^{\phantom{\dagger}}_{r}\rangle$ decays as a power indicating quasi-long range order and a photon SF. At the same time, the photon-qubit Green function, $\langle\sigma^{+}_{0}{\hat{a}}^{\phantom{\dagger}}_{r}\rangle$, also decays as a power law. We see the same behavior in Fig. 10 where several two-point functions are shown on log-log scale. We note, in particular, that $\langle{\hat{a}}^{\dagger}_{0}{\hat{a}}^{\phantom{\dagger}}_{r}\rangle$ decays as a power law, and that, while both $\langle{\hat{a}}_{0}^{\dagger 2}{\hat{a}}_{r}^{2}\rangle$ and $\langle\sigma^{+}_{0}{\hat{a}}^{2}_{r}\rangle$ decay as powers, they are not equal. The leading effect in this phase is then the coherent movement of individual photons.

One expects to encounter similar behavior below the Mott lobe. However, this expectation is not borne out by the numerical results. Figure 11 shows the same correlation functions as Fig. 10 where we see clearly that, although the system is in the SF phase, the quantity $\langle{\hat{a}}_{0}^{\dagger}{\hat{a}}_{r}\rangle$ decays exponentially while it exhibited power law behavior above the MI. The correlation functions $\langle\sigma^{+}_{0}\sigma^{-}_{r}\rangle$ and $\langle{\hat{a}}^{\dagger 2}_{0}{\hat{a}}^{2}_{r}\rangle$ both decay as powers showing the SF nature of this phase. In addition, we find that, unlike the case above the MI (see Fig. 10), the correlation function $\langle\sigma^{+}_{0}{\hat{a}}^{2}_{r}\rangle$ is identical to $\langle{\hat{a}}^{\dagger 2}_{0}{\hat{a}}^{2}_{r}\rangle$, and so we do not show it in Fig.11 to keep the figure clear.

This behavior appears to show that above the MI, we have a photon SF phase (power law decay for $\langle{\hat{a}}^{\dagger}_{0}{\hat{a}}^{\phantom{\dagger}}_{r}\rangle$) whereas below the MI we have a photon pair SF phase (exponential decay of $\langle{\hat{a}}^{\dagger}_{0}{\hat{a}}^{\phantom{\dagger}}_{r}\rangle$ and power law decay of $\langle{\hat{a}}^{\dagger 2}_{0}{\hat{a}}^{2}_{r}\rangle$): Below the Mott insulator the photons appear to form bound pairs which become superfluid.

V The Rabi-Hubbard model

In Appendix A we show the mean field calculation for the Rabi-Hubbard model predicting a phase transition between a disordered (normal) phase and a superradiant phase similar to the Dicke modelgarbe17 . The boundary between the two phases is given by (see Eq. (48))

[TABLE]

where $E_{G}$ is the ground state energy given by (see Eq. (48)

[TABLE]

In addition, the stability condition requires (see Eq. (49),

[TABLE]

The two phases and the unstable region are shown in Fig. 12. This phase diagram is the analog of that of the Dicke model, Fig. 3. However, whereas in the case of the Dicke model, we found excellent agreement between the MF and QMC phase diagrams (Fig. 12), this is not the case for the RH model.

We show in Fig. 13 the average number of excited qubits, $\langle N_{q}\rangle/N$, (top) and the average photon density, $\langle N_{ph}\rangle/N$, (bottom) as functions of the coupling $g$ for several values of $\beta$. We observe a smooth continuous increase of both quantities as $g$ increases and no finite temperature effects. The values are in good agreement with the mean-field approach described in Appendix B, which, for the parameters of Fig. 13, predicts a direct transition from a normal phase to the unbounded region. On the other hand, since the same mean-field approach applied to the JCH Hamiltonian is not able to capture the photon pair SF phase, it is possible that it also misses an SF phase in the RH model for parameters different from shown in the figure. In Fig. 14, we show that for sufficiently large $\beta$ there is no finite size effects. Comparing the top panels of Figs. 13 and 14 with Fig. 1 emphasizes the difference in behavior between the two models and suggests the absence of a phase transition for the RH model, at least for the present parameter values. For the entire range of values of $g$ in Figs.13 and 14, all Green functions, such as $\langle{\hat{a}}^{\dagger}_{0}{\hat{a}}_{r}\rangle$ and $\langle\sigma^{+}_{0}{\hat{a}}^{2}_{r}\rangle$, decay exponentially indicating a normal rather than a superradiant phase.

In short, at least for the parameter values we have considered, the Rabi-Hubbard model seems not to exhibit a phase transition, even though the underlying discrete $Z_{4}$ symmetry allows it. In addition, we have two different mean-field approaches leading to rather different phase diagrams, which seems to indicate that the physics of the Rabi-Hubbard model might be more involved than the one of the Dicke and the JC models. A thorough study is beyond the scope of the present paper and will be presented elsewhere.

VI Conclusions

In this paper we studied the ground state properties of the Dicke, Jaynes-Cummings Hubbard, and Rabi-Hubbard models. We used mean field, perturbation, QMC, and DMRG calculations to elucidate the phase diagrams and the transitions between the various phases. We found that for the Dicke model, exact QMC results agree very well with the mean field phase diagramgarbe17 but that very large $\beta$ was required to probe the ground state properties. For the JCH model, we found that at small hopping, $J$, the system exhibits a single incompressible Mott insulator with two excitons/site. Doping above the MI, the system exhibits a photon SF phase before becoming unstable due to the Hamiltonian becoming unbounded from below. To dope the system above the MI, we add one photon at a time. This is in big contrast to the situation below the MI where we found a surprising photon pair SF: The photons pair up in bound states. Below the MI, one needs to remove two photons at a time to dope. For the JCH model in the SF phase, very large values of $\beta$ are needed. Interestingly, we have not found a phase transition between a normal and a superradiant phase in the case of the RH model; only a normal phase with exponential decay of the Green functions. It would be interesting to study how the three-photon model would differ from what we found here. In particular,would the JCH model still exhibit a MI phase, and would doping below it yield a three-photon bound state superfluid.

Acknowledgements.

S.C and W.G. was supported by the NSFC under Grant No. 11775021 and No. 11734002.

Appendix A Mean field two-photon RH

In this appendix we outline the mean field calculation of the phase diagram of the two-photon Rabi-Hubbard model which proceeds along the same lines as that for the Dicke model in garbe17, . The RH Hamiltonian is:

[TABLE]

We apply the Holstein-Primakoff transformation:

[TABLE]

with $[{\hat{b}}^{\phantom{\dagger}}_{i},{\hat{b}}^{\dagger}_{j}]=\delta_{ij}$, use the mean field approximation,

[TABLE]

and ignore fluctuations in the ${\hat{b}}$ field. This leads to the quadratic mean field Hamiltonian,

[TABLE]

Applying a Fourier transform,

[TABLE]

leads to,

[TABLE]

Now, we use the Bogoliubov transformation:

[TABLE]

To simplify the expressions, we write: $c_{k}\equiv{\rm cos}(k)$, $ch_{k}\equiv{\rm cosh}(X_{k})$, and $sh_{k}\equiv{\rm sinh}(X_{k})$. This yields,

[TABLE]

The first term in Eq. (39) gives the quasi-particle excitations above the ground state; the second term contains the nondiagonal terms which we want to eliminate and the third term gives the ground state energy. When the coefficient of the nondiagonal term is made to vanish, we obtain the condition:

[TABLE]

and using

[TABLE]

we get:

[TABLE]

Writing

[TABLE]

gives:

[TABLE]

and,

[TABLE]

The acceptable values satisfy the condition:

[TABLE]

Then, the energy that must be minimized with respect to ${\bar{b}}$ is

[TABLE]

Note that, putting $J=0$, renders all the hyperbolic functions independent of $k$. Then Eq. (47) agrees with the MF equation in Ref. garbe17, if we take $N=1$ and ${\bar{b}}={\bar{b}}^{*}$.

Equation (47) can be further simplified to the form,

[TABLE]

The requirement that the argument of the square root be positive leads to the important condition:

[TABLE]

This is the stability condition for the system: When it is violated, the Hamiltonian is unbounded from below.

The transition to the superradiant phase is indicated by a nonzero order parameter of the order parameter, ${\bar{b}}\neq 0$. The boundary between the normal and superradiant phases is, therefore, given by:

[TABLE]

This relation cannot be solved analytically, we solve it numerically to get the relation between $g,J,\omega,\omega_{q}$ which gives the phase boundary.

Appendix B Simple mean field for JCH

The mean field calculation in Appendix A applies to the Rabi-Hubbard model and, by ignoring the counter-rotating terms, also to the JCH model. However, we can also perform the mean field calculation on the hopping term of the photon field rather than on the qubit. We assume $\langle a_{i}\rangle=\psi_{i}$, which reduces the Hamiltonian to a single-site form:

[TABLE]

Finding the ground state of this Hamiltonian and a self-consistent solution for the $\psi_{i}$ allows us to obtain the phase diagram. In particular, in the normal phase $\psi_{i}=0$, the ground state of $H_{i}^{MF}$ still depends on the parameters $\omega_{q}$, $\omega$ and $g$, leading, for instance, to the dependence of the qubit and photon densities on the interaction strength observed in Fig.13.

A further approximation amounts to replacing all the photonic operators by the mean field value, $\psi$, thus yielding the leading order contribution which gives the instability line for the JCH model. In that case, assuming that $\psi$ is real, the mean field Hamiltonian becomes

[TABLE]

The ground state energy is found by minimizing

[TABLE]

with respect to $\psi$ which gives

[TABLE]

Near the unbounded region, $\psi\gg 1$ and we can thus ignore the $\omega_{q}$ term and find that for stability we must have,

[TABLE]

The boundary of the unstable region is therefore given by,

[TABLE]

Note the similarity between this equation and the stability condition for the RHM, Eq. (49).

We then define the photon density $n_{ph}=\psi^{2}$ and, using Eq. (54) we obtain,

[TABLE]

Now we consider $\omega$ near the stable-unstable boundary, $\omega_{c}$, so that $\omega-2J\approx g$. We find that the number of photons diverges as the unstable region is approached:

[TABLE]

with $\delta\omega=\omega_{c}-\omega$. This divergence makes the result difficult to demonstrate with QMC and/or DMRG methods.

Appendix C Perturbation calculation of the phase diagram of the JCH model

In this appendix we outline the perturbation calculation of the phase diagram of the two-photon JCH model which proceeds along the same lines as that for the single-photon JCH model in schmidt09, .

The Hamiltonian of the two-photon JCH model is

[TABLE]

which can be split into two parts,

[TABLE]

with the perturbation, $H_{1}$, given by

[TABLE]

and $H_{0}$ is the rest of $H$.

It is convenient to introduce the matrix elements

[TABLE]

with $\sigma,\nu=\pm$. Using Eq. (19), we have

[TABLE]

with

[TABLE]

where

[TABLE]

For the upper boundary of the $n=2$ Mott phase, we compare the ground state energy of the $n=2$ Mott state with that of the state doped by one exciton. However, to determine the lower boundary of the $n=2$ Mott phase, we need to compare the ground state energy with the state doped by 2 holons (with two excitons removed from the system). Similarly, to locate the boundary of vacuum, we need to compare the energy of the state doped by two excitons with vacuum.

We first calculate the energy of the Mott state to second order,

[TABLE]

with

[TABLE]

where the zero-th order wavefunction is

[TABLE]

and the first order wavefunction is

[TABLE]

Now we calculate the lowest energy of the state obtained by doping the $n=2$ Mott state, to second order,

[TABLE]

with

[TABLE]

where

[TABLE]

and

[TABLE]

with

[TABLE]

The equation $E_{M}=E_{d}$ leads to the upper boundary Eq. (21) of the $n=2$ Mott phase. We mention that we have also calculated the ground state energy $E_{2d}$ of the $n=2$ Mott state doped by two excitons and find that the energy at the boundary determined by $E_{M}=E_{d}$ is lower than that at the line determined by $E_{M}=E_{2d}$. The former also matches the DMRG and the SGF results well.

To find the lower boundary of the $n=2$ Mott lobe, we calculate the ground energy of the state with two excitons (holons) removed from (added to) the $n=2$ Mott state:

[TABLE]

with

[TABLE]

where the wave functions are

[TABLE]

The equation $E_{M}=E_{2h}$ leads to the lower boundary, Eq. (23) of the $n=2$ Mott phase. Note that for the lower boundary, it is necessary to remove two excitons to obtain the lowest energy state whereas for the upper boundary, we add only one exciton. This is confirmed by numerical calculations and leads to consequences discussed in the text.

To find the boundary of the vacuum, we need to find the energy of the state with only 2 excitons.

[TABLE]

with

[TABLE]

where the wave functions are

[TABLE]

The vacuum boundary, Eq. (24), is given by the equation $E_{2e}=0$.