Temperature Variation of Ultra Slow Light in a Cold Gas

A model is developed to explain the temperature dependence of the group velocity as observed in the experiments of Hau et al (Nature {\bf397}, 594 (1999)). The group velocity is quite sensitive to the change in the spatial density. The inhomogeneity in the density and its temperature dependence are primarily responsible for the observed behavior.


I. INTRODUCTION
The phenomenon of Bose-Einstein condensation in atomic gases [1,2] lends itself to the study of many fundamental effects. Among them, one aspect presently being investigated both theoretically and experimentally is the interaction of light with atoms in the quantum degeneracy regime [2]. In this context, the propagation of light inside a cold gas is still an open problem. Because of the optical density, it is well known that the transmission of resonant light through a condensate is almost zero [3]. However, electromagnetically induced transparency (EIT) [4] was found to allow the propagation of light by means of quantum coherence between different internal atomic levels [5,6]. In this context, Hau et al discovered a remarkable property of pulse propagation in a Bose condensate. These authors demonstrated the slowing down of the group velocity of the pulse to 17 meters/sec [7]. Furthermore, they have shown a definite dependence of the group velocity on the temperature of the ultra-cold sample. One would like to understand the observed temperature dependence from first principles. For this purpose, it is necessary to extend the standard theory of EIT to a cold gas at finite temperature. However, a theoretical description of this problem is rather complex. Complexities arise when one attempts a systematic treatment of interactions, finite temperature effects and dynamics. Most studies treat these aspects as disjoint: interactions are included in the zero temperature case to study the kinematical aspects [2,8], whereas some dynamical aspects are studied using only the excitations within the electronic ground state [9], and finite temperature effects are usually studied for noninteracting bosons [2,10,11]. A complete theory should study all these aspects together. However, a complete theory of the interaction of light and interacting particles is still unavailable, and a full numerical treatment is a rather hard task. Here, we present approximate but plausible arguments to explain the experimental observations in [7]. The simplicity of our model allows for an analytical expression for the group velocity in the following cases: atoms confined in a box and by a harmonic potential. We obtain results which reproduce the ones in [7] for T > T c . In particular, the treatment brings out the factors playing key roles in the phenomenon. Here, we show that the variation of spatial density of atoms with temperature is the major factor responsible for the temperature dependence of the group velocity. The paper is organized as follows: In Section II the model is introduced. In Section III we derive the group velocity of a pulse propagating in an ideal gas confined inside a box, extend the calculation to the case of an ideal gas in a harmonic oscillator potential, and present and discuss the results in relation to the experiment of Hau et al. In Section IV we present estimates for the group velocity in the interacting case and in the limit of zero temperature.

II. THE MODEL
In this Section we introduce the model used throughout this article. Here, we write the Maxwell-Bloch equations which describe the dynamics of the system consisting of light field and atoms. We derive the linear response of the medium to a weak probe field, taking into account the quantum statistics of the atoms. The group velocity is then defined in the standard manner [12].

Maxwell-Bloch Equations
We consider a gas of N non-interacting bosons. The relevant internal structure corresponds to a three-level atom, with internal levels |g (stable state), |r , (metastable state) and |e excited state, whose energies are ω g , ω r and ω e , respectively (see Fig. 1). The radiative decay rate of the excited state is γ = γ g + γ r , with γ g (γ r ) the rate of decay on the transition |e → |g (|e → |r ). Laser light with frequency ω lg and wave vector k g drives the transition |g → |e , whereas the transition |r → |e is driven by a field of frequency ω lr and wave vector k r . The dynamics of the whole system is given by the Maxwell equation for the electric field vector E and by the optical Bloch equations for the density matrix equations of the N -atom gas. For non-interacting atoms it suffices to consider the equations for the one-atom density matrix ρ, projected on the basis {|j, ǫ } with j = r, g, e and |ǫ the eigenvector of the mechanical motion of one atom at the energy ǫ. They have the form: +i +i where ρ ij (ǫ, ǫ ′ ) = i, ǫ|ρ|j, ǫ ′ (i, j = r, e, g), C j ǫ,ǫ ′ = ǫ| exp(ik j · r)|ǫ ′ ,ρ ej = ρ ej e −iω lj t ,ρ rg = ρ rg e −i(ω lg −ω lr )t and ρ ij = (ρ ji ) * for i = j. Here, c is the speed of light, P is the polarization of the medium, ∆ 0 j = ω e − ω j − ω lj (j = g, r) is the detuning. Rabi couplings are given by g = |d eg · E|/h and Ω = |d er · E|/h, where d ej is the dipole moment of the transition |e → |j . Finally, ρ j ee (ǫ, ǫ ′ ) describes the density matrix after a spontaneous emission event on the transition |e → |j : where α 1 (k) and α 2 (k) form a set of polarization vectors orthogonal tok, andd je = d je /|d je |. In deriving the above equations, we made the rotating wave approximation and transformed to a reference frame rotating at the optical frequency of the laser. Furthermore, in Eqs. (2)-(7) we have introduced the loss-rates Γ ij , which take into account the effects of other mechanisms of decoherence. In the ideal case Γ ge = Γ re = γ/2, whereas Γ gr = 0.

Susceptibility and Group velocity
In an isotropic medium, the linear susceptibility χ is defined by the expression [13] Assuming that two light fields are propagating through the medium along theẑ-direction, we can write the electric field E(r, t) and the atomic polarization P (r, t) as E(r, t) = j=g,r [E j 0 (r, t) exp(+ik j z − iω lj t) + c.c.]/2 and P (r, t) = j=g,r [P j 0 (r, t) exp(+ik j z − iω lj t) + c.c.]/2 where E j 0 , P j 0 are respectively the slowly-varying envelopes of the electric field and atomic polarization at frequency ω lj . The Fourier transform of Eq. (9) gives P 0 (ω) = E 0 (ω)χ(ω). From the relation P(t) = Tr{ρ N (t)d} for the polarization with d the atomic dipole moment operator, we find that the macroscopic polarization of the medium at the position r is where ρ N is the N -atom density matrix, that has the following form in the energy representation.
The N-atom optical coherence density matrix has to be obtained from the solution of Eqs.
(2)-(7) subject to the initial condition: where N (ǫ) is the number of atoms in the ground state with energy ǫ. The Eqs.
(2)-(7) are to be solved to first order in the field E g 0 and to all orders in the field E r 0 . In this work we are interested in the steady state of the atoms with the field, which is assumed to be reached on a time-scale much shorter than the thermalization time-scale of the gas. On the basis of this hypothesis, we assume that the initial condition (12) and Eqs. (2)-(7) determine the steady state solution. Note that the coefficients C ǫ,ǫ ′ determine the one-atom energy states involved in the transition induced by the laser field. For free bosons these coefficients have a simple form, as it can be seen in section III.A. For harmonic oscillator potentials the coefficients C ǫ,ǫ ′ are given in terms of Laguerre polynomials, where the number of vibrational states which are coupled depends on the ratio between the recoil frequency over the trap frequency. In section III.B we use an approximate treatment for this case. Once the susceptibility is known, the dispersion relation of light in the medium is given [13] and we can evaluate the group velocity, defined as v g = ∂ω ∂kg | ω=ωeg . In the limit N χ ≪ 1 the group velocity has the form: where χ ′ = Re(χ).

III. EVALUATION OF THE GROUP VELOCITY
In this section we derive an analytical expression for the group velocity of a laser pulse propagating through an ideal gas of ultracold atoms. We investigate two cases: atoms in a box and atoms confined by a harmonic oscillator potential. Finally, we do the numerical calculations for the case of a gas of sodium atoms and discuss the results in relation with the experimental data of [7].

A. Group velocity in a gas of free non-interacting bosons
We consider a gas of N bosons in a box of volume V . In this case the atomic wave vector eigenstates |k are also energy eigenstates with eigenvalues ǫ =h 2 k 2 2m , where m is the atomic mass. So, we project Eqs. (2)-(7) on the motional basis {|k }. Then the coefficients C j ǫǫ ′ appearing in the density matrix equations have the form C j ǫǫ ′ ≡ C j k,k ′ = δ k,k+k j . We substitute these values into Eqs. (2)-(7) and solve the equations in the steady state limit. To first order in g/Ω and g/Γ, and assuming that at t = 0 the gas is in thermal equilibrium, the steady-state optical coherence ρ eg is found to be where ∆ j is the detuning defined as with∆ 0 j = ∆ 0 j +ω R , and ω R is the recoil frequency defined as ω R =hk 2 /2m. Using (14) in (10), we find the expression for the susceptibility The sum in Eq. (16) is over all the motional states weighted by their statistical occupation N (k), where f is the fugacity, β = 1/K B T and T is the temperature, and k = |k|. Here, χ 0 is the one-atom susceptibility, defined as where λ is the optical wavelength of the transition g → e, λ = 2πc/ω ge , and Γ ge = γ/2. In the following, we assume that |k g − k r | ≪ k g , k r , Γ gr . Therefore, the dependence on k in the denominator of Eq. (16) is mainly due to the first term ∆ g , and we may rewrite Eq. (16) as where k z = k ·ẑ and ζ is a complex number independent of k; After evaluating the susceptibility as given by Eq. (19), the group velocity can be found using Eq. (13).
In the following, we investigate the behaviour of the cloud close to the critical point, dividing our investigation into two regimes: above and below the critical temperature T c .

Above the critical temperature
Above the critical temperature and in the limit of large volumes, one can replace the sum in Eq. (19) with an integral in three-dimensions [14]. The expression to evaluate is now where for convenience we have chosen to integrate in the cartesian coordinates. We write N (k) as and use it in Eq. (21) to obtain This can be written in terms of the standard functions where and where the function w is defined as Given the critical temperature for an ideal Bose gas where n(= N/V ) is the density of atoms, we can rewrite Eq. (24) as with A c = 2 k g Γ geh m n g 3/2 (1) and

Below the critical temperature
For T < T c , we use the expression for the ground state population in the thermodynamic limit [14], to obtain where the second term on the RHS describes the contribution of the condensed phase.

Regime of parameters and approximations
For currently studied optical transitions the argument of the w function in Eq. (28) is y = √ lζ/A ≫ 1 for any value of l ≥ 1. Therefore, the asymptotic expansion of the w function can be applied [16]: We substitute this expansion into Eq. (28), and obtain for T > T c where we have used the relation g 3/2 (f )/g 3/2 (1) = T c /T for T > T c . For T < T c the susceptibility will now have the form Using Eqs. (33) and (34)in the formula (13) we find the group velocity. Note that the dependence on the temperature comes in at higher order in the expansion A/ζ. Clearly a significant temperature dependence for a free gas can come only for narrow optical transitions.

B. Group velocity in a gas of trapped non-interacting bosons
Let us now consider a cloud of atoms trapped in a three-dimensional harmonic potential with cylindrical symmetry, so that the one-atom Hamiltonian describing the mechanical motion has the form where V (r) is the harmonic oscillator potential in cylindrical coordinates with ν r , ν z trap frequencies in the radial and axial directions, respectively. In order to evaluate the susceptibility in the steady state, we solve Eqs. (2)-(7) in the semiclassical limit for the atomic motion, and we sum over the states using the semiclassical statistical distribution [2,15]. This limit is valid when treating the non-condensed fraction of atoms for temperatures T fulfilling the condition K B T ≫hν, and under the condition Γ, ∆ ≫ ν. The hypothesis is justified in the range of parameters of [7] and simplifies considerably the treatment, allowing for an analytical solution of the group velocity. Then the coeffients C ǫ,ǫ ′ simplify to their semiclassical values C ǫ,ǫ ′ ≈ δ p,p ′ +hk δ r,r ′ , where p, r are now the classical canonical coordinates of a harmonic oscillator with energy E = p 2 /2m + V (r). In this limit, the optical coherence ρ eg appearing in Eq. (10) has the form: with ∆ j defined in Eq. (15), and the susceptibility is given by the expression: where the semiclassical statistical distribution is When considering the condensate contribution to the optical susceptibility, one should evaluate ρ eg (ǫ, ǫ ′ ) and sum over the final states with energy ǫ. However, in the regime ω R /ν ≫ 1 we may apply the semiclassical approximation to the final states. The final semiclassical energy is the recoil energy, and we can write the optical coherence for the condensate contribution as .

(40)
The ground state occupation is given in the thermodynamic limit by where T c is the critical temperature of the trapped gas, K B T c =h(ν z ν 2 r ) 1/3 (N/g 3 (1)) 1/3 . Contrary to the case of free bosons, the group velocity is not directly given by the formula (13), because of the spatial variation of the atomic density and therefore of the susceptibility. Here, we evaluate the group velocity using a method equivalent to the experimental one of [7], i.e. we estimate the size D of the cloud and calculate the delay ∆t of a pulse propagating across a selected region of a cold gas with respect to a pulse propagating in the vacuum. The group velocity is then given by the ratio of the size over the delay v exp g = D/∆t. Assuming that the light is propagating along theẑ-axis and cuts a cylinder inside of the volume with section S and centered on theẑ-axis of the cloud, we write the delay ∆t as the average over the section S of all the delays ∆t(r) of pulses propagating at distance r from the axis of the cloud where R is the radius of the illuminated circular section S of the cloud and ∆t(r) is defined as where v g is defined in Eq. (13), and L(r) is half the length of the path along the cloud. Note that R will, in principle, depend on the size of the incoming Gaussian beam. However, in the experiment of [7], R is the radius of a pinhole set before the measuring apparatus. The delay time is experimentally obtained by measuring the difference between the delay time of the pulse propagating across the cloud and the one of a pulse propagating in the vacuum. Assuming that L(r) = L is the distance between a slit before the cloud and the photomultiplier, the final delay will be ∆t = ∆t − L/c. In the following, we evaluate the group velocity as a function of the temperature above and below the critical temperature.

Above the critical temperature
The integral over the momenta in Eq. (38) can be evaluated along the lines of the procedure outlined in Eqs. (21)-(24). One finds Again, the considerations on the w function made in the free case are applicable, and using its asymptotic expansion [Eq.(32)], one gets · g 3/2 f e −βV (r) + g 5/2 f e −βV (r) A 2 ζ 2 .
For L ≫ D z (T ), where D z (T ) is the axial thermal size of the cloud, we can replace L by ∞ in the integral (43). Therefore, the delay of a beam propagating along the z-axis, and entering the cloud at a distance r from the cloud axis is where a 0j = h/mν j is the size of the ground state of the harmonic oscillator in the j direction (j = r, z). In order to evaluate the size of the cloud, we consider that for T < T c , a fraction (T /T c ) 3 of the atoms is outside of the condensate, whereas a fraction 1 − (T /T c ) 3 is in the condensate. Applying the semiclassical approximation to the non-condensate part, we have that which leads us to defining the size of the cloud to be: Dividing D z by ∆t − we find the group velocity below the critical temperature.

C. Numerical Results
In Fig. 2 we plot the group velocity of a gas of sodium atoms in a box (dashed line) and in a harmonic oscillator (dotted line) as a function of temperature T , scaled according to the critical temperature of each case. Density of atoms, number of atoms and trap frequencies have been taken from the data of [7]. The experimental results of [7] are seen to be broadly in agreement with the harmonic oscillator case. The inhomogeneous spatial density of the atoms and its variation with temperature is the key to the understanding of the experimental data. The curve representing the case of free atoms shows that the temperature dependence entering into Eq. (28) as a higher order correction has a negligible effect on the considered scale, and cannot be interpreted as the cause of the behaviour observed in [7]. In Fig.3, we compare the group velocity for two different values of the Rabi frequency Ω coupling |r to |e . The behaviour for T > T c is similar to the corresponding one measured in [7]. The curves we obtain are however steeper, and this can be explained by considering the approximations made in our treatment. In our calculations we have assumed the same number of atoms at every temperature. However, in the experiments lower temperatures are achieved by means of evaporative cooling. This implies that the points of the experimental curve at higher temperatures correspond to larger numbers of atoms, and correspondingly to larger spatial density (for ideal gases). This leads to a smoother gradient of the group velocity versus the temperature than in our case. On the other hand, as the temperature decreases, the effect of the interactions gets stronger causing, among other effects, a lower density of the atoms than in the non-interacting case. Hence, one would expect a group velocity value larger than the evaluated one. Albeit these considerations, the evaluated curve reproduces the experimental one above the critical temperature with some agreement, showing that the ideal gas model provides a qualitative description of the phenomenon. Below the critical temperature the discrepancy between the experimental data and our theoretical predictions is rather dramatic. This is not surprising since the size of the condensate is strongly affected by the effect of the interactions. Already Ketterle and co-workers have reported that the cloud size is much larger in the interacting system compared to the size of the harmonic oscillator ground-state wave function [18]. Therefore, our evaluation can be expected to lead to smaller values of the group velocity than the experimental records. In order to illustrate this point, in the following section, we estimate the group velocity at T = 0 by comparing the ideal case with the Thomas-Fermi case. Finally, we discuss the measure of the group velocity in the two limiting cases for a section S with radius R ≪ D r (T ) and with radius R ≈ D r (T ). This is illustrated in Fig. 4, where the same dependence of the group velocity on the temperature is evident. The orders of magnitude of the pairs of curves corresponding to the same set of parameters are comparable, showing that the behaviour observed in [7] originates mainly from a change in the "average" spatial density of the gas with temperature.

IV. GROUP VELOCITY FOR AN INTERACTING BOSE GAS
In this section, we compare the group velocity value at T = 0 in the two limits: the ideal one, where we consider the particles as non-interacting, and the interacting case, which we treat in the Thomas-Fermi approximation. We estimate the group velocity using the set of parameters of the experiment and the formula (1) of [7]: where n is the density. Therefore, we need to evaluate the group velocity at T = 0 by substituting into Eq. (55) an estimate of the spatial density, which we calculate here as the ratio of the total number of atoms over the volume of the cloud. This evaluation, which corresponds to considering the density as homogeneous, is justified on the basis of the results of Fig. 4, where it is shown that the phenomenon observed in [7] is mainly dependent on the change in the density. For an ideal gas in a harmonic oscillator potential at T = 0, all the atoms are in the ground state, and a rough estimate of the density gives n ≈ N/(4πa 0z a 2 0r /3). Taking N = 10 6 Sodium atoms and ν z = 20 × 2π Hz, ν r = 70 × 2π Hz, Ω = 0.56γ, the ground state dimensions are a 0z ≈ 4.7µ and a 0r ≈ 2.4µ, and we obtain a density n ≈ 8 × 10 15 atoms per cm 3 . Thus, according to (55), the group velocity is v ideal g ≈ 0.03 m/sec. For an interacting gas in the Thomas-Fermi limit, the cloud is an ellipsoid of axes 2R T F r in the radial direction and 2R T F z in the axial direction, where R T F j is the Thomas-Fermi radius: and µ is the chemical potential, defined as with a S scattering length, ν ho = (ν 2 r ν z ) 1/3 geometrical average of the oscillator frequencies and a ho = h/mν ho corresponding oscillator size. Taking a S = 2.75 nm, for the set of parameters of the experiment the Thomas-Fermi dimensions of the cloud are R T F z ≈ 47.4µ and R T F r ≈ 13.6µ . Considering the density of atoms as homogeneous, we obtain n ≈ 3 × 10 13 atoms per cm 3 . From Eq. (55) we find for the group velocity v TF g ≈ 9 m/sec, which is comparable with the value measured in [7] for temperatures below the critical temperature. Therefore, for 10 6 atoms we find a difference of two orders of magnitude in the value of the group velocity between the ideal case and the interacting case. Such difference increases or decreases depending on the total number of atoms in the trap. This estimate substantiates the inference that interactions are responsible for a lower density, and therefore, for a higher average group velocity of the light.

V. CONCLUSIONS
We have derived an approximate analytical expression for the group velocity of a pulse propagating through an ultracold gas which is confined in a box and by a harmonic potential. We have shown that the results reproduce qualitatively the experimental ones presented in [7]. From our analysis it emerges that the definite variation of the group velocity with the temperature of the gas is an effect related to the variation of the spatial density of the gas. We see that the ideal gas model provides a qualitative description of the results for T > T c . However, the behaviour at T < T c can be described in a satisfactory way only by including the interactions and the fact that the cloud is cooled by means of evaporative cooling. The last one has the effect of making the total number of atoms temperature-dependent. Such effects will be the subject of future investigations.

VI. ACKNOWLEDGEMENTS
One of us (GSA) thanks S.E. Harris, M.O. Scully and P. Meystre for interesting discussions on the subject matter of this paper. G.M. wishes to thank J. Schneider for many helpful comments. FIG. 3. Onset: Plot of the group velocity in m/sec in logarithmic scale as a function of the effective temperature θ = T /Tc for a gas of sodium atoms as in [7]. The upper curve corresponds to Ω = 1.2γ, whereas the lower curve corresponds to Ω = 0.56γ. Here, Tc = 432 nK, N = 8.3 × 10 6 , Γgr = 2π × 1000 Hz, νr = 2π × 70 Hz, νz = 2π × 20 Hz and ∆g = ∆r = 0. Inset: Plot of the low temperature behaviour of the corresponding curves in linear scale. The radius of the section S is R=15 µ.