## Abstract

A cold gas of polarizable particles moving in the optical potential of a standing wave high finesse optical resonator acts as a dynamic refractive index. For a sufficiently strong cavity pump the optical forces generated by the intra cavity field perturb the particles phase space distribution, which shifts the optical resonance frequency and induces a nonlinear optical response. By help of the corresponding Vlasov equation we predict that beyond the known phenomenon of optical bi-stability one finds regions in parameter space, where no stable stationary solution exists. The atom field dynamics then exhibits oscillatory solutions converging to stable limit cycles of the system. The linearized analytical predictions agree well with corresponding numerical solutions of the full time dependent equations and first experimental observation in both cases.

©2011 Optical Society of America

## 1. Introduction

Based on the surprisingly fast progress in the experimental technology of laser cooling and manipulation of dilute atomic gases as well as optical resonator and laser stabilization technology, it is now experimentally possible to confine larger and larger numbers of colder and colder atoms within the mode volume of stabilized very high finesse resonators [1,2]. For slow enough particles the optical dipole force induced by the intracavity field significantly influences the atomic motion even in the dispersive limit at very large detunings, where absorption and spontaneous emission only play a minor role. In this limit the dispersive scattering of the particles simply acts as a dynamic refractive index changing the intracavity field evolution [3–5]. The resulting complex coupled atom-field dynamics leads to a wealth of interesting physical phenomena and applications [6, 7], which can be analysed theoretically using a wide range of models of very different complexity.

In previous work we have developed a new approach based on a classical phase space density description of the particles involving a corresponding Vlasov equation together with classical equations for the field mode amplitudes [8]. This approach is particularly suitable to describe very large particle numbers dispersively coupled to a single or a few driven cavity modes, when spontaneous emission and direct interparticle interaction play only a minor role. Related approaches to cold atom dynamics have put forward also by other groups [9, 10].

It has been known theoretically and experimentally for some time now, that the particles in this case act as a nonlinear optical medium through their motional response to the field, even though they are only weakly excited in the regime of linear polarizability at this frequency and intensity [11, 12]. Even for the conceptually most simple case of a single mode standing wave resonator, the coupled atom-field system can exhibit optical bistability. For a BEC in such a cavity self pulsing solutions were experimentally found [13]. They can be qualitatively well described via a fairly simple two-mode expansion for the atomic mean field dynamics as a nonlinearly coupled oscillator model [14]. This model was shown to be mathematically equivalent to a typical optomechanical setup of a high finesse resonator with a single field mode and a movable mirror at one end. Here we study the appearance of a very similar behavior for a general thermal gas in resonator and show that we can treat both limits of very low and high temperatures in a unified mean field model.

The paper is organized as follows. After presenting our general model in section 2 we use solve the Vlasov equation to find stationary states in section 3, which can give several solutions for a single choice of parameters. Using a linear perturbation analysis we then check the stability regions of these solution in section 4. In section 5 we construct and analyze the long time limit cycle solutions in the unstable regime. These are confirmed by numerical solutions of the time dependent Vlasov equation [Eq. (64)].

## 2. Model and basic equations

We consider a large number *N* → ∞ polarizable particles of mass *m* > 0 inside a standing wave cavity interacting with a single mode of the resonator. The cavity photons induce dipole forces on the particles and the distribution of particles in turn shifts the cavity resonance through their dipole moments, thus giving rise to a nonlinearly coupled system. The dynamical quantities of the mean field type model that will be examined in this work are the phase space (quasi-) distribution *f* (*x*,*v*,*t*), which is the Wigner transform of the reduced one particle density matrix *ρ _{p}*

_{,1}(

*x*,

*x*′,

*t*),

*α*(

*t*) = 〈

*â*(

*t*)〉. As detailed in the appendix, in the limit of large particle numbers, the phase space density obeys Wigner’s equation:

*=*

_{c}*ω*−

_{L}*ω*is the detuning between the pump laser frequency

_{c}*ω*and the cavity resonance frequency

_{L}*ω*.

_{c}*η*> 0 is related to the amplitude of the pump laser and will be referred to as pump strength or pump parameter. Moreover,

*U*

_{0}quantifies the polarisability of the particles,

*L*> 0 denotes the cavity length and

*k*> 0 the wave vector of the corresponding resonator mode. The latter decays at a rate

*κ*. To simplify the problem, let us adopt periodic boundary conditions

*f*(0,

*v*,

*t*) =

*f*(

*L*,

*v*,

*t*) and assume that the optical forces exerted on the particles are weak. Since in this case the spatial modulation of the particle distribution will be small, we may write

*f*

_{0}denotes the homogeneous bulk distribution and

*δ f*designates the small deviation therefrom. Note however that due to there being many particles present, this modulation may nevertheless be accompanied by a considerable shift in the cavity resonance. These assumptions imply that

*δ f*approximately satisfies the linearized Wigner equation

*F*

_{0}(

*v*) :=

*L f*

_{0}(

*v*), we are left with

*v*,

*ϕ*(

*v*,.) satisfies the differential equation of a driven harmonic oscillator with a natural frequency given by 2

*k*|

*v*|.

## 3. Steady states

Let’s turn to the determination of the steady states (*ϕ*
_{0},*α*
_{0}) of Eqs. (5)–(6). One easily finds

*δ*:= Δ

*−*

_{c}*NU*

_{0}/2. Whenever the temperature is nonzero, we can define a function

*F*(

*x*), such that or equivalently ${F}_{0}(v)={v}_{T}^{-1}F(v/{v}_{T})$. To shorten the notation, let us introduce the factor

*I*

_{0}= |

*α*

_{0}|

^{2}. From Eq. (9) one gets a cubic equation for the steady state photon number

*ω*:=

_{R}*h̄k*/2

^{2}*m*denotes the recoil frequency and

*ω*=

_{T}*kv*. Clearly, Eq. (12) can have up to three distinct real solutions. A typical response curve is depicted in Fig. (2). In the interval of effective detuning designated by the letter A the system allows for several possible steady states at once. Which one of these (if any) can be attained? A general investigation of the linear stability properties of steady states will shed some light on that question. Let us therefore establish the necessary and sufficient conditions for a given steady state solution (

_{T}*ϕ*

_{0},

*I*

_{0}) to be stable or unstable.

## 4. Stability analysis

#### 4.1. Small Perturbation Analysis

The starting point of our analysis will be Eqs. (5)-(6), linearized around a steady state solution. Hence let us write

Substituting into the equations of motion and keeping only terms that are linear in the deviations*δϕ*and

*δα*, we obtain

*δϕ*,

*δϕ*

^{*},

*δα*and

*δα*

^{*}. To keep the notation to a minimum, we will not introduce new symbols for the transformed quantities. Introducing the “susceptibility”

*I*

_{1},

*I*

_{2}) denotes an initial value contribution. Let

*D*(

*s*) designate the determinant of the matrix on the left hand side of this last equation, i.e., We will henceforth call the complex valued function

*D*(

*s*) dispersion relation. As explained e.g. in [15], instability is equivalent to the existence of zeros of the dispersion relation with positive real part.

#### 4.2. Limit of zero temperature-BEC

Before we direct our attention to the limiting case of a classical gas, let us now look closer at the opposite limit of a condensate at zero temperature. In that case, the velocity distribution is given by

where*δ*(

*v*) is the dirac delta function. The steady-state photon number can be calculated from Eq. (12) with the replacements:

*J*→ − 1 and

*ω*→ 2

_{T}*ω*. In order to make contact with the method of analysis used in other works, we decompose the solution to Eq. (5) according to

_{R}*ϕ*

_{0}is a solution of the homogeneous equation. This last contribution can be neglected for sufficiently smooth initial conditions due to rapid “phase mixing”, i.e.,

*x*(

*t*) := Re {

*ϕ*

_{+}(

*t*) +

*ϕ*

_{−}(

*t*) } one easily checks that Eq. (5) leads to

*δ*> 0, implies instability.

#### 4.3. Thermal gas - classical limit

Let us now investigate the dispersion relation for a nonzero temperature. Note that as it stands, Eq. (19) together with Eq. (38) are defined only for Re(*s*) > 0. To find the boundary that separates unstable from stable states, we need to find the limiting form of the dispersion relation as Re(*s*) → 0^{+}. To this end we recall the Plemjel formula

*s*=

*γ*+

*iω*and writing sloppily

*χ*(

*ω*) = lim

_{γ}_{→0}+

*χ*(

*s*), we find

*F*

_{0}) consequences of Eqs. (19) and (27), for even if we could analytically find the stability boundary in terms of a function ${I}_{0}^{\text{crit}}={I}_{0}^{\text{crit}}$(

*δ*,

*N*,

*U*

_{0},...), we would still have to find the intersection of this curve with the one specifying the actually possible combinations (

*I*

_{0},

*δ*,...) (i.e., Eq. (12)), and this is obviously a hopelessly complicated task. We will therefore concentrate on the physically most relevant velocity distribution, namely the thermal distribution. Well above the condensation threshold, it is given by the Maxwell-Boltzmann distribution given by

*v*≫

_{T}*v*, such that we may replace in all relevant equations, including Eq. (5). Finally, we should state the condition that the deviation of the particle distribution from the homogeneous be small, namely

_{R}

_{γ}_{→0+}

*D*(

*s*) =

*D*(

_{r}*ω*) +

*iD*(

_{i}*ω*), we find

*x*) denotes Dawson’s integral which is defined as Note that the dispersion relation belonging to a thermal velocity distribution can now be seen (upon substituting Δ →

*δ*and ${U}_{0}^{2}{I}_{0}\to 4{\eta}^{2}$) as being identical to the one considered in [8]. Therefore all the results exactly carry over to the present case: A steady state

*I*

_{0}is unstable if and only if we have for all

*ω*

_{0}such that

#### 4.4. Graphical solution of the stability problem

In order to extract the information included in the general criterion
Eq. (34), we will resort to graphical means. To do so, we note that the two conditions *D _{i}* = 0,

*D*= 0 define a function

_{r}*I*

_{0}is unstable if and only if ${I}_{0}>{I}_{0}^{\text{crit}}$. Hence (Δ, ${I}_{0}^{\text{crit}}(\Delta )$) defines the stability boundary in the parameter-space spanned by (Δ,

*I*

_{0}). Using that we can parametrically plot the stability boundary, i.e. the mapping

*δ*,

*I*

_{0}). Together with the response curve defined by Eq. (12), this allows to analyze the stability properties of the system. Consider again the set of parameters used to plot Fig. (2). Figure (3) depicts the response curve and regions of instability defined by Eqs. (12) and (35) for two different values of the pump parameter. Clearly, there exists an interval of effective detunings, where two distinct steady photon numbers are stable (bistability). Perhaps more surprising however is the finding of an interval where there is no stable steady state at all (Fig. (3b))! To summarize the possible behavior of the system when varying the pump strength, we can state that below a certain value, the response curve is single-valued and lies entirely in the stable domain. Above that value, we find multi-valuedness and a region of bistability (Fig. (3a)). Only upon crossing a second threshold does the completely unstable interval appear (Fig. (3b)).

## 5. Long-time behavior

#### 5.1. Classification

Whenever a steady state (two steady states) is (are) available, we expect the system to relax to it (or one of them) starting from arbitrary initial conditions. How does the system behave for parameters that allow for no stable steady state? To answer this question and to see whether the first guess is correct, we discreetized and solved Eqs. (5)-(6) numerically. Figure (4)a depicts the time averaged photon number (divided by the maximally possible) for different effective detunings for the usual set of parameters. It is clearly visible that this averaged photon number deviates from the steady state photon number only for detunings within the unstable interval, confirming our calculations. To our surprise and for generic sets of parameters, the system exhibits **limit cycle oscillations** within the completely unstable interval. Unlike in the case of a collisionless BEC, where oscillatory solutions depend on the initial conditions, these are true limit cycles independent of the initial conditions. This difference can be attributed to the phenomenon of “phase mixing” of the infinite number of oscillators in the classical case. The “observed” frequencies of the latter versus detuning are shown in figure (4)b. For the chosen parameters, these frequencies are of the order of the thermal frequency. For weak instability, the limit cycle oscillation is almost monochromatic. Figure (5) depicts these oscillations of the photon number for two detunings inside the unstable interval. Associated to these oscillations is a standing density wave formed by the particles.

In Fig. (6) we compare the numerical solution of the reduced model with the solution of the Vlasov equation [Eq. (64)] for the same parameters as in Fig. (5). The agreement between the two models is excellent, confirming the validity of the approximations we introduced.

#### 5.2. Limit cycles frequencies close to threshold

As mentioned, one finds numerically that for parameters such that the steady state with Δ > 0 is *weakly* unstable, in its place there exists a periodic oscillation (limit cycle), which is almost monochromatic. This enables us to construct these cycles in a perturbative manner. Without loss of generality, we may assume that the sought for solution satisfies

*χ*

_{0}and

*χ*

_{1}and ${I}_{n}={I}_{-n}^{*}\in \u2102$. Substituting Eq. (38) into Eq. (6) and using

*e*

^{ix}^{cos(}

^{φ}^{)}= Σ

_{n}*J*(

_{n}*x*)

*e*where

^{inφ}*J*are the Bessel functions of the first kind, we obtain the fourier coefficients of the photon number in terms of

_{n}*ω*,

*χ*

_{0},

*χ*

_{1}as

*χ*

_{0},

*χ*

_{1}to the photon number coefficients

*I*self-consistently, we must evaluate the response of the particles to an applied intensity varying like Σ

_{n}*in the form ϕ = Σ*

_{n}I_{n}e^{inωt}*ϕ*

_{n}*(*

_{n}*v*)

*e*, i.e. we have to solve Note carefully that we cannot simply divide this equation by

^{inωt}*nω*+ 2

*kv*to obtain the solution because of the singularity at

*v*= –

*ω*/2

*k*. Instead, we have to regard

*ϕ*(

_{n}*v*) as a generalized function to obtain

*h*(

_{n}*v*) is arbitrary. As a remedy to this uniqueness problem, one may introduce collisions characterized by a collision frequency and let this frequency go to zero in the end. Alternatively, one may replace which means to consider the intensity as turned on infinitely slowly from the infinitely remote past. In either case the arbitrariness is lifted and we get From this solution we find

*ω*,

*χ*

_{0},

*χ*

_{1}characterizing the limit cycle solution we wish to find. To simplify the problem at hand, let us assume that the cycle be such that Then we may use that for $n\ge 0:{J}_{n}(x)\approx \frac{1}{n!}{\left(\frac{x}{2}\right)}^{n}$, if $\left|x\right|\ll \sqrt{1+n}$ and

*J*

_{−}

*= (−1)*

_{n}

^{n}*J*. Then we have

_{n}*ω*, we solve Eq. (47) with Eq. (51) to lowest order in the small quantity. This yields the dispersion relation for limit cycles just above threshold as

*ω*≪

*ω*we find Note that as Δ

_{T}_{0}(if positive) is monotonously increasing with increasing effective detuning

*δ*, so is the limit cycle frequency.

## 6. Conclusions and outlook

The nonlinear dynamical response of a high-Q optical cavity filled with a cold gas of polarizable particles can be treated by an effective Vlasov type mean field approach over a large range of temperatures from close to zero (BEC) to a thermal gas. In both cases density waves of the medium are strongly coupled to the intra-cavity field dynamics. Besides the well known optical bi-stability effect concurrent with a hysteresis effect, we can clearly identify the regions in parameter space, where no stable stationary solutions exists and time evolution converges to a periodic limit cycle. Its frequency is characteristic for the temperature and gas properties and can be explicitly calculated in the BEC limit or close to threshold for a thermal gas. In both cases such oscillations have been reported experimentally and thoroughly studied [13, 17, 18]. The observations agree well with our model. Interestingly, the effect is predicted to persist at much higher temperatures creating a nonlinear response in a linear medium. To observe it in practice at high temperatures would require sufficiently high intensities and densities or very high Q cavities. Here collisions can no longer be ignored, which can obscure this effect or even lead to a qualitative similar behavior based on a quite different physical mechanism. The collisionless model is in general a valid approximation for laser cooled thermal gases but gets more doubtful in the BEC limit. In the current experiments on this system at ETHZ [7] or Tübingen their influence could largely be neglected so far. In the long run the system thus represents a well controllable and experimentally implementable toy system to study generic nonlinear dynamics in the transition range from classical to quantum mechanics. Adding extra modes or atomic species will allow to extend its scope and complexity in many directions.

## 7. Appendix: quantum and classical mean field limit

Here we show, how the mean field Wigner equation can be obtained from the quantum mechanical equations for *N* → ∞ polarizable particles coupled to one or more resonator modes. When spontaneous emission is neglected, the Heisenberg equation of motion for the particle field operator *$\widehat{\psi}$*(*x*,*t*) reads

*V̂*=

*V̂*(

*x*,

*â*

_{1},...,

*â*) =

_{M}*V̂*

^{†}is the potential at

*x*created by the photons of the

*M*modes. Clearly we have [

*V̂*,

*$\widehat{\psi}$*] = 0. Let us introduce the phase space density operator

*p*=

*mv*is the momentum. It may be noted that 〈

*f̂*(

*x*,

*p*)〉 is the Wigner transform of the one body reduced density matrix

*ρ*

_{p}_{,1}(

*x*,

*x*′). It may be shown from Eq. (55) that the phase space density operator satisfies:

*V̂*(

*x*,

*â*) =

*h̄U*

_{0}

*â*

^{†}

*â*cos

^{2}(

*kx*) and with this one readily finds from Eqs. (57) and (58) that the equation for the phase space density operator thus reads:

*â*of the mode:

*N*≫ 1 can be found by decomposing

*â*=

_{j}*α*+

_{j}*δâ*and

_{j}*f̂*=

*f*+

*δf̂*with 〈

*δâ*〉 = 〈

_{j}*δf̂〉*= 0, taking the average of Eq. (57) and neglecting all correlations. This procedure leads to Wigner’s equation (sometimes called the Quantum Vlasov equation) for the expectation value of

*f̂*(

*x*,

*p*)

*h̄k*, we can develop

*k*≫

_{B}T*h̄*

^{2}

*k*

^{2}/2

*m*, i.e. if the thermal energy (or kinetic energy per particle) is much higher than the recoil energy.

## Acknowledgments

This work was supported by the Austrian Science Fund FWF under grant Nr. I119-N16 and SFB F40 FOQUS. We thank Claus Zimmerman for stimulating discussions and Matteo Cristiani and Jürgen Eschner for helpful comments and early communication of their experimental results.

## References and links

**1. **D. Kruse, M. Ruder, J. Benhelm, C. von Cube, C. Zimmermann, P. W. Courteille, T. Elsässer, B. Nagorny, and A. Hemmerich, “Cold atoms in a high-*Q* ring cavity,” Phys. Rev. A **67**, 051802 (2003). [CrossRef]

**2. **M. Khudaverdyan, W. Alt, I. Dotsenko, T. Kampschulte, K. Lenhard, A. Rauschenbeutel, S. Reick, K. Schörner, A. Widera, and D. Meschede, “Controlled insertion and retrieval of atoms coupled to a high-finesse optical resonator,” New J. Phys. **10**, 073023 (2008). [CrossRef]

**3. **P. Domokos and H. Ritsch, “Mechanical effects of light in optical resonators,” J. Opt. Soc. Am. B **20**, 1098–1130 (2003). [CrossRef]

**4. **P. Horak, G. Hechenblaikner, K. M. Gheri, H. Stecher, and H. Ritsch, “Cavity-induced atom cooling in the strong coupling regime,” Phys. Rev. Lett. **79**, 4974–4977 (1997). [CrossRef]

**5. **P. Domokos, P. Horak, and H. Ritsch, “Semiclassical theory of cavity-assisted atom cooling,” J. Phys. B **34**, 187–198 (2001). [CrossRef]

**6. **J. M. Zhang, F. C. Cui, D. L. Zhou, and W. M. Liu, “Nonlinear dynamics of a cigar-shaped bose-einstein condensate in an optical cavity,” Phys. Rev. A **79**, 033401 (2009). [CrossRef]

**7. **S. Ritter, F. Brennecke, K. Baumann, T. Donner, C. Guerlin, and T. Esslinger, “Dynamical coupling between a Bose–Einstein condensate and a cavity optical lattice,” Appl. Phys. B **95**, 213–218 (2009). [CrossRef]

**8. **T. Grießer, H. Ritsch, M. Hemmerling, and G. Robb, “A Vlasov approach to bunching and selfordering of particles in optical resonators,” Eur. Phys. J. D **58**, 349–368 (2010). [CrossRef]

**9. **J. Javaloyes, M. Perrin, G. L. Lippi, and A. Politi, “Self-generated cooperative light emission induced by atomic recoil,” Phys. Rev. A **70**, 023405 (2004). [CrossRef]

**10. **R. Bach, K. Burnett, M. d’Arcy, and S. Gardiner, “Quantum-mechanical cumulant dynamics near stable periodic orbits in phase space: application to the classical-like dynamics of quantum accelerator modes,” Phys. Rev. A **71**, 33417 (2005). [CrossRef]

**11. **S. Gupta, K. L. Moore, K. W. Murch, and D. M. Stamper-Kurn, “Cavity nonlinear optics at low photon numbers from collective atomic motion,” Phys. Rev. Lett. **99**, 213601 (2007). [CrossRef]

**12. **A. Vukics, W. Niedenzu, and H. Ritsch, “Cavity nonlinear optics with few photons and ultracold quantum particles,” Phys. Rev. A **79**, 013828 (2009). [CrossRef]

**13. **F. Brennecke, S. Ritter, T. Donner, and T. Esslinger, “Cavity optomechanics with a Bose–Einstein condensate,” Science **322**, 235–238 (2008). [CrossRef] [PubMed]

**14. **D. Nagy, P. Domokos, A. Vukics, and H. Ritsch, “Nonlinear quantum dynamics of two BEC modes dispersively coupled by an optical cavity,” Eur. Phys. J. D **55**, 659–668 (2009). [CrossRef]

**15. **D. C. Montgomery, *Theory of the unmagnetized plasma* (Gordon & Breach, 1971).

**16. **F. Brennecke, T. Donner, S. Ritter, T. Bourdel, M. Köhl, and T. Esslinger, “Cavity QED with a Bose–Einstein condensate,” Nature **450**, 268–271 (2007). [CrossRef] [PubMed]

**17. **M. Cristiani and J. Eschner (personal communication, 2011).

**18. **T. Valenzuela, M. Cristiani, H. Gothe, and J. Eschner, “Cold Ytterbium atoms in high-finesse optical cavities: cavity cooling and collective interactions,” in “*Lasers and Electro-Optics 2009 and the European Quantum Electronics Conference. CLEO Europe-EQEC 2009. European Conference on*,” (IEEE, 2009), p. 1. [CrossRef]