Quantum stability of self-organized atomic insulator-like states in optical resonators

We investigate a paradigm example of cavity quantum electrodynamics with many body systems: an ultracold atomic gas inside a pumped optical resonator. In particular, we study the stability of atomic insulator-like states, confined by the mechanical potential emerging from the cavity field spatial mode structure. As in open space, when the optical potential is sufficiently deep, the atomic gas is in the Mott-like state. Inside the cavity, however, the potential depends on the atomic distribution, which determines the refractive index of the medium, thus altering the intracavity field amplitude. We derive the effective Bose-Hubbard model describing the physics of the system in one dimension and study the crossover between the superfluid -- Mott insulator quantum states. We determine the regions of parameters where the atomic insulator states are stable, and predict the existence of overlapping stability regions corresponding to competing insulator-like states. Bistable behavior, controlled by the pump intensity, is encountered in the vicinity of the shifted cavity resonance.

In recent years considerable attention has been paid to a new regime of CQED, which we term CQED with many-body systems. These studies focus onto the mechanical effects of the resonator field on the atomic motion, and on the non-linearity arising from the interdependence between the cavity field and the atoms dynamics. Following the theoretical prediction of Ref. [10], signatures of self-organization have been measured in the light scattered by laser-cooled atoms in a transversallypumped cavity [11]. These structures and their properties have been theoretically studied in detail in Refs. [12,13]. In different setups, Bragg scattering of atomic structures inside optical resonators has been experimentally investigated in Ref. [14].
While in all the cases mentioned so far the atomic motion is essentially classical, the stability and properties of these structures in the quantum regime are still largely unexplored. This question acquires a special relevance in view of the recent experimental progress of CQED with ultracold atoms. In fact, strong atom-field coupling between Bose-Einstein condensed atoms (BEC) and the mode of a high-finesse optical cavity has been realized in the experiments reported in Refs. [15,16,17,18]. Moreover, CQED techniques were used to measure pair correlations in the atom laser [19], and have been proposed for characterizing quantum states of ultracold matter [20,21,22].
In [23] we investigated the ground state of ultracold atoms in the optical lattice formed by the interaction with the cavity mode. This system combines CQED with the many-body physics of strongly correlated ultracold atoms. In particular, the non-linear dependence of the cavity field on the atomic motion opens novel perspectives to the rich scenario of ultracold atomic gases in optical lattices. In open space, in fact, these systems offer the possibility to realize paradigmatic systems of quantum many-body physics [24,25], such as various versions of Hubbard models [26]. A prominent example is the Bose-Hubbard model [27], exhibiting the superfluid (SF) -Mott insulator (MI) quantum phase transition [28], whose realization with ultracold atoms was proposed in Ref. [29], and demonstrated in Ref. [30]. In [23] we addressed the question whether and how this transition is modified when the atoms are inside a resonator, where the optical lattice due to the intracavity-field depends on the atomic density.
In this article we report the details of the derivations presented in [23] and extend them to novel regimes. The system we consider consists of ultracold atoms inside a resonator, which is driven by a laser. Due to the strong coupling between cavity and atomic degrees of freedom, the atoms shift the cavity resonance, thus modifying the intracavity field intensity. This in turn determines the depth of the cavity potential. At ultralow temperatures we assume that the atoms occupy only states of the lowest band of the periodic optical potential. In this regime, we present the detailed derivation of the Bose-Hubbard model for atoms in the one-dimensional potential of an optical resonator, which complements and extends the derivation for few atoms by Maschler and Ritsch in Ref. [31,32] to large numbers of atoms and which is valid in an appropriate thermodynamic limit. Using this model we study the SF -MI crossover as a function of the system parameters: the chemical potential µ, the pump strength and frequency, and the atomic density. Assuming the tight-binding regime, we may describe the MI states using Wannier functions [33], whose form is determined by the optical potential. The Wannier functions are used then to calculate the coefficients in the Bose-Hubbard model, as in standard textbooks. We determine the boundaries of the MI states, using the strong coupling expansion of Ref. [34], which is a quite accurate method for the calculation of the phase-diagram of the BH model in 1D [25].
It must be stressed that the derivation of the Bose-Hubbard model in the cavity does involve certain novel aspects. Namely, the periodic optical potential depends functionally on the atomic density, and hence on the Wannier functions. The problem is hence highly nonlinear: the coefficients of the equations determining the Wannier functions depend functionally on the Wannier functions themselves, and the latter have thus to be determined self-consistently. This property has also physical implications. In fact, since in our system the coefficients of the Hubbard model depend on the atomic density, consequently, the diagrams in the µ-t plane, where t is the tunneling energy, exhibit overlapping, competing Mott states, that may even consist of disconnected regions for a wide range of parameters. In the vicinity of the shifted cavity resonance in the strong coupling regime, the situation is even more complex: one encounters there also dispersive bistable behavior [3]. Thus, from the quantum optical perspective, this paper investigates stability of Mott-like phases, i.e., insulator-like states, in an optical resonator.
This article is organized as follows. In Sec. II we present our model that describes a system of two level atoms confined along the axis of an optical cavity. In Sec. II A we introduce the single-atom Hamiltonian. The many-body dynamics including the quantum noise is introduced in Sec. II B. Section II C is devoted to the physical discussion of the role of the various physical parameters on the system dynamics. The effective Bose-Hubbard Hamiltonian is derived in Sec. III. In Sec. IV we discuss the ground-state properties of our model. We use the Gaussian approximation for the Wannier wave functions, checking carefully its validity. Within the Gaussian approximation and the strong coupling expansion method of H. Monien [34] the stability regions for the Mott states are obtained analytically, up to the solution of the non-linear self-consistent equation for the width of the Wannier functions. Numerical results are reported and their physical meaning is discussed in Sec. IV B. The validity of the approximations is addressed in Sec. IV C. We conclude in Sec. V, while in the appendices the details of the derivation of the Bose-Hubbard Hamiltonian and of the strong coupling expansion method are reported.

II. THE MODEL
In this section we generalize the quantum optical model of a single atom inside a cavity to the many-body case, considering particle collisions at ultralow temperatures and quantum statistics. At this purpose, we first intro-duce the single-atom dynamics, write then the Hamiltonian in second quantization, introducing atom-atom collisions, and discuss the basic properties.

A. Single-particle dynamics
We consider a single atom of mass M inside a cavity. The atomic dipole transition at frequency ω 0 , between the ground state |g and the excited state |e , couples quasi-resonantly with an optical mode of the resonator at frequency ω c , wave vector k and position-dependent coupling strength g(x) = g 0 cos(kx), g 0 being the vacuum Rabi frequency. The resonator is driven by a classical field of amplitude η and oscillating at frequency ω p . We consider the atomic motion along the cavity axis, which coincides here with the x-axis, and assume tight confinement along the radial plane so that the transverse motion can be considered frozen out. Atomic center-of-mass position and momentum operators are x and p, fulfilling the uncertainty relation [x, p] = ih. In the reference frame oscillating at the frequency ω p of the pump field, the normally-ordered Hamiltonian describing the coherent dynamics of the atomic and cavity-mode states reads where ∆ a = ω p − ω a and ∆ c = ω p − ω c are the pumpatom and pump-cavity detunings, a (a † ) the annihilation (creation) operator of a cavity photon at frequency ω c , fulfilling the commutation relation a, a † = 1, and σ = |g e|, σ † = |e g| are the dipole lowering and raising operators. Spontaneous emission of the atomic dipole at rate γ and cavity losses at rate κ are described within the quantum Langevin equations formalism, such that the quantum Heisenberg-Langevin equations for the dipole and cavity operators read [1] where σ z = σ † σ − σσ † and a in , f in are the input noise operators, whose mean value vanishes and which are δ−correlated in time, namely, At large atom-pump detuning the adiabatic elimination of the excited atomic state can be performed. Assuming that the changes in the atomic position are negligible during the time scale in which the atom reaches its internal steady state, namely when k B T ≪h|∆ a |, we solve the Heisenberg-Langevin equations at a fixed value of the atomic position x. Hence, for |∆ a | ≫ g 0 n , γ, |∆ c |, we set σ z (t) ≈ −1 in the equations, and obtain σ † ≈ ig(x)a † /∆ a . After tracing out the internal degrees of freedom, the single-particle Hamiltonian for cavity and atomic center-of-mass degrees of freedom reads where we have used the explicit form of the cavity spatial mode function, and is the depth of the single-photon dipole potential.

B. Many-body dynamics
We now extend the previous model and derive the corresponding effective Hamiltonian for a gas of N bosons at ultralow temperatures. The particle interactions are modeled by s-wave scattering. We introduce the field operators Ψ j (x), Ψ † j (x), with j = g, e labeling the internal ground state, such that In second quantization the Hamiltonian (1) becomes H and is decomposed according to H = H 0 + H 1 , where with u j the strength of the onsite interaction depending on the atomic state, and In the above description we omit to write the Hamiltonian term describing the collisions between atoms in different internal states, as we will consider that the excited state is essentially empty in the parameter regime we choose. The quantum Heisenberg-Langevin equations for atomic and field operators reaḋ where f in (t) and a in (t) are the noise operators defined in the previous section. Solving the equation for Ψ e (x, t) in the limit of large detuning, |∆ a | ≫ γ, g 0 n , we find where the adiabatic approximation lies on the assumption thath|∆ a | ≫ k B T , as in the single-particle case, and we have neglected the input noise term, assuming the decay rate γ ≪ |∆|. Substituting this value into Eq. (15), the Heisenberg-Langevin equation for the field is given byȧ is the integral of the density of atoms in the electronic ground state, weighted by the cavity spatial-mode function squared. We now assume the bad-cavity limit, namely the cavity field relaxes to the steady state on a much faster time scale than the one in which the density of the atomic medium varies. This limit implies κ ≫ k B T /h, and consistency with the previous assumption imposes |∆ a | ≫ κ ≫ k B T /h. In this limit the dependence of the field on the initial condition is negligible, and its solution is essentially the inhomogeneous one that can be written as Here, we have discarded the input-noise terms, as they are at higher order in the perturbative expansion and we will be dealing with normally-ordered equations, so that two-time correlations of the noise operators vanish, see Eq. (5). We also introduced the operator which is a function of the atom operators in the ground state. Substituting Eq. (18) into the equation for the ground-state field operator we obtaiṅ where C. Discussion Equation (20) shows explicitly the effect of the coupling with the resonator on the atom dynamics: The coupling to the common cavity mode induces a non-linear interaction, which enters in the equation through operator (19). It is useful to consider the average number of photons at steady state n ph = a † a St , which we obtain from Eq. (18) and reads The average number of photons n ph hence depends on the atomic density distribution. On the other hand, it determines the depth of the confining potential, V ≈h|U 0 |n ph , and thus the atomic density distribution. In particular, the confining potential reaches a maximum for the values at which the denominator of Eq. (22) is minimum. From the form of operator (18) one infers that n ph can reach the maximum value when the parameters ∆ c and U 0 have the same sign (the operator Y is positive valued). From Eq. (8) this requires that the detunings ∆ c and ∆ a have equal signs. This property highlights the role of the detuning in the dynamics as control parameters. We now comment on the parameters required for accessing the regime in which the effect of the non-linearity will be important, and its consistency with the derivation we performed. We first review the important assumptions, on which our model is based. Spontaneous decay is neglected over the typical time scales of the system. This imposes that the effective spontaneous scattering rate γ ′ , due to off-resonant excitation of the dipole transition, fulfills the inequality γ ′ ≪ κ. Using that γ ′ ∼ n ph g 2 0 γ/∆ 2 a , where n ph is the mean value of intracavity photons, Eq. (22), then spontaneous emission can be neglected provided that As our model is based on a single-mode cavity, we also require that the detuning between atom and cavity mode is smaller than the free spectral range δω. This reduces to the condition A further important assumption relies on the relaxation time of the cavity field, which has to be much faster than the typical time scale of atomic motion. This can be estimated as κ B T ≪hκ. Finally, in order to have that the non-linear effect on n ph is sufficiently large for a small number of atoms, we have required that U 0 ∼ κ. This condition is however not strictly necessary: strong nonlinear effects can be observed for smaller values of U 0 when the number of atoms is sufficiently large [3].
Let us now estimate the number of intracavity photons which are usually needed, in order to find the atoms in the Mott-insulator state. We consider specifically the case in which overlap (bistability) regions between different Mott zone can be observed. In Sec. IV B we find that this occurs at values of the pump amplitude η ∼ 20κ. This value was evaluated for 50 to 100 atoms in the resonator. Correspondingly, the number of intracavity photons is n ph ∼ 100. From condition (23) we find that γ ≪ 2|∆| a /n ph , where we used U 0 = κ, and which is fulfilled for γ = 2π × 3 MHz and |∆| a = 2π × 10 GHz. Condition (24) is then satisfied when the free spectral range δω is of the order of THz. Once |∆| a is fixed, we find that g 0 ∼ 2 × 0.1 κ/2π MHz. Using the value κ = 2π × 53 MHz [16], this requires g 0 ∼ 2π × 700 MHz, which is presently at the border of experimental reach. However, for smaller values of U 0 , say U 0 ∼ 0.1κ, and for larger numbers of atoms, say N ∼ 1000, the peculiar CQED effects on ultracold atoms we predict in this work could be well observed for parameter regimes of present experiments, see for instance [16].

A. Derivation of the Bose-Hubbard Model
We now derive a Bose-Hubbard type of model for the dynamics of the atoms in their self-sustained potential when the atoms are well localized in the minima of the potential itself. Starting from the assumption that the atoms are in a Mott-insulator state, we decompose the atomic field operator into the operators b † i and b i , which create and annihilate, respectively, atoms at the lowest band of the potential site centered at x = x i , according toΨ wherebyw(x−x i ) are Wannier functions, which are to be determined by solving self-consistently the equations of motion. The commutation rules of operators b i , b † i obey the bosonic commutation relations in the regular Bose-Hubbard model, where the potential is independent of the state of the atoms. We will show that in our case this is not apriori warranted, due to the non-linear dependence of the potential on the atomic density distribution, which gives rise to a non-linear equation for the atomic wave function. However, the bosonic commutation relations are still recovered in a properly defined thermodynamic limit, which we will identify.
We rewrite now Eq. (20) within this Wannier decomposition,ḃ where H (BH) 0 and C are obtained from H 0 and C, respectively, using the Bose-Hubbard decomposition. They read and The coupling matrix elements in Eqs. (27)-(28) read with l = 0, 1 as we keep only on-site and nearestneighbour couplings. In Eq. (28) we introduced the operator where operatorŶ is the Bose-Hubbard decomposition of Y, Eq. (17), after neglecting couplings beyond the nearest neighbors, and takes the formŶ = J 0N + J 1B . In order to determine the Bose-Hubbard Hamiltonian, we now derive an effective Hamiltonian H BH such that C = [b ℓ , H BH ]/h. This is performed in the limit in which we can expand operator F in Eq. (32) in the small quantity J 1 . The details of the derivation are reported in App. A. The Bose-Hubbard model is recovered for a large number of atoms, according to a properly defined thermodynamic limit. We define the thermodynamic limit by letting N and the cavity volume to infinity, keeping finite the number of atoms per potential site. This implies the scaling U 0 ∼ 1/N . Additionally, we impose the scaling η ∼ √ N , which corresponds to keeping the potential depth constant as N increases. This scaling corresponds to ramping up the pump intensity with √ N . The Bose-Hubbard type of Hamiltonian We notice that the coefficients of Hamiltonian (33) are operator-valued, hence imposing a Wannier expansion such that the coefficients depend on the operatorN , namelyw Hence, the commutation relations between the operators b i are not the ones of bosonic operators as in the typical Bose-Hubbard model. Nevertheless, in the thermodynamic limit one finds We therefore perform the Wannier expansion in this thermodynamic limit, consistently with the assumptions made in order to obtain Hamiltonian (33).

B. The Bose-Hubbard Hamiltonian
We rescale now Hamiltonian (33) in units of the strength of the on-site interaction U , which is defined in Eq. (31). The rescaled HamiltonianĤ =Ĥ/U readŝ contains a rescaled chemical potential, while the tunnel parametert is expressed in terms of the coefficient The higher order terms in J 1B , describing long-range interaction, have been neglected. Note that the number of particles is conserved since [N,Ĥ] = 0. We also remark that the term f (N )/N tends to a constant and finite value in the thermodynamic limit. An important physical quantity, which will be useful for the following study, is the depth V of the cavity potential, V =hU 0 n ph with n ph the number of photons in the Bose-Hubbard expansion, Eq. (22). At leading order in the expansion in J 1 it takes the form Hamiltonian (37) and potential (41) are the starting points of our analysis for the determination of the system ground state.
Let us now make some considerations about the system for a fixed number of atoms N . From the form of the potential (41), and in particular from the form of the coefficient ζ, Eq. (40), we observe that for equal signs of the detunings ∆ a and ∆ c one can have that ζ vanishes. This case corresponds to driving the system on resonance, and gives a maximum of the cavity mode potential. This resonance situation occurs for atom numbers N that maximize the photon number, and gives rise to bistability [3], which modifies substantially the properties of the model with respect to the standard Bose-Hubbard one.

IV. DETERMINATION OF THE GROUND-STATE
In this section we determine the ground state of the system for a fixed number of particles. Moreover, we discuss the situation when the number of atoms is fluctuating. Our purpose is to identify the parameter regime in which the atoms are in the Mott-insulator state.
Starting from the assumption that the system is in the Mott-insulator state, we use the strong coupling expansion [34] to verify its validity. In particular, we apply a standard degenerate perturbative calculation in the parametert = t/U , and determine the ground-state energy E M (n 0 ,μ,t) for the Mott state with n 0 particles per site, and the ground-state energies E ± (n 0 ,μ,t) when one particle is added or subtracted to the n 0 -th Mott state. The condition determines the boundariesμ ± (n 0 ) of the n 0 -th Mott phase as a function of the coupling parameter. For µ + (n 0 ) >μ − (n 0 ) the region between the two chemical potentials determines the Mott zone. The Mott state gets unstable as the parameters are varied such that µ + (n 0 ) =μ − (n 0 ) and finallyμ + (n 0 ) <μ − (n 0 ). In this section we determine the boundaries of the Mott state in a diagram, in which we plotμ as a function of relevant parameters. We remark that, in the typical Bose-Hubbard model, when the system exits the Mott phase, then it is in a superfluid state. In our case, this is probably verified in most cases, which we will discuss individually. Finally, the parametert can be controlled by varying the pump amplitude η, which is straightforwardly related to the number of of photons inside the cavity and hence to the height of the potential. Alternatively it can be changed by varying the atom-pump detuning ∆ a and the cavity-pump detuning ∆ c , which enter in the dynamics through the coefficient (40) in the denominator of Eq. (39), and correspond to changing the refractive index of the atomic medium.
In the following we first study the functional dependence of the integrals on the system parameters using the Gaussian ansatz. We then determine numerically the regions of the Mott-insulator state in the diagram where the chemical potentialμ is studied as a function of the pump intensity η.

A. Coefficients in the Gaussian approximation
We determine the boundaries of the Mott-insulator regions using the Gaussian approximation, hence replacing the Wannier functions by Gaussian functions in the integrals (29)- (31). Thus, the Wannier functions are replaced with Gaussian functions such that where σ is the width to be determined. This treatment allows us to identify the dependence of the coefficients on the physical parameters, reproducing with good approximation the results obtained with the Wannier functions in the parameter regimes we discuss in Sec. IV C. In particular, we modify the Gaussian functions in order to fulfill the orthogonality condition, In this way we avoid small, but unphysical contributions. Let K be the number of lattice sites. The width σ of the Gaussian functions is found from the depth V of the cavity-mode potential, Eq. (41). In particular, σ 2 =h/ 2m|V |k. In order to determine the boundaries of the Mott states in the diagram ofμ as a function of η, we determine the coefficients for the three cases (1) N = Kn 0 + 1, (2) N = Kn 0 and (3) N = Kn 0 − 1, and introduce the subscript (i) with i = 1, 2, 3 for the corresponding coefficient. We evaluate the integrals in Eqs. (29)-(31) for these three cases and express them as a function of the dimensionless parameter where E R is the recoil energy. In term of y (i) , they read where a s is the scattering length, ∆ yz is the atomic wave packet transverse width and the sign ± depends on the sign of ∆ a . In the limit J ± 0(i) ≫ |J ± 1(i) | the potential amplitude according to (41) is given by As J ± 0(i) depends on V (i) which, on the other hand, depends itself on J ± 0(i) , the above equations must be solved self-consistently. This is a consequence of the atomdensity dependence on the coupling parameters. In particular, for ∆ a > 0 (atoms at the nodes), J 0(i) → 0 in the strong pumping limit, η → ∞, and the results become independent of the number of atoms. On the other hand, if ∆ a < 0 (atoms at the antinodes) the parameter J 0(i) → 1 for sufficiently large pumping, and the nonlinearity is strongest.

B. Numerical results
In this section we study the regions of the Mottinsulator state as a function of the chemical potential and of the inverse pump amplitude η −1 . The boundaries are determined by numerical evaluation of Eqs. (29)-(31) using the modified Gaussian functions. The atomic parameters we choose correspond to 87 Rb atoms with scattering length a s = 5.77 nm and atomic transition wavelength λ = 830 nm. The optical potential has K lattice size and the transverse width of the atomic wave packet is ∆ y = ∆ z = ∆ yz = 30 nm. We evaluate the "phase diagrams" for K = 50 − 10000 at fixed number of atoms N , scaling N so to keep the atomic density constant. The results for the Mott zones agree over the whole range of values, so in the figures we report the ones obtained for K = 50 for different values of the detunings. a function of the dimensionless parameter κ/η. Interestingly, the extension of the Mott zones seems to decrease roughly as n −1 0 in both cases. We first analyse the case displayed in Fig. 1(a). For ∆ a < 0 the atoms are trapped at the maxima of the intracavity field. Hence, the coupling with the cavity mode is maximum when the confinement is very tight. Here, for large values of the pump intensity (i.e., for small values of κ/η) the Mott zones at different values of n 0 show some overlap. This overlap is a cavity QED effect, in fact Mott states with larger number of atoms per site are favored as they increase the coupling strength to the cavity mode, and thus the depth of the potential. The overlap is only at the border of the boundaries, as atom-atom collisions compete with this effect. In Fig. 1(b) the detuning ∆ a > 0, and the atoms are hence trapped at the nodes (the zeroes) of the intracavity field. Hence, the coupling with the cavity mode is minimum when η → ∞. Indeed, here we observe that for large values of η (small values of κ/η), the Mott zones almost do not overlap. However, for smaller values of η they exhibit an "exotic" behavior: overlap, disappear and reappear.
Further insight is gained in Fig. 2, where we study the depth of the cavity potential as a function of the pump parameters. The curves displayed in Fig. 2(a) correspond to the parameters of the phase diagram in Fig. 1(a). Here, one observes that the potential amplitude increases monotonically as V ∼ η 2 in the parameter regime where the non-linearity is weak. Correspondingly, the width σ (i) of the Wannier functions, giving atomic localization at the minima of the potential, decreases smoothly as σ (i) ∝ |V (i) | −1/4 ∼ 1/ √ η, see Eq. (44). For larger pump strengths, when the non-linearity becomes important, the behaviour is slightly changed. The curves in Fig. 2(b) correspond to the parameters of the phase diagram in Fig. 1(b). Here, one finds that the potential depth increases rapidly where the corresponding Mott zones exhibit a minimum in the value ofμ. Correspondingly, the width σ (i) ∝ |V (i) | −1/4 diminishes rapidly. This behaviour changes at the value of η where the potential gradient increases abruptly. In this regime σ i varies very slowly. This can be understood as a competition between the cavity field, which tends to localize the atoms at the minima, and the atomic quantum fluctuations: When the potential is sufficiently high to trap the atoms within a small fraction of the wavelength, the cavity field is pumped more effectively.  1(b), respectively. When the non-linearity is strong the potential amplitudes differ from the linear situation where |V | ∼ η 2 . The average number of cavity photons is found by multiplying the rescaled potential depths in the plots by the factor fn = Er/|hU0|, which here is fn ≈ 0.006.
We now consider the situation in which the detunings ∆ a and ∆ c have the same sign. In this case the parameter ζ(N ) in Eq. (40) vanishes when the condition J ± 0(i) = ∆ c /U 0 N is fulfilled, whereby 0 < J + 0(i) < 1/2 and 1/2 < J − 0(i) < 1, see Eq. (46). This resonance condition gives rise to bistability, leading to an abrupt change of the potential depth. As a consequence, the Mott state may become unstable. The upper plot of Fig. 3 displays the potential amplitudes V (i) /E r for one atom per site as function of η/κ for U 0 = −κ and ∆ c = −45κ, exhibiting the typical functional behavior of bistability. The lower plot displays the corresponding phase diagram. For the case of one atom per site, the Mott ground-state will suddenly disappear for η ∼ κ when the pumping is adiabatically lowered. Clearly the first "jump" occurs in the potential V (i) (corresponding to the lowest atom density), and comparing the two plots one finds that this takes place exactly when the first Mott zone suddenly ends. The system most likely jumps into a state where higher Bloch bands are populated. In this case the single-band and Gaussian approximations break down.

FIG. 3:
The bistability behaviour of the potential amplitudes V (i) as a function of κ/η (upper plot) and corresponding phase diagramμ − η −1 (lower plot) for ∆c = −45κ, U0 = −κ. At n0 = 1 the Mott region suddenly ends for η ∼ κ, where the corresponding potential V (i) jumps to a lower value. Here, the system most likely is in a state where higher Bloch bands are populated, due to non-adiabatic effects and the small potential depth of this solution.
The overlapping of the Mott zones and the bistability, which we observe in the phase diagram, are novel features when compared with the typical scenario of cold atomic gases trapped by an external potential. Let us first dis-cuss on the existence and uniqueness of the ground-state. When the Mott-insulator state is stable, given the number of atoms N, the ground state is fully determined once the atomic density ρ = N/K is fixed. Outside of the Mott zones we expect superfluidity in the parameter regimes in which there cannot be optical bistability (detunings with opposite signs). In the situation of multiple solutions of the Eqs. (45)-(49), the system will most likely be found in the one solution which minimizes the energy.  fig. 1 (a), where the atomic density ρ has been included as a third axis. Here, the contour lines correspond to a fixed atom density ρ, such that for a given ρ the scaled chemical potential depends on the pumping strength η according to this particular contour line. The dotted contour curves indicate the lines with exactly n0 atoms per site. The projection of the Mott-zones onto thẽ µ − η −1 plane is shown.
A more complete picture of the phase diagram can be obtained by considering the dependence on the atomic density. The strong coupling method for higher orders is cumbersome once the number of added/subtracted particles to the Mott states becomes larger than one. However, the first-order corrections are still easily obtainable for any atom number. In fig. 4 we present schematically the extended phase diagram of fig. 1 (a) where the atomic density has been included as a third axis. This diagram has been obtained by fitting the intermediate lines between the Mott zones, and by verifying that it reproduces the first order calculations for small values of κ/η. We remark that the "overhangs", corresponding to the overlapping zones in Fig. 1(a), constitute the novel feature, which we encounter in this model as compared to the standard Bose-Hubbard model.
When the number of atoms is not fixed [35], the atomic density may take multiple values where the phase diagram exhibits "overhangs". For sufficiently long times, we expect that the system will be found in the number of atoms such that the energy is minimized. This implies also that there may be a competition between a Mott and a superfluid state at two different values of the density, which happen to be at similar energies. Keeping this situation in mind, we restrict our analysis to different and overlapping Mott states, and compare their energy. Figure 5 displays a phase diagram on theμ − η −1 plane, whereby the Mott states with higher energy are plotted on top of the ones with lower energy. We observe that for large pumping strengths the Mott states with a higher number of atoms n 0 have in general a greater energy, while for lower or moderate pumping strengths this is not necessary true. For example, the end of the Mott zone with n 0 = 4 has smaller energy than the corresponding one for the n 0 = 3. This is a pure CQED effect.

C. Validity of the approximations
We now discuss the regime of validity of the calculations, from which we extracted the phase diagrams presented in this section. The derivation of the system coupling parameters relies on the assumption J 0 ≫ |J ± 1 |. The maximum value of the ratio |J ± 1 |/J 0 ≈ 0.056 occurs at |V | ≈ E r /2, hence the nearest-neighbor coupling is at least 17 times smaller than the on site coupling. The expansion to first order in J 1 B of Eq. (32) is motivated for any number of atoms since the perturbative parameter λ ≡ J 1B /J 0N ∼ J 1 /J 0 is strictly smaller than unity.
The values of the chemical potential, as in Eq. (52), are derived from a third-order perturbation expansion in the parametert = t/U and it is expected to break down for larget. We verified that in generalt < 0.25. Moreover, we compared the phase diagrams with the ones obtained by truncating at the second order int, and could verify that they do not differ substantially one from the others. We remark that the perturbation calculations are carried out assuming periodic boundary conditions, while the system here studied has a fixed number of sites, K = 50. We checked the validity of the assumption by comparing the results obtained for different lattice sites, up to 10000, keeping the density fixed.
As it concerns the tight-binding approximation (i.e., only including nearest neighbor couplings), the singleband approximation (i.e., expanding the field operators Ψ(x) and Ψ † (x ′ ) using only the lowest band Wannier functions), these are both related to the regime of validity of the Gaussian approximation. Within this approximation one finds |J 1 /J n | = exp (n 2 − 1) π 2 4y ≫ 1, also indicating validity of the tight-binding approximation in this regime. Figure 6(a) displays the difference ∆−∆ T BA between the width ∆ of the first Bloch band, obtained from diagonalization of the single-particle Hamiltonian in Eq. (7), and the width ∆ T BA evaluated with the tightbinding approximation, as a function of y −1 . Figure 6(b) displays the difference between the coupling parameters obtained by using the Wannier functions and the modified Gaussian functions as a function of y −1 . We note that for values y −1 < 1 the validity of both the tight binding and Gaussian approximation visibly breaks down. This has also been verified by recalculating some of the above phase-diagram using the Wannier functions.  (7), where, in scaled units, there is only a single parameter of interest, namely the dimensionless potential amplitude or equivalently y.

V. CONCLUSIONS
We have shown that ultracold bosonic atoms inside a resonator may form stable insulator-like states, and thus enter the Mott phase, which is sustained by and sustains the cavity potential. The low temperature properties of the system are determined by the competition between the quantum electrodynamic effects and the quantum fluctuations of the atomic matter waves. This competition gives rise to a non-trivial dependence of the regions of stability and of the collective atomic states on the system parameters. Since the cavity potential depends on the state of the atoms, the behavior of the ultracold atomic gas in the cavity differs hence significantly from the one which is encountered in open space. We have derived the Bose-Hubbard Hamiltonian for the cavity confined system, and have shown that the coefficients of this Hamiltonian depend explicitly on the number of atoms. We have determined regions of parameters where the atomic insulator states are stable, predicting the existence of overlapping stability regions for competing Mott states. Bistable behavior is encountered in the vicinity of the shifted cavity resonance, controlled by the pump parameters.
Our theory allows us to determine the state of the atoms when their number is fixed, while for fluctuating, non-fixed atom number, in general, the system will choose the state of minimum energy. This will also take place when an external inhomogeneous potential, such as a harmonic trap potential, is additionally applied to the atoms. In such case we envision a possibility of hysteresis effects in the harmonic potential, when the frequency of this potential is slowly increased and, subsequently, slowly decreased. However, the question how the presence of an inhomogeneous potential will modify the insulator-like states requires further careful studies, since the state of the system depends in a highly nontrivial way on global parameters, which in turn determine the local density of the atoms and the intracavity field. The condition that atoms may affect locally the potential, hence giving rise to phonon-like feature [? ], may be reached in multi-mode resonators, allowing for localized polaritonic excitations [37,38]. Further novel features are expected when fermions are considered instead of bosons. These questions will be tackled in future works. lowship; Consolider Ingenio 2010 "QOIT").

APPENDIX A: DERIVATION OF THE EFFECTIVE HAMILTONIAN
We consider Eq. (26), and rewrite it aṡ where C is defined in Eq. (28) andŶ = J 0N + J 1B . We aim at finding an effective Hamiltonian H BH of the Bose-Hubbard form, such that in some thermodynamic limit to be identified.
We expand now operator C at first order in J 1 , assuming J 1 ≪ J 0 , as it is verified in the Mott-insulator state, using [N ,B] = 0 and where we have introduced the notation At first order in J 1 , we find Let us now consider the commutation relations between the various operators entering this expression. We note that and it is hence of order 1/N . Similarly, the commutator [b ℓ , B] = b ℓ+1 + b ℓ−1 is at higher order in the expansion in 1/N . Henceforth, we can rewrite which is valid at the considered order in the expansion in 1/N . At leading order in 1/N , Eq. (A6) is a differential equation, such that G ′ (J 0N ) = −U 0 η 2 F † (J 0N )F (J 0N ). Using the explicit form of operator F (x), Eq. (32), we find G ′ (x) = −U 0 η 2 /(κ 2 + (∆ c − U 0 x) 2 ) which gives and finally the effective Hamiltonian in Eq. (33).

APPENDIX B: PERTURBATIVE DERIVATION OF THE ZONE BOUNDARIES
We consider Hamiltonian as given in Eq. (37). This Hamiltonian differs from the standard Bose-Hubbard Hamiltonian, as the coefficients depend on the operatorN . We apply now to Eq. (B1) the method of Ref. [34], which allows to determine the region of stability of the Mott-insulator states. The method consists in a perturbative expansion in the parameter t, which is assumed to be small within the parameter regime of interest. In this limit, for large onsite interaction strength U (hard core limit), in the optical lattice the configuration which is energetically favorable has the smallest number of atoms per site. For a lattice of K sites and N = Kn 0 + j atoms, with j < K, there will be either n 0 or n 0 + 1 atoms per site. Clearly, when N = Kn 0 atoms (j = 0), there exists only one possible ground-state, while for N > Kn 0 several ground state configurations exists, and one has to apply degenerate perturbation theory. The ground-state of Hamiltonian (B1) is found after imposing periodic boundary conditions, and diagonalizing operatorB in the momentum representation. At t = 0 the ground-state for N = Kn 0 is given by |Ψ 0 (n 0 ) = |n 0 , n 0 , ..., n 0 , corresponding to n 0 atoms per site, while for N = Kn 0 + j, with j > 0, they are defined by the relation whereÂ † kj = 1 √ K K n=1 e inkj ab † n √n n + 1 (B4) creates one particle in a site starting from the lowest energy states. There is an analogous state for one hole. Here, a = π/k is the distance between neighboring sites, and the wave vector k j = 2πj/Ka, with j = − K 2 , − K 2 + 1, ..., K 2 − 1 (assuming K even for simplicity). The ground state energy is calculated applying perturbation theory in third-order int to this unperturbed basis. Due to symmetry, only zeroth and second-order in the perturbation of t(N )B contribute to the ground-state energies of the Mott-insulator state. For N = Kn 0 one finds U (2) 2Kn 0 (n 0 + 1) while the SF energies for the added particle/hole energies are n 0 (n 0 + 1)(n 0 − 1).

(B6)
Here we have used the subscript (i) corresponding to the three different cases, N = Kn 0 and N = Kn 0 ± 1, see Sec. IV. The limit of stability of the Mott-insulator state is found when the states |Ψ M (n 0 ) and |Ψ ± (n 0 ) are degenerate. The conditions E M (n 0 ) − E + (n 0 ) = 0 and E M (n 0 ) − E − (n 0 ) = 0 determine the boundaries µ ± (n 0 ) of the Mott states in the phase diagramμ −t, thus obtaining the results in Eqs. (52).