3 respect to the relaxation of the condensate. One can assume for any given atomic configuration that the res- onator field has already reached its steady state value, which, by Eq. (4a), is atomic density Ψ † ( x, t )Ψ( x, t ). However, this frequency shift is an operator and couples the equations of motion of the radiation field to those of the atomic field opera- tors in a nonlinear way. The equation of motion of the atomic field operator reads In the atomic part of the equation of motion (2c), in addition to the inert trap potential V$_{ext}$ , the atom-light interaction creates a λ/ 2 periodic optical potential for the atoms, proportional to the dynamical photon number operator a † a . This steady-state mean field provides the optical poten- tial in Eq. (4b). So in the end one only integrates solely Eq. (4b), with α ( t ) adiabatically eliminated and inserted from Eq. (5), instead of the coupled Eqs. (4a) and (4b). It is worth noting that if the time scales of the resonator field and BEC dynamics do not differ that much, some complex coupled solutions can exist [30], which need the simultaneous integration of Eqs. (4a) and (4b). After the adiabatic elimination of the photon field α , the most direct method to numerically calculate the steady state of the BEC wavefunction ϕ ( x ) is the one based on the imaginary-time propagation of Eq. (4b). In real-time the steady-state solution has the time depen- dence: To solve the coupled nonlinear operator equations (2) simultaneously is a hard task. The most convenient ap- proximation is the mean-field approximation, when one separates first the operators into a mean value and to fluctuations around it: α$_{ss}$ = iη Δ$_{C}$ − N 〈 U 〉 + iκ . (5) i ℏ ∂ ∂t Ψ( x, t ) = [ − ℏ $^{2}$Δ 2 m + V$_{ext}$ ( x ) + ℏ a † ( t ) a ( t ) U ( x ) + g Ψ † ( x, t )Ψ( x, t ) $^{]}$Ψ( x, t ) . (2c) ϕ ( x, t ) = ϕ ( x ) e − iµt/ $^{ℏ}$, (6) Ψ( x, t ) = √$_{Nϕ}$$_{(}$$_{x, t}$$_{) +}$ δ Ψ( x, t ) , (3a) a ( t ) = α ( t ) + δa ( t ) . (3b) i d d t α ( t ) = [ − Δ$_{C}$ + N 〈 U 〉 − iκ $^{]}$α ( t ) + iη, (4a) i ℏ ∂ ∂t ϕ ( x, t ) = [ − ℏ $^{2}$Δ 2 m + V$_{ext}$ ( x ) + ℏ | α ( t ) | $^{2}$U ( x ) + g N | ϕ ( x, t ) | 2 ] ϕ ( x, t ) , (4b) C. Mean-field solution with µ/ ℏ the lowest frequency eigenvalue of the nonlinear problem (4b). In imaginary time all fluctuations around the steady-state die out, since they propagate with higher frequencies in real time and consequently vanish faster in imaginary time than the steady-state solution. One just has to renormalize the solution ϕ ( x, t ) after some finite propagation time. Note that, since all quantities on the left hand side of Eq. (4b) are real, the condensate wave- function ϕ ( x ) can also be chosen real. We also note at this point that, as we will see later, in the case of an ef- fective blue detuning Δ$_{C}$ − N 〈 U 〉 > 0, the resonator field heats the atomic motion (some excitations have positive imaginary parts) and there is no steady-state condensate wavefunction at all. However, due to the method of imag- inary time propagation, one can find a BEC wavefunction even in this case, corresponding to a dynamically unsta- ble equilibrium situation. The mean values are c-numbers, defined as ϕ ( x, t ) = N − 1 / $^{2}$〈 Ψ( x, t ) 〉 the so-called condensate wavefunction, which is normalized to unity; and α ( t ) = 〈 a ( t ) 〉 the co- herent part of the cavity field. Consequently, the fluctu- ations have zero mean. The time evolution of the mean values is obtained by substituting (3) into the equations of motion (2) and neglecting all terms containing fluctu- ations. By this way one arrives at a Gross-Pitaevskii like set of equations of the coupled dynamics with the notation 〈 U 〉 ≡ ∫ ϕ ∗ ( x, t ) U ( x ) ϕ ( x, t ) dx . There is a simplification of the numerical problem due to the possible separation of time scales. The time evolu- tion of Eq. (4a) is governed by two characteristic frequen- cies, namely the detuning Δ$_{C}$ and the photon loss rate κ . In Eq. (4b) the characteristic frequency is set by the recoil frequency ω$_{R}$ = ℏ k $^{2}$/ (2 m ). In experiments the latter one is usually several orders of magnitude smaller than the former ones. For example in the Esslinger group experi- ments [1, 3, 4], the parameters | Δ$_{C}$ | ∼ κ ≈ 2 π × 1 MHz, while the recoil frequency is about ω$_{R}$ ≈ 2 π × 4 kHz . In this situation the dynamics of the resonator field relaxes very fast compared to the dynamics of the atomic mo- tion and therefore can be considered instantaneous with Having obtained the steady-state values of the BEC wavefunction and the resonator field amplitude, one can look for the fluctuations of the annihilation (and cre- ation) operators δa ( t ), δ Ψ( x, t ) ( δa † ( t ), δ Ψ † ( x, t )) in lin- ear order. This linear stability analysis corresponds to the Bogoliubov theory of the BEC system [31], and also have an analogy in optomechanics, and also in other non- linear systems, especially in hydrodynamics. Inserting the separation of the field operators (3) into Eqs. (2a) and (2c), and neglecting fluctuations higher than first order, one arrives to D. Fluctuations around the mean-field solution