Band structure loops and multistability in cavity-QED
BBand structure loops and multistability in cavity-QED
B. Prasanna Venkatesh, J. Larson, and D. H. J. O’Dell Department of Physics and Astronomy, McMaster University,1280 Main St. W., Hamilton, ON, L8S 4M1, Canada Department of Physics, Stockholm University, Se-106 91, Stockholm, Sweden
We calculate the band structure of ultracold atoms located inside a laser-driven optical cavity. Forparameters where the atom-cavity system exhibits bistability, the atomic band structure developsloop structures akin to the ones predicted for Bose-Einstein condensates in ordinary (non-cavity)optical lattices. However, in our case the nonlinearity derives from the cavity back-action ratherthan from direct interatomic interactions. We find both bi- and tri-stable regimes associated withthe lowest band, and show that the multistability we observe can be analyzed in terms of swallowtailcatastrophes. Dynamic and energetic stability of the mean-field solutions is also discussed, and weshow that the bistable solutions have, as expected, one unstable and two stable branches. Thepresence of loops in the atomic band structure has important implications for proposals concerningBloch oscillations of atoms inside optical cavities [Peden et al ., Phys. Rev. A , 043803 (2009),Prasanna Venkatesh et al ., Phys. Rev. A , 063834 (2009)]. PACS numbers: 42.50.Pq, 37.10.Jk, 37.10.Vz, 37.30.+i, 06.20.-f
I. INTRODUCTION
Optical bistability is a manifestation of nonlinearity inoptical systems which is well known in the laser physicscommunity [1, 2] (see also [3] and [4] and referencestherein). It describes a situation in which there are twopossible stable output light intensities for a single inputintensity, and occurs when an optical medium with a non-linear refractive index is placed inside an optical cavityformed from two mirrors. The bistable behaviour resultsfrom the combination of the nonlinearity of the mediumwith the action of the feedback provided by the mirrors.A new addition to the family of systems displaying op-tical bistability has recently been demonstrated in ex-periments performed by the ultracold atom groups atBerkeley [5] and the ETH [6] who found optical bista-bility in systems comprising of vapors of ultracold atomstrapped inside optical cavities which are driven by laserlight. The atomic vapor acts as a dielectric medium and,despite being tenuous, can significantly perturb the lightin a cavity with a small mode volume and high finesse ifthe cooperativity C N is in the regime C N ≡ N g / κγ > N is the number of atoms, g is the single photonRabi frequency, 2 γ is the atomic spontaneous emissionrate in free space, and 2 κ is the cavity energy dampingrate. The perturbation of the light by the atoms is non-linear and distorts the cavity lineshape away from beinga lorentzian which is symmetric about the resonance fre-quency into one with an asymmetric shape. For largeenough cooperativity the lineshape becomes folded over(see Fig. 2 below), so that for a certain range of frequen-cies there are three possible output light intensities (twostable, one unstable) from the cavity for a single inputintensity. The experiments [5] and [6] exhibited this op-tical bistability as a hysteresis effect seen by chirping thelaser frequency through the cavity resonance from aboveand below the resonance: a sudden jump in the inten-sity of light transmitted through the cavity was observed which occurred at two different frequencies, dependingupon the direction of the chirp.An important difference between traditional laser sys-tems and the ultracold atom experiments [5, 6] is theorigin of the nonlinearity. In the former case the nonlin-earity of the medium occurs in its polarization response,i.e. it arises from the internal degrees of freedom of theatoms. By contrast, in the ultracold atom experimentsthe detuning of the cavity from atomic resonance waslarge enough that the polarization response was in thelinear regime. The nonlinearity was instead due to the re-sponse of the centre-of-mass wave function of the atoms:the atoms re-arrange their position distribution accord-ing to the balance between the dipole force applied bythe intracavity light field (which forms a periodic lattice)and their zero-point energy. As a consequence, the depthof the optical lattice that forms inside the cavity in ex-periments like [5] and [6] is not fixed purely by the drivelaser intensity, as is the case in standard optical latticesmade by interfering laser beams in free space. Rather,when C N > a r X i v : . [ qu a n t - ph ] J un tices [16]. Of special interest to us here are the curiousswallowtail loops that occur in the band structure (en-ergy versus quasi-momentum) of atomic Bose-Einsteincondensates (BECs) in one-dimensional optical lattices.These loops have been studied theoretically by a num-ber of groups [17–21] in order to explain the breakdownin superfluidity observed in experiments where a BECflows through a lattice [22, 23]. The loops correspond tomultiple solutions for the atomic wave function within asingle band for a certain range of quasi-momenta. Theycan occur either around the boundaries of the Brillouinzone or the center, depending upon the band and the signof the interactions. They manifest themselves physicallyvia a dynamical instability that destroys the superflow.However, the nonlinearity responsible for the swallow-tail loops in the free-space lattices is provided by theinteratomic interactions, which become important at thedensities required for Bose-Einstein condensation. Theloops occur when the strength of the interactions is abovea critical value [17, 18], and therefore non-interactingatoms in an optical lattice do not display these instabili-ties. Our purpose in this paper is to investigate whetherthe cavity feedback nonlinearity associated with opticalbistability in ultracold atoms can also lead to loops in theatomic band structure. As we shall see, the answer to thisquestion is in general affirmative, and so band structureloops appear to be a robust phenomenon which appearwhatever the source nonlinearity, although the structureand location of the loops does depend on the details ofthe nonlinearity.One consequence of loops in the atomic band structureis a hysteresis effect [19] if the quasi-momentum is sweptthrough the band, and a consequent loss of adiabatic-ity even if the quasi-momentum is swept infinitely slowly[24]. Recent experiments on a two-site lattice have con-firmed this scenario [25]. These effects have implicationsfor experiments that use Bloch oscillations of cold atomsin optical lattices for high precision measurements, for ex-ample to determine the fine structure constant α [26], tomeasure gravity [27, 28], or to investigate Casimir-Polderforces [29]. The band structure hysteresis is reminiscentof the hysteresis described above in the context of theoptical bistability experiments [5] and [6]. Indeed, as weshall show in this paper, for atoms in an optical cavitythe two effects are different sides of the same coin, onebeing seen in the light and the other in the atoms. Ofparticular relevance, therefore, are two recent proposals[30, 31] concerning Bloch oscillations of atoms inside op-tical cavities that rely upon the modification of the trans-mitted/reflected light arising from the feedback nonlin-earity as a non-destructive measurement of the Bloch os-cillations. The presence of loops will severely disrupt theBloch oscillations. In the case where the nonlinearityhas its origin in atom-atom interactions the loops canbe eliminated by, e.g., using spin polarized fermions [29]or by reducing the interactions via a Feshbach resonance[32]. However, when the nonlinearity is due to the feed-back nonlinearity these methods do not apply, and one of our aims here is see if there are regimes where the feed-back nonlinearity is present and leads to a modificationof the light, but remains below a critical value neededto form loops in the atomic band structure. In order toallow the origin of the nonlinearity to be clearly identi-fied, we shall only consider non-interacting atoms in thispaper, but our calculations can be easily extended to in-clude interactions.The plan for the rest of the paper is as follows. In Sec-tion II we give a brief description of the system and in-troduce the mean-field hamiltonian describing cold atomsdispersively interacting with the single mode of a cavity.In Section III we derive a reduced hamiltonian describ-ing the nonlinear evolution for the atomic field after theadiabatic elimination of the light field. In Section IV wecalculate the band structure by two different methodsthat solve for the spatially extended eigenstates (Blochstates) of the reduced hamiltonian, and show that thetwo methods give the same results. We find that forcertain parameter regimes, the energy as a function ofquasi-momentum develops loop structures. In Section Vwe recall the optical bistability in atom-cavity systemsdiscussed by [5, 6] and make the connection between theloop dispersions and bistability. In Section VI we developan analytical estimate for the critical pump strength nec-essary to generate loops and in Section VII we test thisresult by illustrating how the band structure depends onlaser detuning and laser intensity, i.e. the birth and deathof loops as these parameters are varied. In Section VIII,we show that for quasi-momentum q (cid:54) = 0 the system canexhibit tristability (five solutions for the steady state cav-ity photon number).A convenient mathematical framework for describingbifurcations of a system whereby there is a sudden changein the number of solutions as a parameter is smoothlychanged is catastrophe theory [33–37]. Catastrophe the-ory was applied to the problem of optical bistability inthe late 1970s by a number of authors including Gilmoreand Narducci [38] and Agarwal and Carmichael [ ? ]. InSection IX we relate the problem of atom-cavity multi-stability to catastrophe theory and show that the sys-tem under consideration can be described by swallow-tail catastrophes organized by an (unobserved) butter-fly catastrophe and use this to understand multistability.This is followed in Section X by a discussion of the stabil-ity of these solutions and finally we conclude in SectionXI. There is also an appendix which summarizes someconcepts of catastrophe theory. II. THE HAMILTONIAN
Consider a gas of N ultracold atoms inside a Fabry-Perot optical cavity. The atoms interact quasi-resonantlywith a single cavity mode of the electromagnetic fieldof frequency ω c , and it varies along the cavity axis ascos( k c z ), where k c = ω c /c . This cavity mode is exter-nally pumped by a laser of frequency ω p . The relevantfrequency relations can be characterized with the two de-tunings ∆ c ≡ ω p − ω c , (1)∆ a ≡ ω p − ω a , (2)where ω a is the atomic transition frequency. In the dis-persive regime, the occupation of the excited atomic stateis vanishingly small and it can be adiabatically elim-inated. A one-dimensional hamiltonian for the atom-cavity system in the dispersive regime can then be writ-ten as [9, 10] H = − (cid:126) ∆ c ˆ a † ˆ a + i (cid:126) η (cid:0) ˆ a † − ˆ a (cid:1) + (cid:90) d z ˆΨ † (cid:20) − (cid:126) M ∂ ∂z + (cid:126) g ∆ a ˆ a † ˆ a cos ( k c z ) (cid:21) ˆΨ , (3)where ˆΨ( z, t ) and ˆ a ( t ) are the field operators for theatoms and the cavity photons, respectively. We havewritten the hamiltonian in a frame rotating with thepump laser frequency ω p , which leads to the appearanceof the two detunings. The first term is just the free fieldevolution of the cavity mode. The second term representsthe laser coherently pumping the cavity at rate η , and thethird term describes the atomic part of the hamiltonian.It consists of a kinetic energy part and a potential energypart. The potential energy term can either be understoodas the atom moving in a periodic potential with ampli-tude ( (cid:126) g / ∆ a ) (cid:104) ˆ a † ˆ a (cid:105) , or, if combined with the first termin the hamiltonian, as a shift in the resonance frequencyof the cavity due to the coupling between the atom andthe field (see also Eq. (5) below).The Heisenberg equations of motion for the atom andphoton field operators can be obtained from the abovehamiltonian. In this paper we are interested in proper-ties at a mean-field level, where operators are replaced bytheir expectation values and correlations between prod-ucts of operators are neglected. In other words, thelight is assumed to be in a coherent state with amplitude α ( t ) = (cid:104) ˆ a (cid:105) , and the atoms are assumed to all share thesame single-particle wave function ψ ( x, t ) = (cid:104) ˆΨ (cid:105) / √ N .From the Heisenberg equations these amplitudes obeythe following coupled equations of motion [40] i ∂ψ∂t = (cid:18) − ∂ ∂x + U n ph cos ( x ) (cid:19) ψ, (4)d α d t = (cid:18) i ∆ c − iN U (cid:90) d x cos ( x ) | ψ ( x ) | − κ (cid:19) α + η, (5)where we have scaled energies by the recoil energy E R = (cid:126) k / M , and time by (cid:126) /E R . The depth of the peri-odic potential seen by the atoms is then given by U n ph ,where n ph ≡ | α | is the mean number of photons in thecavity, and U is defined as U ≡ (cid:126) g / (∆ a E R ) . (6)The damping rate κ (in units of E R ) of the amplitudeof the light field in the cavity has been added into the equation of motion in order to account for a markoviancoupling between the cavity mode and a zero tempera-ture bath. We have also introduced the dimensionlesslength x ≡ k c z . In the above equations, any fluctuationsinduced by the reservoirs have been neglected. This isjustified when considering relatively large quantum num-bers, for corrections see reference [12].In this paper we are interested in the band structureand this requires the quasi-momentum to be a good quan-tum number. Physically, this implies that we are study-ing delocalized atomic wave functions which cover manylattice sites (this is closer to the situation realized in theexperiment [6] than that realized in [5] where the atomswere localized to just a few sites). We shall thereforeexpand the wave function ψ ( x, t ) in Bloch waves, as willbe detailed in subsequent sections. It thus makes senseto normalize ψ ( x, t ) over one period of the potential as (cid:82) π | ψ | dx = π , and also evaluate the integrals appearingin the above equations over one period. In particular, theintegral in Eq. (5) which determines the coupling betweenthe atoms and the light will be defined by (cid:104) cos ( x ) (cid:105) ≡ π (cid:90) π | ψ ( x ) | cos ( x )d x. (7) III. THE REDUCED HAMILTONIAN
In this Section we shall eliminate the optical degreesof freedom from the hamiltonian in order to obtain adescription only in terms of the atoms. This results in anonlinear Schr¨odinger equation and an energy functionalwe shall refer to as the reduced hamiltonian.Setting d α/ d t = 0 in Eq. (5) gives the steady statelight amplitude in the cavity, α = ηκ
11 + i ∆ c − NU (cid:104) cos ( x ) (cid:105) κ (8)which leads to the following equation for the steady statephoton number n ph = η κ + (∆ c − N U (cid:104) cos ( x ) (cid:105) ) . (9)In fact, this expression also holds to high accuracy inmany time-dependent situations because κ is typically fargreater than any frequency associated with the evolutionof the external degrees of freedom of ultracold atoms.The light field is then slaved to the atoms and “instan-taneously” adjusts itself to any change in ψ ( x, t ). Thesteady state solution can then be substituted back intoEq. (4) to give us a single, closed, nonlinear Schr¨odingerequation for the atomic wave function i ∂ψ∂t = (cid:32) − ∂ ∂x + U η cos ( x ) κ + (∆ c − N U (cid:104) cos ( x ) (cid:105) ) (cid:33) ψ . (10)The stationary solution ψ ( x, t ) = ψ ( x ) e − iµt of thisequation leads to an expression for the eigenvalue µ ofEq. (10), µ [ ψ, ψ ∗ ] = 1 π (cid:90) π dx (cid:32)(cid:12)(cid:12)(cid:12)(cid:12) dψdx (cid:12)(cid:12)(cid:12)(cid:12) + U n ph cos ( x ) | ψ ( x ) | (cid:33) . (11)If the Schr¨odinger equation (10) were linear, then theeigenvalue µ would be the energy of the atoms in state ψ . Furthermore, the functional (11) would serve as theenergy functional from which this Schr¨odinger equationcould be obtained as an equation of motion using Hamil-ton’s equation i (cid:126) ∂ψ∂t = δE [ ψ, ψ ∗ ] δψ ∗ . (12)i.e. E [ ψ, ψ ∗ ] = µ [ ψ, ψ ∗ ] for a linear equation. However,this is not the case here. Instead, the eigenvalue µ is thechemical potential which is related to the true energy E via µ = ∂E/∂N (a prominent example that illustratesthis distinction is the Gross-Pitaevskii equation and itsenergy functional [21]). Using this fact, one can showthat the true energy functional from which the equationof motion (10) can be derived is in fact [10, 13] E [ ψ ] = Nπ (cid:90) π dx (cid:12)(cid:12)(cid:12)(cid:12) dψdx (cid:12)(cid:12)(cid:12)(cid:12) − η κ arctan (cid:32) ∆ c − NU π (cid:82) π d x | ψ | cos ( x ) κ (cid:33) . (13)as can be verified by applying Hamilton’s equation (12).We shall refer to this functional as the reduced hamilto-nian. The first term represents the kinetic energy. Thesecond term is an atom-light coupling dependent termthat can be interpreted as follows. The phase shift of thesteady state light field inside the cavity relative to thepump laser is, from Eq. (8), φ = arctan Im( α )Re ( α ) = arctan (cid:32) ∆ c − NU π (cid:104) cos ( x ) (cid:105) κ (cid:33) . (14)This allows us to rewrite the reduced hamiltonian as E [ ψ ] = Nπ (cid:90) π d x (cid:12)(cid:12)(cid:12)(cid:12) d ψ d x (cid:12)(cid:12)(cid:12)(cid:12) − I ph φ , (15)where I ph = η /κ is the incident photon current fromthe pump laser. Note that this hamiltonian looks similarin form to the effective quantum two-mode hamiltoniansobtained in the optomechanical limit in [12] and [13], andalso the ones describing bistability in [10]. IV. BAND STRUCTURE
From now on we specialize to Bloch wave solutions.We begin by describing two methods for calculating theBloch waves and their energies. Agreement between thetwo methods will be demonstrated.
A. Method 1: Energy extremization
The first method, which adapts that detailed in [18] fora static optical lattice, hinges on the observation that thepotential in the Schr¨odinger Eq. (10) is periodic with pe-riod π (in dimensionless units). Despite the nonlinearity,the Bloch theorem [31, 42, 43] applies so that the eigen-functions can be written as the product of a plane wavewith wavenumber q , called the quasi-momentum, and anamplitude U q ( x ) which is periodic with the period of thelattice ψ q ( x ) = e iqx U q ( x ) , (16) U q ( x + π ) = U q ( x ) . For the linear problem, the energies E ( q ) are arrangedinto bands separated by gaps. In the so-called reducedzone scheme for the band structure, q lies in the firstBrillouin zone − ≤ q <
1, and the band energies arefolded back into the first Brillouin zone.The periodicity of U q ( x ) implies that it can be ex-panded as a Fourier series. The Bloch wave can thereforebe written ψ q ( x ) = e iqx (cid:88) n a q,n e i nx . (17)The notation for the n th expansion coefficient a q,n re-flects the fact that it depends on the chosen value of q .This expansion can now be substituted into the reducedhamiltonian Eq. (13), and the resulting function numer-ically extremized with respect to the parameters a q,n ,maintaining the normalization of ψ q ( x ) throughout theprocedure. We take the parameters a q,n to be real forthe same reasoning as given in [18]. The Fourier series isterminated at some n = R , which is determined by theconvergence of the energy eigenvalues as R is varied.In Fig. 1 we plot the low lying part of energy spec-trum E [ ψ q ] as a function of quasi-momentum resultingfrom the extremization. The values of κ and N are veryclose to the values used in the experiment with rubid-ium atoms described in [6], and the other parametersare chosen so as to exhibit the interesting behavior tobe discussed in the rest of the paper. The two panels ofFig. 1 differ only in the value of the pump-cavity detuning∆ c . Fig. 1a shows a result familiar from linear problemsinvolving quantum particles in periodic potentials, butFig. 1b shows a very different story: there are two otherbranches that together form a looped structure that is amanifestation of the nonlinearity of the reduced hamilto-nian. As will be discussed more below, the loop shown in −1 −0.5 0 0.5 10.20.40.60.811.21.41.6 quasi-momentum q (units of k c ) E [ ψ ] / N ( un i t s o f E R ) −1 −0.5 0 0.5 10.20.30.40.50.60.70.80.91 quasi-momentum q (units of k c ) E [ ψ ] / N ( un i t s o f E R ) (a)(b) FIG. 1. Energy loops in the first band. The curves were obtained by extremizing the reduced hamiltonian (13). For both (a)and (b) the laser pumping η = 909 . ω R , the cavity decay κ = 350 ω R , the atom-light coupling U = ω R , and number of atoms N = 10 . In (a) the pump-cavity detuning was ∆ c = 1350 ω R , which gives no loops, and in (b) it was ∆ c = 3140 ω R whichgives a loop symmetric about the band center as shown. As explained in the text, the number of photons in the cavity andhence the lattice depth change with q . For example, in (a) at q = 0 we have n ph = 0 .
06 and at q = 1 we have n ph = 0 .
68. In(b) at q = 0 we have for the lowest branch n ph = 4 .
13, for the middle branch n ph = 0 .
28, and for the upper branch n ph = 2 . q = 1 we have n ph = 1 .
08. At the point where the middle and upper branches meet we have n ph = 0 . Fig. 1b belongs to the first band (in particular, it does notcorrespond to the second band because it only extendsover part of the first Brillouin zone). Looped structureshave been found before for BECs in static optical lat-tices [17, 18]. We will come back to the similarities and differences between our results and [17, 18] in the nextsection.It is important to appreciate that, by virtue of thenonlinearity of the problem, the lattice depth n ph U isdifferent at each value of the quasi-momentum shown inFig. 1. Furthermore, the lattice depth is different foreach of the curves even at the same values of q (except atdegeneracies). In this sense, the band structure plots wedisplay in this paper are non-standard because for staticlattices the depth of the lattice is fixed for all values of q . In order to obtain the lattice depth at each point ona curve E [ ψ q ] in Fig. 1, one should take the wave func-tion that extremizes the energy functional (13) at thatpoint and enter it into the integral (cid:104) cos ( x ) (cid:105) appearingon the righthand side of Eq. (9) for the photon number.This change in lattice depth with detuning is reported inFig. 2, but for the reader’s convenience we have includedsome values at particular points q in the caption of Fig. 1.The fact that, depending on the detuning ∆ c , thereare either one or more steady state photon numbers forthe cavity reminds us of the optical bistability observedin [5, 6]. There, as the cavity driving laser detuning wasswept through the resonance, bistability was observed forcertain pump strengths η . The bistability derives fromquantum effects [41] in the sense that it is due to changesin the atomic center-of-mass wave function. It is distinctfrom standard semi-classical optical bistability [4]. Theconnection between the loops in the atomic band struc-ture and optical bistability will be examined in detail inSection V. To complete this section we look at an alter-native method for determining the eigenfunctions of theeffective hamiltonian which makes the connection withbistability more transparent. B. Method 2: Self-consistency equation
In the second method, the strategy is to solve Eq. (9)directly for the steady-state photon number. This equa-tion contains n ph both explicitly on the left hand sideand implicitly on the right hand side through the atomicwave function ψ q ( x, n ph ) appearing in the integrand ofthe integral (cid:104) cos ( x ) (cid:105) = 1 π (cid:90) π cos ( x ) | ψ q ( x, n ph ) | d x. (18)The photon number is not a parameter set from the out-side (e.g. by the pump laser) but must be solved forself-consistently. In our notation for the wave functionwe have therefore included n ph in the list of argumentsrather than the list of parameters. The dependence of ψ q ( x, n ph ) upon the number of photons is because thedepth of the lattice in which the atoms sit is given by U n ph , as can be seen directly from the Schr¨odingerequation (4). Therefore, Eq. (9) must be solved self-consistently for n ph , and we will often refer to it as the“self-consistency equation”. As mentioned in the intro-duction, the physical mechanism that gives rise to thefeedback between the atoms and the light stems for theatom-light coupling, Eq. (7) [or Eq. (18)]. The atoms’spatial distribution controls the phase shift suffered bythe light on traversing the cavity, and hence the cavityresonance condition, and therefore the amplitude of the light in the cavity. At the same time, the amplitude ofthe light determines the depth of the lattice which influ-ences the atomic wave function.One class of solutions to the self-consistency problemis provided by the Mathieu functions [44]. In fact, theseprovide exact solutions because the Schr¨odinger equation(4) for a particle in a sinusoidal potential is nothing butthe Mathieu equation. Despite the fact that the ampli-tude of the sinusoidal potential in (4) itself depends onthe solution ψ q ( x, n ph ) of the equation, this amplitudeevaluates to a constant because ψ q ( x, n ph ) appears un-der the integral given by Eq. (7). All that is necessary isthat the wave function that enters into the integral is thesame as the wave function appearing in the rest of theSchr¨odinger equation. This is the case precisely when theself-consistency equation (9) is satisfied. This leads us toa very important point: by virtue of the self-consistencyequation (9), the photon number n ph is completely equiv-alent , in the sense of the information it carries, to thewave function ψ . Said another way, the Mathieu func-tions are specified by only two quantities: the value of thequasi-momentum and the depth of the potential. Thus,once we set the quasi-momentum, which is a parameter,the wave function is entirely determined by U n ph , where U is also a parameter. We shall frequently take advan-tage of the equivalence between ψ and n ph in the restof this paper because it allows us to replace the wavefunction by a single number n ph .We shall denote the Mathieu functions by χ m,q,n ph ( x ),where m is the band index. They satisfy Mathieu’s equa-tion which in our problem takes the form (cid:32) ∂ ∂x + U n ph cos ( x ) (cid:33) χ m,q,n ph ( x )= µ m,q,n ph χ m,q,n ph ( x ) . (19)Our notation for the Mathieu functions indicates the fullparametric dependence with the exception of the atom-light coupling strength U . Note that in the Mathieufunctions we list n ph as a parameter, like q , because thatis the role it plays in the Mathieu equation. We thereforehave that ψ m,q ( x, n ph ) = χ m,q,n ph ( x ) . (20)Mathieu’s functions can, of course, be written in Blochform. In this paper we focus on the first band and so weshall suppress the band index from now on. We thereforehave χ q,n ph ( x ) = exp( iqx ) U q,n ph ( x ). Substituting thisBloch form into the time-dependent Schr¨odinger equa-tion (4) as ψ q ( x, n ph , t ) = χ q,n ph ( x ) exp( − iµt ) one ob-tains (cid:18) ∂∂x + iq (cid:19) U q,n ph ( x ) + U n ph cos ( x ) U q,n ph ( x )= µ q,n ph U q,n ph ( x ) . (21)This equation can either be solved numerically fromscratch, or a package such as Mathematica can be usedwhich gives the Mathieu functions for each value of q ,and U n ph . For a particular choice of q , these Mathieufunctions can now be used in (9) to find the value(s) of n ph that give self-consistency.There are two main differences between method 1 andmethod 2 described above. Firstly, method 1 is a vari-ational calculation, whereas method 2 exploits the def-inition of the steady state photon number to obtain asingle nonlinear integral equation (9) which must be sat-isfied by the atomic wave function. Secondly, in method2 we can explicitly set the band index whereas the vari-ational wave function used in method 1 is rather moregeneral. In spite of these differences, we find that thetwo methods are in complete agreement (to within nu-merical accuracy) for all the parameter ranges we wereable to test (for very deep lattices method 1 becomes un-workable because a very large number of terms in theFourier expansion are required). We have also checkedthat the two methods agree for higher bands.It may at first seem surprising that two such seem-ingly different methods are equivalent. We emphasizethat both stem from the non-linear Schr¨odinger equa-tion (10), which is just Mathieu’s equation with a po-tential depth which is set self-consistently. Method 1minimizes an energy functional that is derived from thisnon-linear Schr¨odinger equation. In principle, one couldminimize ansatze other than the Bloch functions we usehere (e.g. localized functions), but these would not thensatisfy the original non-linear Schr¨odinger equation (10)exactly. Method 2 is based upon the observation thatwave functions that satisfy the self-consistency equation(9) are precisely the required solutions of the non-linearSchr¨odinger equation (10) providing we restrict ourselvesto solutions which are Mathieu functions. Again, onecould find other types of solution, but these would notsatisfy Eq. (10).An interesting question to ask is whether the nonlin-earity of our problem mixes the linear bands, so that, forexample the self-consistent first band is a superpositionof Mathieu functions from different bands. This is notwhat we find. Instead, the solutions we have found allcorrespond to being from one or other band, but not asuperposition. Method 2, in particular, allows us to ex-plicitly set the band index and we are therefore able toidentify all three curves shown in Fig. 1 as belonging tothe first band (we have also checked that the Mathieufunctions corresponding to all three curves are orthog-onal to the Mathieu functions for the second band forthe same three lattice depths). Actually, we do not findloops in higher bands for the parameter values consid-ered in Fig. 1. Although in this paper we restrict ourattention to the first band, we have found that we canhave loops in the higher bands as well. The calculation ofenergy dispersions using the self-consistent method is nu-merically less demanding and so we will continue to usethe latter for the remainder of this paper. In the nextsection we discuss optical bistability and its relation toloops in the band structure. −10 −5 0 5 10024681012 (∆ c - N U / /κ P h o t o nnu m b e r n ph η = 0.5 η cr η = 1.5 η cr η = 2.5 η cr η = 3.5 η cr FIG. 2. (Color online) The steady state photon number insidethe cavity as a function of detuning ∆ c for the parameters κ =350 ω R , U = ω R , N = 10 . Each curve is for a different valueof the pump strength η : thick blue η = 0 . η cr , dashed brown η = 1 . η cr , thin magenta η = 2 . η cr , and dash-dotted red η =3 . η cr . As can be seen, as η increases the lineshapes becomemore and more asymmetric and fold over at the critical pumpstrength η cr ( q = 0) ≡ η = 325 ω R . The atomic wave functioncorresponds to the first band with q = 0. V. BISTABILITY AND LOOPS
As mentioned in the Introduction, optical bistabilitywas discovered in the ultracold atom experiments [5] and[6] via a hysteresis effect as the detuning of the pumplaser was swept from above and below the cavity reso-nance. This effect can be described theoretically by usingthe self-consistency equation to calculate the number ofphotons n ph in the cavity at each value of the detuning(the number of cavity photons can be measured directlyfrom the photon current transmitted through the cav-ity which is given by n ph κ ). The results are plotted inFig. 2 for different values of the pump strength η . Inthe absence of atoms the cavity lineshape is a lorentziancentered at ∆ c = 0 with a full width at half maximum2 κ . At small pump intensity, the presence of the atomsshifts the center of the resonance by N U (cid:104) cos ( x ) (cid:105) whilethe shape is largely unaltered. But as η is increased, thelineshape curve becomes asymmetric and eventually foldsover when η is above a critical value η cr ( q = 0) ≡ η , in-dicating multiple solutions and hence bistability. η cr ( q )depends on the quasi-momentum and η is defined as its −10 −5 0 5 10−0.4−0.200.20.40.60.811.2 (∆ c - N U / /κ E [ ψ ] / N ( un i t s o f E R ) η = 0.5 η cr η = 1.5 η cr η = 2.5 η cr η = 3.5 η cr FIG. 3. (Color online) The energy given by the reduced hamil-tonian (13) as a function of detuning ∆ c for the parameters q = 0, κ = 350 ω R , U = ω R , N = 10 . Each curve is for adifferent value of the pump strength η : thick blue η = 0 . η cr ,dashed brown η = 1 . η cr , thin magenta η = 2 . η cr , anddash-dotted red η = 3 . η cr . The critical pump strength is η = 325 ω R . For η > η swallowtail loops develop corre-sponding to bistability. The loops grow in size as η increases. value at q = 0.In the folded over region, only one the solutions (corre-sponding either to the highest or the lowest photon num-ber) is seen at a time, depending upon the direction ofthe sweep. The middle branch cannot be accessed usingthis experimental protocol and is in any case unstable,as will be discussed at more length in Section X.Note that in Fig. 2 we have chosen the quantum stateof the atoms to be the q = 0 Bloch state. In fact, this isthe case for all the plots in this section because that isclosest to the situation in the experiments that have beenconducted so far. One of the main points of this paper isessentially to examine the extra degree of freedom con-ferred by the quasi-momentum. For atoms in ordinaryoptical lattices the quasi-momentum can be controlledby accelerating the lattice (rather than the atoms) via aphase chirp [45]. An accelerated lattice is hard to achieveinside a Fabry-Perot cavity, but if the atoms have a mag-netic moment one can instead subject them to an inho-mogeneous magnetic field so that a force is applied (or,of course, in a vertically oriented cavity the atoms will beaccelerated by gravity). As is well known from the theoryof Bloch oscillations, under a constant force F the quasi-momentum evolves according to the Bloch acceleration n ph E [ ψ ] / N ( un i t s o f E R ) increasing ∆ c FIG. 4. (Color online) The double well structure of theenergy of the reduced hamiltonian (13) as a function ofcavity photon number n ph . Each curve is for a differ-ent value of the detuning ∆ c . The values are ∆ c =1600 , , , , , ω R . The arrow indicateshow the curves evolve as ∆ c increases. The top most curve(blue), and the bottom most curve (yellow), have only oneminimum, whereas the rest of the curves have two minima(the inset shows a zoom-in of the curves close to n ph = 0)indicating bistability. Consequently, bistability only occursfor a certain limited range of ∆ c . The other parameters are q = 0, κ = 350 ω R , U = ω R , N = 10 , η = 3 . η , where η = 325 ω R . theorem [42, 43] q ( t ) = q + F tk c E R (22)where q and t are expressed in the dimensionless unitsgiven in Section II, and q is the quasi-momentum attime t = 0. The Bloch acceleration theorem requiresthat the evolution be adiabatic, and the implications ofthis for intra-cavity optical lattices have been discussedin [31], albeit without loops. The effect of a constantforce is thus to sweep the system through the band ata constant rate and, in principle, any quasi-momentumcan be achieved by switching off the magnetic field afteran appropriate time delay.Figure 3 depicts the energy versus detuning curves cor-responding to the photon number versus detuning curvesshown in Fig. 2. In the bistable regime, the energy curvesin Fig. 3 develop swallowtail loops. This can be under-stood in terms of the familiar connection between bista-bility, hysteresis, and the change in the energy manifolddescribed in detail by [19]. Consider one of the curvesin Fig. 3 where there is bistability, e.g. the curve with η = 3 η . For ∆ c values to the left and right of the swal-lowtail loop, the energy functional Eq. (13) has a singleextremum corresponding to a particular wave function ψ q ( x, n ph ). In the bistable region, the energy functionalhas the structure of a double-well, furnishing three ex-trema: two minima and one maximum, that give thethree branches of the loop corresponding to three differ-ent wave functions. This double-well structure is shownin Fig. 4 as a function of n ph for different values of de-tuning ∆ c . Note that, as already observed above, theself-consistency equation (9) provides a direct mappingbetween the wave function and n ph .Figures 2, 3, and 4 all show different aspects of hystere-sis as ∆ c is swept either from above or below the cavityresonance. It is enlightening to see how it arises in Fig. 4.If the detuning is swept from below the resonance theninitially there is a single solution for n ph , given by theminimum in the reduced energy functional which occursat the very left hand side of Fig. 4, as best seen in theinset. When ∆ c is increased another solution appears ata larger value of n ph . However, this state of the systemis not realized (for this direction of the detuning sweep),even when it becomes the global minimum, until the en-ergy barrier between the two solutions vanishes and theleft hand minimum disappears. The system then jumpsto the new minimum at a larger value of n ph . The reversehappens when ∆ c is swept in the other direction. Thishysteretic behavior is corroborated by Figs 2 and 3.In the Section VI, a method for determining η cr ( q ) isdescribed. Generally, this requires a numerical computa-tion, but for small values of the intracavity lattice depthan analytical expression can be worked out. It turnsout that the dependence of η cr on the atomic state quasi-momentum can be used to explain the loops in the energyquasi-momentum plots (Fig. 1). VI. CRITICAL PUMP STRENGTH FORBISTABILITYA. Conditions for bistability
Returning to the cavity lineshape shown in Fig. 2, werecall that as the pump strength η is increased the steadystate photon number in the cavity can exhibit bistabilityfor a certain range of detuning ∆ c . Bistability first devel-ops at a single value of the detuning, which we denote by∆ c = ∆ . The critical pump strength at which this bista-bility at ∆ occurs is η cr ( q ), and in this section we wantto calculate it. Let us first re-write the self-consistencyequation (23) as n ph n max = 11 + (cid:16) ∆ c − NU f ( U n ph ,q ) κ (cid:17) , (23)where n max ≡ η /κ is the maximum number of the pho-tons that can be in the cavity at steady state. In order toreinforce the idea that the wave function and the number n ph h c o s ( x ) i FIG. 5. (Color online) Plot of the atom-light overlap integral f ( n ph , q ) = (cid:104) cos ( x ) (cid:105) , first defined in Eq. (7), as a functionof the cavity photon number n ph . The atomic wave functionis taken to be the q = 0 Bloch wave of the first band, andthe atom-light interaction is set at U = 5 ω R . Note that themaximum value f ( n ph , q ) can take is one half, irrespective ofthe values of U and q . As n ph → ∞ we find that f ( n ph , q ) → of photons are really equivalent quantities, we have re-placed the notation for the integral (cid:104) cos ( x ) (cid:105) appearingin (23) by f ( U n ph , q ) ≡ (cid:104) cos ( x ) (cid:105) . (24)This function is plotted in Fig. 5 for blue detuning(∆ a >
0) where we see that as the lattice gets deeper theatomic wave function adjusts to become more localizedon the lattice nodes and reduce the overlap between thelight and the atoms. Furthermore, the steep gradient atshallow lattices implies that the system is more sensitive,and so more nonlinear, at small photon numbers.It is instructive to solve Eq. (23) graphically as theintersection between two functions of n ph , as shown inFig. 6. The left hand side is a straight line whose gradientis 1 /n max . For very small n max the gradient is very largeand there is only one solution close to the origin. As n max is increased the gradient is reduced and the straight linejust grazes the curve at a critical value of n max at whichthere are now two solutions. Increasing n max further,there is then a range of values of n max for which thereare three solutions. Finally, for large n max there is onlyone solution again. When three solutions exist for certainvalues of n max , the system becomes bistable and a plotof the input intensity (proportional to n max ) versus theoutput intensity (proportional to n ph ) has the classic s-0 n ph n ph / n m a x FIG. 6. (Color online) A graphical solution of the self-consistency equation in the form (23). The blue curve rep-resents the right hand side of Eq. (23) for typical values ofthe cavity parameters ( q = 0). The red dash-dotted, greendashed and black straight lines represent the left hand side ofthe equation plotted for different values of n max ; they inter-sect the blue curve at one, three, and one points, respectively.The blue curve tends to a finite value at n ph = 0 which is setby the fact that for n ph = 0 we have f = 1 / shaped form shown in Fig. 7. This picture suggests aconvenient way to determine the conditions for bistabilitybecause the two points where the curve turns over delimitthe bistable region. These points satisfy ∂n max /∂n ph =0, giving κ + (∆ c − N U f ) − n ph (∆ c − N U f ) N U ∂f∂n ph = 0 . (25)This equation can be solved numerically for n ph for dif-ferent values of ∆ c , assuming that κ , N , and U are fixed.As expected from Fig. 7 (see also Fig. 14), depending onthe value of ∆ c , there are either zero, one, or two valuesof n ph that satisfy Eq. (25). For large values of ∆ c thereare no solutions. As ∆ c is reduced, a single solution for n ph suddenly appears at ∆ c = ∆ , which, by substitut-ing this value of n ph into the self-consistency equation(23), gives us η cr ( q ). Reducing ∆ c further, this single so-lution immediately branches into two solutions for n ph .Referring to Fig. 7, we see that these two solutions for n ph define a range of values of n max , and hence also of η , where bistability occurs. We can find this range ofvalues of η by inserting the two solutions for n ph intothe self-consistency equation to give us η and η . When
20 40 60 80 100051015202530 n max n ph FIG. 7. (Color online) Input intensity vs output intensity fora bistable cavity system. In this example the atomic wavefunction in the cavity is in the q = 0 state and κ = 350 ω R ,∆ c = 1500 ω R , and U = ω R . The points where the curvefolds over are given by the solution of Eq. (25). η < η < η bistability occurs. Note that the values of η cr ( q ), η , and η , all depend on the state of the atomicwave function and hence are different for different quasi-momenta. This dependence on quasi-momentum is whatlies behind the existence of loops in the band structure. B. Critical pump strength in shallow lattices
In the regime of shallow lattices, the bistability condi-tion given by Eq. (25) can be solved analytically. Firstconsider the critical pump strength of the q = 0 case, forwhich we use the notation η cr ( q = 0) ≡ η . As describedin [6], and as can be seen from Fig. 5, for small latticedepths we can linearize the atom-light overlap integral as f ( U n ph , q = 0) = 12 − U n ph . (26)Substituting this into the self-consistency equation (9),we obtain a cubic equation in n ph , which is reminiscentof the classical Kerr nonlinearity in a medium with anintensity dependent refractive index [4] (note that whenwe come to numerically solve Eq. (25) below in this sec-tion and in the rest of the paper, we are going beyond theclassical Kerr effect). The condition (25) for bistabilityin this limit then reduces to the solution of a quadratic1equation in n ph ,3 n N U − n ph N U ( N U − c )+ 64(( N U − c ) + 4 κ ) = 0 . (27)The vanishing of the discriminant of the above equationrequires that ∆ = N U / − √ κ , η = (cid:113) κ √ NU , and n = NU √ κ . In the last expression, n is the number ofphotons in the cavity at the critical point ∆ c = ∆ . C. Critical pump strength as a function ofquasi-momentum
We now extend the above analysis for shallow latticesto include non-zero quasi-momentum. Expanding theatomic Bloch state in a Fourier series, and assuming shal-low lattice depths, we can truncate the series after threeterms ψ q ( x, n ph , t ) = e iqx ( c ( t ) + c ( t ) e i x + c ( t ) e − i x ) . (28)In this state one can explicitly calculate (cid:104) cos ( x ) (cid:105) as (cid:104) cos ( x ) (cid:105) = 12 + 12 (Re c c ∗ + Re c c ∗ ) ≡
12 + 12 ( X ( t ) + Y ( t )) . (29)Rewriting the equations of motion (4) and (5) for thenewly defined variables one findsd X d t + (4 q + 4) X + ( q + 1) U α ∗ α | c | = 0 , d Y d t + (4 q − Y − ( q − U α ∗ α | c | = 0 , d α d t = (cid:18) i ∆ c − i N U − i N U X + Y ) − κ (cid:19) α + η. The atomic state has therefore been mapped onto twocoupled oscillators X and Y . The oscillators are not cou-pled to each other directly, but do interact through thelight field α which acts as a driving term for both ofthem. The above equations resemble Eqns (3) and (4)of [6], and introduce an analogy to optomechanics [46].Solving these equations at steady state gives α = ηκ − i ( ˜∆ c − N U / X + Y )) , (4 q + 4) X + ( q + 1) U η κ + (cid:16) ˜∆ c − N U / X + Y ) (cid:17) | c | = 0 , (4 q − Y − ( q − U η κ + (cid:16) ˜∆ c − N U / X + Y ) (cid:17) | c | = 0 , where ˜∆ c = ∆ c − N U /
2. Combining the steady statesolutions into a single equation for the variable p = X + quasi-momentum q (units of k c ) η c r / κ FIG. 8. (Color online) Comparison between exact numericalcalculation (dots) and analytical estimate (line) for the criticalpump strength η cr at which loops appear as a function of thequasi-momentum. The values of the parameters are U = ω R , N = 10 , and κ = 350 ω R . The analytical estimate is fromEq. (31) which is accurate for small lattice depths. Note thatthe agreement is good for quasi-momentum close to q = 0. Y , and assuming | c | ≈ p = ¯ n max (cid:16) ∆ c κ − NU κ − NU p κ (cid:17) , (30)where ¯ n max = U η q − κ . Comparing Eqns (30) and (23)in the limit when f = 1 / − U n ph /
16, we finally obtainan expression for the critical pump strength above whichbistability occurs as a function of qη cr ( q ) = (cid:115) κ (1 − q )3 √ N U (31)where we remind the reader that the frequencies in thisexpression are in units of the recoil frequency ω R . Thisestimate for η cr ( q ) is compared to the full numerical so-lution of Eq. (25) in Fig. 8. The parameters are suchthat the maximum intracavity depth is only of the orderof one atomic recoil energy E R , and hence the approxi-mation agrees well with the numerical calculation. Theagreement deteriorates as q →
1. This is due to the factthat the above two mode approximation fails at q = 1because the coefficient c in Eq. (28) is equal to zero atthe edge of the Brillouin zone.Let us now connect the above results to the phe-nomenon of loops in the band structure. Because η cr ,2and also η , and η depend on the value of q , as we vary q we expect that the conditions required to have bista-bility won’t necessarily be met over the entire Brillouinzone. That is, as q is varied, η may no longer lie in therange η ( q ) < η < η ( q ). In that case, we expect anyadditional solutions to form closed loops extending onlyover part of the Brillouin zone, rather than entire bandscovering the whole Brillouin zone. The dependence of η cr upon q given by Eq. (31) suggests that for shallowlattices the loops will form first at the edge of the Bril-louin zone and then propagate inwards as η is increased.Looking back at Fig. 1, which shows the energy plottedas a function of quasi-momentum at a fixed value of η and ∆ c , we see a loop centered at q = 0, but which doesnot extend out to q = 1, in apparent contradiction towhat is predicted by Eq. (31). This is because the latticein Fig. 1 is too deep for Eq. (31) to apply. In the nextsection we shall examine how loops appear and disappearand, in particular, we will see that loops are indeed bornat the edges of the Brillouin zone. VII. THE BIRTH AND DEATH OF LOOPS
In this Section we examine how loops appear and dis-appear in the band structure as the detuning and pump-ing are varied. We have already seen in Sections IV andV how multiple solutions and swallowtail loops developas the detuning from the cavity resonance is varied, butthis was for fixed quasi-momentum q . Here we includethe quasi-momentum dependence. In Fig. 9 we plot theevolution of the loop structures that appear in the bandstructure as ∆ c is varied. The detuning increases fromthe top to the bottom panel. In the plots, the pumpstrength is fixed at η = 2 . η and we see that for smalldetunings, when bistability initially sets in, swallowtailshaped loops appear at the outer edges of the Brillouinzone. As the detuning is increased the swallowtail loopsfrom the two edges move closer and merge. Initially, themerged loop lies partly below the lowest band, but as thedetuning is further increased it moves up and separatesfrom the lowest band. The loop then shrinks in size andvanishes.One important point to notice is that the swallowtailloop in plot (b) of Fig. 9 is qualitatively different fromthe ones obtained in [17, 18] for an interacting BEC inan optical lattice because in our case the energy disper-sion continues to have zero slope at the band edge evenwhen the loops have been formed. The nonlinearity in[17, 18], which is due to interatomic interactions, hasquite a different form to that considered here. For exam-ple, the nonlinearity arising from interactions is spatiallydependent due to the variation in density of a BEC in anoptical lattice, whereas the nonlinearity considered heredoes not have a spatial dependence because it appearsunder an integral over space, see Eq. (7).In Fig. 10 we plot the same thing as Fig. 9, except thatwe have reversed the sign of the atom-light coupling U E [ ψ ] / N ( un i t s o f E R ) quasi-momentum q (units of k c ) (a)(b)(c)(d)(e) FIG. 9. (Color online) The birth and death of band struc-ture loops as the laser-cavity detuning ∆ c is varied, for thecase when the laser is blue-detuned from atomic resonance(∆ a > c increases as one goes from (a) to (e) as follows:1500 ω R , 2100 ω R , 2600 ω R , 3100 ω R and 3600 ω R . The rest ofthe parameters are κ = 350 ω R , U = ω R , N = 10 , η = 2 . η and η = 325 ω R . E [ ψ ] / N ( un i t s o f E R ) quasi-momentum q (units of k c ) (a)(b)(c)(d)(e) FIG. 10. (Color online) The birth and death of band struc-ture loops as the laser-cavity detuning ∆ c is varied, for thecase when the laser is red-detuned from atomic resonance(∆ a < c increases as one goes from (a) to (e) as follows: − ω R , − ω R , − ω R , − ω R and − ω R .The rest of the parameters are κ = 350 ω R , U = − ω R , N = 10 , η = 2 . η and η = 325 ω R . so that it is negative. Experimentally, this is the casewhen ∆ a <
0, i.e. the pump laser is red detuned fromthe atomic resonance. Note that the effect of the signflip upon the potential term U n ph cos ( x ) occuring inthe atomic Schr¨odinger equation is equivalent to a spa-tial translation of π/
2. This transforms the atom-lightoverlap integral (7) as (cid:104) cos ( x ) (cid:105) → − (cid:104) cos ( x ) (cid:105) (32)where the left hand side refers to the U < U > n ph = η κ + (∆ c + N | U | − N | U | f ( | U | n ph , q )) . (33)Figure 11 plots the evolution of the band structure asthe external laser pumping η is varied, for fixed detuning∆ c . We see that as η increased the loops first appear asswallowtails at the edges of the Brillouin zone (see insetin Fig. 11a), in agreement with the predictions of Eq.(31).Figures 12 and 13 are both 3D plots of the loops asfunctions of ∆ c and q , but each one covers a differentrange of ∆ c . In the first plot the merging of the swal-lowtail loops into the band center loops is shown. In thesecond plot, as ∆ c increases the band center loops moveupwards in energy and separate from the ground band,shrink in size and eventually disappear.As will be demonstrated in the next section, in cer-tain parameter regimes the spectrum may qualitativelychange compared to the above. In particular, we showthat for some q (cid:54) = 0, and for sufficiently larger pumpstrength η , we can achieve tristability. Moreover, weshow how this multistable behaviour can be understoodin terms of catastrophe theory. VIII. TRISTABILITY
Thus far we have seen that there are regions of theparameter space { ∆ c , η, q, U } where the self-consistencyequation (9) applied to the first band admits either oneor three solutions. In the latter case we have bistability.It is natural to ask whether there are other regions ofparameter space where even more simultaneous solutionscan occur? A recent paper discussing two-componentBECs in a cavity has found regions of parameter spacethat support tristability [47] and in this section we wantto examine whether tristability is possible in the ordinaryone-component case but with finite values of the quasi-momentum.In Fig. 14 we show the region of { ∆ c , η } space wherebistability occurs. This plot is for fixed values of U and q . In particular, for q = 0 we at most find bistabil-ity. The crescent shape in Fig. 14 demarks the regionwith three solutions: the number of solutions changesby two as one crosses its boundary. The location of η ,4 E [ ψ ] / N ( un i t s o f E R ) quasi-momentum q (units of k c ) (b)(c)(d)(e)(a) FIG. 11. (Color online) The birth and death of band structure loops as the pumping rate η is varied. In this figure the detuningis held constant at ∆ c = 2900 ω R . The value of η increases from (a) to (e) as follows: 0 . η , η , η , η , η , where η = 325 ω R as usual. The inset shows a zoom-in for η = 0 . η , illustrating that as η is increased, the loops are born at the edges of theBrillouin zone. FIG. 12. (Color online) Energy as a function of quasi-momentum q and detuning ∆ c . At smaller values of ∆ c , theswallowtail loops occur in pairs close to the edges of the Bril-louin zone at q = ±
1, and as the detuning is increased theypropagate inwards and merge as shown in this plot. ∆ c in-creases out of the page. Parameters are κ = 350 ω R , U = ω R , N = 10 , η = 2 . η , and η = 325 ω R .FIG. 13. (Color online) Energy as a function of quasi-momentum q and detuning ∆ c . ∆ c increases out of the page.For larger values of ∆ c , the band center loops move up in en-ergy and do not touch the lower band (shaded red in the plot).Eventually, for still higher values of detuning, they shrink anddisappear as shown in this plot. Parameters are κ = 350 ω R , U = ω R , N = 10 , η = 2 . η , and η = 325 ω R . the smallest pump strength for which bistability occurswhen q = 0, is indicated by an arrow on the vertical axisin order to make contact with the discussion of Sec. V.However, when we allow non-zero values of q we find thatwe can indeed have regions of tristability, due to the pres-ence of five solutions, which occur inside the swallowtailshaped region in Fig. 15, which is plotted for q = 0 . −15 −10 −50 2.9 5.7 8.6 11.414.317.120 (∆ c - N U / /κ η / κ η FIG. 14. (Color online) Bifurcation structure of the solutionsto the self-consistency equation Eq. (9) in the { η, ∆ c } planewith q = 0, U = ω R , κ = 350 ω R , N = 10 . The numberson the plot indicate the number of solutions that exist for thesteady state photon number in the cavity. The critical valueof the pumping η = η cr ( q = 0) for bistability for these pa-rameters is indicated by the arrow. Inside the crescent shapedregion the system supports three solutions (one unstable), i.e.it is bistable. Fig. 16 plots the corresponding photon number versuslaser pumping curve for a fixed value of ∆ c , and illus-trates how, as the laser pumping strength is changed,the system goes from one solution, to three solutions, tofive solutions, back to three solutions and finally back toone solution again. This curve is calculated for a verti-cal slice through the parameter space shown in Fig. 15.We give an example of the band structure when there istristability in Fig. 17.So far we have conducted a rather ad hoc explo-ration of the four-dimensional parameter space given by { ∆ c , η, q, U } . Furthermore, in the two-dimensional slicesshown in Figures 14 and 15 we have glimpsed snap shotsof a rather complex looking structure of solutions. Weare therefore led to ask whether there is any order in thiscomplexity? Are the geometric structures seen in theplots random, or is there in fact an underlying structurethat organizes them? In order to make progress with un-derstanding multistability it would be useful to have amore systematic framework for analyzing our solutionsand just such a framework is provided by catastrophe6 −14 −12 −10 −8 −6 −4 −20 2.9 5.7 8.6 11.414.3 (∆ c - N U / /κ η / κ
13 53
FIG. 15. (Color online) Bifurcation structure of the solutionsto the self-consistency equation Eq. (9) in the { η, ∆ c } planewith q = 0 . U = ω R , κ = 350 ω R , N = 10 . The numberson the plot indicate the number of solutions that exist forthe steady state photon number in the cavity. Inside theswallowtail shaped curve there are five solutions and hencemultistability, see Fig. 16. theory, which we now discuss. As will become clear, thestructures we see in parameter space are not only generic,but are organized in a very particular, and therefore pre-dictable, way. IX. CATASTROPHE THEORY ANALYSISA. Overview of catastrophe theory
Catastrophe theory is a branch of bifurcation theorythat concerns the study of singularities of gradient maps[35, 37]. Examples of gradient maps abound in physics,for example Hamilton’s principle of least action in me-chanics, and Fermat’s principle of least time in optics.In both theories the physical paths (rays) are given bythe stationary points of a generating function Φ( s ; C ),which in mechanics is the action and in optics is the op-tical distance. In both cases the gradient map is obtainedfrom a variational principle that takes the form ∂ Φ( s ; C ) ∂s = 0 . (34) n max n ph FIG. 16. Plot of multistable steady state photon number n ph versus n max (equal to η /κ ) for ∆ c = 1630 ω R and q = 0 . U = ω R , κ = 350 ω R , N = 10 . The ∆ c value is chosen fromthe region which supports five solutions in Fig. 15. −1 −0.5 0 0.5 10.40.50.60.70.80.911.11.21.31.4 quasi-momentum q (units of k c ) E [ ψ ] / N ( un i t s o f E R ) FIG. 17. (Color online) Plot of tristable band structure. Theparameters are given by η = 980 ω R , ∆ c = 1640 ω R , κ =350 ω R , N = 10 , and U = ω R . (cid:45) C (cid:45) (cid:45) C (cid:45) (cid:45) FIG. 18. (Color online) The cusp catastrophe is generatedby the quartic potential function Φ given in Table I, whichcan be viewed as representing a double-well potential. Thered folded-over surface plotted in this figure obeys the cubicstate equation ∂ Φ /∂s = s + C s + C = 0, which gives thestationary solutions s i of Φ. When C < s , s and s , for each value of C and C . These points are the two minima and single maximumof the double-well potential. When C > s corresponding to the minimum of a singlewell. A vertical slice through the figure such that C = 0gives a pitchfork bifurcation. The { C , C } plane forms thetwo dimensional control space where the cusp catastrophe it-self lives, and this is shown at the bottom of the figure. Thecusp catastrophe is formed of two fold curves joined at a sin-gular cusp point. The cusp catastrophe demarks the regionof control space that sustains three solutions for s , and so itis given by the projection of the folded-over part of the statesurface onto the { C , C } plane. Crossing the fold lines frominside the cusp to outside it, two of the solutions (the maxi-mum and one of the minima) annihilate. Mathematically, thisis described by the potential function being stationary to thenext higher order [35], i.e. ∂ Φ /∂s = 3 s + C = 0. Eliminat-ing s by combining this equation with the state equation givesthe equation for the cusp catastrophe as C = ± (cid:112) − C / C = C = 0, all three stationary pointscoalesce simultaneously to leave a single solution. In catastrophe theory, this equation is sometimes referredto as the state equation , the generating function Φ( s ; C )is called the potential function and the variables appear-ing in the potential function are considered to be of twobasic types: the state variables s = { s , s , s . . . } andthe control parameters C = { C , C , C . . . } . The statevariables parameterize all possible paths (not just thepaths corresponding to the stationary points) and the control parameters determine the conditions imposed onthe paths. For example, if we are interested in the path(s)which pass through the point { X, Y, Z } in three dimen-sional space, then the coordinates { X, Y, Z } form thecontrol parameters. For a fixed set of control parame-ters, the potential function defines the height of a land-scape with coordinates s . The classical paths (rays) arethen the stationary points s i (labelled by the index i ifthere are more than one) of this landscape, namely themountain peaks, valley bottoms and saddles [35].The gradient map becomes singular when there is morethan one stationary point for a given set of control pa-rameters. In optics this is the phenomenon of focusing,because more than one ray passes through the same phys-ical point C = { X, Y, Z } , leading to a caustic. The caus-tic, or catastrophe, as it is known in catastrophe theory,lives in control space, which in the standard optics caseis physical 3-dimensional space. As C is varied one canexplore the structure of the caustic. This is what wasshown above in Figures 14 and 15, which are two dimen-sional slices through the parameter space { ∆ c , η, q, U } .The crescent and swallowtail shapes are the catastrophes,whose full structure only becomes apparent when viewedin four-dimensional parameter space.Catastrophes are points, lines, surfaces, and hypersur-faces in control space across which the number of solu-tions to the problem changes. Catastrophe theory clas-sifies these catastrophes in terms of their codimension K , which is the difference between the dimension of thecontrol space and the dimension of the catastrophe. Forexample, if we consider the two dimensional space shownin Figures 14 and 15, we find “fold” curves ( K = 1) and“cusp” points ( K = 2). If we add a third dimension thenwe would find fold surfaces ( K = 1), cusp edges ( K = 2),and “swallowtail” points ( K = 3).In order to make the foregoing discussion more con-crete, consider the structure shown in Fig. 18 which il-lustrates the cusp catastrophe. The surface shown inthe figure is the state surface ∂ Φ /∂s = 0 plotted in acomposite of control and state space. The control spaceis two dimensional and is given by the C , C plane atthe bottom of the figure. As listed in Table I, the cuspcatastrophe is described by a quartic potential functionwhich, by varying the control parameters C and C ,can be tuned between being a double or a single wellpotential. A prominent physical example of such quar-tic potential function is the thermodynamic potential inLandau’s phenomenological theory for continuous phasetransitions [48]Φ( s ; P, T, h ) = Φ + Bs + As + hs. (35)The order parameter s for the phase transition can beidentified as the state variable. The parameter h de-scribes an external field, and A and B are functions ofpressure P and temperature T . At first sight, it appearsas though the Landau potential function contains threecontrol parameters, and so does not correspond to thecusp catastrophe. However, it is easy to see that one of8 TABLE I. Standard forms of the cuspoid catastrophesName Φ( s ; C ) KFold s / Cs s / C s / C s s / C s / C s / C s s / C s / C s / C s / C s the parameters is redundant because the state equationcan be written in terms of only two control parameters: C = h/B and C = A/B . The Landau thermodynamicpotential can therefore be seen to correspond to the cusppotential function.In the present case of ultracold atoms in an opticalcavity, the potential function is the reduced hamiltonian(13), the state variable is the wave function ψ , and thecontrol parameters are { ∆ c , η, q, U } . We assume thatthe cavity decay rate κ and number of atoms N areconstants that are unchanged throughout the analysis.The stationary Schr¨odinger equation, obtained from thetime-dependent Schr¨odinger equation (10), is thereforethe state equation which determines the allowed clas-sical “paths” (rays) ψ . However, because we are in-terested here in solutions of the Bloch wave form, our“paths” ψ are Mathieu functions labelled by the quasi-momentum and the depth of the optical lattice. Thequasi-momentum q is one of the control parameters, butthe depth of the optical lattice is determined uniquely, viathe self-consistency equation (9), by the number of pho-tons in the cavity n ph . Therefore, we can choose to workwith n ph rather than ψ as already discussed in SectionIV B. The state equation which determines its stationaryvalues is the self-consistency equation (9).The purpose of formulating our problem in terms ofcatastrophe theory is that we can now take advantage ofa very powerful theorem [33]. This states that there areonly a strictly limited number of different forms whichthe potential function Φ( s ; C ) can take in the neighbour-hood of a catastrophe. The first four are listed in Table I.Note that each of the standard forms is a polynomial in s , but is linear in the control parameters. In fact, for con-trol spaces of dimension four or less there are only sevendistinct structurally stable potential functions. Three ofthese require two state variables, whereas we only requireone state variable, so that leaves the four so-called cus-poid catastrophes, which are the ones listed in Table I.This remarkable result allows us to predict, at leastqualitatively, the structures seen in Figures 14 and 15given only very rudimentary information such as thenumber of control parameters. However, the sting in thetail is that it is rare for the potential function that ap-pears in any particular problem to already be in one ofthe standard forms shown in Table I. Rather, it is gener-ally necessary to perform various transformations uponthe variables in order to manipulate the raw potentialfunction into one of the standard forms. We already sawthis for the rather simple case of the Landau theory dis- cussed above, and we shall see below that this is also truefor the problem of multistability in atom-cavity systems.From Table I, we see that the butterfly potential func-tion is the only one which gives up to five stationary so-lutions and has a four dimensional control space. We cantherefore immediately say that our problem correspondsat least to a butterfly catastrophe because we have al-ready found parameter regimes which give five solutions.This does not rule out the possibility of a higher catas-trophe (the higher catastrophes contain the lower ones,as can be seen in Fig. 18 for the case of the cusp andthe fold), but until we find regimes with a higher num-ber of solutions (and, indeed, we have not) we shall workwith the hypothesis that we are dealing with a butterflycatastrophe. We might therefore expect to find a veryspecial point in parameter space where all five solutionsmerge into one (the “butterfly point”). However, thisrequires us to be able to maneuver through parameterspace in order to find this point. The delicate issue ofwhether the four experimental parameters { ∆ c , η, q, U } can be transformed into the butterfly’s four linearly in-dependent control parameters { C , C , C , C } , and soallow us to fully explore the butterfly catastrophe, willbe studied in Section IX C below. First, we begin withthe simpler case of shallow lattices. B. Application of catastrophe theory to shallowlattices
The starting point for our analysis will not be the po-tential function Φ, which in our case is the energy func-tional (13) that explicitly depends on the atomic wavefunction ψ q ( x, n ph ) and implicitly on the photon num-ber n ph . Instead, we begin one step further along withthe state equation which, as explained above, is providedby the self-consistency equation Eq. (9) for the photonnumber. This was also the approach adopted in [ ? ]in order to tackle bistability in traditional laser systems.The photon numbers that satisfy the state equation forgiven values of the control parameters form the set ofstationary points of the potential function. For our statevariable we shall actually choose v ≡ U n ph . (36)In terms of v the state equation can be written v + v (∆ c − N U f ( v, q )) − η U = 0 . (37)where (cid:104) cos ( x ) (cid:105) = f ( v, q ) is evaluated in a Bloch state ψ q ( x, v ), and the choice of v as the state variable is moti-vated by the dependence of the Bloch state on the prod-uct U n ph , which is the depth of the optical lattice, seeEq. (4). In the above equation we have also rescaled allfrequencies by κ (i.e. we have divided throughout by κ and set κ = 1).Equation (37) is not straightforward to analyze be-cause we do not have a closed-form analytical expression9for f ( v, q ). Thus, we find ourselves in the common situ-ation, as mentioned above, that it is not obvious whichstandard potential function, and hence which standardstate equation, corresponds to our problem. However,we have already seen in Sec. V that when the depth ofthe optical lattice is small, we can perform a series ex-pansion and obtain an approximate analytical expressionfor f ( v, q ). This is the approach we shall follow first andwe will take up the general case in the next subsection.Specializing to the case of q = 0, we use Eq. (26) to write f ( v, q = 0) = 1 / − v/
16, and upon substitution intoEq. (37) this gives v + b v + b v + b = 0 , (38) b = 32∆ c N U − ,b = 64 (cid:0) N U − c ) (cid:1) N U ,b = − η N U . The leading term in the above equation is cubic in thestate variable v (similar in that respect to the classicalKerr non-linearity [3, 4]). It is therefore close to, butnot yet identical with, the standard form of the stateequation for a cusp s + C s + C = 0 . (39)Complete equivalence to the cusp can be achieved byremoving the quadratic term from (38) via the transfor-mation s = v + b /
3, leading to the following values for C and C C = b − b , (40) C = b − b b + 227 b . (41)Thus, we see that the canonical control parameters inthe final mapping to the standard form are complicatedfunctions of the physical control parameters η, ∆ c , and U .Let us compare the above prediction to the case illus-trated in Fig. 14. The figure shows the number of steadystate solutions for the photon number in the { η, ∆ c } plane, for fixed values of U and q . The figure was cal-culated for q = 0, and the part of it close to the horizon-tal axis corresponds to low photon number, and so theshallow lattice theory outlined above applies in that re-gion. We indeed find that the first derivative of the stateequation vanishes at all points along the curves separat-ing regions with one and three solutions (this is how thecurves were computed), which are therefore fold curves,while at the point with η = η we find the second deriva-tive also vanishes, identifying it as a cusp point where allthree solutions coalesce into a single solution. We there-fore find that catastrophe theory correctly accounts forthe structure seen in Fig. 14. A key point to note from the above analysis is that theunderlying catastrophe that we identified had only twocontrol parameters even though there were three “exper-imental” parameters that could be varied in the originalstatement of the physical problem (we set q = 0). Wemet a similar situation for the Landau theory of continu-ous phase transitions discussed above. The question thenarises, how do we identify the underlying catastrophe incases where the transformation to standard form is hardto find? One way to proceed is via the defining characterof each potential function in Table I, which is the highestderivative that vanishes at the most singular point. Forthe cusp the most singular point is s = C = C = 0where the first, second, and third derivatives of the po-tential function with respect to s vanish. We shall takeadvantage of this defining character in order to tacklethe general case of a lattice of arbitrary depth in thenext subsection.It is worth commenting on the role of the extra di-mension in control space that was present in the originalstatement of the shallow lattice problem, as given by Eq.(38). It is easy to imagine that, given a basic catastro-phe, we can always embed it in a control space of higherdimension without fundamentally changing the catastro-phe providing we extend it into the extra dimension in atrivial way. For example, given the cusp catastrophe thathas a minimal control space which is two-dimensional,as shown in Fig. 18, we can always add a third dimen-sion such that the cusp becomes a structure rather likea tent, with the ridge pole being a cusp edge (a continu-ous line of cusps). The existence of the line of cusps canbe inferred from the shift of variables s = v + b / b which is the extracontrol parameter. C. Application of catastrophe theory to lattices ofarbitrary depth
Lifting the restriction of small photon number, we de-fine the notation G ( v ; ∆ c , η, q, U ) for the left hand sideof the state equation (37) G ( v ; ∆ c , η, q, U ) ≡ v + v (∆ c − N U f ( v, q )) − η U = 0 . (42)The fact that we have four experimental parametersholds out the possibility that the above state equationis that of a butterfly catastrophe (see Table I). Further-more, because we have already discovered regimes withfive solutions, we must have at least a butterfly. How-ever, from our experience with the shallow lattice case,we know that we may not be able to fully explore thecatastrophe if some of the control parameters are trivial.We shall therefore investigate whether we can find a sin-gular point (the butterfly point) where the fifth and alllower derivatives of the potential function with respect0 quasi-momentum q (units of k c ) / ( N U ) q sw FIG. 19. A plot showing the values of 1 / ( NU ) obtainedfrom the simultaneous solution of Eqns (43)–(45) for differ-ent values of quasi-momentum. These equations give thefirst three derivatives of the state equation (42), and hencecorrespond to swallowtail points. Only the crosses satisfy1 / ( NU ) > q > q sw . to v simultaneously vanish (i.e. the fourth and lower or-der derivatives of G ( v ; ∆ c , η, q, U ) simultaneously van-ish) and Eq. (42) is also satisfied.Taking the derivatives of Eq. (42), and simplifyingslightly the resulting equations, we find (cid:18) ∆ c N U − f (cid:19) − vf f i (cid:18) ∆ c N U − f (cid:19) = − N U (43)2 f f i + v (cid:0) f f ii + ( f i ) (cid:1) f i + vf ii = ∆ c N U (44)3 f f ii + 3( f i ) + v (cid:0) f i f ii + f f iii (cid:1) vf iii + 3 f ii = ∆ c N U (45)12 f i f ii + 4 f f iii + v (cid:0) f i f iii + f f iv + 3( f ii ) (cid:1) vf iv + 4 f iii = ∆ c N U (46)where the first, second, third, and fourth derivatives ofthe function f ( v, q ) with respect to v are denoted by f i , f ii , f iii , f iv , respectively. In fact, we have seen Eq.(43) before, as it is the same as Eq. (25) that we used asa condition for bistability. Our strategy will be to findsolutions to Eqns (43–46) numerically. In order to fa-cilitate this, observe that Eq. (44) and Eq. (45) can be quasi-momentum q (units of k c ) v FIG. 20. (Color online) A plot showing the values of the statevariable v = U n ph for which the first three derivatives of thestate function G vanish simultaneously (swallowtail points).For approximately 0 . < q < . q . combined into a single equation:2 f f i + v (cid:0) f f ii + ( f i ) (cid:1) f i + vf ii − f f ii + 3( f i ) + v (cid:0) f i f ii + f f iii (cid:1) vf iii + 3 f ii = 0 . (47)We solve this equation for v at different values of q bynumerically finding the zeros of the left hand side. Oncewe find the zeros for a particular q , we can use Eq. (44) tocalculate ∆ c /N U at these values. Next, these values of v and ∆ c /N U are used in Eq. (43) to calculate 1 / ( N U ) .In this last step we find that we obtain the unphysicalresult 1 / ( N U ) < q > q sw , where q sw = 0 . / ( N U ) computed using the above method for differ-ent values of the quasi-momentum q in the neighborhoodof q sw . Only the crosses satisfy 1 / ( N U ) >
0. Wedrop the solutions corresponding to the dots, for which1 / ( N U ) <
0. For q > q sw there are always values ofq at which 1 / ( N U ) >
0. Note that q sw = 0 .
545 is auniversal result since it does not depend on any otherparameter values.The final step is to check if Eq. (46) for the fifth deriva-tive is satisfied at the values of ∆ c , N U , and v that wecomputed using the lower derivatives. We did not findany value of q > q sw where this was the case. The im-portant part of our numerical computation involves the1calculation of the derivative of the function f ( v, q ) forwhich we do not have an analytical expression. We usedthe MATLAB R (cid:13) [49] routine DERIVSUITE [50] to cal-culate the derivatives. This routine also provides errorson the derivatives and we can compound errors and findvalues for expressions like the left hand side of Eq. (47)with error. We use this error as the tolerance in our zerofinding. Some more details of this procedure are providedin the Appendix.The fact that we did not find a point where the fourhigher derivatives of the function G simultaneously van-ish in the range 0 < q < q because f ( v, q ) is symmetric under q → − q ) means thatalthough the underlying catastrophe that organizes thesolutions is at least a butterfly (because we find five solu-tions), the four experimental parameters at our disposal { ∆ c , η, q, U } do not translate into four linearly indepen-dent coordinates in control space { C , C , C , C } [51].We are therefore not able to navigate freely through thefour-dimensional control space and locate the butterflypoint at C = C = C = C = 0. This is the extensioninto four dimensions of the situation we already found inSection IX B for shallow lattices.The identification of places where three derivatives of G simultaneously vanish means that the highest singu-larities we have in the parameter space { ∆ c , η, q, U } areswallowtail points (swallowtails are contained within agreater butterfly catastrophe). In the Appendix we out-line a proof that in the neighbourhood of a point wherethree derivatives of the state equation vanish the poten-tial function must be equivalent to that of a swallowtailcatastrophe. Note that these swallowtails occur entirelyin control space and are thus true swallow tail catas-trophes in the sense of Table I. The swallowtails shownpreviously in Figures 3, 9, 10, 11, 12 and 13 are not swal-lowtail catastrophes because these figures show a com-bination of state and control spaces. By contrast, Fig-ures 14, 15 and 21 are solely in control space. In par-ticular, Fig. 21 shows a two-dimensional slice throughthe three-dimensional swallowtail catastrophe in controlspace. Unlike Fig. 15, this slice includes the swallowtailpoint, which is the highest singularity on the swallowtailcatastrophe, and is the point where four solutions of thestate equation (42) simultaneously coalesce so that num-ber of solutions changes by 4 (i.e. the point where threederivatives of G simultaneously vanish). We emphasizethat this can only occur when q > q sw . We also note fromFig. 20 that between 0 . < q < . q , which is again an indica-tion of the presence of a higher underlying catastrophe.As described in reference [52], swallowtails contain twocusps, and butterflies contain two swallowtails. X. STABILITY ANALYSIS
The stability of cold atoms in an optical cavity hasbeen treated in Refs. [12, 53]. In this case we follow an −15 −10 −5 00 2.9 5.7 8.6 11.414.3 (∆ c - N U / /κ η / κ −11.5−11−10.5−10−9.51.12.3 FIG. 21. (Color online) Bifurcation structure of the solutionsto the self-consistency equation Eq. (9) in the { η, ∆ c } planewith q = 0 . U = 1 . ω R , κ = 350 ω R , N = 10 , andhence NU / . κ . The numbers on the plot indicate thenumber of solutions for n ph in that region of the parameterspace. The inset shows a swallowtail singularity point wherefive solutions coalesce into a single solution. The coordinatesof the swallowtail point shown in this figure are v = 0 . η = 1 . κ , ∆ c = 6 . κ . approach more in line with [18], where the energy func-tional, Eq. (13), and the nonlinear equation of motion,Eq. (10), for the atomic wave function serve as the start-ing points for an examination of energetic, and dynamicstability, respectively. Hence, we examine the stability ofthe Bloch states at different values of quasi-momentumat fixed values of η and ∆ c . Before going into the detailsof the calculation we note that it is well known from thestudy of bistability in classical nonlinear cavity systems,that the back-bent branch of the lineshape profile shownin Fig. 2 is unstable, see, for example, [3]. From the threedimensional plots Fig. 12 and Fig. 13, we see that thisback-bent branch corresponds to the upper branch of theloop in energy-quasi-momentum space. Thus, we expectto find this branch unstable.We first consider energetic stability in the spirit of [18].The grand canonical potential per unit length is G [ ψ ] =2 E [ ψ ] − µNG [ ψ ] N = 1 π (cid:90) π dx | d ψ dx | − η κ N arctan (cid:16) ∆ c − NU π (cid:82) π dx | ψ | cos ( x ) (cid:17) κ − µ (cid:90) dx | ψ ( x ) | . (48)We perturb the wavefunction as ψ ( x ) = ψ ( x ) + δψ ( x ),where ψ extremizes G, i.e. one of the solutions thatwe obtained in Sec. V. Since ψ is an extremum, thefirst order variation in G vanishes and the second ordercontribution can be written as δG N = (cid:104) δψ | H | δψ (cid:105) + ρ (cid:104) δψ | cos ( z ) | ψ (cid:105) + ρ (cid:104) ψ | cos ( z ) | δψ (cid:105) + 2 ρ (cid:104) δψ | cos ( z ) | ψ (cid:105)(cid:104) ψ | cos ( z ) | δψ (cid:105) , (49)where H = − d dz + U n ph cos ( z ) (50)and ρ = η N U κ c − NU (cid:104) ψ | cos ( z ) | ψ (cid:105) κ (cid:16) ∆ c − NU (cid:104) ψ | cos ( z ) | ψ (cid:105) κ (cid:17) . (51)Equation (49) can be cast into a simple matrix form [18] δG N = 12 (cid:90) dx Ψ † ( x ) A Ψ( x ) . (52)Here Ψ( x ) = (cid:32) δψ ( x ) δψ ∗ ( x ) (cid:33) and A = (cid:16) H +2 ρ cos ( x ) ψ ( x ) I ∗ [ ... ] 2 ρ cos ( x ) ψ ( x ) I [ ... ]2 ρ cos ( x ) ψ ∗ ( x ) I ∗ [ ... ] H +2 ρ cos ( x ) ψ ∗ ( x ) I [ ... ] (cid:17) , where I [ . . . ] is an integral operator defined by I [ δψ ∗ ( x )] ≡ (cid:90) d x cos ( x ) ψ ( x ) δψ ∗ ( x ) . (53)The eigenvalues of the matrix A decide the energetic sta-bility. If A is positive definite, the solution ψ is energet-ically stable.In order to examine dynamic stability, we linearize theequation of motion Eq. (10) by writing ψ ( x, t ) = [ ψ ( x )+ δψ ( x, t )] e − iµt . This leads to i d δψ d t = (cid:20) − d d x + U | α | cos ( x ) − µ (cid:21) δψ ( x ) (54)+ 2 ρ cos ( x ) ψ ( x ) ( I ∗ [ δψ ( x )] + I [ δψ ∗ ( x )]) . −1 −0.5 0 0.5 10.20.30.40.50.60.70.80.91 quasi-momentum q (units of k c ) E [ ψ ] / N ( un i t s o f E R ) FIG. 22. (Color online) Energetic and dynamical stability ofthe band structure loops. The upper branch of the loop (blackdashed line) is energetically and dynamically unstable. Theother branches (red lines) are energetically and dynamicallystable. Parameters are, ∆ c = 3140 ω R , κ = 350 ω R , U = ω R , N = 10 , η = 2 . η , and η = 325 ω R . One can write a similar equation for δψ ∗ and combinethe two into a matrix equation similar to Eq. (52) i d δ Ψd t = σ z Aδ Ψ , (55)where σ z is the Pauli z -matrix. The solution ψ is dy-namically stable if all the eigenvalues of σ z A are real.Thus, the occurrence of complex eigenvalues of σ z A sig-nals dynamical instability. Before we quote the results,a comment is in order about the form of the pertur-bations δψ . The integral operator in Eq. (53) cou-ples the perturbation and ψ . If ψ = e iqx U q ( x ), with U q ( x ) = U q ( x + π ), i.e. a Bloch function with quasi-momentum q , the form of δψ that leads to non-zero cou-pling is δψ ( x ) = e iqz (cid:88) j b j e i jx . (56)That is, the perturbation should be a Bloch wave with thesame quasi-momentum as ψ . In Eq. (49) we consider thechange in the grand canonical potential per unit length,but the above choice is made to satisfy the requirementthat the integral operator I over the system size gives anon-zero answer. A physical way to motivate the abovechoice goes as follows: in the absence of interatomic inter-actions the allowed excitations of the quasi-momentumstate q have to be in multiples of the crystal momen-tum 2 k c (simply 2 in dimensionless units) since the only3source of perturbation is the interaction with the cavityfield. The above form respects this requirement. Also,the number of terms in the Fourier expansion Eq. (56)has to be less than the number of terms in the originalexpansion for ψ in Eq. (16) to avoid spurious instabili-ties [18]. Using Eq. (56), we find that the upper branchof the looped dispersions is both dynamically and ener-getically unstable as expected. The other two branchesare stable. This is shown in Fig. 22 for one particularcase. We shall not perform the stability analysis for thecase when there are five solutions, but anticipate by anextension of the case for the bistability scenario, that twoof the solutions will be dynamically unstable and threewill be stable. XI. SUMMARY AND CONCLUSIONS
In this paper we have analyzed bistability in atom-cavity systems in situations where the atoms are in spa-tially extended states (Bloch waves) with non-zero quasi-momentum q . We find that bistability in the number ofphotons in the cavity goes hand-in-hand with the emer-gence of loops in the band structure. Both are manifes-tations of a bifurcation in the stationary solutions of thecoupled atom-light equations of motion.We have studied how the loops appear and disappearas the laser detuning and the laser pumping rate arechanged. In particular, Eq. (31) provides an analyticalestimate of the critical pump strength η cr ( q ) at whichbistability sets in. It depends on the quasi-momentum ofthe atomic state, and predicts that loops first appear atthe edges of the first Brillouin zone ( q = ±
1) and thenmove inwards. This is indeed what we find upon solvingthe coupled atom-light equations numerically: swallow-tail loops appear at the edges of the first Brillouin zoneas the pump strength η is increased above η cr ( q = 1). As η is increased further the swallowtails extend inwards,merge, and detach from the rest of the band to form aseparate loop centred at q = 0 which ultimately closesup and vanishes. A rather similar behaviour is observedas the pump-cavity detuning ∆ c is swept from below thecavity resonance to above it.The loops we find are qualitatively different from thosethat occur for BECs in static optical lattices in the pres-ence of repulsive interatomic interactions [17–19]. There,the loops are centered at the edge of Brillouin zone andcause the dispersion to have a finite slope at that point.By contrast, the band structure we find always has zeroslope at the edge of the Brillouin zone. Nevertheless,there are also many similarities, including the stabilityof the various branches of the loops. We find that theupper branch of the loops are energetically and dynam-ically unstable, as expected from optical bistability con-siderations.The extra degree of freedom afforded by the quasi-momentum (over considering only q = 0) results in thepossibility of tristability, namely regions of parameter space where there are five solutions, three stable and twounstable. The complexity of the solutions in parame-ter space led us to perform an analysis of the problemin terms of catastrophe theory which is a useful mathe-matical tool for understanding the organization of bifur-cations of solutions. The key to our treatment was therecognition that, because exact solutions for the atomicwave functions are Mathieu functions which are speci-fied only by the lattice depth (once one chooses a quasi-momentum), the photon number n ph which determinesthe cavity depth provides a completely equivalent de-scription to the wave function.In the case of shallow lattices we were able to proceedanalytically and found that the structure of the solutionsin parameter space when q = 0 corresponds to a cuspcatastrophe, at the most singular point of which threesolutions (two stable and one unstable) form a pitch-fork bifurcation, and this describes the onset of bista-bility as the laser pumping is increased. Interestingly,the three experimental parameters { ∆ c , η, U } reducedto just two effective control parameters. In the generalcase of arbitrary lattice depth and 0 ≤ q ≤
1, the high-est singularities we found were swallowtail catastropheswhere four solutions simultaneously merge. The swallow-tails only exist when q > q sw = 0 . { ∆ c , η, q, U } reduced to three effective control pa-rameters meaning that generically one is unable to locatethe butterfly point (where five solutions simultaneouslymerge).The band structure loops found here have importantimplications for Bloch oscillations of atoms in cavities[30, 31]. Bloch oscillations are essentially an adiabatic ef-fect where, as a force is applied to the atoms they remainin the same band but their quasi-momentum evolves lin-early in time, as shown in Eq. (22). Swallowtail loops inthe band structure will have a deleterious effect on Blochoscillations because, as the quasi-momentum evolves, theatoms will reach the edge of a loop where the branchthey are following vanishes. This will lead to a sudden,non-adiabatic, disruption in the state of the atoms asthey are forced to jump to another branch or even an-other band. For BECs in ordinary static optical latticesthese non-adiabatic jumps are thought to be the cause ofthe destruction of superfluidity during Bloch oscillations[22, 23]. We have not included interactions in our treat-ment (interactions are necessary for superfluidity), butrelated effects will likely occur, especially considering theadded heating due to the fact that the lattice depth willalso abruptly change at the same point. However, whenthe loop detaches from the main band it will no longer af-fect Bloch oscillations. Furthermore, loops only occur forcertain limited regions of parameter space, i.e. inside thecusp and swallowtail catastrophes shown in Figs 14, 15,and 21. For experiments involving Bloch oscillations wetherefore recommend that parameter regimes are chosenwhich lie outside to these regions.4Finally, we add that although we have only consideredBloch waves in this paper, localized states (for exam-ple Wannier functions) can be formed from superposi-tions of Bloch waves with different values of the quasi-momentum. In this sense, localized states therefore con-tain all values of the quasi-momentum and so might beexpected to display tristability too. However, it shouldbe borne in mind that the nonlinearity of the systemmeans that superpositions of Bloch states of different q but the same lattice depth will not in general obey theeffective Schr¨odinger equation (10). ACKNOWLEDGMENTS
We gratefully acknowledge E. A. Hinds, D. Peli-novsky and M. Trupke for discussions. For funding,DHJO and BPV thank the Natural Sciences and Engi-neering Research Council of Canada, and the OntarioMinistry of Research and Innovation, and JL thanksVR/Vetenskapsr˚adet.
APPENDIX
We shall now sketch a proof showing that the function G defined in Eq. (42) produces swallowtail catastrophesbetween 0 . ≤ q ≤ C = { ∆ c = 0 . , η = 14 . , q = 0 . , U =0 . } and v = 7 . , n ph = 50 . κ and the number of atoms is set at N = 10 ). The numerical package [50] we used for cal-culating the derivatives gives error bounds allowing us toestimate the accuracy of our calculations. For the point { C , v } we found that the right hand side of Eq. (47) wasequal to 8 . × − with an error of 4 . × − . Thismeans that the third derivative of the state function van-ishes within error, indicating a swallowtail point. How-ever, the smallest value we found anywhere in parameterspace for the difference between the left hand side andthe right hand side of Eq. (46) was − . × − witherror 6 . × − . This means that the fourth deriva-tive of the state function does not vanish within error,suggesting there is no butterfly point. The value of thequasi-momentum at the point where we found this min-imal fourth derivative was q = q sw = 0 . n is the number of state variables, thenconsider a function p : R n → R• j k p is the Taylor expansion of p to order k • J k p is j k p minus its constant term • p is k -deteminate at 0 if any smooth function p + q ,where q is of order k + 1 (leading order term of the Taylor expansion is of order k + 1), can be locallyexpressed as p ( y ( s )) with y : R n → R n being asmooth reversible change of co-ordinates. • E kn is the vector space of polynomials in s , . . . , s n of degree ≤ k . • J kn is the subspace of E kn with zero constant term • ∆ k ( p ) is the subspace of J kn spanned by all Qj k (cid:16) ∂p∂s i (cid:17) k , where 1 ≤ i ≤ n , Q ∈ E kn , and thebar symbol represents the restriction to the k -thorder of the expansion. • The codimension of a function p is the codimen-sion of ∆ k ( p ) in J kn for any k for which p is k -determinate. • An r − unfolding of p at 0 is a function: P : R n + r → R , ( s , . . . , s n , t , . . . , t r ) (cid:55)→ P ( s, t ) = P t ( s ) , such that P ,... ( s ) = p ( s ).More informally, the term “unfolding” refers to how thecatastrophe unfolds as one moves away from the originin control space. At the origin in control space the catas-trophe reduces to its most singular part, known as its germ . For example, from Table I, the germ of the swal-lowtail catastrophe is given by s . The terms in the po-tential function which depend on the control parametersare called the unfolding terms, and the number of themis equal to the codimension. If P is an r − unfolding of p , set (cid:8) w k ( P ) , . . . , w kr ( P ) (cid:9) = (cid:26) ∂∂t (cid:0) J k ( P t , ,..., ) (cid:1) , . . . , ∂∂t r (cid:0) J k ( P ,...,t r ) (cid:1)(cid:27) . (57) W k ( P ) is the subspace of J kn spanned by { w k ( P ) , . . . , w kr ( P ) } .Referring to the standard forms given in Table I, thepotential function and state equation for the swallowtailare given byΦ( s ; C ) = 15 s + C s + C s + C s (58) G ( s ; C ) ≡ s + C s + C s + C = 0 . (59)Notice that the state equation for a swallowtail catastro-phe is similar to the potential function for a cusp catas-trophe up to a constant C . Since the state function G is the central object in the treatment given in SectionIX rather than the potential function, instead of provingthat the underlying potential function is equivalent tothat of a swallowtail catastrophe, we will prove that thestate function G around the singular point v and C isequivalent to the potential function of a cusp catastrophe5(note that this is different from the small photon numbercase we studied in Section IX B where we showed thatthe underlying potential was equivalent to the potentialfor a cusp catastrophe). To that end, first notice thatthe role of the constant term C in Eq. (59) is played by − η U in Eq. (42). Subtracting this function we havea modified form of G (where we have also dropped thedependence on q since we are focusing on a particularquasi-momentum): F ( v ; { ∆ c , U } ) = F ( v ; C ) = v + v (∆ c − N U f ( v, q )) (60)which satisfies F (cid:48) ( v ; C ) = F (cid:48)(cid:48) ( v ; C ) = F (cid:48)(cid:48)(cid:48) ( v ; C ) =0. In order to have a function defined in the neighbor-hood of v and C , let us set the origin of v at v and theorigin of C at C and define F ( v ; { ∆ c , U } ) ≡ F ( v + v , C + C ) − F ( v , C ) . (61)Thus, we have F (0 ,
0) = 0 and for the function g ( v ) ≡ F ( v, { , , } ) the most singular point is at v = 0 where g (cid:48) , g (cid:48)(cid:48) , g (cid:48)(cid:48)(cid:48) vanish. The function g is the germ which wedescribed above, and is the key feature which identifiesthe catastrophe. When g is Taylor expanded around 0one has g ( v ) = g ( iv ) v + g ( v ) v + O ( v ) (62)where g ( iv ) (0) is the first non-zero Taylor coefficient. Thismeans that g is 4-determined around 0 and we say that g ∼ v around 0. According to Table I, the canonicalunfolding of the 4-determined germ around 0 is the cuspcatastrophe Φ( s ; C ) = s / C s / C s , where v and s are related via a differomorphism (smooth transforma-tion of coordinates).Next we calculate the codimension of g . The Jacobianideal for g is in this case ∆ ( g ) = { v , v + g v g iv | v =0 v } .Hence, the codimension of g is dim( J ) − dim(∆ ( g )) =4 − F is thus a 2 parameter unfoldingof the germ g . In order to prove that the function F canbe described by a cusp catastrophe, we need to provethat F is isomorphic as an unfolding to the canonical form Φ( s ; C ) = s / C s / C s . In order to do thiswe need to invoke the idea of transversality.Transversality generalizes what we know of two inter-secting lines in a two dimensional plane to multidimen-sional manifolds. Two subspaces of a manifold are trans-verse if they meet in a subspace that is as small in di-mension as possible. If X (dim r ) and X (dim t ) aresubspaces of X (dim n ), X and X are transverse iftheir intersection is empty or if it is of the dimensionmax(0 , r + t − n ). Our first aim is to prove that the2 − unfolding of F is a versal unfolding. To do this weuse a defining theorem for versality from [36] which statesthat: an r − unfolding P of p , where p is k -determinateis versal if and only if W k ( P ) and ∆ k ( p ) (defined above)are transverse subspaces of J kn . We have already found∆ ( g ), the polynomial space W ( F ) is spanned by thevectors: w ( F ) = ∂∂U (cid:0) J ( F (0 , U , , )) (cid:1) ,w ( F ) = ∂∂ ∆ c (cid:0) J ( F (0 , , ∆ c )) (cid:1) , The expressions depend on the derivatives of the couplingfunction f ( v, q ) and the value of the parameters at thesingular point { v , C } . They are too cumbersome tostate here but their general forms are given by w i ( F ) = (cid:88) j =1 .. z ij v j , which we determined numerically and all of the z ij ’sare non-zero. The polynomials w i are linearly indepen-dent which gives the dimensionality dim( W ( F )) = 2.Furthermore, we have verified that the rank of the ma-trix formed by the polynomial coefficieints of ∆ ( g )and W ( F ) is 4 and this combined with the fact thatdim(∆ ( g )) + dim( W ( F )) = 2 + 2 = dim( J ) provesthat ∆ ( g ) and W ( F ) are transverse. Thus, by the the-orem stated above F is a versal unfolding of the germ g and since it is a 2 − unfolding (codimension of g = 2) itis also universal [36]. This proves the equivalence of theunfolding of F to the cusp catastrophe. [1] H. M. Gibbs, S. M. McCall, and T. N. C. Venkatesan,Phys. Rev. Lett. , 1135 (1976).[2] R. Bonifacio and L. A. Lugiato, Phys. Rev. Lett. , 1023(1978); R. Bonifacio and L. A. Lugiato, Phys. Rev. A ,1129 (1978).[3] P. Meystre and M. Sargent, Elements of Quantum Optics ,3rd edition, (SpringerVerlag, Berlin, 1997).[4] C.W. Gardiner and P. Zoller,
Quantum Noise , (Springer,Berlin, 2004).[5] S. Gupta, K. L. Moore, K. W. Murch, and D. M.Stamper-Kurn, Phys. Rev. Lett. , 213601 (2007). [6] F. Brennecke, S. Ritter, T. Donner, and T. Esslinger,Science , 235 (2008); S. Ritter, F. Brennecke, K.Baumann, T. Donner, C. Guerlin and T. Esslinger, App.Phys. B , 213 (2009).[7] C.J. Hood, M. S. Chapman, T. W. Lynn, and H. J. Kim-ble, Phys. Rev. Lett. , 4157 (1998); P. M¨unstermann,T. Fischer, P. W. H. Pinkse, and G. Rempe, Opt. Comm. , 63 (1999); M. Trupke, J. Goldwin, B. Darqui´e, G.Dutier, S. Eriksson, J. Ashmore, and E. A. Hinds, Phys.Rev. Lett. , 063601 (2007).[8] J. Ye, D. W. Vernooy, and H. J. Kimble, Phys. Rev. Lett. , 4987 (1999); C. J. Hood, T. W. Lynn, A. C.Doherty, A. S. Parkins, and H. J. Kimble, Science ,1447 (2000); P. W. H. Pinkse, T. Fischer, P. Maunz, G.Rempe, Nature , 365 (2000).[9] C. Maschler and H. Ritsch, Phys. Rev. Lett. , 260401(2005); I. B. Mekhov, C. Maschler and H. Ritsch, Nat.Phys. , 319 (2007); C. Maschler, I. B. Mekhov, and H.Ritsch, Euro. Phys. J. D , 545 (2008).[10] J. Larson, B. Damski, G. Morigi, and M. Lewenstein,Phys. Rev. Lett. , 050401 (2008); J. Larson, S.Fernandez-Vidal, G. Morigi, and M. Lewenstein, New J.Phys , 045002 (2008).[11] J. M. Zhang, F. C. Cui, D. L. Zhou, and W. M. Liu,Phys. Rev. A , 033401 (2009).[12] D. Nagy, P. Domokos, A. Vukics, and H. Ritsch, Eur.Phys. J. D, , 659 (2009).[13] L. Zhou, H. Pu, H. Y. Ling, and W. Zhang, Phys. Rev.Lett., , 160403 (2009).[14] D. Nagy, G. Szirmai, and P. Domokos, Euro. Phys. J.D , 127 (2008); S. Fernandez-Vidal, G. De Chiara, J.Larson, and G. Morigi, Phys. Rev. A , 043407 (2010);D. Nagy, G. Konya, G. Szirmai, and P. Domokos, Phys.Rev. Lett , 130401 (2010).[15] K. Baumann, C. Guerlin, F. Brennecke, and T. Esslinger,Nature , 1301 (2010).[16] O. Morsch and M. Oberthaler, Rev. Mod. Phys. , 179(2006).[17] B. Wu and Q. Niu, Phys. Rev. A , 061603(R) (2001);B. Wu, R. B. Diener, and Q. Niu, Phys. Rev. A ,025601 (2002); B. Wu and Q. Niu, N. J. Phys. , 104(2003).[18] Dmitri Diakonov, L. M. Jensen, C. J. Pethick, and H.Smith, Phys. Rev. A , 013604 (2002); M. Machholm,C. J. Pethick, and H. Smith, Phys. Rev. A , 053613(2003); M. Machholm, A. Nicolin, C. J. Pethick, and H.Smith, Phys. Rev A , 043604 (2004).[19] E. J. Mueller, Phys. Rev. A , 063603 (2002).[20] A. Smerzi, A. Trombettoni, P. G. Kevrekidis, and A. R.Bishop, Phys. Rev. Lett. , 170402 (2002).[21] C. J. Pethick and H. Smith, Bose Einstein Condensa-tion in Dilute Gases , 2nd edition (Cambridge UniversityPress, Cambridge, 2008).[22] S. Burger, F. S. Cataliotti, C. Fort, F. Minardi, M. In-guscio, M.L. Chiofalo and M. P. Tosi, Phys. Rev. Lett. , 4447 (2001).[23] S. Cataliotti, L. Fallani, F. Ferlaino, C. Fort, P Mad-daloni, and M. Inguscio, New J. Phys. , 71 (2003).[24] B. Wu and Q. Niu, Phys. Rev. A , 023402 (2000).[25] Y.-A. Chen, S. D. Huber, S. Trotzky, I. Bloch, and E.Altman, Nat. Phys. , 61 (2011).[26] R. Battesti, P. Clad´e, S. Guellati-Kh´elifa, C. Schwob,B. Gr´emaud, F. Nez, L. Julien, and F. Biraben, Phys.Rev. Lett. , 253001 (2004); P. Clad´e, E. de Miran-des, M. Cadoret, S. Guellati-Kh´elifa, C. Schwob, F. Nez,L. Julien, and F. Biraben, Phys. Rev. Lett. , 033001(2006).[27] G. Roati, E. de Mirandes, F. Ferlaino, H. Ott, G. Mod-ugno, and M. Inguscio, Phys. Rev. Lett. , 230402(2004).[28] G. Ferrari, N. Poli, F. Sorrentino, and G. M. Tino, Phys.Rev. Lett. , 060402 (2006).[29] I. Carusotto, L. P. Pitaevskii, S. Stringari, G. Modugno,and M. Inguscio, Phys. Rev. Lett. , 093202 (2005).[30] B. M. Peden, D. Meiser, M. L. Chiofalo, and M. J. Hol- land, Phys. Rev. A , 043803 (2009).[31] B. Prasanna Venkatesh, M. Trupke, E. A. Hinds, and D.H. J. O’Dell, Phys. Rev. A , 063834 (2009).[32] M. Gustavsson, E. Haller, M. J. Mark, J. G. Danzl, G.Rojas-Kopeinig, and H.-C. Nagerl, Phys. Rev. Lett. ,080404 (2008).[33] R. Thom, Structural Stability and Morphogenesis , (Ben-jamin, London, 1975).[34] V. I. Arnold, Russ. Math. Surveys , 1 (1975).[35] Michael Berry Singularities in Wave and Rays in LesHouches, Session XXXV, 1980
Physics of Defects , editedby R. Balian et al. (North-Holland Publishing Company,Amsterdam, 1981).[36] T. Poston and I. Stewart,
Catastrophe Theory and ItsApplications , (Dover Publications, NewYork, 1996).[37] J. F. Nye
Natural Focusing and Fine Structure of Light ,(Institute of Physics, Bristol, 1999).[38] R. Gilmore and L. M. Narducci, Phys. Rev. A ,1747(1978).[39] G. P. Agrawal and H. J. Carmichael, Phys. Rev. A ,2074 (1979).[40] Langevin fluctuation terms do not appear in Eq. (5)if the bath is at thermal equilibrium. See p 343 ofC. Cohen-Tannoudji, J. Dupont-Roc and G. Grynberg, Atom-Photon Interactions (Wiley, New York, 1992).[41] J. Larson, G. Morigi, and M. Lewenstein, Phys. Rev. A , 023815 (2008).[42] F. Bloch, Z. Phys. , 555 (1928).[43] C. Zener, Proc. R. Soc. London, Ser. A , 523 (1934).[44] M. Abramowitz and I. Stegun, Handbook of MathematicalFunctions (National Bureau of Standards, Washington,1964).[45] M. Ben Dahan, E. Peik, J. Reichel, Y. Castin, and C.Salomon, Phys. Rev. Lett. , 4508 (1996); O. Morsch, J.H. M¨uller, M. Cristiani, D. Ciampini, and E. Arimondo,Phys. Rev. Lett. , 140402 (2001).[46] T. J. Kippenberg and K. J. Vahala, Science , 1172(2008); F. Marquardt and S. M. Girvin, Physics , 40(2009).[47] Y. Dong, J. Ye, and H. Pu, Phys. Rev. A , 031608(R)(2011).[48] L. D. Landau and E. M. Lifshitz, Statistical Physics / ( NU ) = 0, and q can take anyvalue. This is a rather unphysical situation because, as-suming a finite number of atoms, it can only occur when U → ∞ and this then means v → ∞ . From Fig. 5, wesee that when v → ∞ the function f ( v, q ) becomes veryflat and so all its derivatives vanish allowing Eqns (43)-(46) to be satisfied in a rather trivial way. Therefore, weignore this solution.[52] A. E. R. Woodcock and T. Poston, A Geometrical Studyof the Elementary Catastrophes , (Springer-Verlag, Berlin,1974).[53] P. Horak and H. Ritsch, Phys. Rev. A , 023603 (2001).[54] J. Gaite, J. Margalef-Roig, and S. Miret-Art´es, Phys.Rev. B57