8 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 20000 24000 28000 32000 36000 | α ss | 2 (a) full GPE 2 mode model 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 20000 24000 28000 32000 36000 | α ss | 2 (b) 0 0.5 1 1.5 2 2.5 20000 24000 28000 32000 36000 | α ss | 2 Δ$_{C}$ in units of ω$_{R}$ (c) dependence of Eq. (24), with β$_{k}$ chosen real. Eq. (29) can also be solved either by imaginary time propagation or by direct algebraic means using Eq. (24) and the nor- malization condition: β 2 0 + β 2 1 = 1. The equations of motion for the fluctuation operators can be derived analogously to the way as Eqs. (25) was obtained. It is appropriate to separate the trivial time dependence due to the chemical potential via the defini- tion δ ˜ b$_{k}$ = N − 1 / 2 e iµt/ ℏ δb$_{k}$ . With this The matrix K can be thought of being the grand canon- ical Hamiltonian of the system, and by virtue of Eqs. (29) and (24) K$_{k,l}$ β$_{l}$ = 0. Equations (31a) and (31b) form a linear eigenvalue problem for the fluctuations. It is easy to check that the trial function of δa = 0 and δ ˜ b$_{k}$ = β$_{k}$δ ˜ b is a constant solution, i.e. it is a zero mode. It also follows from the normalization of β$_{k}$ that δ ˜ b is an- tihermitian. With the decomposition δ ˜ b$_{k}$ = β$_{k}$ δ ˜ b + γ$_{k}$ δ ˜ c one can arrive at a closed set of equations between the fluctuations δa and δ ˜ c and their hermitian adjoints. with δ = Δ$_{C}$ − N 〈 U 〉 , and i d d t δ ˜ a = ( − δ − iκ ) δa + N α$_{ss}$β$_{k}$U$_{k,l}$ ( δ ˜ b$_{l}$ + δ ˜ $_{b}$† $_{l}$) , (31a) i d d t δ ˜ b$_{k}$ = K$_{k,l}$δ ˜ b$_{l}$ + ( α ∗ $_{ss}$δa + α$_{ss}$δa † ) U$_{k,l}$β$_{l}$, (31b) K$_{k,l}$ = 4 ω$_{R}$ $^{(}$P$_{1}$) k,l − µ $_{ℏ}$δ$_{k,l}$ + | α$_{ss}$ | $^{2}$U$_{k,l}$ . (31c) i d d t δ ˜ a = ( − δ − iκ ) δa + N α$_{ss}$ ( β$_{k}$U$_{k,l}$γ$_{l}$ )( δ ˜ c + δ ˜ c † ) , (32a) i d d t δ ˜ c = ( γ$_{k}$K$_{k,l}$γ$_{l}$ ) δ ˜ c + ( γ$_{k}$U$_{k,l}$β$_{l}$ ) ( α ∗ $_{ss}$δa + α$_{ss}$δa † ) . (32b) FIG. 2: (color online) The mean cavity photon number | α$_{ss}$ | 2 as a function of the cavity detuning Δ$_{C}$. The pa- rameters are: N = 6 × 10 $^{4}$, U$_{0}$ = 0 . 96 ω$_{R}$ , κ = 363 . 9 ω$_{R}$ . The pumping strength is different on the 3 panels: η = (80 . 06 , 283 . 8 , 549 . 5) ω$_{R}$ for panels a, b, c, respectively. A. The mean-field solution For a comparison between the full GPE solution and the optomechanical model we numerically solve the Gross-Pitaevskii Equation (4b) on a 200 point grid with imaginary time propagation with the steady-state value of the mean radiation field amplitude (5) and compare it to the solution of (29) (also in imaginary time). The results for the mean photon number are plotted in Fig. 2. The bistable behavior is nicely exhibited in these fig- ures. In panel a) the pumping strength η = 80 . 06 ω$_{R}$ is below a threshold value, the resonance curve gives a unique solution for all detunings Δ$_{C}$. Notice, that in Eqs. (4a), or (5) the effective detuning of the resonator field is δ = Δ$_{C}$ − N 〈 U 〉 , i.e., the light shift further de- tunes the cavity. On the left side of the resonance, the effective detuning δ is negative, the cavity is detuned to the red. We will see in the next subsection that this situ- ation corresponds to cavity cooling, where the imaginary part of the fluctuation spectrum of the combined atom- photon system is negative (or nonpositive for the many mode system). Fluctuations decay to zero. On the right side of the resonance, where δ is positive, the imaginary Note, that in the optomechanical model s-wave scatter- ing can not be included in general, because it populates states with q = 0, which are disregarded in this model. It is a possibility to include some remaining effects of the atom-atom collision by coupling only the two states under consideration, however, it would result in a nonlo- cal and unphysical interaction in coordinate space. Or, in order to extend the optomechanical model, although still lacking the full consideration of s-wave scattering ex- plicitly, the higher quasi-momentum states can be taken into account either by considering a (close to, but non-) plane wave condensate wavefunction [25], or in a two- fluid model [38]. In the following we summarize and discuss our results based on the numerical solutions of Eqs. (4) and (7), and compare them with those of Eqs. (29) and (32) and to the findings of the Esslinger group [3, 4]. III. RESULTS