Nonlinear two-photon Rabi-Hubbard model: superradiance and photon/photon-pair Bose-Einstein condensate

We study the ground state phase diagram of a nonlinear two-photon Rabi-Hubbard (RH) model in one dimension using quantum Monte Carlo (QMC) simulations and density matrix renormalization group (DMRG) calculations. Our model includes a nonlinear photon-photon interaction term. Absent this term, the RH model has only one phase, the normal disordered phase, and suffers from spectral collapse at larger values of the photon-qubit interaction or inter-cavity photon hopping. The photon-photon interaction, no matter how small, stabilizes the system which now exhibits {\it two} quantum phase transitions: Normal phase to {\it photon pair} superfluid (PSF) transition and PSF to single particle superfluid (SPSF). The discrete $Z_4$ symmetry of the Hamiltonian spontaneously breaks in two stages: First it breaks partially as the system enters the PSF and then completely breaks when the system finally enters the SPSF phase. We show detailed numerical results supporting this, and map out the ground state phase diagram.


I. INTRODUCTION
Light-matter interaction is ubiquitous in nature and is, therefore, the focus of much theoretical and experimental work. The relatively simple Rabi model [1,2] gives a good description of the interaction of photons with a two level system (qubit) and recent advances in manipulating the interactions of photons with single atoms have resulted in its experimental realization using twolevel atoms in cavities [3] (cavity quantum electrodynamics, QED) or Josephson junctions on solid state chips [4][5][6][7][8] (circuit QED). In the strong coupling limit, g/ω 0.1, one can exploit the random wave approximation (RWA) to justify ignoring the counter-rotating (CR) terms resulting in a Hamiltonian which conserves the exciton number (number of photons plus excited qubits) due to its U (1) symmetry. This results in the Jaynes-Cummings (JC) model [9] which can be solved exactly. Such cavities can be arranged in a chain where the coupling between near neighbors can be controlled via a tunable tunneling of the photon mode thus resulting in the Jaynes-Cummings-Hubbard model [10] which describes itinerant bosons (the photons) interacting with localized qubits. This model has been shown to behave like the one dimensional Bose-Hubbard model (BHM) [11] with a superfluid phase (of excitons) and incompressible Mott insulator (MI) lobes [12][13][14][15][16][17][18][19]. The MI is approximately a product of single site states obtained from a superposition of photons and excited atoms [18] and exhibits behavior similar * waguo@bnu.edu.cn † george.batrouni@inphyni.cnrs.fr to photon blockade [20] where there is a finite energy gap opposing the addition of a photon.
In the experimentally attainable [21][22][23][24][25][26], ultra-strong coupling regime, g/ω ∼ 1, the CR terms can no longer be ignored. Their restoration to the Hamiltonian reduces the U (1) symmetry to Z 2 . Consequently, the number of excitons is no longer conserved thus removing the possibility of a MI phase. The phase diagram now consists of a disordered phase and an ordered (coherent) one separated by a quantum phase transition in the universality class of the two-dimensional Ising model [10,[27][28][29]. The ordered phase, therefore, exhibits a Bose-Einstein condensate (BEC) of photon. This transition resembles the incoherent/coherent (normal/superradiant) phase transition in the Dicke model [30,31].
The two-photon Rabi model has also attracted much attention as new systems are realized where multi-photon processes come into play. For example, it has been used to describe second order processes in Rydberg atoms in cavities [32] and quantum dots [33,34], and mechanisms have been proposed to realize it in circuit QED [35]. The two-photon model undergoes spectral collapse where the Hamiltonian is no longer bounded when the coupling exceeds a certain value [36][37][38][39][40][41][42][43]. In the strong coupling regime, the CR terms can be ignored, as in the onephoton case, and result in the two-photon JC model with U (1) symmetry and conserved number of excitons. In the ultra-strong regime, the CR terms are restored and the symmetry is reduced to Z 4 . The ground state of the many-body two-photon model was studied in the context of the Dicke model [44] using mean field [45] and QMC [46]. Excellent agreement was found between the two methods which showed the system to exhibit a quantum phase transition between a normal (disordered) and superradiant phase. Here, however, the term superradiant is used to indicate a macroscopic change in the average number of photons but which remains relatively small. This is in contrast with the one-photon case where there is a very large number of photons in the superradiant phase which form a BEC [10,[27][28][29]. QMC simulations were also used to examine the two-photon JCH and Rabi Hubbard (RH) models [46]. It was found that the JCH model has only one MI lobe with two excitons/cavity, unlike the onephoton case where there is a succession of MI lobes [12][13][14][15][16][17][18][19]. Furthermore, it was found that there are two superfluid (SF) phases; a single photon SF phase (SPSF) and a photon pair SF phase (PSF). In the former, the single particle Green function decays as a power indicating the quasi-long range order for the photons. In the latter case, however, the single particle Green function decays exponentially while the photon pair Green function decays as a power indicating quasi long range order for bound photon pairs but not for simgle photons. When the CR terms are restored, the symmetry of the Hamiltonian is reduced from U (1) to Z 4 and it was found that the two-photon RH model does not exhibit any quantum phase transitions; it has only the disordered phase and the spectral collapse region [46].
In this paper, we examine the possibility that nonlinear photon terms in the Hamiltonian could stabilize the system and allow the appearance of quantum phase transitions in the RH model. Non-linear terms, i.e. photon-photon interactions, have attracted experimental and theoretical interest as a means to generate topological photon pairs [47][48][49] which have robust transport properties. To this end we exploit the density matrix renormalization group [50,51] (DMRG) with open boundary conditions (OBC) using the ALPS library [52], and the stochastic Green function [53,54] (SGF) quantum Monte Carlo method, with periodic boundary conditions (PBC), to study the phase diagram of the two-photon RH model with a nonlinear photon term. We show that the photon-photon interaction term, stabilizes the system, eliminating spectral collapse, and leads to the appearance of two quantum phase transitions, the first is the transition from the disordered phase to PSF phase where the Z 4 symmetry patially breaks, the second is the transition between the PSF and the SPSF phases where the Z 4 symmetry is completely broken. We show that in the PSF phase, we have a BEC of photon pairs whereas in the SPSF phase we have a BEC of photons.
The paper is organized as follows. In section II we present the model and discuss briefly the methods used to perform the numerical calculations and simulations. In section III we present and sicuss our results, and in section IV we discuss our conclusions.

II. MODEL AND METHOD
The one-dimensional nonlinear two-photon Rabi-Hubbard (RH) model we study is governed by the Hamiltonian, where N is the number of sites (or cavities) andâ i (â † i ) is the photon destruction (creation) operator in the ith cavity. Photon tunneling between sites is governed by the hopping parameter J. The photon frequency is ω and the qubit energy spacing is ω q ; σ z i and σ x i = σ + i + σ − i are the Pauli matrices acting on the ith qubit, σ + i (σ − i ) is the corresponding raising (lowering) operator. The last term in Eq.(1) describes the onsite photon-photon interaction. Ignoring the CR terms in Eq.
i ), leads to the Jaynes-Cummings-Hubbard (JCH) model where the system is invariant under the generalized rotation operator, for any θ, thus exhibiting U (1) symmetry and conservation of the number of excitons N exc = N ph + 2N q where N q is the number of qubits in the excited state and N ph is the number pf photons. However, the CR terms break the U (1) symmetry, they pick up a phase exp(i4θ) due to the action of the operator R(θ); this reduces the symmetry to Z 4 . Now the system is left invariant by this rotation only for θ = n2π/4, n = 0, 1, 2, 3 and N exc is no longer conserved. The discrete Z 4 symmetry can break spontaneously in the ground state leading to an ordered BEC phase.
To characterize the various possible phases, we calculate several Green functions, where α and β denote creation and annihilation operators of the photons (â † i andâ i ) or the qubits (σ + j and σ − j ). O denotes the ground state expectation value, GS|O|GS , for DMRG, and the statistical average for QMC, where Z = Tr e −βH O and β = 1/T . For example, the one-body photon Green function at equal time is given by, The qubit Green function is, and the following two functions will be particularly useful: and The photon pair Green function is given by, The Green functions described by Eqs. (3)(4)(5)(6)(7)(8)(9) are defined on a lattice with perdiodic boundary conditions used in the QMC simulations. When using DMRG with open boundaries, we choose the origin of the correlation function in the center of the lattice and calculate the correlations from that point. Power law decay of one of these Green functions would indicate quasi-long range order for the corresponding quantity. If the Green function decays to a finite constant, it signals long range order and the spontaneous breaking of the corresponding symmetry. The Fourier transform of the single particle Green function, Eq.(5), gives the photon momentum distribution, n ph (k); the Fourier transform of the pair Green function gives the pair momentum distribution, n pair (k); the Fourier transform of the qubit Green function, Eq. (6), gives the qubit momentum distribution, n q (k) . These quantities will indicate whether a condensate (i.e. long range order) is present. We also measure the average number of excitons, We calculate these quantities using the ALPS DMRG package [52] and the SGF QMC algorithm [53,54]. For DMRG, we verified that, in all cases, the number of states we kept and sweeps we performed were sufficient for proper convergence (up to 240 states and 300 sweeps for large systems).

III. RESULTS
The results we present here were obtained using DMRG and SGF QMC depending on the physical quantities being studied and the boundary conditions: Open boundary conditions were typically used with DMRG and periodic ones with SGF. We begin, therefore, by comparing the two methods to ensure that results yielded by one can be reliably compared with the other. Figure 1 compares DMRG and SGF results, both with periodic boundary conditions, for several Green functions in the PSF phase (see below) versus distance. The finite temperature SGF QMC simulations were performed at very low temperature, β = 360, and show excellent agreement with the zero-temperature DMRG calculations. In what follows, we take U = 1 to fix the energy scale. We also choose the photon frequency ω = 1 and the qubit energy spacing ω q = 0.1 and determine the phase diagram of the system governed by Eq.(1) in the (J, g) plane. It was seen in Refs. [45,46] that the average number of excited qubits per site, n q = N q /L, is a good probe for the phase transition between the normal (incoherent) phase and the superradiant phase. Here, we will use both n q and n ph = N ph /L (the average number of photons per site) as probes for possible phase transitions. We use Padé approximants [46] to fit our numerical results for n q and n ph as functions of g or J. Differentiating the Padé forms, and also numerically differentiating the DMRG data, reveal the presence of peaks indicating the location of the phase transitions. The analytical Padé and the numerical derivatives are in excellent agreement. Figure 2 shows DMRG results for n ph and n q and the corresponding numerical and Padé derivatives, dn ph /dg and dn q /dg, versus g for constant J = 0.4 and for L = 12, 16, 20, 24. In Fig.2(a) we show n ph versus g at fixed J = 0.4, and in (b) we show dn ph /dg which, remarkably, displays two peaks indicating the presence of two quantum phase transitions. n q and its derivative are shown in Fig.2(c) and (d) respectively. It is seen that dn q /dg exhibits only one peak which matches the first peak observed for dn ph /dg. We conclude, therefore, that the qubits have one type of potential order while the photons have two. This will be examined below. It is seen that the positions of the peaks shift to smaller g values as L increases; this shift is used to extrapolate the location of the transition to the thermodynamic limit and construct the phase diagram. Similar behavior is observed when g is fixed and and J is changed, as shown in Figs.3. We conclude, therefore, that typically, for this system, two quantum phase transitions are encountered as g (J) is fixed and J (g) is varied. We recall that for the two-photon Rabi model in the absence of the quartic term, no phase transition is present [46]. The addition of the quartic term, U i (a † i a i ) 2 , stabilized the system and allowed the appearance of quantum phase transitions.
To determine the nature of the three phases present, we study the behavior of the Green functions, Eqs. (3)(4)(5)(6)(7)(8)(9). Figure 4 shows four of these functions, obtained with DMRG, on log-log scale clearly exhibiting decay faster than power (exponential) with distance and, therefore, corresponding to the normal (incoherent) phase. This behavior is characteristic of the region in parameter space, (J, g), where J and/or g are small, i.e. before the first peak in Figs.2 and 3 leading to the conclusion that this region corresponds to a disordered (incoherent) phase.
Choosing, for example, (J, g) = (0.6, 1.02) puts the system between the two peaks observed in Figs.2 or 3. We see in Fig.5 that, in this case, the one-body photon Green function, Eq.(5), still decays exponentially whereas all other Green functions shown, rapidly saturate to constant values. The exponential decay of G a † ,a (r) indicates the absence of off-diagonal long range or quasi-long range photon order and, consequently, no photon SF, and â = 0. However, the saturation (constant value) of the other Green functions as the distance increases, indicates the establishment of true off-diagonal long range order (ODLRO) for photon pairs. The saturation value is the square of the order parameter, e.g. G saturated a † 2 ,a 2 ∼ | a 2 | 2 . Therefore, between the two peaks of Figs.2 and 3, the system is in a photon pair SF (PSF) phase with a pair BEC. In terms of symmetry breaking, this means that the Z 4 symmetry of the Hamiltonian, Eq.(1), generated by the operator Eq.(2), is only partially broken: â 2 = 0, and â = 0. Increasing g and/or J further, puts the system on the right of the second peaks in Figs.2 and 3. We see in Fig.6 that now all Green functions saturate to constant values as the distance increases. This shows that now the system is in a phase with true off-diagonal long range photon order indicating the complete breaking of the Z 4 symmetry and the establishment of a single photon SF. This shows that the Z4 symmetry is completely broken and the system is in the single particle SF phase with a photon BEC.
To elucidate further the nature of the two SF phases, we examine the behavior of the qubit, the one-body and the pair momentum distributions. In Fig.7 we show the normalized photon momentum distribution, n ph (k), in the three phases we have identified and for several system sizes. We see that, as the system size increases, n ph (k = 0) decreases when the system is in the normal or PSF phases indicating the absence of a photon BEC and, therefore, the absence of (ODLRO), â i = 0. However, n ph (k = 0) remains constant in the SPSF phase indicating that here the photons have formed a BEC, â i = 0, and the Z 4 symmetry is broken. On the other hand, Fig.8 shows the pair momentum distribution, n pair (k), for the same parameters as in Fig.7. Here we see that in the normal phase, n pair (k = 0) decreases as L increases but remains constant in the PSF phase. This means that in this phase, bound photon pairs have formed a BEC even though the photons themselves have not, and, therefore, â 2 i = 0 but â i = 0. In the SPSF phase, we already saw, Fig.7, that n ph (k = 0) = 0, and, therefore, it is not surprising that n pair (k = 0) is also (trivially) nonzero in this phase. Figure 9 shows the momentum distribution of the qubits. In the disordered phase, n q (k = 0) decreases as L increases showing that there is no order. In both the PSF and SPSF, n q (k) behaves in the same way: n q (k = 0) remains constant as L increases showing that the qubits are ordered. In fact, when the qubits order at the transition between the normal and PSF phases, they remain ordered in the same manner as the system transitions from PSF to SPSF. This is manifested by appearance of only one peak in the derivatives of n q as seen above.
To map out the phase diagram, we perform simulations along several lines of constant g while varying J and also constant J while varying g. As shown above, the two peaks exhibited by dn ph /dg and by dn ph /dJ indicate the location of the quantum phase transition for the system size being studied. We do the calculations for several system sizes (L = 12, 16,20,24) and extrapolate is the critical g for the transition between the disordered and PSF phases, and g (2) c is the critical g for the transition between the PSF and SPSF. δg → 0 as J increases. Similar behavior is seen for dn ph /dJ as g increases. The data were obtained with DMRG.
the critical values to the thermodynamic limit. The separation between the two peaks in the derivatives of n ph is not constant; it decreases as g or J get larger. Figure  10 illustrates this for several values of J. The inset shows the separation between the two peaks, δg ≡ g (2) c − g (1) c , where g (1) c is the critical g for the transition between the disordered and PSF phases, and g the transition between the PSF and SPSF. We see that δg gets smaller as J increases. We observe similar behavior for dn ph /dJ as J increases. Putting all these results together leads to the phase diagram we show in Fig.11. The figure shows the three phases discussed above: Normal incoherent phase, superradiant photon condensate phase (single particle SF, SPSF) and sandwiched in between is the pair SF phase (PSF). As the two transitions approach each other, for example at large J, they become hard to distinguish and appear to merge eventually into one single transition leading to a direct passage from the normal phase to the SPSF phase without passing first through the PSF phase. We also show in Fig.11 the Gutzwiller mean field result for the boundary between the SPSF and the other phases (see appendix for details).

IV. CONCLUSIONS
A lot of attention has been given to the two-photon Rabi model, its applications [32][33][34][35] and the spectral collapse [36][37][38][39][40][41][42][43] it undergoes in certain regions of its parameter space. Interestingly, in contradistinction with the one-photon Rabi-Hubbard model which, in its ground state, exhibits a quantum phase transtion from a disordered phase to a superradiant one, chracterized by a spontaneously broken symmetry (Z 2 ) and a photon BEC, the two-photon RH model only exhibits the disordered phase. In this work, we used exact computational methods (DMRG and QMC) and showed that the system can be stabilized by a nonlinear (quartic) term which models effective photon-photon interactions [47][48][49]. This sta-bilization eliminates the spectral collapse of the model and, in fact, exposes two quantum phase transtions in the ground state. The first transition, a consquence of partial spontaneous symmetry breaking, takes the system from the disordered phase to the pair superfluid (PSF) phase which is characterized by a nonvanishing pair condensate order parameter, a 2 = 0, while at the same time a = 0. The second transition is from the PSF to the single particle SF (SPSF) phase and completes the symmetry breaking with the photon condensate order parameterm acquiring a nonvanishing expectation value, a = 0. An interesting question to ask is what happens in higher order photon processes, for example the threephoton model? Without stabilization, the three-photon model undergoes spectral collapse but, with a nonlinear photon-photon term like we used here, it will be stabilized. It would be interesting to study this model both in its Jaynes-Cummings limit (ignoring the counter-rotating terms) where the symmetry is U (1) and will not break in one dimension, and also in its full Rabi form where the CR terms are kept and where the symmetry is now Z 6 . Specifically, will there be two or more kinds of SF phases? For example, will there be photon triplet SF and BEC where three photons act as a single boson that condenses? Will the Z 6 symmetry break down in more stages than the smaller Z 4 and, consequently, result in more phases? photons in the cavity; g (e) denotes a qubit in the ground (exicted) state. The wave function can then be written in terms of this basis, |Ψ = nmax i=0, k=(g,e) c k,i |i, k . (A3) The order parameter, ψ = Ψ| a|Ψ becomes, ψ = c g,0 * c g,1 + √ 2c g,1 * c g,2 . . . + √ n max c g,nmax−1 * c g,nmax + c e,0 * c e,1 + √ 2c e,1 * c e,2 + . . . + √ n max c e,nmax−1 * c e,nmax (A4) The coefficients are first chosen randomly, the Hamiltonian matrix is calculated and diagonalized. This gives a new estimate for the ground state wavefunction which is then used iteratively to calculate an improved ground state wavefunction and so on until the process converges. This way, we calculate ψ for a chosen fixed value of U and many values of ω, ω q , J, g, and determine the region in phase space where ψ = 0. For example, for U = 1, we obtain the dashed red line in Fig.11.