Quantum phase transitions of light in a dissipative Dicke-Bose-Hubbard model
Ren-Cun Wu Lei Tan Wen-Xuan Zhang Wu-Ming Liu





Institute of Theoretical Physics, Lanzhou University, Lanzhou 730000, China

tanlei@lzu.edu.cn
Institute of Theoretical Physics, Lanzhou University, Lanzhou 730000, China


Institute of Theoretical Physics, Lanzhou University, Lanzhou 730000, China

Beijing National Laboratory for Condensed Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190,China


The impacts that the environment has on the quantum phase transition of light in the Dicke-Bose-Hubbard model are investigated. Based on the quasibosonic approach, mean field theory and the perturbation theory, the formulation of the Hamiltonian, the eigenenergies and the superfluid order parameter are obtained analytically. Compared with the ideal cases, the order parameter of the system evolves with time as the photons naturally decay in their environment. When the system starts with the superfluid state, the dissipation makes the photons tend to localize, and a greater hopping energy of photon is required to restore the long-range phase coherence of the localized state of the system. Furthermore, the Mott lobes disappears and the system tends to be classical with the number of atoms increasing; however, the atomic number is far lower than that expected under ideal circumstances. Therefore, our theoretical results offer valuable insight into the quantum phase transition of a dissipative system.


42.50.Nn, 42.50.Pq, 05.70.Fh 




Quantum phase transitions of light in a dissipative Dicke-Bose-Hubbard model
    Wu-Ming Liu
    December 30, 2023
============================================================================







§ INTRODUCTION

Quantum simulation has become a research frontier and an indispensable tool in quantum information science<cit.>,  it's remarkable development in experiment realization has led to incredible advances in the field of quantum optics and atomic physics<cit.>. Among the recent developments, the system of coupled cavity arrays embedded with cold atoms has been intensively investigated as a platform to realize and simulate quantum many body phenomena because of its extremely high tunability, individual addressability and flexibility in its geometric design<cit.>. A wide range of condensed matter system have been theoretically investigated and  many proposals for probing them have been proposed including the quantum phase transition<cit.>, spin glasses<cit.>, photon crystals<cit.>, the emergence of gauge fields<cit.>, the quantum Hall effects<cit.>, the pfaffian-like topological state<cit.> and the supersolid <cit.>.

The simplest physical model of light-matter coupling in a coupled cavity array system is the Jaynes-Cummings Hubbard model, which presents an array of optical cavities that each contain a single two-level atom(TLA) in the photon-blockade regime<cit.>. A modified Jaynes-Cummings Hubbard model based on each cavity embedded a three-level atom has been proposed recently; this model circumvents the drawbacks of the excited-state spontaneous emission and provides a tunable extension of two-polariton bound states of the standard Jaynes-Cummings Hubbard <cit.>. As the number of atoms in each cavity increases, the collective effects due to atomic interactions among themselves give rise to intriguing many-body phenomena. In quantum optics the Dicke model is a paradigm of collective behavior<cit.> that describes the interaction of ensembles of TLAs that are collectively coupled to the single mode of radiation of a cavity<cit.>. Numerous investigations of interesting physical effects and their experiment realization<cit.>, such as the super-radiation phase<cit.>, the superradiant Mott Insulator<cit.> and the dynamical phase transition<cit.>, are discussed. Thus, a Dicke-Bose-Hubbard (DBH) model that includes more than one identical coupled cavities and  N identical TLAs in each cavity has been conducted to study the quantum phase transitions of light without considering the counter-rotating terms<cit.>. The transfer of excitations under a large range of operative conditions is also demonstrated and explored by tuning the controlling parameters in the DBH model<cit.>. Both the emergence of a polaritonic glassy phase<cit.> and the quantum phase transitions from the superfluid to the Bose-glass and the Mott-insulator states<cit.> have also been studied. Most recently, the localization-delocalization quantum phase transition of photons of the DBH model including counter-rotating terms has been presented<cit.>. The model shows that under the influence of the counter-rotating terms, the Mott lobes are fully suppressed.

As is well known, a realistic quantum optical system can rarely be isolated from its surrounding completely, particularly in an experiment; rather it is usually coupled to its external environment with an infinite number of degrees of freedom. To date, an investigation of quantum phase transition of photons in an dissipative DBH model is still lacking. To treat the interplay between the coupled cavity arrays and its environment in a more general setting, we developed a quasi-bosonic approach to describe the quantum phase transition and photon transport in an open quantum optical systems <cit.>. Without the requirement of considering the finite environment's degrees of freedom, the quasibosonic method is a great concept, that has a computational advantage. In the present paper, we use the quasibosonic approach to obtain an effective Hamiltonian of the dissipative DBH model. The coordinates of a bath can be eliminated and the system can be considered an ensemble of quasi-bosons in less time than its decay rate. Next, the eigenenergies and the superfluid order parameter of the system are also derived analytically for two TLAs on resonance, and we numerically demonstrate the phase diagram of an arbitrary number of TLAs. The theoretical analysis presented here will be an essential reference for future experiments to explore the quantum effects for multi-atom system.

The paper is organized as follows. In section 2, the dissipative DBH model is introduced based on the quasibosonic approach. Section 3 is devoted to deriving the eigenvalues and eigenstates for two atoms in each cavity. The analytical solutions of the superfluid order parameter for dressed states are given, and the properties of quantum phase transition are discussed in section 4. The extension to
an arbitrary number of TLAs is also given in this section. Section 5 gives the conclusion.





§ THE DISSIPATIVE DICKE-BOSE-HUBBARD MODEL


The system considered is depicted in Fig. 1. The Hamiltonian of the DBH model considering the coupling to its environment is give by (with ħ=1)<cit.>

    H   =   ∑_iH^DM_i-κ∑_ija^†_ia_j-μ∑_iN_i+H_R
    
    H^DM_i   =   ω_aJ^+_iJ^-_i+ω_ca^†_ia_i+β(a_iJ^+_i+a^†_iJ^-_i)
    
    H_R   =   ∑_i∑_ω_kω_kr_k^†r_k+H_CR+H_AR
    
    H_CR   =   ∑_i∑_ω_k[η_c(ω_k) r_k^†a+h.c] 
    
    H_AR   =   ∑_i∑_ω_k[η_a(ω_k) r_k^†J^-+h.c]

where the indices i and j are individual cavities range over all sites. a^†_i and a_i are the photon creation and annihilation operators, respectively. J^±_i=∑_jδ^±_j are the atomic collective raising and lowering angular momentum operators and the total number of excitations is N_i=a^†_ia_i+J^+_iJ^-_i. The transition energy of TLA is ω_a and ω_c is the frequency of the cavity field. All the atoms couple to cavities with the same coupling β<cit.>. We assume that the hopping energy of photons κ_ij=κ between sites i and j and the chemical potential in the grand canonical ensemble μ_i=μ are the same for all cavities. The coupling Hamiltonian of the system with the environment and the Hamiltonian of the environment are described as H_R. ω_k is a model of bath, and r_k^†(r_k) is the creation and annihilation operators of environment in the kth model. H_CR is the interaction of cavity with environment. The interaction of the atoms with the environment is denoted as H_AR.

Considering the influence of the environment, the decoherence of every cavity and the two-level atom would result in the incoherent or dissipative propagation of the incident photon, thus  nonequilibrium dynamics for the open quantum many-body system will arise.
In general, simulations of nonequilibrium many-body effects for a finite freedom of the system can be performed using the master equation and the mean-field decoupling approximation<cit.>. However it is a formidable task to solve a fairly large parameter space because of the infinite freedom of the environment. To address this problem, our group proposed a quasibosonic approach to eliminated the infinite freedom of the environment, in which the operators of the environment can be treated as a c-number and then the dissipative system can be solved easily<cit.>. One can obtain an effective Hamiltonian for the system based on the quasibosonic approach.

    H    =   ∑_iH^DM_i-κ∑_ijã^†_iã_j-μ∑_iÑ_̃ĩ
    
    H^DM_i   =   ω̃_aJ̃^+_iJ̃^-_i+ω̃_cã^†_iã_i+β(ã_iJ̃^+_i+ã^†_iJ̃^-_i)

where ω̃_a=ω_a-iγ_a, ω̃_c=ω_c-iγ_c. γ_a and γ_c are decay rates of atoms and cavities, respectively. ã_i^†(ã_i) is a quasiboson creation(annihilation) operator. J̃^+_i(J̃^-_i) is the dressed atomics raising(lowering) angular momentum operator. The dissipation becomes an inherent property for the DBH model considered here.

A superfluid order parameter ψ, with the mean field assumption ψ≡⟨ã_i⟩, is usually introduced to gain insight into the role of dissipation in the quantum phase transition. For ψ≠0, the system is in superfluid phase. When ψ=0, the system is in the Mott-insulator phase. In the present case, the expected value of ã_i is in general complex with the formation ⟨ã_i⟩=ψ-iψ_γ. ψ_γ is a solvable small quantity
as a function of decay rates of the system, and vanishes in the limit of ideal cases. Using the decoupling approximation ã^†_iã_j=⟨ã^†_i⟩ã_j+⟨ã_j⟩ã^†_i-⟨ã^†_i⟩⟨ã_j⟩. The mean-field Hamiltonian of Eq.(2.2) can be written as

    H^MF   =   ∑_iH_i^MF
    
    H^MF_i   =    H_i^DM-κψ(ã^†_i+ã_i)+κ|ψ|^2-μ∑_iÑ_̃ĩ

This mean-field Hamiltonian is assumed to be the same for every site.




§ EIGENVALUES AND EIGENSTATES OF THE DISSIPATIVE DICKE-BOSE-HUBBARD MODEL

In the following, the case of two TLAs in each cavity is investigated as an example to provide a detailed illustration. The extension for an arbitrary number of two-level atoms is given in 4, which can be easily calculated by using the
same approach. The bare states of system are |0,e^⊗2⟩|n-2⟩, |g,e⟩|n-1⟩, and |g^⊗2,0⟩|n⟩ with photon number n running from 0, 1, 2, 3 to ∞<cit.>. For two TLAs system, the case in which two atoms are in the excited state can be denoted as |0,e^⊗2⟩, only one atom in the excited state is denoted by |g,e⟩, and |g^⊗2,0⟩ is for the case that the two atoms are in the ground state. Here the total of 3n bare state bases form a group for the whole Hilbert space. Based on these states, the matrix elements for H^MF_n can be obtained.
 
 
    H^MF_n = |
    [     2ω̃_a+(n-2)ω̃_c-nμ             √(2(n-1)β)                      0;             √(2(n-1)β) 2ω̃_a+(n-1)ω̃_c-(n+1)μ                 √(2nβ);                      0                 √(2nβ)               nω̃_c-nμ;                        ]|+κ|ψ|^2


with ω̃_c=ω̃_a=ω̃(ω̃=ω-iγ), γ=γ_a+γ_c. The eigenvalues can be obtained by diagonalizing the matrix in Eq.(3.1), and the corresponding eigenstates can be found.

    E^(0)_|0,n⟩=nω̃


    E^(0)_|±,n⟩   =   (2n+1)ω̃±β R(n,ω̃/β)/2



    |0,n⟩=-√(n-1)|0,e^⊗2⟩|n-2⟩+√(n)|g^⊗2,0⟩|n⟩/√(2n-1)


    |±,n⟩   =   √(n)|0,e^⊗2⟩|n-2⟩+1/2√(2)[ω̃/β± R(n,ω̃/β)]|g,e⟩|n-1⟩+√(n-1)|g^⊗2,0⟩|n⟩/√(2n-1+{1/2√(2)[ω̃/β± R(n,ω̃/β)]}^2)


Here R(n,ω̃/β)=√(8(2n-1)+(ω̃/β)^2) is the effective Rabi frequency. The energy levels split into three branches corresponding to the upper branch E^(0)_|+,n⟩, center branch E^(0)_|0,n⟩ and the lower branch E^(0)_|-,n⟩ as shown in Fig. 1(b).


§  THE QUANTUM PHASE TRANSITION

In this section, we use the perturbation theory to obtain the superfluid order parameter and study the quantum phase transition by changing the controlling parameters. We have assumed that cavities are coupled weakly to each other, thus, the interaction term between cavities can be considered as a perturbation term when the two-level atoms system is coupled strongly to the cavity field. The effective Hamiltonian Eq.(2.3) thus reads

    H^MF_i   =   H^DM_i+H^'_i
    
    H^'_i   =   -κψ(ã^†_i+ã_i)+χ|ψ|^2-μÑ_̃ĩ

which is valid on each site, we therefore drop the subscript i in the following.
Considering the analogy of transition from the Mott-insulator to superfluid state between the Jaynes-Cummings model and the Bose-Hubbard model and the fact that the analytical results obtained by the second and fourth-order perturbations are in good agreement with the exact diagonalization numerical calculation<cit.>, we derive the analytically solution of the system in terms of the second-order perturbation for simplicity.  Eq.(3.2) and Fig. 1(b) show that a center energy level E_|0,n⟩ is required to perform the translation, thus, the on-site repulsion U based on the center branch state |0,n⟩ is independent of the atom-cavity coupling β,  which is different from one defined by the state |±,n⟩. To study the quantum phase transition in detail, the superfluid order parameter must calculated separately for different cases.

Preparing in the center branch of the dressed state: According to the definition of the superfluid order parameter ψ=⟨Φ_n(t)|ã_i|Φ_n(t)⟩, |Φ_n(t)⟩ can be obtained based on the second-order perturbation theory. We first obtain the second-order corrections of energy eigenvalues E^(2)_|0,n⟩ and (normalized) eigenstates ϕ̃^(2)_|0,n⟩ with respect to the dressed basis Eq.(3.4).

    E^(2)_|0,n⟩=(n-1)(2n-2)^2κ^2ψ^2/(2n-1)(2n-3)(ϵ-iγ)
     +4n^3κ^2ψ^2/(2n-1)(2n+1)(-ϵ+iγ)


    ϕ̃^(2)_|0,n⟩=√(n-1)(2n-2)(-κψ)/√((2n-1)(2n-3))(ϵ-iγ)|n-1⟩
    +2n√(n)(-κψ)/√((2n-1)(2n+1))(ϵ-iγ)|n+1⟩

where ϵ=ω-μ. Therefore, the eigenvalue of the dissipative system based on the second-order perturbation theory is

    E_|0,n⟩≡ E_s+iE_γ

with

    E_s   =   nϵ+κ|ψ|^2+(-8n^3+12n^2+4n-4)κ^2ψ^2ϵ/(2n-1)(2n-3)(2n+1)(ϵ^2+γ^2)
    
    E_γ   =   nγ+(-8n^3+12n^2+4n-4)iκ^2ψ^2γ/(2n-1)(2n-3)(2n+1)(ϵ^2+γ^2)

When the system is in the Mott-insulator state, ψ=0, we have E_γ= nγ. When ψ≠ 0, one can take E_γ≈ nγ because we assume that the coupling strength κ between cavities is weak.  Up to second order, the expression for the (normalized) eigenstates is

    ϕ_|0,n⟩=1/√(Ñ)ϕ̃_|0,n⟩


    ϕ̃_|0,n⟩=√(n-1)(2n-2)(-κψ)/√((2n-1)(2n-3))(ϵ-iγ)|n-1⟩+|n⟩
    +2n√(n)(-κψ)/√((2n-1)(2n+1))(ϵ-iγ)|n+1⟩


    Ñ=1+(n-1)(2n-2)^2κ^2ψ^2/(2n-1)(2n-3)(ϵ^2+γ^2)
    +4n^3κ^2ψ^2/(2n-1)(2n+1)(ϵ^2+γ^2)

Ñ is the normalized constant. For the open system considered here, the superfluid order parameter ψ is time-dependent. According to Eq.(4.5), the (normalized) eigenstates is a function of the time, however, its time derivative can be ignored because of the second-order correction. Thus, the approximative time-dependent wave function of the system can be written as

    Φ_n(t)=f(t)ϕ_|0,n⟩

Using the Schrödinger equation, one can find

    Φ_n(t)=ϕ_|0,n⟩e^-iE_|0,n⟩t

Therefore, the superfluid order parameter ψ for the state |0,n⟩ can be obtained


    ψ_1=e^-nγ t√((8n^3-12n^2-4n+4)ϵ/(16n^4-32n^3+12n^2+4n-4)κ-(2n-1)(2n+1)(2n-3)(ϵ^2+γ^2)/(16n^4-32n^3+12n^2+4n-4)κ^2e^-2nγ t)



Eq. (4.7) shows that ψ_1 is a function of the parameters κ, γ, t and μ (In present case, μ is a constant). The superfluid order parameter evolves and decays with time with a decay rate proportional to the number of photons n.

Preparing in the negative branch of the dressed state: Assume that each site is prepared in the
negative branch of the dressed state |-,n⟩. We can find the second-order deviations using a similar procedure, although the calculations become quite tedious when using our current formulation. The superfluid order parameter ψ_2 can be obtained by solving the following equation.


    ψ_2=Re{e^-2γ nt/Ñ^̃'̃[2[2√(n(n-1)(n-2))+√(n-1)/8(ω+iγ/β-R_n-1^†)(ω-iγ/β-R_n)]^2(-κψ_2)/[2ϵ+2iγ-β(R^†_n-R_n-1^†)][2n-1+1/8(ω-iγ/β-R_n-1)^2][2n-3+1/8(ω+iγ/β-R^†_n-1)^2]+
    2[2√(n(n-1)(n+1))+√(n)/8(ω+iγ/β-R_n^†)(ω-iγ/β-R_n+1)]^2(-κψ_2)/[-2ϵ+2iγ-β(R_n-R_n+1)][2n-1+1/8(ω+iγ/β-R^†_n)^2][2n+1+1/8(ω-iγ/β-R_n+1)^2]]}


Ñ^̃'̃ (Appendix) is the normalized constant. In what follows, we use Eqs. (4.7) and (4.8)  to numerically investigate the features of quantum phase transition arising from the competition between the on-site repulsion U_n and the hopping rate under the influences of the environment.



Analyses. As illustrated in Fig. 2. First, we start with a superfluid phase. The time evolution of the superfluid order parameter ψ for different hopping rates κ and decay rates γ are shown. Comparing Figs. 2a(b) with 2c(d), a clear quantum phase transition is found for different initial states. The ideal cases are also given in Fig. 2 for comparison, which shows that the system is still in the coherent state that was prepared initially.  The evolution of the dissipative system clearly reflects the expected decay of the coherence, which is the most obvious characteristic different from ideal cases. For a small t, although ψ decreases slightly, the system remains in a superfluid state. At a sufficiently large t, the effects of the environment become large, and the coherence of the system is initially destroyed in a pronounced manner and is then gradually reduced. Thus, ψ decays rapidly and the system undergoes a phase transition into a Mott-insulating phase. With the increase of the photon number n, the long-range order parameter will decrease rapidly, as shown in Fig. 3(b) and (d), because of the decay time is  proportional to n.  The critical point t_c is a function of controlling parameters and can be found by setting ψ=0 in Eqs. (4.7) and (4.8), which yields t_c=1/2nγln(4n^2-4n-4)κ/(2n+1)(2n-3)ϵ. It follows that, for certain cavity decay rates κ, one may change the other controlling parameters according to t_c to enable the dissipative system to maintain coherence for a relatively long time.  Obviously, when the
external environment is considered, the decoherence
of every resonator and the TLAs would result in the decay of the superfluid order parameter. In the experiment, dynamical decoupling<cit.> and feedback control<cit.> have been proposed to hamper the decay of the cavity field and TLAs and thus improved the coherence time.




In contrast, as shown in Fig. 3, we seek to determine how a Mott-insulator state in the beginning can restore the coherence by changing the intercavity hopping rate κ for a dissipative system.
For a small κ, there are no enough excitations for hopping between cavities. By rasing the hopping rate to a certain value κ_c=(8n^3-12n^2-2n+3)ϵ e^2nγ t/(8n^3-12n^2-4n+4), the system will restore its long-range coherence and a phase transition from Mott-insulator to the superfluid phase appears. According to Eqs.(4.7) and (4.8), the photon hopping rate is also found to decrease because of the effect of the environment, thus the long-range coherence can only occur when the increase of the photon hopping rate is faster than its decay. Figs. 3(a) and 3(c) also demonstrate that the influence of environment accumulate over time. With an increase in time, a large hopping rate is required to restore the coherence. Because the system and the environment have been recognized as a whole system in the effective Hamiltonian, then the dissipation is the inherent nature of the system. Therefore, for t=0, the system is also dissipative, and the hopping rate required for the phase transition to occur is higher than the rate expected in the ideal case.  In addition, increasing the number of photons to n, the dissipation of the system is also enhanced correspondingly, a higher hopping energy is thus required to induce a phase transition, as shown in Figs. 3(b) and 3(d).

In what follows, we extend the model to an arbitrary number of TLAs cases. The Dressed-state basis can be written by the general method to  diagonalize the effective DBH Hamiltonian (1.1) by numerical computation. The phase diagrams of the dissipative DBH model are plotted in Fig.4. For comparison, we also show the ideal cases. In dissipative cases, we choose t=0, which implies that the dissipative system is nearly equilibrium. As shown in Fig. 4, as interaction with the environment destroys the coherence of the system, the Mott lobes becomes smaller and the area of the coherent phase decreases. Next, the realization of the superfluid state requires a large hopping rate to derive the localized photons in each cavity. It can also be found that, in a regime with a small hopping rate κ, fewer TLAs could cause the system to become a localized phase compared with the ideal cases. With an increase in the number of TLAs, the coherent state may disappear rapidly for the dissipative system.




§ CONCLUSION

Based on the quasibosonic approach, a realistic situation of a DBH  model coupled to its environment was considered. The analytical solution of the superfluid order parameter for two TLAs per cavity was derived. The transition behaviors of the superfluid to Mott-insulation phase and the restoring coherence were discussed. The phase diagram for an arbitrary number of TLAs was also investigated. As the number of TLAs increases, Mott lobes may disappear and such a system tends to be classical. Most importantly, the atomic number is far lower than that under ideal circumstances. This work can provide parameters for reference to simulate strongly correlated many body systems in the actual operation.


This work was supported by the National Natural Science
Foundation of China under Grant No.11274148.


99
 M M. J. Hartmann, J. Opt. 18, 104005 (2016).
 Hur K. L. Hur, L. Henriet, A. Petrescu, K. Plrkhanov, G. Roux, and M. Schiró, C. R. Physique 17, 808-835 (2016).
 Noh C. Noh, and D. G. Angelakis, Rep. Prog. Phys. 80, 016401 (2017).
 Michael M. Knap, E. Arrigoni, and W. V. D. Linden, Phys. Rev. B 82, 045126 (2010).
 Srinivasan K. Srinivasan, and O. Painter, Nature(London) 450, 862 (2007).
 Aoki T. Aoki, B. Dayan, E. Wilcut, W. P. Bowen, A. S. Parkins, T. J. Kippenberg, K. J. Vahala, and H. J. Kimbe, Nature(London) 443, 671 (2006).
 Birnbaum K. M. Birnbaum, A. Boca, R. Miller, A. D. Boozer, T. E. Northup, and H. J. Kimble, Nature(London) 436, 87 (2005).
 Mabuchi H. Mabuchi, and A. C. Doherty, Science 298, 1372 (2002).
 Raimond J. M. Raimond, M. Brune, and S. Haroche, Rev. Mod. Phys. 73,565 (2001).
 Toyoda K. Toyoda, Y. Matsuno, A. Noguchi, S. Haze, and S. Urabe, Phys. Rev. Lett. 111, 160501 (2013).
 Underwood D. L. Underwood, W. E. Shanks, J. Koch and, A. A. Houch, Phys. Rev. A 86, 023837 (2012).
 Hartmann M. J. Hartmann, F. G. S. L. Brandão, and M. B. Plenio, Nat. Phys. 2, 849 (2006).
 Greentree A. D. Greentree, C. Tahan, J. H. Cole, and L. C. L. Hollenberg, Nat. Phys. 2, 856 (2006).
 Angelakis D. G. Angelakis, M. F. Santos, and S. Bose, Phys. Rev. A 76, 031805(R) (2007).
 Koch J. Koch, and K. L. Hur, Phys. Rev. A 80, 023811 (2009).
 Brandao M. J. Hartmann, F. G. S. I. Brandão, and M. B. Plenio, Laser Photon. Rev. 2, 527 (2008).
 Strack P. Strack, and S. Sachdev, Phys. Rev. Lett. 107, 277202 (2011).
 Li X. P. Li, and W. V. Liu, Phys. Rev. A 87, 063622 (2013).
 Jin J. Jin, D. Rossini, and R. Fazio, Phys. Rev. Lett. 110, 163605 (2013).
 Umucal1lar R. O. Umucal1lar, and Z. Carusotto, Phys. Rev. A 84, 043804 (2011).
 Carusotto R. O. Umucal1lar, and Z. Carusotto, Phys. Rev. Lett. 108, 206809 (2012).
 Andrew A. L. C. Hayward, and A. M. Martin, Phys. Rev. A 93, 053614 (2016).
 Bujnowski B. Bujnowski, J. K. Corso, A. L. C. Hayward, J. H. Cole, and A. M. Martin, Phys. Rev. A 90, 043801 (2014).
 Guo L. J. Guo, S. Greschner, S. Y. Zhu, and W. Z. Zhang, arXiv: 1611.06404.
 Minar J. Minář, S. Günes, Söyler, and Z. Lesanovsky, New J. Phys. 18, 053035 (2016).
 Maggitti A. Maggitti, M. Radonjié, and B. M. Jelenkovié, Phys. Rev. A 93, 013835 (2016).
 Dicke R. H. Dicke, Phys. Rev. 93, 99 (1954).
 Badshah F. Badshah, S. Qamar, and M. Paternostro, Phys. Rev. A 90, 033813 (2014).
 Chitra R. Chitra and O. Zilberberg, Phys. Rev. A 92, 023815 (2015).
 Mlynek J. A. Mlynek, A. A. Abdumalikov, C. Eichler, and A. Wallraff, Nature Commun. 5, 5186 (2014).
 Baumann K. Baumann, C. Guerlin, F. Brennecke, and T. Esslinger, Nature(London) 464, 1301 (2010).
 Hepp K. Hepp and E. H. Lieb, Phys. Rev. A 8, 2517 (1973).
 Klinder J. Klinder, H. Keßler, M. R. Bakhtiari, M. Thorwart, and A. Hemmerich, Phys. Rev. Lett. 115, 230403 (2015).
 Kebler J. Klinder, H. Keßler, M. R. Bakhtiari, M. Thorwart, and A. Hemmerich, PNAS 112, 11 (2015).
 Lei S. C. Lei and R. K. Lee, Phys. Rev. A 77, 033827 (2008).
 Rossini D. Rossini, and R. Fazio, phys. Rev. Lett. 99, 186401 (2007).
 Na N. Na, S. Utsunomiya, L. Tian, Y. Yamamoto, Phys. Rev. A 77, 031803(R) (2008).
 Lu Y. C. Lu, and C. Wang, Quantum Inf. Process 15, 4347-4359 (2016).
 Tan L. Tan and L. Hai, J. Phy. B 45, 035504 (2012).
 Liu K. L. Liu, L. Tan, C. H. Lv, and W. M. Liu, Phys. Rev. A 83, 063840 (2011).
 Buzek V. Bužek, M. Orszag, and M. Roško, phys. Rev. Lett. 94, 163601 (2005).
 Nissen F. Nissen, S. Schmidt, M. Biondi, G. Blatter, H. E. Türeci, and J. Keeling, Phys. Rev. Lett. 108, 233603 (2012).
 J  M. J. Hartmann, Phys. Rev. Lett. 104, 113601 (2010).
 Gerace I. Carusotto, D. Gerace, H. E. Tureci, S. D. Liberato, C. Ciuti, and A. Imamoǧly, Phys. Rev. Lett. 103, 033601 (2009).
 Aron C. Aron, M. Kulkarni, H. E. Türeci, Phys. Rev. X 6, 011032 (2016).
 Oosten D. V. Oosten, P. V. D. Straten, and H. T. C. Stoof, Phys. Rev. A 63, 053601 (2001).
 Bylander J. Bylander, S. Gustausson, F. Yan, F. Yoshihara, K. Harrabi, G. Fitch, D. G. Cory, Y. Nakamura, T. S. Tsai, and W. D. Oliver, Nat. Phys. 7, 565 (2011).
 Xue S. B. Xue, R. B. Wu, W. M. Zhang, J. Zhang, C. W. Li, and T. J. Tarn, Phys. Rev. A 86, 052304 (2012).





§ APPENDIX

The normalized constant of superfluid order parameter ψ_2.


    Ñ^'   =   4κ^2ψ^2AA^†+4κ^2ψ^2BB^†



    A=2√(n(n-1)(n-2))+√(n-1)/8(ω+iγ/β-R_n^†)(ω-iγ/β-R_n-1)/[2ϵ-2iγ-β(R_n-R_n-1)]√([2n-1+1/8(ω+iγ/β-R^†_n)^2][2n-3+1/8(ω-iγ/β-R_n-1)^2])



    B=2√(n(n-1)(n+1))+√(n)/8(ω+iγ/β-R_n^†)(ω-iγ/β-R_n+1)/[-2ϵ+2iγ-β(R_n-R_n+1)]√([2n-1+1/8(ω+iγ/β-R^†_n)^2][2n+1+1/8(ω-iγ/β-R_n+1)^2])



 The conjugates of A and B are A^† and B^†, respectively.
