Superradiant phase transition in the ultrastrong coupling regimeof the two-photon Dicke model

# Superradiant phase transition in the ultrastrong coupling regime of the two-photon Dicke model

L. Garbe Laboratoire Matériaux et Phénomènes Quantiques, Université Paris Diderot, CNRS UMR 7162, Sorbonne Paris Cité, France    I. L. Egusquiza Department of Theoretical Physics and History of Science, University of the Basque Country UPV/EHU, Apartado 644, E-48080 Bilbao, Spain    E. Solano Department of Physical Chemistry, University of the Basque Country UPV/EHU, Apartado 644, E-48080 Bilbao, Spain IKERBASQUE, Basque Foundation for Science, Maria Diaz de Haro 3, 48013 Bilbao, Spain    C. Ciuti    T. Coudreau    P. Milman    S. Felicetti Laboratoire Matériaux et Phénomènes Quantiques, Université Paris Diderot, CNRS UMR 7162, Sorbonne Paris Cité, France
July 15, 2019
###### Abstract

The controllability of current quantum technologies allows to implement spin-boson models where two-photon couplings are the dominating terms of light-matter interaction. In this case, when the coupling strength becomes comparable with the characteristic frequencies, a spectral collapse can take place, i.e. the discrete system spectrum can collapse into a continuous band. Here, we analyze the thermodynamic limit of the two-photon Dicke model, which describes the interaction of an ensemble of qubits with a single bosonic mode. We find that there exists a parameter regime where two-photon interactions induce a superradiant phase transition, before the spectral collapse occurs. Furthermore, we extend the mean-field analysis by considering second-order quantum fluctuations terms, in order to analyze the low-energy spectrum and compare the critical behavior with the one-photon case.

preprint: APS/123-QED

## I Introduction

In quantum optics, the superradiant phase transition Mallory69 (); Hepp75 () (SPT) is the abrupt change in the behavior of the ground state properties of a quantum many-body system, while a physical parameter is continuously varied. It is a quantum phase transition, i.e. it can be accessed at zero temperature and it is due to quantum fluctuations. The archetypal model known to display such a transition is the Dicke model Gelfand54 (); Brandes05 () (DM), which describes the interaction of an ensemble of two-level quantum systems, or qubits, with a single bosonic mode. In the limit of a large number of qubits, the DM undergoes a SPT Emary03 () in the ultrastrong coupling (USC) regime, where the collective light-matter coupling becomes comparable to the qubit and field bare frequencies Ciuti05 (). Although the DM is commonly used to describe atomic and solid-state systems, whether it provides a reliable description of the system ground state when approaching the critical coupling is still the subject of debate Rzazewski75 (); Knight78 (); Gawedzki81 (); Keeling07 (); Deliberato14 (); Vukics14 (); Griesser16 (). In particular, the presence of the so-called diamagnatic term is expected to prevent the SPT. The debate has been recently extended to the framework of circuit QED Chen07 (); Lambert09 (); Nataf10 (); Viehmann11 (); Jaako16 (); Bamba16 (), where the USC regime has been experimentally achieved Niemczyk10 (); Forn-Diaz10 (); Forn-Diaz16 (); Chen16 (); Yoshihara16 ().

However, a compelling way to circumvent no-go theorems consists in using driven systems to engineer effective Hamiltonians. Indeed, the SPT has been observed in driven atomic systems which effectively reproduce the DM Baumann10 (); Baden14 (); Klinder15 (). In general, driven atomic or solid state systems represent a powerful tool to access the USC regime of quantum optical models, both in few Crespi12 (); Ballester12 (); Pedernales2015 (); Langford16 (); Felicetti16 () and many-body physics Dimer07 (); Carusotto13 (). In the USC regime, even apparently simple models entail a very complex physics. This is the case for the quantum Rabi model Braak11 (); Rossatto16 (), which corresponds to a single-qubit DM. Furthermore, the engineering of effective Hamiltonians allows to implement generalized quantum optical models Genway14 (); Zou14 (), including anisotropic couplings or two-photon interactions. In the case of anisotropic couplings, reaching the USC regime leads to parity-symmetry breaking Xie14 (); Cui15 () and to a rich phase diagram in the many-body limit Baksic14 ().

Similarly, the two-photon Rabi model has highly counter-intuitive spectral and dynamical features. Its spectrum collapses into a continuous band for a specific value of the coupling strength Emary02 (); Dolya09 (); Travenec12 (); Peng16 (). In the transition from the strong to the USC regime of the two-photon Rabi model, a continuous symmetry breaks down into a four-folded discrete symmetry, identified by a generalized-parity operator Felicetti15 (). However, so far there are no known results on the ground state of two-photon models in the many-body limit.

In this work, we perform first a mean-field analysis of the two-photon Dicke model and we find that the system exhibits a phase transition in the thermodynamic limit. This transition is superradiant in the sense that it is characterized by a macroscopic change in the average photon number. The boundary of the phase transition is set by the critical value of the coupling strength, which depends on the qubit and field energies. Interestingly, the coupling strength for which the spectrum collapses depends only on the field frequency. For larger values of the coupling strength, the Hamiltonian is not bounded from below and the model is not well defined. We define the parameter regime where the SPT could be accessed within the validity region of the model, that is, where the critical coupling is smaller than the collapse coupling strength. Finally, we go beyond mean-field by including second-order quantum fluctuations. This lets us characterize the system phases and analyze the differences with the SPT of the standard DM. We find fundamental differences in the critical scaling of the bosonic field.

The two-photon DM could be implemented using trapped ions Felicetti15 (), which have been used to realize spin systems composed of hundreds of qubits Britton12 (). Similar schemes can be conceived for other atomic or solid state systems. Particularly promising are superconducting devices, where bosonic modes have been coupled to increasingly large spin ensembles Macha14 (); Kakuyanagi16 (). In any implementation, the critical issue would be the number of qubits that can be effectively coupled to a single bosonic mode. In our case, we have considered the thermodynamic limit , with the number of qubits. We show that, in this limit, the superradiant phase transition of the two-photon DM entails squeezing and spin-squeezing properties. These features provide a signature of the phase transition that could be observed Sorensen99 (); Vitagliano14 (); Saideh16 () for a smaller number of qubits.

## Ii Mean-Field Analysis

We consider an ensemble of N qubits interacting with a bosonic mode via two-photon interaction, as sketched in Fig.1. The system Hamiltonian is given by:

 ^H=ω^a†^a+ωq2N∑j=1^σjz+gNN∑j=1^σjx(^a2+^a†2), (1)

with , and are Pauli operators describing the j-th ion, and and are bosonic ladder operators. In Ref.Felicetti15 (), Felicetti et al. have proposed to implement this model using a chain of ions in a trap illuminated by two lasers, with the motional degree of freedom of the chain playing the role of the bosonic field. The coupling strength , qubit energy spacing and bosonic frequency are then effective and tunable parameters which depend on the frequencies and amplitudes of the two lasers. In principle, it is then be possible to reach the ultrastrong coupling (USC) regime: . In the following, we will consider this regime of parameters, as well as the thermodynamic limit . In the USC regime, the two-photon Dicke model exhibits a spectral collapse Travenec12 (); Felicetti15 (): for , the energy levels of the system collapse into a continuum. Beyond this limit, the ground state of (1) is no longer defined, which renders the notion of phase transition meaningless. Thus, the goal of this work is to study the existence of a phase transition for . For this purpose, it is convenient to describe the ensemble of qubits by pseudospin operators: we define , , , which gives us:

 ^H=ω^a†^a+ωq^Jz+gN(^J++^J−)(^a2+^a†2). (2)

In the following, we will speak only in terms of the fluctuations and polarization of this "spin", which are physically interpreted in terms of population and coherence of the states of the qubits. We study the phase diagram of our model by using a mean-field approach, inspired by the analysis performed in Emary03 () and Baksic14 (). The point here is to determine how the properties of the ground state evolve when increases. First of all, we use the Holstein-Primakoff (HP) transformation to turn our pseudospin operators into bosonic operators:

 ^J+=^b†√N−^b†^b,^J−=√N−^b†^b^b,^Jz=^b†^b−N2, (3)

with . Here and in the following, we restrict ourselves to the eigenspace of associated with eigenvalue (that is, the maximal angular momentum eigenspace ) note subspace ().

Next, we shift and with respect to their mean values in the ground state of (2):

 (4)

As a zeroth order approximation, we neglect altogether the spin fluctuations and , which gives:

 ^H=ω^a†^a+gβ(^a2+^a†2)+ωq|β|2−ωqN2, (5a) gβ=gN(β+β∗)√N−|β|2. (5b)

In other terms, we replace the spin operators by their classical mean values. This Hamiltonian is quadratic in , and can therefore be diagonalized by Bogoliubov transformation. The ground state is a squeezed vacuum state with squeezing parameter:

 r(MF)a=12arctanh(2gβω). (6)

The corresponding ground state energy is given by:

 EG=12cosh(2r(MF)a)ω(ω2−4g2β)+ωq|β|2−ωqN2−ω2. (7)

The final step consists in minimizing this energy in order to determine the value of the system adopts in its ground state. then plays the role of order parameter for our system: a change in its behavior leads to a modification of the qualitative properties of the ground state, which indicates a phase transition. We find that, for , is minimal for . For , we have two degenerate minima:

 β=β∗=±√N2(1−√1−μ4μ2λ2−μ)1/2=±β0, (8)

where we have defined:

 λ=ω2ωqN≥0;μ=4g2ω2≥0. (9)

Thus, our system exhibits two phases: in the first, and the squeezing parameter are zero for all , meaning that the field is in the vacuum state and the pseudospin is polarized along the -axis: , . In the second phase, the ground state is twice degenerate. This degeneracy comes from the fact that can be positive or negative. The field is in a squeezed vacuum state, the direction of squeezing depending on the sign of : (i) For , the quadrature is squeezed, the fluctuations of are amplified. (ii) For , the squeezed quadrature is . In both cases, when increases, the squeezing parameter increases and the pseudospin polarization evolves towards the -axis: decreases until it hits zero for , and increases at the same time. This phase diagram is reminiscent of the one-photon Dicke model and its superradiant phase transition Emary03 (); Baksic14 (). Here, only the pseudospin acquires macroscopic mean value in the second phase; the mean value of remains zero. However, the average number of photons becomes nonzero at the transition; thus, we argue that our transition may still be dubbed "superradiant". Note also that the zero mean value of the bosonic field comes from the fact that we have considered only quadratic terms for in the Hamiltonian. If a linear term is present, will acquire nonzero mean value at the mean-field level (See Appendix A).

Let us finally notice that, for the phase transition to occur before the spectral collapse, the following condition must be satisfied: . Since and are adjustable effective parameters, it is possible in principle to meet this condition. From this point on, we will consider that the order of magnitude of the parameter is (). We can notice that the approximation that we have made to obtain (5) amounts simply to keeping only terms of order in the Hamiltonian, and to neglecting higher-order terms note zero omegaq ().

## Iii Beyond Mean-Field

Let us now take the fluctuations of into account, in addition to the mean value. We will expand the Hamiltonian keeping terms of order , and . Then, we will make use of a technique that Hwang & al. applied to the Rabi and Jaynes-Cumming models article Hwang JC (); article Hwang Rabi (). This method is inspired by the Schrieffer-Wolff transformation ref Schrieffer-Wolff (); ref adiabatic elimination (), and can be summarized as follows: one starts with an Hamiltonian describing a spin operator and another degree of freedom, that can be written as

 ^H=^H0+ϵ^Hc, (10)

with a small parameter , an operator that does not couple the eigenstates, and that does. In article Hwang JC () for instance, this method was applied to a Jaynes-Cummings Hamiltonian:

 ^HJC=ω0^a†^a+Ωspin2^σz−ϵg(^a^σ++^a†^σ−). (11)

The method consists in decoupling the spin eigenspaces up to a certain order of . This is done by finding a transformation , with an anti-hermitian operator such that commutes with up to a certain order in . It is then possible to project in one of the spin eigenspaces, which allows to effectively suppress the spin degree of freedom and to diagonalize the Hamiltonian more easily.

In our case, it would seem natural to decouple the eigenspaces of the pseudospin . Instead, however, we are going to apply the HP transformation once more, and shift the operators with respect to their mean-field expectation value .

We define the following operators:

 ^K0=12(^a†^a+12), (12) ^K+=12^a†2, (13) ^K−=12^a2, (14)

which obey spin-like commutation relations , and . Note that the commutation algebra here is SU(1,1) instead of SU(2), hence these are not spin operators even if the commutation relations are the same. In order to apply the method described above, we take profit of these commutation relations to decouple the eigenspaces of these pseudospin operators up to a certain order of the small parameter . Then, we project out the degree of freedom. This gives us an effective Hamiltonian describing the low-energy fluctuations of above the ground state; the detailed calculations can be found in Appendix B, C and D.

While not very intuitive at first glance, this manipulation can be justified by two arguments. Firstly, it is interesting to use the HP transformation on again in order to make a link to our previous study and to re-use our results. Hence, only is available to play the role of the pseudospin whose eigenspaces we seek to decouple. Next, as we mentioned previously, we consider the following regime of parameters: . Since , this means that the typical energy separation between the eigenstates, given by , is much bigger than the typical energy separation for the (or ) eigenstates, given by . As a consequence, the eigenspaces associated with the pseudospin (12) are well separated in energy, but the eigenspaces of are not. A schematic representation of the energy levels is shown in Fig.2. Thus, it is possible to project the Hamiltonian into the lowest-energy eigenspace of and to study the fluctuations of while staying inside this subspace, but not the opposite.

In the first phase, after having decoupled the eigenspaces and projected the system into the lowest-energy one, we obtain the following Hamiltonian:

 ^H=−ωqN2+ωq^d†^d−g2Nω(^d+^d†)2+O(ωN√N), (15)

where has been defined in Eq.(4). The associated ground state for is a squeezed vacuum state, with squeezing parameter (let us note that this parameter is negative, meaning that the squeezed quadrature is instead of ). This constitutes a piece of information about the fluctuations of , while only mean values were accessible in our first analysis. The field is found to be in a coherent vacuum state, with no modifications with respect to the mean-field analysis.

In the second phase, we also find squeezing properties for the ground state, with squeezing parameter:

with . The results for the field are identical to the mean-field case with a slight correction, i.e. we have a squeezed vacuum state with squeezing parameter which differs from by a correction of order (see Appendix D). Going back to the definition of the field (3), the behavior for the operators can be summarized as follows: in the first phase, is polarized along the -axis: . When increases, the fluctuations of are damped, while the fluctuations of are amplified proportionally. The amplification factor goes to infinity when approaching the transition: even though the approximations we have used break down near the critical point, the divergence of the fluctuations in our model may have observable consequences in an actual experiment. In the second phase, the polarization will gradually evolve from the -axis to the -axis as increases. Near the spectral collapse , will be polarized along the -axis, and we will have squeezing properties for the fluctuations in the and directions. As of the field, we have a coherent vacuum state in the first phase; then, at the transition, the field acquires squeezing properties. In the second phase, the squeezing parameter increases with and diverges at the spectral collapse. The behavior of spin squeezing fluctuations is schematically depicted in Fig.3. We can also characterize the ground state energy and the excitation energy in both phases: the results for the excitation energy are displayed on Fig.4.

Finally, using the effective Hamiltonian in both phases, we can compute the critical exponents of several observables that exhibit critical behavior at the transition, that is,

 A(g→gt)∝|g−gtgt|γA,

being the critical exponent of . We can compare those results to the critical exponents in the one-photon case Emary03 (); Vidal06 (), with an interacting term of the form (note that this convention differs slightly from the one usually found in the literature, ). In the one-photon model, the atomic and photonic excitations are hybridized to form polaritons. The energy spectrum exhibits two polaritonic branches, one of which goes to zero at the transition. In the limit , however, the excitations cease to hybridize; the lowest- and highest-energy polaritons are purely atomic and photonic, respectively, which corresponds to our low-energy excitations and high-energy () excitations. Notice that this situation is similar to what is considered in Ref.Mottl12 (), where an atomic ensemble interacts with a strongly detuned optical cavity. We compare the critical exponents of the one- and two-photon Dicke models in the regime , in a such a way that the polaritonic branches have the same atomic and photonic components in both cases. The results (computed for ) are summarized in Table 1.

Notice that the variance of remains constant near the transition in the two-photon case, but diverges in the one-photon case, which is a marked difference. However, the diverging terms in the one-photon case are high-order terms; near the transition, . Thus, as long as one does not get too close to the transition, the behavior of remains identical in the one-photon and the two-photon cases. When , though, the results begin to differ.

On the other hand, the critical exponents are the same in both cases for and . Nevertheless, for these quantities, some of the higher-order terms we have neglected in our analysis may diverge faster near the transition that the terms we took into account, leading to the breakdown of our analysis when comes close enough to the critical value . More precisely, the scaling parameter for is no longer valid when becomes comparable with .

## Iv Summary and Outlook

As a conclusion, we have shown the presence of a superradiant phase transition of the two-photon Dicke model in the ultrastrong coupling regime. We have characterized the behavior of the qubits and bosonic field in both phases, and we have studied some of their critical properties near the transition. With respect to the one-photon case, fundamental differences arise in the behavior of the bosonic field , which does not acquire macroscopic occupation in the second phase, and whose fluctuations do not diverge at the critical point, at least at the order considered. The pseudospin , on the other hand, does exhibit diverging fluctuations at the transition that could lead to observable phenomena which would mark the transition in experimentally accessible situations. Indeed, as an extension of this work, it would be interesting to analyze the spin-squeezing properties of the ground state of the system for a finite number of qubits. Quantum-phase transitions and the ground-state properties could be analyzed also modifying the symmetries of the model Baksic14 (), considering anisotropic couplings or including multiple bosonic fields. Finally, it would be interesting to study dynamical features and the emergence of quantum chaos Emary03 (); Alvermann12 (); Bakemeier13 () in the two-photon Dicke model.

## Acknowledgments

This work was supported by the French Agence Nationale de la Recherche (SemiQuantRoom project ANR-14-CE26-0029-01) and from University Sorbonne Paris Cite EQDOL contract. S.F. acknowledges funding from the PRESTIGE program, under the Marie Curie Actions-COFUND of the FP7. E.S. and I.L.E. acknowledge funding from Spanish MINECO/FEDER FIS2015-69983-P, Basque Government IT986-16, UPV/EHU UFI 11/55.

## Appendix A Linear term expansion of the Hamiltonian

Let us consider an extension of our model, obtained by adding a linear term in :

 ^H=ω^a†^a+ωq^Jz+g2N^Jx(^a2+^a†2)+g1N^Jx(^a+^a†). (16)

This gives, in the mean-field approximation:

 ^H=ω^a†^a+gβ2(^a2+^a†2)+gβ1(^a+^a†)+ωq|β|2−ωqN2, (17a) gβ1=g1N(β+β∗)√N−|β|2, (17b) gβ2=g2N(β+β∗)√N−|β|2. (17c)

We can shift the operator: . A proper choice of , namely, , allows to suppress the linear term in . This gives us:

 ^H=ω^c†^c+gβ2(^c2+^c†2)+ωq|β|2−ωqN2−(gβ1)2ω+2gβ2, (18)

which is just Hamiltonian (5) with an additional constant term. Hence, in the ground state is in a displaced squeezed state, with squeezing factor and . There is also a correction on the values of and . Importantly, the presence of the term in (18), which is not invariant under the transformation , lift the degeneracy between and .

## Appendix B Description of the method used to go beyond mean-field

We are now going to detail the method used in article Hwang JC () and article Hwang Rabi (). We begin with (10), and we assume that we have a two-dimensional pseudospin operator . is then a block diagonal operator with respect to the eigenbasis, and is block anti-diagonal with respect to this eigenbasis (in the following, we will simply talk about block diagonal and block anti-diagonal operators, and about pseudospin eigenstates and eigenspaces: the reference to is implicit). We apply the unitary transformation to , which yields:

 ^H′=e−^S^He^S=∞∑k=01k![^H,^S](k), (19)

with and . We impose to be block anti-diagonal. If we split into its block diagonal and block anti-diagonal parts, we get:

 ^H′d=∞∑k=012k![^H0,^S](2k)+∞∑k=012k+1![ϵ^Hc,^S](2k+1), (20a) ^H′a=∞∑k=012k+1![^H0,S](2k+1)+∞∑k=012k![ϵ^Hc,^S](2k). (20b)

These equations stem from the following relations: consider and two arbitrary block diagonal operators, and and two arbitrary block anti-diagonal operators. Then, and are block diagonal and is block anti-diagonal.

At this point, the idea is to expand with respect to : , with block anti diagonal for all . By properly choosing the for from to with arbitrary, it is possible to cancel up to order , thus decoupling the pseudospin eigenspaces up to this order.

Our case, however, is slightly different, because we have an infinite number of eigenstates instead of only two. When SU(1,1) is presented in terms of quadratic combinations of creation/annihilation operators the relevant representations are those with Bargmann parameter and . The first one corresponds to even number of excitations, while the second one is associated with an odd number article Duan Bargmann representation and Rabi analytic solution (). They are not connected by any of the operators of the algebra. Since we will focus on the ground state, we concentrate on the case. Then, we have an infinity of eigenstates , with . This means that we have to modify slightly the method above if we want to apply it to bigger Hilbert spaces.

For this, let us consider the operators diagonal in the basis, which can be written as with arbitrary . Let us call the ensemble of these operators ; , for instance, belongs to this ensemble. Then, let us call the ensemble of all operators of the form , with arbitrary coefficients; is an element of . contains the operators of the form ; and in a general way, we define the ensemble that contains operators that can be written as:

 ^F(j)+^F(j−2)+^F(j−4)+...=∑p

Let us note here that there is a redundancy in this definition: an element of also belongs to , for all . We will use the following property: for all belonging to and belonging to , belongs to , which we denote symbolically by

 [M(i),M(j)]op⊆M(i+j). (22)

The idea in the following will be to isolate the elements of the various and to cancel those that couple the eigenstates of to a certain order.

## Appendix C Study of the first phase

As indicated in the main text, we perform the Holstein-Primakoff transformation, as well as the transformation (12-14), and we consider the case : . This gives us for the Hamiltonian:

 ^H=−12ω−ωqN2+2ω^K0+ωq^d†^d+2g√N(^K++^K−)(^d+^d†). (23)

With a redefinition of parameters and , we get:

 ^H1 = (^H+12ω+ωqN2)2ω (24) = ^K0+ω1N^d†^d+ω2√N(^d+^d†)(^K++^K−),

which is explicitly of the form (10); does not couple the eigenstates, but does. In our case, the small parameter is .

The generalized method goes as follows: all the operators implied in our calculations can be (not uniquely) decomposed into elements belonging to the different : , where can take an arbitrary high value. We then propose the following decompositions:

 ^H1=^H(0)1+^H(1)1, (25a) ^H(0)1=^K0+ω1N^d†^d, (25b) ^H(1)1=ω2√N(^d+^d†)(^K++^K−)=1√N^V(1), (25c)
 (26)

with and that belong to and that belongs to . Using (22) and (19), we then obtain the following decomposition for :

 ^H′(0)1=^H(0)1=^K0+ω1N^d†^d, (27a) ^H′(1)1=1√N^V(1)+1√N[^K0,^P]+1N[^K0,^Q]+O(1N√N), (27b) ^H′(2)1=1N[^V(1),^P]+12N[[^K0,^P],^P]+1N[^K0,^R]+O(1N√N), (27c) ^H′(3)1=O(1N√N), (27d)

since couples the eigenvalues of , we are going to cancel it at order . For this purpose, we need and . We note that these conditions can be met by the following choice of operators:

 ^P=−ω2(^d+^d†)(^K+−^K−), (28a) ^Q=0. (28b)

As of the term, we find, using the above expressions for and :

 ^H′(2)1=−4ω222N(^d+^d†)2^K0+1N[^K0,^R]+O(1N√N). (29)

In general, we should choose so as to cancel the non-diagonal terms of this operator. Here, though, is already diagonal; it is thus sufficient to set . We speculate that, for higher order expansions, if can be made diagonal at order , we can make it diagonal at order by adding terms of the form to the expansion of , with that belongs to and . It could be interesting in a future study to consider higher-order terms in a more systematic way to confirm or to invalidate this conjecture.

In the end, all those manipulations yield:

 ^H′1=e−^S^H1e^S=^K0+ω1N^d†^d−2ω22N(^d+^d†)2^K0+O(1N√N). (30)

This Hamiltonian does commute with ; projection in the ground state gives (according to the definition (12)). Restoring the constants in the Hamiltonian yields the operator described in the main text:

 ^H=−ωqN2+ωq^d†^d−g2Nω(^d+^d†)2+O(ωN√N). (31)

This effective Hamiltonian describes the fluctuations near the ground state. Being quadratic in , it can be diagonalized by using a Bogoliubov transformation of parameter , which corresponds to a squeezed vacuum state as described in the main text. The state of is given by the ground state of ; here it is just the coherent vacuum state, as in the mean-field scenario. Finally, the ground state energy and the excitation energy are computed: we obtain

 E(1)exc=ωq√1−4g2Nωωq, (32a) E(1)G=−ωqN2+E(1)exc−ωq2. (32b)

## Appendix D Study of the second phase

We proceed in a similar fashion, using HP transformation once more, and setting ; but this time . As a starting point, we use the value of obtained by our mean-field analysis. For ease of notation, we define the following parameters:

 α=β√N=O(1), (33a) χ=√1−β2N=√1−α2=O(1), (33b) δ=1−β2N−β2=O(1). (33c)

Expanding the Hamiltonian gives:

 \vspace15pt^H2 = (^H+12ω+ωqN2−ωqNα2)/2ω = λ0^K′0+λ1√N(d+^d†)+λ2√N(d+^d†)(^K′++^K′−)+λ3√N(d+^d†)^K′0 +λ4N^d†d−2N^K′0^V1(^d)+1N(^K′++^K′−)^V2(^d)+O(1N√N).

Where we have defined new operators and parameters:

 ^K′0=cosh(2r(2)a)^K0+12sinh(2r(2)a)(^K++^K−), (35a) ^K′++^K′−=cosh(2r(2)a)(^K++^K−)+2sinh(2r(2)a)^K0, (35b) r(2)a= 12arctanh(4gαχω+gαωχN) (35c) =12arctanh(2gβω+gαωχN),
 λ0=cosh(2r(2)a)−(4gαχω+gαωχN)sinh(2r(2)a), (36a) λ1=ωqNα2ω, (36b) λ2=gχδωcosh(2r(2)a), (36c) λ3=−2sinh(2r(2)a)gχδω, (36d) λ4=ωqN2ω, (36e) ^V1(^d)=sinh(2r(2)a)(−gαχω^d†^d−gω(α2χ+α34χ3)(^d