text:
# Non-Gaussian Velocity Distributions in Optical Lattices

## I Introduction

Laser cooling is now a well established technique to produce narrow velocity distributions for dilute samples of atomic gases (see e.g. nobel97 ). The interaction between the atoms and the radiation modes removes kinetic energy from the atoms, and extremely cold samples can be obtained. In the standard context of Doppler or sub-Doppler laser cooling, atom-atom interactions are neglected and hence a thermodynamic temperature cannot be defined. Nevertheless, measured velocity distributions are generally very well fitted by a Gaussian function, and assigning a ‘kinetic temperature’ to the distribution is a useful way to characterize a laser cooled atomic sample.

One of the simplest theoretical models of laser cooling assumes a moving two-level atom interacting with counter-propagating pairs of laser beams, tuned slightly below the atomic resonance (Doppler cooling doppler ). This will yield Doppler shifts, asymmetric with regards to velocity, and thus a damping force (friction). Doppler cooling is counteracted by momentum diffusion due to absorption and emission of photons. If a spatial average is taken of diffusion as well as friction, one obtains a stationary Gaussian velocity distribution. This is valid since, in steady-state, most atoms have velocities well above spatial modulations in the light shift potential (caused by the interaction between the induced dipole moment and the light), and thus the dynamics can be described in terms of a Fokker-Planck equation with constant friction and diffusion coefficients. High irradiance results in light shifts of the involved energy levels that can be comparable to the kinetic energy, and one can no longer assume a constant velocity as atoms travel over a wavelength. Spatial averaging can still be performed, but one does not obtain the standard description of laser cooling in terms of competition between a friction force and a diffusion effect, since these are not simply functions of velocity. The resulting velocity distribution will in this case not be Gaussian and different distributions have been proposed Kazantsev . However, for practical Doppler cooling configurations, this effect is negligible, and there are no known observations of clearly non-Gaussian distributions.

For a multilevel atom, population transfer and coherences between degenerate levels open up the possibility for more subtle cooling mechanisms. These are not limited by the radiative lifetimes of the upper levels, and can therefore lead to narrower distributions. In particular, Sisyphus cooling Dalibard89 ; Chu891 ; Chu892 ; Phillips89  is based on a laser beam configuration that results in a periodic modulation of the polarization of the light, and thus spatially modulated optical pumping and steady-state population distribution between different degenerate substates. The light shift will also be periodic, and will differ for different substates. The combination of hamiltonian motion and optical pumping cycles transfers atomic energy to the vacuum modes Dalibard89 ; Chu891 ; cohen90 . A rule of the thumb for Sisyphus cooling tells us that the ‘temperatures’ obtained correspond to kinetic energies that are of the order of the light shift. This behavior has been experimentally verified Salomon90 ; Gatzke ; Jersblad ; Ellmann  down to kinetic temperatures of a few recoil energies. A seminal analysis of Sisyphus cooling, by Dalibard and Cohen-Tannoudji Dalibard89 , is again based on spatially averaged friction and diffusion coefficients. Even though the final regime corresponds to a situation where one can no longer assume atoms moving at constant velocity over many wavelengths, the scaling law obtained by this approach appears to be excellent.

In more rigorous full quantum mechanical analyses, Castin et al. Castin ; Castin91  find that Sisyphus cooling ought to lead to non-Gaussian distributions. In particular, for irradiances close to the lower limit for efficient laser cooling, the effects of recoils due to absorbed and emitted photons become prominent. Then, atomic trajectories become very irregular and the velocity cannot be assumed to be constant. Therefore one cannot compute a spatially averaged velocity dependent force. Also, the atoms will be trapped in microscopic potential minima (forming optical lattices mennerat02 ; jessen ), and the ensemble should be characterized by a distribution of vibrational modes and unbound modes, rather than by a velocity distribution.

Essentially all experimental investigations of Sisyphus cooling result in distributions that are well fitted by Gaussians. The reason for this is probably a combination of several facts. Many experiments are done in a regime where an average friction coefficient seems adequate (sufficiently large light shift). The deviations from Gaussian distributions are subtle and are mainly hidden in the noisy wings of the recorded distribution. Furthermore, it is difficult to set-up an experimental velocity probe with the required resolution. Nevertheless, deviations from Gaussian velocity distributions for laser cooled atoms have been reported in one recent paper Salomon . However, to our knowledge, there has been no systematic experimental study of the non-Gaussian distributions, nor any attempts to approach the observed distributions with more precise functions.

In this work, we report a detailed study of velocity distributions, as a function of the irradiance (and thus the light shift) for a three dimensional Sisyphus cooling configuration. We also perform a one-dimensional numerical simulation of velocity distributions, based on a full-quantum Monte-Carlo wave function technique. This is applied for the atomic angular momentum which is relevant in our experiment. We fit the recorded data, the experimental as well as the numerical, to different functions and compare the outcomes.

## II Fitting functions and motivations

The main purpose of this paper is to present more details about the velocity distributions of atomic samples cooled and trapped in optical lattices, where the Sisyphus cooling theory is expected to apply. A further step is to provide a function that gives a good approximation of the velocity distribution. The choice of a fitting function is made difficult by the complex dynamics of the atoms in the lattice. Indeed, even if the seminal process described in Dalibard89  gives very good insights in the dynamical behavior of the atoms, it is not sufficient in regimes relevant for typical experimental situations, where the intercombination of hamiltonian motion in the modulated potentials and optical pumping cycles, with time scales of the same order, makes it difficult to perform analytical calculations Castin . Along the following lines we justify a priori the choice of three types of functions (Gaussian, Tsallis and double Gaussian) that we used to fit the experimental and the numerical recorded data. As we will see, these choices are based on simple considerations about well-known generalizations of the model presented in Dalibard89 .

#### Gaussian function

In the standard description of 1D-Sisyphus cooling, the internal atomic state is adiabatically eliminated in such a way that the atomic dynamics is described in simple terms of a force $`F\left(v\right)`$ and fluctuating forces of momentum diffusion coefficient $`D_v\left(v\right)`$. $`F\left(v\right)`$ accounts for the optical pumping-assisted Sisyphus cycles and $`D_v\left(v\right)`$ corresponds, on the one hand, to the random recoils due to absorption and emission of photons, and on the other hand, to changes of potential curves. The velocity distribution, $`W(v)`$, is thus governed by a Fokker-Planck equation (FPE) Risken89 ; VanKampen92 :

$$_tW=_v\left(\frac{1}{M}F\left(v\right)W\right)+_v\left(D_v\left(v\right)_vW\right);$$
(1)
with $`M`$ being the atomic mass. In the linear regime for the atomic velocity, one finds cohen90 :

$`F\left(v\right)=\alpha v`$ (2)
$`D_v=D_v^{(1)}+D_v^{(2)}.`$
In this context, $`\alpha `$ and $`D_v`$ depend on the lattice parameters and are independent of the velocity. $`D_v^{(1)}`$ corresponds to the random absorption and emission of photons while $`D_v^{(2)}`$ represents the fluctuations of the light-shift induced force cohen90 . The steady-state solution of Eq. (1) with vanishing probability current ($`F\left(v\right)W+MD_v\left(v\right)_vW=0`$) is thus a Gaussian function with rms width $`\sigma _v=\sqrt{MD_v/\alpha }`$:

$$W\left(v\right)=W_0\mathrm{exp}\left(\frac{\alpha v^2}{2MD_v}\right).$$
(3)
#### Tsallis function

Beyond the linear regime for atomic velocity, the friction force and the velocity diffusion coefficients have to be refined into Castin ; Hodapp95 :

$`F\left(v\right)={\displaystyle \frac{\alpha v}{1+\left(v/v_\text{c}\right)^2}}`$ (4)
$`D_v\left(v\right)=D_v^{(1)}+{\displaystyle \frac{D_v^{(2)}}{1+\left(v/v_\text{c}\right)^2}},`$
where $`v_\text{c}`$ is the capture velocity which corresponds to the typical atomic velocity above which the Sisyphus process breaks down. Now, it is straightforward to show that the steady-state solution with vanishing probability current of Eq. (1) reads Lutz

$$W(v)=W_0\left[1\beta \left(1q\right)v^2\right]^{\frac{1}{1q}}$$
(5)
$`q=1+{\displaystyle \frac{2MD_v^{(1)}}{\alpha v_\text{c}^2}}`$ and $`\beta ={\displaystyle \frac{\alpha /2M}{D_v^{(1)}+D_v^{(2)}}}.`$ (6)
The function in Eq. (5) is the so-called Tsallis function and is in fact very general. It particularly provides a broad class of fitting functions including Gaussian functions ($`q`$ approaching one), Lorentzian functions ($`q=2`$) and inverted parabolas ($`q=0`$). At this stage, it is interesting to note that the Tsallis function has been introduced in the context of non-extensive thermodynamics Tsallis1 ; Tsallis1bis . The large amount of literature in this context allows one to find many papers dealing with problems already addressed in laser cooling; in particular anomalous diffusion in the presence of external forces Quarati ; Tsallis2 ; Plastino95 , multiplicative noise problems, and the relation to the edge of chaos in mixed phase space dynamics Lyra98 ; Weinstein02 . It is known that Sisyphus cooling can give rise to anomalous diffusion Katori97 ; Zoller96 , in particular for shallow optical potentials, where an atom can travel over many wavelengths before being trapped again. Even though we do not have a detailed analysis of the dynamics of the atoms in an optical lattice, for parameters corresponding to our situation, one cannot rule out anomalous diffusion and/or chaotic behavior.

#### Double Gaussian function

As Sisyphus cooling results in a situation where the kinetic energies of the atoms are of the order of the light shift potential, one can neither neglect atoms with lower energy (‘trapped’ in the potential wells) nor those moving more or less freely above the potential modulation (as in a ‘conduction band’). This leads to a description of the atomic sample in terms of a bimodal dynamics. Note that such a bimodal description has been shown to be relevant for the prediction of the diffusive properties of atoms in an optical lattice palencia . The kinetic equation of the ‘high energy’ atoms might very well be described by spatially averaged friction and diffusion coefficients resulting in a Gaussian distribution as shown previously. The ‘low energy’ atoms will be trapped, and subject to a different kinetic equation, and we assume that their velocity distribution is again a Gaussian. Our trial function is thus the sum of two Gaussian distributions with different widths (double Gaussian). One corresponding to ‘trapped’ atoms and the other one to ‘high energy’ atoms.

## III Experiments

### III.1 Experimental setup

The experimental setup has been described in detail previously (see e.g. Jersblad ; Ellmann ). Briefly, we first accumulate 133Cs atoms in a magneto-optic trap (MOT). We adjust the irradiance and the detuning, then we turn off the magnetic field and leave the atoms in an optical molasses with even further reduced irradiance. Thus we cool the atoms to a temperature of 3-4 $`\mu `$K. The atoms are transfered to a three-dimensional optical lattice, which is based on four laser beams of equal irradiance and detuning (for a review of optical lattice set-ups, see e.g. mennerat02  or jessen ). The detuning is a few tens of $`\mathrm{\Gamma }`$ below the ($`F_\text{g}=4F_\text{e}`$ = 5) resonance for the 133Cs D2 line at 852 nm ($`\mathrm{\Gamma }=2\pi 5.21`$ MHz is the linewidth of the excited state). The detuning ($`\mathrm{\Delta }`$) and irradiance ($`I`$) of the beams can be easily changed in order to control the depth of the light shift potential $`U_0I/|\mathrm{\Delta }|`$. The beams are aligned as in Fig. 1: two laser beams are linearly polarized along the $`x`$-axis and propagate in the $`yz`$-plane symmetrically with respect to the $`z`$-axis, whereas the other two beams are polarized along the $`y`$-axis and propagate in the $`xz`$-plane symmetrically with respect to $`z`$. This yields a tetragonal pattern of points with pure circular polarization, alternately $`\sigma ^+`$ and $`\sigma ^{}`$. These points correspond to potential wells where the atoms are trapped and optically pumped into the extreme $`m_F`$-levels (+4 and -4 respectively in $`\sigma ^+`$\- and $`\sigma ^{}`$-wells).

For high atomic velocities, this configuration will correspond to a three-dimensional version of the Sisyphus cooling model. As the atoms approach equilibrium, their kinetic energies will get lower than the modulation depth of the optical potential, and thus atoms become trapped in lattice sites. They will get distributed in bound states, where the lowest states closely resemble harmonic oscillator states.

In two different sets of runs, we let the atoms equilibrate in the optical lattice for 25 ms and 50 ms respectively. The velocity distribution is then recorded with a standard time-of-flight method (TOF) Phillips89 . After the lattice period the trapping field is turned off, and the atoms are released in the gravitational field; approximately 5 cm below the trap region a thin sheet of resonant laser light crosses the vertical axis along which the atoms fall, and the induced fluorescence is recorded with a photo-diode. Each vertical velocity component at the time of release will correspond to a specific arrival time at the probe beam. The probe beam is carefully spatially filtered and focussed by a cylindrical lens. The interaction region is less than 50 $`\mu `$m thick, and the trapped cloud of atoms is approximately 400 $`\mu `$m in diameter. This gives a velocity resolution of 0.05 mm/s, or 0.015 $`v_\text{R}`$ (where $`v_\text{R}`$ = 3.5 mm/s is the velocity corresponding to the recoil from one absorbed photon resonant with the D2-line). Our statistics is good enough not to contribute to this resolution. The optical lattice beams are turned off, by switching an acousto-optic modulator, faster than a microsecond. This is fast enough to avoid adiabatic release of the atoms in the lattice, which could greatly influence the velocity distribution, in particular in the high velocity tails 111 If the optical lattice beams are turned off too slowly, the atoms may partially equilibrate in the gradually decreasing potential. There may also be adiabatic cooling Kastberg . In both these cases, the cooling during a slow turn off can greatly influence the velocity distribution, in particular in the high velocity tails. Such adiabatic switching is often used in order to achieve lower ‘temperatures’..

### III.2 Experimental results

We recorded the velocity distributions for several modulation depths and we fitted them with the functions introduced in section II with a slight modification that accounts for atomic losses. During the long optical lattice phase, we have a constant loss of atoms, probably due to spatial diffusion. Therefore, the baseline is higher for atoms with a downward velocity (short times, $`v<0`$) than it is for atoms with a upward one ($`v>0`$). We compensate for this by adding a sharp step function to the fit, with the amplitude of the step as a free parameter. The amplitude of this step function is found to increase sharply for decreasing potential depths between $`U_\text{0}=200E_\text{R}`$ and $`100E_\text{R}`$. A probable reason is that spatial diffusion increases rapidly when the potential depth falls below some threshold, which takes place for higher potential depths than the threshold for cooling (usually called ‘décrochage’) Hodapp95 . This is consistent with previous studies carminati . In principle, we could have used a linearly decreasing function instead of the step function, but then this would have had be terminated by a sharp step. We avoied this in order to minimize the number of free parameters and also because we wanted to simplify as much as possible in the absence of detailed knowledge of the loss of atoms.

In Fig. 2, we show the rms width of the distributions, $`\sigma _v`$, as a function of the depth of the optical potential $`U_0`$, as derived from the fits to single Gaussian functions. The width, which is normally associated with a kinetic temperature, increases for deeper potential depths as usual.

In Figs. 3 and 4, typical recorded velocity distributions, together with Gaussian fits, are shown for low and high modulation depths. Figure 3 shows data taken with an equilibration time of 25 ms, and for Fig. 4 the equilibration time was 50 ms. This corresponds to typically $`10^6`$ radiative lifetimes. The plots with low irradiance are averages of twenty measurements and those of high irradiance of five measurements. For high values of the irradiance, a Gaussian function fits the velocity distribution extremely well. However, for low irradiance, it is clear that the wings of the distribution is not so well fitted. For the short equilibration time, this is more pronounced.

For all data, even below ‘décrochage’, the attempt with Lorentzian fits worked very poorly. Fits to double Gaussians and Tsallis functions, however, reproduced recorded distributions better than single Gaussians. In insets in figs. 3a and 4b we show fits to double Gaussians for shallow potentials. In Fig. 5, we compare the errors from the fits for these three types of functions. When the irradiance is varied, the signal-to-noise changes substantially, and so does the magnitude of the loss pedestal at short times, and the width and shape of the distribution. This makes it very hard to achieve a consistent normalization of the quality of the fits. The value of $`\chi ^2`$ ($`\chi ^2=\mathrm{\Sigma }(y_\text{i}x_\text{i})^2`$, where $`y_\text{i}`$ is the measured and $`x_\text{i}`$ the fitted value) for an individual fit includes information about both noise and systematic deviation from the fit function, which are difficult to separate. The data displayed in Fig. 5 are ratios between unnormalized values of $`\chi ^2`$ for the different fit functions. The displayed data are for the equilibration time of 25 ms. The other data set has the same features. For deep potentials, all fits are essentially equally good. At more shallow potentials, a Tsallis function reproduces the data better than a Gaussian. For the whole range, a double Gaussian gives the best fit. For the most shallow potentials, the fitted step becomes too important for $`\chi ^2`$ in order to draw any major conclusion from this analysis.

The parameter $`q`$ in Eq. (5) can be regarded as a measure of the shape of the distribution. A $`q`$ approaching 1 will be identical to a Gaussian distribution, whereas $`q=2`$ corresponds to a Lorentzian distribution. In Fig.6, we show a plot of the fitted value $`q`$, for 25 ms equilibration time. For decreasing irradiances, $`q`$ increases smoothly from one, and eventually reaches a value higher than $`q=1.6`$. For the longer equilibration time, the same trend is evident, but it is much less pronounced, and $`q`$ dose not reach higher than $`q=1.3`$.

The good fit to a double Gaussian can be interpreted as a sign of a bimodal velocity distribution. In Fig. 7a, we show the fitted widths of the two Gaussians for both data sets. This should correspond to the ‘temperatures’ of the two modes. Both these ‘temperatures’ increase linearly with potential depths. The areas of the two Gaussians should be a measure of the fraction of atoms being in one or the other of the modes. In Fig. 7b is the calculated relative populations.

The ‘cold mode’ with narrow velocity distribution always contains most of the atoms, but the relative number of atoms in the ‘hot mode’ gets larger for decreasing potential depths. For potentials deeper than $`U_\text{0}=250E_\text{R}`$ there is no measurable portion of atoms in the ‘hot modes’. The thermal energy of the ‘hot mode’ is of the same order (whithin the large uncertainties) as the energy barrier of the optical potential, i.e. the modulation depth (shown in the dashed line in Fig. 7a).

## IV Numerical simulations

In order to further analyze the results of our experimental data, we performed numerical simulations for the quantum dynamics of atoms undergoing Sisyphus cooling. We consider the case of a $`J=4J=5`$ transition, as for the 133Cs atoms used in the experiments, and for the sake of simplicity we restrict the motion of the atoms into one dimension (1D). The laser configuration is the well-known 1D-lin$``$lin configuration Dalibard89  which in fact corresponds to the $`z`$-direction in our three dimensional (3D) experimental setup (Fig. 1) with a different lattice spacing. This restriction is legitimate because, first, the temperature has been shown to be independent of the lattice spacing palencia ; carminati  and, second, the temperature is very similar for both 1D- and 3D- configurations (See the comparison between 3D-experimental and 1D-numerical results in Fig. 2. The slight deviation can partly be attributed to the difference in dimensionality). We first describe the numerical method for the integration of the dynamics equations (section IV.1) and then we present the results of the simulations (section IV.2).

### IV.1 Integration of the dynamics equations

Consider a two level atom, with Zeeman degeneracy, interacting with a laser field

$$\stackrel{}{E}_\text{L}(z,t)=\frac{1}{2}\left\{\stackrel{}{E}_\text{L}^+\left(z\right)\text{e}^{\text{i}\omega t}+\stackrel{}{E}_\text{L}^{}\left(z\right)\text{e}^{+\text{i}\omega t}\right\}.$$
(7)
The laser light is strong enough to be treated as a classical field so that we can separate the coupling between the atom and the electromagnetic field into a coupling to the laser field and a coupling to vacuum. The coupling to the laser $`\widehat{V}_{\text{AL}}`$ induces a hamiltonian evolution for the atom. On the contrary, because of the coupling to the vacuum modes, $`\widehat{V}_{\text{AV}}`$, the atom is an open quantum system for which the evolution has to be treated in the density matrix formalism. The evolution of the atom is thus governed by the generalized optical Bloch equations (OBE) CCT-JDR-GG88 ; cohen75 . In the regime of low saturation, where the experiments are performed, the excited state relaxes much faster than the typical evolution time of the ground state and thus it can be adiabatically eliminated from the OBE. The evolution of the system then reduces to a master equation for $`\varrho `$, the atomic density matrix restricted to the ground state cohen90 :

$`\text{i}\mathrm{}{\displaystyle \frac{\text{d}\widehat{\varrho }}{\text{d}t}}=\widehat{H}_{\text{eff}}\widehat{\varrho }\widehat{\varrho }\widehat{H}_{\text{eff}}^{}+_{\text{relax}}\left(\widehat{\varrho }\right)`$ (8)
$`\text{where }\text{aa}\widehat{H}_{\text{eff}}={\displaystyle \frac{\widehat{p}^2}{2M}}+{\displaystyle \frac{\widehat{V}_{\text{AL}}^{}\widehat{V}_{\text{AL}}^+}{\mathrm{}\left(\mathrm{\Delta }+\text{i}\mathrm{\Gamma }/2\right)}}.`$ (9)
Here, $`\widehat{p}`$ is the momentum operator, $`\mathrm{\Delta }`$ is the detuning between the laser and the atomic transition, $`M`$ is the mass of one atom, and $`\widehat{V}_{\text{AL}}^\pm =\widehat{\stackrel{}{D}}^\pm \stackrel{}{E}_\text{L}^\pm `$ are the raising and lowering parts of the dipole interaction operator respectively. In Eq. (8), $`\widehat{H}_{\text{eff}}`$ is a non-hermitian operator describing the atom-laser interaction222 The non-hermitian part of $`\widehat{H}_{\text{eff}}`$ comes from the relaxation of the excited state. and $`_{\text{relax}}`$ is an operator describing the coupling to the vacuum field, i.e. spontaneous emission of photons. The integration of the master equation is performed via a full quantum Monte-Carlo wave function method dalibard92 ; dum92  in which $`\varrho `$ is substituted with a set of stochastic wave functions. The pseudo-hamiltonian evolution (first term in Eq. (8)) of each wave function $`|\psi `$ is governed by a Schrödinger-like equation involving the non-hermitian hamiltonian $`\widehat{H}_{\text{eff}}`$:

$$\text{i}\mathrm{}\frac{\text{d}|\psi }{\text{d}t}=\widehat{H}_{\text{eff}}|\psi .$$
(10)
Since equation (10) does not include the filling terms of the ground state from the excited state due to spontaneous emission, $`|\psi `$ is not normalized and the instantaneous spontaneous emission rate is given by: $`\frac{d\psi |\psi /dt}{\psi |\psi }`$. To take into account the emission of photons, the pseudo-hamiltonian evolution (Eq. (10)) is interrupted by quantum jumps, whose repetition rate is determined with accordance to the spontaneous emission rate. It follows from the emission of a photon of wave vector $`\stackrel{}{\kappa }`$ and polarization $`\stackrel{}{ϵ}`$ that the wave function is instantaneously changed into

$$|\psi |\psi _{\stackrel{}{\kappa },\stackrel{}{ϵ}}=1_{\stackrel{}{\kappa },\stackrel{}{ϵ}}|\widehat{V}_{\text{AV}}\left(|\psi _\text{e}|0\right)$$
(11)
with relative probabilities $`\left||\psi _{\stackrel{}{\kappa },\stackrel{}{ϵ}}\right|^2`$. Here the excited state wave function $`|\psi _\text{e}=\frac{\widehat{V}_{\text{AL}}^+|\psi }{\mathrm{}\left(\mathrm{\Delta }+\text{i}\mathrm{\Gamma }/2\right)}`$ is determined by the adiabatic elimination procedure of the excited state and $`|0`$ and $`|1_{\stackrel{}{\kappa },\stackrel{}{ϵ}}`$ represent the electromagnetic field states respectively without any photon, and with one photon of wave vector $`\stackrel{}{\kappa }`$ and polarization $`\stackrel{}{ϵ}`$. The Monte-Carlo integration then provides a set of time dependent stochastic wave functions $`|\psi `$, which represent the atomic state through the average, $`\overline{\sigma }`$, of the density matrices associated to the wave functions, $`\sigma =|\psi \psi |`$. It is then straightforward to show that the quantum master equation for $`\overline{\sigma }`$ is the same as the master equation for the actual density matrix $`\varrho `$ (Eq. 8). Hence, the value of any observable $`\widehat{𝒪}`$ for the quantum system represented by $`\varrho `$ is equal to the ensemble average of the value of the same observable for each stochastic wave function represented by $`|\psi `$ moelmer93 : at any time,

$$\overline{\psi \left|\widehat{𝒪}\right|\psi }=\text{Tr}\left(\widehat{\varrho }\widehat{𝒪}\right).$$
(12)
### IV.2 Results of the simulations

In this work, we are interested in the particular observable that represents the momentum distribution of the atoms. We have performed the full quantum Monte-Carlo integration of the dynamics equations for a set of $`200`$ wave functions for various lattice parameters (detuning and modulation depth). Because the width of the momentum distributions are typically broader than several $`\mathrm{}k`$, the spontaneous emission pattern can be approximated by photons emitted along the 3D coordinate axes $`x`$, $`y`$ or $`z`$. With such an approximation, all operators in Eqs. (9) and (11) couple states of the form $`|m,p`$ to states of the form $`|m^{},p\pm \mathrm{}k`$ (where $`m`$ and $`m^{}`$ represent the internal sub-level of the atomic ground state). It is then convenient to perform the integration in the $`|p`$-representation. The state $`|\psi `$ is decomposed onto the basis of the $`|p`$ states (with $`p=n\mathrm{}k`$, with $`n`$ an integer positive or negative). Finally, for usual situations considered in this work, the typical momenta are smaller than 20 $`\mathrm{}k`$, so that we take $`|n|100`$. From the simulations, we determined the mean kinetic energy as a function of time. After a thermalization period, the energy reaches a steady-state during which the momentum distribution was recorded and averaged. The thermalization period was chosen to be $`1/\mathrm{\Gamma }`$, corresponding to a time in the order of a millisecond. Since the calculation is performed in 1D, this time cannot be directly compared to the thermalization times in the 3D experiment.

In order to identify whether the momentum distribution is compatible with a Gaussian curve or not, we first compare the root-mean-square (rms) momentum $`p_{\text{rms}}`$ defined by $`E_\text{K}=p_{\text{rms}}^2/2M`$ (where $`E_\text{K}`$ is the mean kinetic energy of the atomic sample) and $`p_\text{e}`$ which represents half the width at $`1/\sqrt{\text{e}}`$ of the stationary momentum distribution. For a Gaussian distribution, those two values are equal.

We plot in Fig. 8, the numerical results for $`p_{\text{rms}}`$ and $`p_\text{e}`$ as a function of the potential depth $`U_0`$ for three different detunings $`\mathrm{\Delta }`$.

We find that these values are independent of the detuning within the numerical errors. Several points for lower values of $`U_0`$ have also been calculated but the atomic cloud was found not to thermalize. For those cases, the temperature increases more or less linearly and the velocity distribution becomes almost flat. It is also clear in Fig. 8 that $`p_{\text{rms}}`$ and $`p_\text{e}`$ have different behaviors. $`p_{\text{rms}}`$ reproduces the well known dependence of the kinetic energy versus the modulation depth: $`p_{\text{rms}}`$ scales as $`\sqrt{U_0}`$ for high values of $`U_0`$ and abruptly increases as $`U_0`$ reaches very low values, typically lower than 150 $`E_\text{R}`$ (the point of décrochage). The minimum value of $`p_{\text{rms}}`$ is found to be of the order of $`\left(p_{\text{rms}}\right)_{\text{min}}4.1\mathrm{}k`$. On the contrary, we find that $`p_\text{e}`$ increases monotonically versus $`U_0`$ for low values as well as for high values of $`U_0`$. The minimum value of $`p_\text{e}`$ is obtained for the minimum value of $`U_0`$ for which a steady-state velocity distribution can be obtained ($`U_078E_\text{R}`$) and is found to be of the order of $`\left(p_\text{e}\right)_{\text{min}}3.4\mathrm{}k`$. We identify two different regimes that can be distinguished: For $`U_0`$ above décrochage ($`U_0150E_\text{R}`$), both $`p_\text{e}`$ and $`p_{\text{rms}}`$ increase and $`p_\text{e}`$ is slightly larger than $`p_{\text{rms}}`$ that is to say that the momentum distribution is wider than a Gaussian distribution with the same $`p_{\text{rms}}`$. For $`U_0`$ below décrochage ($`U_0150E_\text{R}`$), $`p_\text{e}`$ decreases while $`p_{\text{rms}}`$ increases rapidly as $`U_0`$ decreases; the momentum distribution has large wings and becomes narrower than a Gaussian distribution. These different characteristics are illustrated in Fig. 9 where we plot the simulated velocity distributions together with Gaussian fits in the two regimes $`U_0150E_\text{R}`$ and $`U_0150E_\text{R}`$.

One should note that this result is in disagreement with earlier calculations performed for atoms with a theoretical $`J=1/2J=3/2`$ transition for which Castin et al. find that $`p_{\text{rms}}>p_\text{e}`$ for any value of the potential depth $`U_0`$ Castin . In fact, when running our simulation for the $`J=1/2J=3/2`$ transition, we were able to reproduce the results of Castin  and we thus conclude that the discrepancy is due to the different atomic transitions considered in Castin  and in the present work. We finally conclude that in general, the momentum distribution significantly differs from a Gaussian distribution. Moreover, we find that the threshold for $`p_{\text{rms}}`$ at low values of $`U_0`$ do not affect $`p_\text{e}`$.

We now turn to a more detailed analysis of the momentum distributions. We first fit the velocity distributions to Tsallis functions. The dependence of the Tsallis parameter $`q`$ on the modulation depth is also shown in Fig. 10 and show a linear dependence of $`q`$ versus $`U_0`$.

For all numerical data $`q`$ differs from $`1`$ only by less than $`5\%`$ and is less than $`1`$. Moreover, $`q`$ is found to tend to $`1`$ for shallow potentials indicating that the best Tsallis fit is close to a Gaussian curve in disagreement with the previous discussion. The discrepancy between numerical simulations and experimental measurements may be caused by the different dimensionality considered in the experiments and in the simulations.

Consider now fits to double Gaussians. We plot in Fig. 11 numerically recorded velocity distributions in logarithmic scale for potential depths in the two regimes corresponding to shallow and deep potentials, together with fits to double Gaussians. For the deep potential ($`U_0=293E_\text{R}`$), the profile is essentially parabolic and thus corresponds to a Gaussian distribution. For the shallow potential ($`U_0=78E_\text{R}`$), we clearly identify two contributions: in addition to a narrow parabolic profile (corresponding to low energetic atoms), a broad one (corresponding to high energetic atoms) appears.

This supports the interpretation of the dynamics in terms of a bimodal atomic distribution, with each mode corresponding to ‘trapped’ atoms and to nearly ‘free’ atoms. The whole distribution is well fitted by a double Gaussian function. We plot in Fig. 12a the widths of both the modes as functions of $`U_0`$.

The numerical results are in good agreement with experimental ones (see Fig. 7). For shallow potentials, we find two Gaussian components with widths that both increase with the potential depth $`U_0`$, whereas for deep potentials, the ‘hot component’ is almost undetectable.Thus the route to ‘décrochage’ for shallower potentials can be interpreted as a transfer from the cold mode to the hot mode. This is supported by the results for the populations of the two Gaussian contributions to the velocity distribution plotted in Fig. 12b. We actually find that the cold mode is largely dominant even for very shallow potentials close to ‘décrochage’.

Finally, we compare the numerical and experimental results. A direct quantitative comparison is not adequate, since the simulations are done in 1D. However, qualitatively, the experimental data are reproduced excellently. Figures 9-12 show numerical data corresponding to the experimental ones in Figs. 3-7. The single Gaussian works for high irradiance but fails to fit the wings of the distribution for low irradiances. A Tsallis function does not fit the distribution any better than a single Gaussian for the numerical data. Again, the distribution is best fitted by a double Gaussian and this is particularly pronounced for shallow potentials. The fits to double Gaussians also reproduce the signature of one ‘hot’ and one ‘cold’ mode for shallow potentials. This strongly supports assumptions of a bimodal distribution.

We find that the ‘cold’ mode is largely dominant even for very shallow potentials close to ‘décrochage’in both the experimental ($`70\%`$) and the numerical ($`90\%`$) results (see Fig. 7b and 12b). The quantitative discrepancy between the limit population of the ‘hot’ and ‘cold’ modes results from the different dimensionnalities (3D experiments versus 1D calculations). Indeed, the fraction of bound atoms in a $`d`$D dimension situation can be estimated by

$$n_{\text{bound}}^{(d\text{D})}=\alpha \underset{j}{}p_{d\text{D}}(E_j)\underset{1\mu d}{}\mathrm{𝟏}_{E_{j\mu }<E_{\text{max}\mu }}$$
(13)
where $`p_{d\text{D}}(E_j)`$ represents the population of a state of energy $`E_j`$, $`\alpha `$ in a normalization factor and $`E_{\text{max}\mu }`$ stands for the maximum energy of bound states in the potential wells along the direction $`\mu `$ ($`E_{\text{max}\mu }`$ is of the order of the potential depth). Now, if the space directions are separated (which is the case in the harmonic approximation that one can assume to be valid for bound states), then

$$n_{\text{bound}}^{(d\text{D})}=\underset{1\mu d}{}\alpha _\mu \underset{j}{}p_{1\text{D}}(E_{j\mu })\mathrm{𝟏}_{E_{j\mu }<E_{\text{max}\mu }}.$$
(14)
Hence, assuming that all directions are equivalent, the fraction of atoms in non-bound states reads

$`n_{\text{non bound}}^{(d\text{D})}`$ $`=`$ $`1{\displaystyle \underset{1\mu d}{}}\left(1n_{\text{non bound}}^{(1\text{D}\mu )}\right)`$ (15)
$``$ $`dn_{\text{non bound}}^{(1\text{D})}.`$
Therefore, the limit populations of the ‘hot’ mode at ‘décrochage’ are consistent in the experiments ($`30\%`$ and $`d=3`$) and in the simulations ($`10\%`$ and $`d=1`$).

## V Conclusions

In this work, we have studied the velocity distributions of cold atomic samples obtained by Sisyphus cooling both in experiments with 133Cs and in full quantum numerical simulations performed for the actual $`45`$ transition of 133Cs. We stressed in particular the deviation from a Gaussian distribution. This has already been forecasted via semi-classical as well as quantum simulations for a simplified $`1/23/2`$ transition showing the difference of the rms velocity $`v_{\text{rms}}`$ and the velocity $`v_\text{e}`$ corresponding to half the width at $`1/\sqrt{\text{e}}`$ of the distribution Castin . We recovered such a property but with a significantly different behavior of the ratio $`v_{\text{rms}}/v_\text{e}`$. This shows that the non-Gaussian behavior of the velocity distributions is certainly not a trivial effect in Sisyphus cooling.

### V.1 Summary of our results

Our results (experiments as well as numerical simulations) show that the velocity distributions are compatible with Gaussian functions for deep enough potentials (typically for $`U_0`$ larger than a hundred recoil energies). Note that in this case, the atoms are trapped in the potential wells (i.e. the kinetic energy of the atomic cloud is significantly smaller than the potential depth). The deviation of the velocity distribution from Gaussian functions become more prominent for shallow light shift potentials. We tested several types of functions to better fit the shape of the velocity distributions in the range of parameters corresponding to deeper potentials. We found that a better fit (corresponding to smaller $`\chi ^2`$) can be obtained by using a Tsallis function or a double Gaussian.

#### Tsallis functions -

The use of a Tsallis function is related to the details of the dynamics of atoms cooled by the Sisyphus mechanism which is known to be slightly more complicate than a Brownian motion. The Tsallis function introduces a new parameter $`q`$ which deviation from $`1`$ measures the deviation of the velocity distribution from a Gaussian function. The parameter $`q`$ can be calculated in the ‘jumping regime’ mennerat98  and it is straightforward to show that $`q`$ tends to $`1`$ for high values of the potential depth (thus corresponding to a weak deviation from a Gaussian), and increase for shallow potentials. An ab initio calculation of $`q`$ is more tricky in the ‘oscillating regime’ which correspond to the domain of parameters for shallow potentials, near the point of ‘décrochage’ mennerat98 . Nevertheless we can plot the value of $`q`$ corresponding to the best fit of the measured velocity distributions as a function of modulation depth. For large modulation depths, we find that $`q`$ approaches $`1`$, which corresponds to a Gaussian distribution, in agreement with the analytical calculation (see section II and Lutz ). When reducing the potential depth, we clearly observed an increase in $`q`$ and this corresponds to a velocity distrubution with wings larger than in a Gaussian function. In our case the maximum $`q`$ is close to 1.6 and this corresponds in our experiments to a potential depth $`U_060E_\text{R}`$. For $`U_0<60E_\text{R}`$, the atomic cloud does not reach a steady state and the optical lattice disintegrates. It is interesting to note that the rms velocity of Tsallis distributions with $`q`$ above $`q_{\text{cr}}`$ = 5/3 diverge Tsallis99 . If one would plot rms ‘temperatures’ of the atoms using the rms velocity, this would correspond to a diverging temperature. As one is often limited by noise in the wings of the velocity distribution, one has a tendency to restrict the analysis to atoms with velocities several times below the 1/e value of the distribution. Any divergence is hence avoided. Note also that such divergences are very familiar: the wings of a Lorentz distribution are also known to cause a divergence of the rms value of the distribution. One can also recall that in the case of narrow line cooling, the rms velocity diverges Wallis89 ; Katori99  when one approaches the atomic resonance, and that for very small detunings one can no longer even have a normalized distribution function Wallis89 .

#### Double Gaussian functions -

Fitting the recorded (experimental or numerical) distribution functions to double Gaussians works even better than the Tsallis function. On the one hand, it is not surprising that a fitting procedure with more free parameters gives better fits. On the other hand, the velocity distribution in logarithmic scale in Fig. 11a clearly exhibits two components with very different widths. For deep potentials, one recovers a Gaussian distribution of ‘cold atoms’ bound in in the potential wells as expected from the above discussion. When decreasing the potential depth, a fraction of ‘hot atoms’ grows up (for $`U_0<120E_\text{R}`$). These atoms have an energy larger than the potential depth and are not trapped in the potential wells. We found that the fraction of ‘hot atoms’ can be significant for very shallow potentials. It reaches $`30\%`$ in 3D experiments and $`10\%`$ in 1D numerical simulations just above ‘décrochage’ at $`U_060E_\text{R}`$ (the discrepancy between the experiments and the simulations is due to the different dimensions as shown at the end of section IV.2). This result strongly supports assumptions that an optical lattice has a bimodal velocity distribution. A straightforward interpretation would be that some atoms are bound at lattice sites, whereas others have enough energy to move around on top of the modulated potential. An interesting results of our work is that the phenomenon of ‘décrochage’ does not correspond to a sharp increase of the width of the velocity distributions corresponding to each mode but to a continuous transfer from the ‘cold mode’ to the ‘hot mode’. We found that when décrochage occurs, the fraction of atoms in the ‘hot mode’ does not exceed a few tens percent.

### V.2 Perspectives

The results shown in this paper stronlgy suggest that the simple picture for Sisyphus cooling, based on a competition between a diffusion and a friction (see Eq. 4), is not adequate to describe the ‘coldest’ velocity distributions. Even though one has to be careful before generalizing the conclusions of this paper to other situations of laser cooling and/or trapping, the existence of two velocity modes might provide a useful guide to understand the dynamics and limits of laser cooling. One can note e.g. that for shallow potentials, one has fewer bound states, and the fraction of atoms in the conduction band gets more prominent, as shown in Figs. 7 and 12.

These atoms will experience a friction force corresponding to the classical Sisyphus cooling model. The route to equilibrium for the bound atoms is less clear. One hypothesis castpriv  is that bound levels are uniformly ‘watered’ from the conduction band, whereas high lying levels are more likely to escape. Thus, the route to equilibrium is not quite a competition between cooling and heating. A drawback with this theory is that it would not yield Gaussian velocity distributions. However, this theory has the advantage that the rate of equilibration should depend linearly on irradiance, which is consistent with previous experiments raithel ; palencia2 . In contrast, the standard Sisyphus cooling theory predicts a cooling rate independent of irradiance Dalibard89 . An interesting experiment would be to measure the velocity distribution as a function of time after a sudden change of the light shift potential, and see if the two populations would evolve differently.

Figures 3 and 4 seem to indicate a time dependence of the experimentally recorded velocity distribution. However, with the current data set (using only the two cooling times 25 ms and 50 ms), and with the current experimental uncertainties, we cannot draw any quantitative conclusion for the time dependence of the velocity distributrion. In future work, we will study the velocity distribution as a function of time.

It would also be interesting to extend the test functions used in this paper to a narrow-line cooling scheme, which become more and more used with the laser cooling of earth-alkaline atoms. At this stage, one can however note, that a non normalized distribution function will have as an effect that there is no steady state distribution and that in this case atoms will diffuse to large velocities. This will appear in an experiment as a leakage rate of the atoms from the optical lattice. The background observed in our experiment become more and more dominant for shallow potential wells. One might expect this to have a contribution from a diffusion of the atoms beyond the capture range of the optical lattice corresponding in practice to a non-normalized distribution function. A detailed analysis of the velocity distribution of atoms in optical lattices thus appears as a promising tool to study new statistical effects.

Experiments as well as full quantum simulations (in 1D and 3D) should allow one to get new insights in the dynamics of such systems. Apart from the suggestions above, future work could e.g. focus on the phase space dynamics of atoms in optical lattices and of quantum transport properties of ultra-cold atoms or even Bose-Einstein condensates.

###### Acknowledgements.

LSP thanks the swedish group for warm hospitality during the period when a part of this work was achieved. He also acknowledges financial support from the Swedish Foundation for International Cooperation in Research and Higher Education (STINT). RK thanks the Wenner-Gren foundation for a travel grant. We would like to thank Dr. Peter Olsson at Umeå University for letting us use the LINUX cluster and also for support during the simulations at the theoretical physics department at Umeå University. We also thank Eric Lutz and Philippe Verkerk for fruitful discussions. This work has been supported by the Carl Trygger foundation and the Knut & Alice Wallenberg foundation.


format: text

text:
# Rapid convergence of time-averaged frequency in phase synchronized systems

## I Introduction

The rich collective behavior, including mutual entrainment and self-synchronization, in systems of coupled oscillators has been considered by several investigators in the past few years (see, for example, Refs. Pikovsky et al. (2001); Strogatz (2000); Kuramoto (1984) and references therein). Recently, a considerable amount of research has been devoted to the study of coupled *chaotic* oscillators and, in particular, to the phenomenon of phase synchronization. Provided that the phase can be defined Pikovsky et al. (1997); Josić and Mar (2001), two coupled nonidentical chaotic oscillators are said to be phase synchronized if their frequencies are locked but their amplitudes are not Rosenblum et al. (1996); Pikovsky et al. (2001). This appears to be a general phenomenon and it has been observed in systems as diverse as electrically coupled neurons Elson et al. (1998); Makarenko and Llinás (1998), biomedical systems Schäfer et al. (1998), chemical systems Kiss et al. (2002), and spatially extended ecological systems Blasius et al. (1999). Moreover, the potential role of phase synchronization in brain functions has been explored Tass et al. (2003); Fitzgerald (1999).

The most common theoretical approach to phase synchronization of chaotic oscillators is based on an analogy with the evolution of the phase of a periodic oscillator in the presence of external noise Pikovsky et al. (2001). This approach leads to the conclusion that the dynamics of the phase is generally diffusive and the phase performs a random walk. However, the effective ”noise” in such a description cannot be considered as a Gaussian $`\delta `$-correlated noise in all circumstances. In this paper we show that the effective ”noise” exhibits strong temporal (anti-)correlations for a general class of chaotic attractors. In particular, we present evidence from simulations and experiments that many phase synchronized systems of chaotic oscillators where the chaotic state is reached through a period-doubling cascade show a rapid convergence of the time-averaged frequency. The speed of convergence toward the natural frequency scales as the inverse of the measurement period. This implies that short measurement times may suffice for reliable determination of frequencies in those systems.

The outline of the paper is as follows: Section II presents a sketch of the theoretical background of phase dynamics in chaotic systems and the expected scaling of the time-averaged frequency with the observation time. The results of numerical simulations of globally coupled arrays of Rössler oscillators are presented in Sec. III where features of the chaotic attractor leading to rapid convergence of the time-averaged frequency are identified. Section IV contains a description and analysis of experimental results on a globally coupled array of electrochemical oscillators where rapid frequency convergence is observed.

## II Theoretical Background

An appropriate definition of the phase for chaotic self-sustained oscillators can be obtained from the Poincaré map of the flow. Such a map can be constructed if a surface of section in the phase space of the autonomous continuous-time chaotic system exists which is crossed transversally by all trajectories of the chaotic attractor. Then, for each piece of a trajectory between two crossings of this surface, we define the phase as a linear function of time:

$$\varphi (t)=2\pi \frac{tt_n}{t_{n+1}t_n}+2\pi n,$$
(1)
for $`t_nt<t_{n+1}`$. Here, $`t_n`$ is the time of the $`n`$th crossing of the surface of section. The definition is ambiguous because it depends on the choice of the Poincaré surface. Yet, any choice of a phase variable for chaotic oscillators investigated in this paper leads to the same macroscopic behavior Josić and Mar (2001).

With this phase definition, the phase dynamics can be described by

$`A_{n+1}`$ $`=`$ $`(A_n),`$ (2)
$`d\varphi /dt`$ $`=`$ $`\omega (A_n)\omega _0+\eta (A_n),`$ (3)
where the amplitude $`A_n`$ is the set of coordinates of the phase point on the Poincaré surface at the nth intersection and $``$ defines the Poincaré map that takes $`A_n`$ to $`A_{n+1}`$. The “instantaneous” frequency $`\omega (A_n)=2\pi /T_n`$ is determined by the Poincaré return time $`T_n=t_{n+1}t_n`$ and depends in general on the amplitude. Assuming chaotic behavior of the amplitudes, we can consider the term $`\omega (A_n)`$ to be the sum of the average (natural) frequency $`\omega _0`$ and some effective “noise” $`\eta (A_n)`$ with zero mean, although this stochastic term has a purely deterministic origin Pikovsky et al. (2001). Thus, Eq.(3) has the solution

$$\varphi (t)=\varphi _0+\omega _0t+_0^t\eta (\tau )𝑑\tau .$$
(4)
The variance of the integral $`_0^t\eta (\tau )𝑑\tau `$ is given by

$`\left({\displaystyle _0^t}\eta (\tau )𝑑\tau \right)^2`$ $`=`$ $`2{\displaystyle _0^t}𝑑\tau ^{}(t\tau ^{})K(\tau ^{})`$ (5)
$``$ $`2t{\displaystyle _0^{\mathrm{}}}𝑑\tau ^{}K(\tau ^{})2tD_\eta .`$
Here, $`K(\tau ^{})=\eta (\tau )\eta (\tau +\tau ^{})`$ and the average $`\mathrm{}`$ is taken over the invariant measure of the chaotic attractor. The approximate equality holds provided that $`\eta (\tau )`$ has a finite correlation time $`t_c`$ such that $`K(\tau ^{})0`$ for $`\tau ^{}>t_c`$ and one considers times $`tt_c`$. This shows that the dynamics of the phase is diffusive and the phase performs a random walk as long as $`t_c`$ is finite and $`D_\eta 0`$. However, as pointed out by Pikovsky and co-workers Pikovsky et al. (2001), the dynamics of the phase generally differs from stochastic Brownian motion because the effective “noise” cannot be considered as a Gaussian white noise process. This was recently confirmed for a system of locally coupled Rössler oscillators where temporal correlations were shown to exist Davidsen and Kapral (2002). These correlations determine the speed of convergence of the time average $`\overline{\omega }(T)=T^1(\varphi (T)\varphi (0))`$ toward $`\omega _0`$. The speed of convergence — as measured by the standard deviation of the ensemble distribution of $`\overline{\omega }(T)`$ — scales as $`1/T`$. This is in contrast to what one would expect if $`\eta `$ were a white noise process where one obtains $`1/\sqrt{T}`$ scaling since $`\overline{\omega }(T)\omega _0=1/T_0^T𝑑t\eta (t)`$ and the standard deviation of $`_0^T𝑑t\eta (t)`$ scales with $`\sqrt{T}`$. Here, we provide further numerical and experimental evidence that such a scaling generally occurs for oscillators where the chaotic state is reached through a period-doubling cascade. In particular, the type of coupling between the chaotic oscillators seems to be irrelevant. Moreover, we show that the temporal correlations induce an alternating behavior of $`K(\tau ^{})`$. This leads to an extremely small value of $`D_\eta `$ for these systems and suggests that they can generally be phase synchronized.

## III Globally Coupled Rössler Oscillators

We consider a $`L\times L`$ array of *globally* coupled, nonidentical chaotic Rössler oscillators:

$$\frac{𝐱(𝐫,t)}{t}=𝐑(𝐱(𝐫,t))+K\underset{\widehat{𝐫}𝒩}{}\left(𝐱(\widehat{𝐫},t)𝐱(𝐫,t)\right),$$
(6)
where $`R_1=\omega (𝐫)x_2x_3`$, $`R_2=\omega (𝐫)x_1+0.2x_2`$, $`R_3=x_1x_35.9x_3+0.2`$. The sites of the lattice are labelled by $`𝐫`$, $`K=K^{}/L^2`$, $`K^{}`$ is the coupling constant and $`𝒩`$ is the set of all sites on the lattice. We take $`L=64`$ and choose the $`\omega (𝐫)`$’s randomly from a uniform distribution on the interval $`[0.99,1.01]`$. This ensures that all oscillators are chaotic. We find that the speed of convergence of $`\overline{\omega }(T)`$ toward $`\omega _0`$ in the phase synchronized state scales as $`1/T`$ (see Fig. 1). This is also true for higher $`\omega (𝐫)`$ dispersion.

The $`1/T`$ scaling observed in the globally coupled Rössler system follows from the fact that the local chaotic attractors, although one-banded in the phase synchronized state (see Fig. 2), have internal structure.

In the absence of coupling ($`K=0`$) the local non-identical Rössler oscillators display a variety of banded chaotic attractors because of period-doubling cascades and period-3 windows in their vicinity. In the phase synchronized state for $`K=0.06`$, the local chaotic attractors are very similar to one another and have an internal structure similar to that of a merged two-banded attractor. Although this structure is not obvious in the left panel of Fig. 2, it can be seen clearly in the local next amplitude maps (Fig. 3) and in the density distribution $`\rho (l)`$ (inset in Fig. 3). The intersection point of the map with the bisectrix defines four regions in the map such that the two-bandedness of the attractor can be characterized by the percentage $`p`$ of points in the upper left and lower right quarter-planes. Indeed, only very few points $`(1p=3\%)`$ — typically those for small $`A_{n,1}`$’s which are in the lower left quarter-plane of the amplitude map — stay in the same region of the bimodal density distribution, implying that the majority of the iterates map regions of high density on the left to regions of high density on the right and vice versa.

This nearly two-banded character is also seen in the series of $`T_n`$’s (see Fig. 4). An alternation of long and short return times can be observed with rare exceptions. This shows the close relation between topological characteristics (next amplitude map) and temporal evolution ($`T_n,T_{n+1}`$ map) and suggests the following explanation for the observed scaling:

Assuming that each internal ”band” of the attractor is associated with a distinct mean frequency $`\omega _1`$ or $`\omega _2`$ such that $`\omega _0=2\omega _1\omega _2/(\omega _1+\omega _2)`$, and the system alternates between bands, the ”random” force in the evolution equation (3) for the phase of a single oscillator takes the form

$$\eta (t)=(1)^{\kappa +1}\frac{\omega _\kappa (\omega _1\omega _2)}{(\omega _1+\omega _2)}+\stackrel{~}{\eta }(t),$$
(7)
where $`\kappa =1(2)`$ when the system is on band 1(2). The noise term $`\stackrel{~}{\eta }(t)`$ accounts for deviations of the period within a band and deviations from the strict alternations between bands. If $`T`$ is the time needed for $`n`$ oscillations, it follows from Eq.(4) that

$`\overline{\omega }(T)\omega _0=\{\begin{array}{cc}_0^T𝑑t\frac{\stackrel{~}{\eta }(t)}{T}\hfill & n\text{ even},\hfill \\ (1)^{\kappa +1}\frac{2\pi (\omega _1\omega _2)}{(\omega _1+\omega _2)T}+_0^T𝑑t\frac{\stackrel{~}{\eta }(t)}{T}\hfill & n\text{ odd},\hfill \end{array}`$ (8)
where $`\kappa =1(2)`$ corresponds to band 1(2) appearing first in the series. The term $`\frac{2\pi (\omega _i\omega _j)}{(\omega _1+\omega _2)T}`$ scales with $`1/T`$, while the term $`_0^T𝑑t\frac{\stackrel{~}{\eta }(t)}{T}`$ scales with $`1/\sqrt{T}`$ provided that $`\stackrel{~}{\eta }(t)`$ is $`\delta `$-correlated noise. Since the amplitude of the $`1/\sqrt{T}`$ term is very small, two scaling regimes can be identified for the maximal deviation from $`\omega _0`$: A $`1/T`$ scaling on intermediate time scales which is dominated by the switching between the ”bands” and a $`1/\sqrt{T}`$ scaling on long time scales. The intermediate time scale scaling behavior can indeed be observed for each oscillator: The envelope of $`\overline{\omega }(T)\omega _0`$ decays as $`1/T`$. Note that the standard deviation of the ensemble distribution of $`\overline{\omega }(T)`$ decays as $`1/\sqrt{T}`$ for identical chaotic Rössler oscillators and $`K=0`$ — as expected from the law of large numbers. In the case of nonidentical oscillators, phase synchronization implies that the phases are synchronized such that $`\sqrt{(\overline{\omega }(T)\omega _0)^2}`$, where the average is taken over the ensemble, adopts the scaling of a single oscillator reproducing the observed scaling and explaining the results. This result does not depend on the type of coupling of the oscillators as confirmed by the numerical results presented here local  and in Ref. Davidsen and Kapral (2002). This explanation is based exclusively on the fact that the series of the instantaneous frequencies alternates with the exception of rare events and that the system is phase synchronized.

The above argument easily generalizes to attractors with multiple internal “bands” which we observe in the phase synchronized state for higher coupling. Figure 5 shows that the deviations from a three-banded structure are rather small for $`K=0.2`$. Three bands A, B and C can be identified such that A is mapped to B and C to A. However, band B is mapped to C and Z, and Z back to B such that a pure three-banded structure does not exist. The invariant measure of the map in area Z is very low (approximately $`10\%`$).

The alternation of positive and negative values for $`\eta (t)`$ implies that $`K(\tau ^{})`$ alternates as well. Hence, $`D_\eta `$ is very small, i.e., a high degree of phase coherence persists. Our findings suggest that this should be generally expected for oscillators where the chaotic state is reached through a period-doubling cascade. Thus, these oscillators would be phase coherent suggesting that they can generically be phase synchronized.

The change in the structure of the local attractors as the coupling strength is changed is the analog of the bifurcation structure seen in homogeneous locally coupled Rössler oscillators where one finds period doubling different from the isolated Rössler oscillator Goryachev and Kapral (1999).

## IV Globally Coupled Electrochemical Oscillators

We now study the convergence of the time-averaged frequency in an array of globally coupled electrochemical oscillators. A standard three electrode electrochemical cell consisting of a nickel working electrode array (64 1-mm diameter electrodes in an $`8\times 8`$ geometry with 2 mm spacing), a $`\mathrm{Hg}/\mathrm{Hg}_2\mathrm{SO}_4/\mathrm{K}_2\mathrm{SO}_4`$ reference electrode and a Pt mesh counter electrode were used. The potentials of all of the electrodes in the array were held at the same value ($`V=1.310\mathrm{V}`$) with a potentiostat (EG&G Princeton Applied Research, model 273). Experiments were carried out in 4.5 M $`\mathrm{H}_2\mathrm{SO}_4`$ solution at a temperature of $`11^\mathrm{o}\mathrm{C}`$. The working electrodes were embedded in epoxy and reaction takes place only at the ends. The currents of the electrodes were measured independently at a sampling rate of 100 Hz and thus the rate of reaction as a function of position and time was obtained.

The electrodes were connected to the potentiostat through 64 individual parallel resistors ($`R_\mathrm{p}`$) and through one series collective resistor ($`R_\mathrm{s}`$). We employed a method of altering the strength of global coupling while holding all other parameters constant. The total external resistance was held constant while the fraction dedicated to individual currents, as opposed to the total current, was varied. A total resistance can be defined as

$$R_{\mathrm{tot}}=R_\mathrm{s}+R_\mathrm{p}/64.$$
(9)
In these experiments $`R_{\mathrm{tot}}=14.2\mathrm{\Omega }`$. The series resistor couples the electrodes globally. The parameter $`\epsilon `$, the ratio of series to total resistance, is a measure of the global coupling

$$\epsilon =\frac{R_\mathrm{s}}{R_{\mathrm{tot}}}.$$
(10)
For $`\epsilon =0`$, the external resistance furnishes no additional global coupling; for $`\epsilon =1`$, maximal external global coupling is achieved.

The population of chaotic oscillators is characterized by certain amount of heterogeneity Kiss et al. (2002). Without added coupling there is a distribution of frequencies of the oscillators with a mean of $`1.219\mathrm{Hz}`$ and a standard deviation of $`18\mathrm{mHz}`$. With weak added global coupling ($`\epsilon =0.1`$) a nearly phase synchronized state occurs in which 63 of 64 oscillators have a frequency of $`1.230\mathrm{Hz}`$, and that of the remaining element has a frequency of $`1.237\mathrm{Hz}`$. The results are shown for this region of high phase synchrony.

### Experimental Results

The chaotic attractor of a representative single element of the coupled system is shown in Fig. 6.

The chaotic state is adjacent to a period-three window and a period doubling sequence. The reconstructed attractor (Fig. 6a) is low-dimensional (correlation dimension $`2.3\pm 0.1`$). To obtain the phase, the Hilbert transform approach is applied Pikovsky et al. (2001). The two-dimensional embedding using Hilbert transform is shown in Fig. 6b. As in the case of globally coupled Rössler oscillators, it is difficult to distinguish the banded character of the attractor.

The merged banded structure is more clearly seen in return maps obtained from the series of return times ($`T_n`$). The series of the return time clearly shows an alternation of long and short oscillations (Fig. 7a). The map constructed using $`T_n`$ exhibits an approximately one-dimensional character (Fig. 7b).

The oscillators in the population are slightly different. Figure 8a shows a superposition of the 1D maps of all the oscillators.

Clearly, the majority of the phase points lie in the upper left and lower right boxes indicating two-banded character. The probability of two-banded oscillations varies from one site to the other (Fig. 8b) with a mean of $`p_{\mathrm{m}ean}=0.85`$ and a standard deviation of $`0.05`$. We note that the two-banded character of the uncoupled system of uncoupled oscillators ($`\epsilon =0`$, $`p_{\mathrm{m}ean}=0.7`$, Figs. 8c and d) is smaller than that of the phase synchronized region. Therefore, during the transition to phase synchronization with increasing the coupling strength the two-bandedness observed in the time series of oscillations is more pronounced.

As expected from the numerical simulations and theoretical considerations, the strong (anti)correlations in the phase dynamics have pronounced effects on the speed of convergence of frequencies.

Figure 9 shows that the standard deviation of frequencies is proportional to $`1/T`$ at $`\epsilon =0.1`$. (The fluctuations around the fitted line are caused by slowing down and speeding up within a cycle.) Similar scaling results are obtained for other (larger) coupling strengths.

## V Conclusion

For many systems of coupled oscillators, the chaotic state is reached through a period-doubling cascade. This implies that a Poincaré surface exists such that the next-amplitude map is similar to that in Figs. 3 or 5. As long as the correspondence between the next-amplitude map and the relation of subsequent return times holds, this should lead to a small value of $`D_\eta `$ and, in particular, to a $`1/T`$ scaling in the phase synchronized state. The experimental results show that such a structure of the series of return times can exists even if no well-behaved next-amplitude map can be identified. Our findings further suggest that chaotic oscillators where the chaotic state is reached through a period-doubling cascade can generally be phase synchronized due to the small value of $`D_\eta `$, i.e., high degree of phase coherence. This is in accord with findings in Ref. Pikovsky et al. (2001) where different chaotic systems were analyzed with respect to their ability to be phase synchronized. There it was found that the Lorenz system (where the chaotic attractor is not reached through a period-doubling cascade) cannot be phase synchronized. In particular, the next amplitude map has a very different structure, i.e., a logarithmic singularity.

Finally, in many biological and physical systems only short time series are available. Our results show that the determination of frequencies in weakly coupled oscillatory systems may be obtainable from such short time series.

###### Acknowledgements.

This work (JLH and IZK) was supported in part by the National Science Foundation (CTS-0000483) and the Office of Naval Research (N00014-01-1-0603). RK and JD were supported in part by a grant from the Natural Sciences and Engineering Research Council of Canada. We thank Yumei Zhai for help with the experiments.