Doppler cooling of a Coulomb crystal

We study theoretically Doppler laser-cooling of a cluster of 2-level atoms confined in a linear ion trap. Using several consecutive steps of averaging we derive, from the full quantum mechanical master equation, an equation for the total mechanical energy of the one dimensional crystal, defined on a coarse-grained energy scale whose grid size is smaller than the linewidth of the electronic transition. This equation describes the cooling dynamics for an arbitrary number of ions and in the quantum regime. We discuss the validity of the ergodic assumption (i.e. that the phase space distribution is only a function of energy). From our equation we derive the semiclassical limit (i.e. when the mechanical motion can be treated classically) and the Lamb-Dicke limit (i.e. when the size of the mechanical wave function is much smaller than the laser wavelength). We find a Fokker-Planck equation for the total mechanical energy of the system, whose solution is in agreement with previous analytical calculations which were based on different assumptions and valid only in their specific regimes. Finally, in the classical limit we derive an analytic expression for the average coupling, by light scattering, between motional states at different energies.


I. INTRODUCTION
The development of laser cooling and trapping techniques in the last decades has allowed for many spectacular experimental achievements [1]. Hand in hand with the experimental work, a theory of laser cooling has been widely developed, providing a precise description of many experimental situations [2,3]. Yet, the great majority of the theoretical studies have dealt with the interaction of laser light with single particles (being the same as many noninteracting particles). The treatment of a many-body system coupled to light is considerably more complex, due to the large number of degrees of freedom which results from the interaction among the constituents of the system.
In this work we investigate laser cooling of a many-body system, taking as a representative -and experimentally relevant -example a one-dimensional Coulomb cluster, i.e. a crystallized structure of ions which are trapped in a Paul or Penning trap [4][5][6][7]. We develop a model for the dynamical behaviour of the crystal's mechanical energy. In this system the physical processes in play are the trapping potential which confines the ions, their mutual Coulomb interaction, and their interaction with laser light. For sufficiently small kinetic energy the motion of the ions is properly described by the collective excitations of the cluster, the so-called harmonic modes of the crystal, due to the interplay of trapping potential and Coulomb repulsion. The problem is then characterized by three main frequency scales: (i) the oscillation frequencies of these modes, which are determined by the trap frequency ν and the number of particle N (for ions of equal mass and charge); (ii) the recoil frequency ω R characterizing the exchange of mechanical energy between radiation and atoms, ω R = 2π 2h /mλ 2 , with m mass of the ions and λ wavelength of the light coupling quasi-resonantly to the electronic transition, and (iii) the linewidth of the electronic transition γ, which characterizes the rate at which photons are scattered by the atoms. Doppler cooling of particles in a harmonic trap focuses on the regime γ > ν.
The theory is well developed for single trapped particles in specific regimes: in the semiclassical limit γ ≫ ν, ω R , when the mechanical motion can be treated classically, and in the Lamb-Dicke limit ν ≫ ω R , when the size of the motional wave function is much smaller than the laser wavelength [3]. In the semiclassical case, the standard procedure consists in writing the equations for the system dynamics in the Wigner representation, adiabatically eliminating the excited electronic state and then expanding inh up to the second order. In this way a Fokker-Planck equation in the position and the velocity of the ion is derived [9]. In the Lamb-Dicke limit, corresponding to an expansion in ω R /ν up to the first order, the cooling dynamics are described by a set of rate equations projected on the electronic ground state and on the eigenstates of the motion, which can be analytically handled [3,8]. Both treatments yield in the limit γ ≫ ν the same final distribution of cooled atoms in energy space, and the same dependence of the final temperature on the cooling parameters. Our approach shows how they are special cases of the same equation that we derive for the mechanical energy of the crystal.
The presence of many mutually interacting particles complicates the treatment considerably, because the mechanical Hamiltonian, having an increased number of degrees of freedom, often does not allow for a simple and transparent solution. An immediately visible effect arises in the spectrum of the mechanical energy, where the increased number of degrees of freedom is in general connected with the appearance of quasi-degeneracies and with a dense distribution of the energy levels [10]. The situation is facilitated again if one restricts the treatment to the Lamb-Dicke regime [11,12] but in the more general case we want to discuss, different steps of simplification can be made as will be shown.
The treatment in this paper focuses on Doppler cooling of an N -body one-dimensional Coulomb crystal, in particular on the type of linear ion chains obtained in linear Paul traps [7,13]. Here, the many-body mechanical Hamiltonian can be approximated by a set of harmonic oscillator modes, all having different frequencies ν 1 , ..., ν N , where each mode corresponds to a collective excitation of the crystal. Since for Doppler cooling the rate of photon scattering, and thus the cooling dynamics is determined by γ, we define a grid in energy space ∆E ≪hγ and study the cooling dynamics among the energy shells defined by the grid. From this starting point, with a procedure of successive averaging we derive a rate equation for the population at each motional energy on the grid. From this result, the semiclassical limit is recovered by expanding in ω R /γ. In particular, we get a Fokker-Planck equation for the motional energy describing cooling of N ions, whose solution agrees with the ones of standard semiclassical treatments [9,12,14,15]. In the Lamb-Dicke regime, starting from the full master equation, we discuss the conditions under which a compact equation for cooling can be found. Finally, we derive an explicit form of the rate equation in the limit γ ≫ ν by evaluating an analytic expression for the average coupling by light scattering between motional states at different energies.
Although we restrict the investigation to a one-dimensional Coulomb crystal, the theory can be directly extended to three dimensions, and the results are applicable to clusters in both Paul and Penning traps. More in general, the idea behind our treatment is that -in incoherent processes -if the rate determining the dynamics of interest of the system can be singled out, the contribution of processes occurring on faster time scales is often well represented by the average value on the slower time scale, characterizing the incoherent process. This work is organized as follows. In Section II the model from which our derivation starts is introduced and discussed. In Section III we discuss and apply a series of approximation, from which we obtain an equation for the cluster's mechanical energy describing cooling of the crystal. In Section IV we study the equation in the Lamb-Dicke and Semiclassical limits, and we compare our results with those obtained in the same limits with other treatments. In Section V we derive an explicit functional form of the rate equation in the limit γ ≫ ν. Finally, we draw our conclusions.

II. MODEL
In this section we introduce the model which we study throughout the paper. We first discuss the mechanical properties of a one-dimensional Coulomb crystal, then the interaction of laser light with the internal electronic transition of each ion, and finally the mechanical action of the light on the collective mechanical degrees of freedom of the crystal.

Mechanical properties
The mechanical potential on an ion cloud in a Paul trap is the sum of the potential exerted by the trap on each ion and of the Coulomb repulsion among the ions. Sufficiently far away from the trap electrodes, the potential of a Paul trap can be considered harmonic and the total potential for N ions of charge e has the form: where r j = (x j , y j , z j ) is the position of ion j and u j 0x , u j 0y , u j 0z depend on the trap parameters and on the mass of the jth ion. At sufficiently low temperatures the ions crystallize around the classical equilibrium positions r (0) i which are solutions of the set of equations ∂V /∂ r i | r (0) i = 0 [4]. Sufficiently below the crystallization temperature the motion of the ions around these equilibrium points is harmonic to good approximation. In this limit, the total mechanical potential can be described by its Taylor series truncated at the second order in the expansion around { r (0) i }. In a Paul trap with cylindrical geometry and very steep potential in the radial (x andŷ) direction the ions crystallize along the trap (ẑ-) axis. In this regime the amplitude of the radial oscillations is much smaller than the axial ones. Here, we assume the radial degrees of freedom to be frozen out such that the motion is one-dimensional along theẑ-axis, and we take the axial potential to be electrostatic, as it is the case in a linear Paul trap [13]. Then the truncated mechanical potential has the form: where V jk is a real, symmetric and non-negative matrix (the explicit form of its elements can be found, for example, in [16]), while q j = z j − z j . Given m mass of each ion, the secular equation for the harmonic motion of the crystal has the form: where ν α are the eigenvalues and b α j the associated eigenvectors, which are complete and orthonormal. In the basis q ′ α = i b α i q i the motion is described by the Lagrangian for N independent harmonic oscillators of frequency ν α . Given the canonical momentum p ′ α =q ′ α conjugate to q ′ α , the motion is quantized by associating a quantum mechanical oscillator with each mode. Denoting by a α and a † α the annihilation and creation operators for the mode α, respectively, the coordinates q i are now written as and the Hamiltonian for the mechanical motion has the form: The energy eigenstates of each mode α are the number states |n α with eigenvalues (energies) E nα =h(n α + 1 2 )ν α . In the following we use states |n = |n 1 , n 2 , ...n N to describe the eigenstates of the motion of the crystal corresponding to the eigenvalues E n = α E nα .
In summary, N trapped ions crystallized in one dimension can be described by N harmonic oscillators with frequencies ν 1 , ..., ν N . These frequencies are solutions of (3) and, what is important, they are incommensurate. Therefore the spectrum of mechanical energies of the system does not exhibit the discreteness and equispacing property of the single harmonic oscillator spectrum, but shows a dense distribution of levels, and at sufficiently high energies and large N it assumes a quasi-continuum character. In this limit we can define the density of states of the system, g(E), which is a smooth function of the energy and is defined on a grid of energies ∆E, such that the number of states D(E) contained in the interval F (E) = [E − ∆E/2, E + ∆E/2] is D(E) ≫ 1, and D(E) ≃ g(E)∆E. For N modes (ions) g(E) can be evaluated by solving the following integral: i.e. by taking the derivative of the ratio between the volume in phase space of energy ≤ E and the volume occupied by a single state, (hν 1 )(hν 2 )...(hν N ). In Fig. 1 we plot the number of states as a function of the total mechanical energy for a chain of three ions, taking a grid ∆E =hν 1 /5, and compare the exact value with the smoothed function g(E)∆E, with g(E) given by (6).

Interaction with light
We consider that laser light with frequency ω L and wave vector k = (k x , k y , k z ) drives some ions of the crystal, coupling to their internal two-level transition with electronic ground state |g , excited state |e and resonance frequency ω 0 . In the Rotating Wave Approximation and in the frame rotating with the laser frequency, the total Hamiltonian has the form: Here H at is the Hamiltionian for the internal degrees of freedom, defined as where δ = ω L − ω 0 is the detuning of the laser from the atomic resonance and j labels the ions. H mec is the mechanical Hamiltonian defined in (5), and H AL describes the interaction between laser and atoms, where Ω(z j ) is the Rabi frequency at the position z j , σ + j , σ − j are the raising and lowering dipole operators, respectively, for the internal state of the j th ion, while {j} is the set of driven ions. Assuming that the light intensity does not vary rapidly in the vicinity of z The operator e −ikz zj in (9) represents the mechanical effect of the light interaction, i.e. a shift in momentum space by one photon recoil which goes along with the excitation. Using (4), its explicit form is: where η α j is the Lamb-Dicke parameter for mode α and ion j, defined as [17]: For some ions we can have Ω j = 0, which means that not all ions of the chain are driven. Such condition can be achieved either when the ions are sufficiently spaced to allow for their individual addressing [18], or by introducing a different type of ion into the crystal, such as an isotope of the trapped element [19], or a different species [17,20] whose transition frequency is not resonant with the laser. In this case, some differences in the mechanical behaviour arise, which, however, are not relevant for the results that we derive below.

Master Equation
We describe the dynamics of the driven crystal through the master equation for the density matrix ρ of the N -ion system: where H has been defined in (7), and L is the Liouvillian describing the incoherent evolution of the system due to its coupling to the modes of the electromagnetic field: Here γ is the decay rate of the internal excited states |e j , andρ j is the density matrix describing the feeding for the ions j, shifted in momentum space by the recoil of the spontaneously emitted photon, with N (cos θ) being the dipole pattern of the spontaneous decay, and k = | k|.

Low saturation limit
In the limit of low saturation, Ω ≪ γ, the excited states |e j can be eliminated from (12) in second order perturbation theory [21]. Thereby we obtain a closed equation for the internal ground state |g = N j=1 |g j , which we project on the basis of mechanical states |n . The equations describing the evolution of this system are [10]: Here, the coefficients describing the coupling among the populations, n|ρ|n , and the coherences, n|ρ|m , are proportional to the Franck-Condon coefficients, which have the form: where r α = min(l α , n α ), and L |lα−nα| rα (x) is a generalized Laguerre polynomial. They represent the probability amplitude of a transition from the initial motional state |n to the final motional state |l by absorption or emission of a photon.

III. AN ERGODIC EQUATION FOR LASER COOLING OF THE CRYSTAL
In the following, starting from the full quantum mechanical Eqs. (15) and using consecutive steps of averaging, we derive an equation for Doppler cooling of an ion crystal. This equation will be ergodic, in the sense that it describes the dynamics of the system in terms of the population P (E) at the crystal's mechanical energy E.
In the Doppler cooling limit (which implies γ > ν α for all α = 1, ..., N , and δ = O(γ)) we define a grid of energies ∆E such that in an interval F (E) = [E − ∆E/2, E + ∆E/2] the number of states D(E) ≫ 1, and the Lorentzian describing the resonant response of the atom varies infinitesimally on each interval. The resulting shells of energies can be visualized as shown in Fig. 2 for the case of a two-ion crystal, where each axis represents the energy of one mode, and each point is a state |n = |n 1 , n 2 . Here, the states with total energy E fall into the shell of width ∆E along the line E =hν 1 (n 1 + 1 2 ) +hν 2 (n 2 + 1 2 ).
Coarse grained energy space for the case of two modes of frequency ν1, ν2. The points are the states with energy En = En 1 + En 2 . The broad lines represent two energy shells F(E) and F(E ′ ).
We can now decompose the sums over the motional states appearing in (15) as follows: where now E k is one energy value on the grid, and F (E k ) = [E k − ∆E/2, E k + ∆E/2] is the energy interval centered in E k with width ∆E. Using these definitions we rewrite the equations for the populations in (15) and sum over all states that lie within the same energy shell F (E k ): where we have separated the terms involving the populations, (18) and (20), from the ones which involve the coherences, (19) and (21). We now compare the term (18), the loss rate of the population, with the coupling to the coherences in (19): Line (18) is proportional to the sum of moduli square of Franck-Condon-coefficients, thus it is the sum of positive terms. The sum in line (19) adds up coefficients with alternating signs, as can be verified from Eq. (16). For D(E) ≫ 1 we can therefore safely assume that This is equivalent to a Random Phase Approximation. On the basis of this consideration the coupling between populations and coherences in (19) may be neglected in comparison to the loss rate described by the term in (18). Analogously we can neglect the terms in (21) in comparison with (20), and we obtain the following set of rate equations: +γ where all Rabi frequencies have been set equal, Ω j = Ω. This simplifies the discussion in the next subsection, but the more general case of different Ω j can be treated with the same methods under conditions which will be discussed below.

A. Ergodic hypothesis
Equation (23) describes how laser light couples to the internal electronic transition and exchanges mechanical energy with the crystal. The probability amplitude for the radiative event to occur is weighted by a Lorentzian distribution, which below saturation has width γ, and which is a function of the difference between the mechanical energies of the crystal before and after the scattering. The mechanical effect of light on the motion of the crystal is described by the operator (10), and the probability for the crystal initially in state |n to be scattered into state |k by the absorption or emission of a laser photon is given by the modulus square of the Franck-Condon coefficient, Eq. (16). When we regard the modulus square of Eq. (16) as a distribution over the states |k after an absorption or emission event, given the initial state |n , then the average motional energy transferred to the ion cluster and its variance are calculated as (see Appendix A): where ω R =hk 2 /2m is the recoil energy. Thus, the first and second moments of this distribution depend on the (single-ion) recoil energy ω R and on the energy per ion of the initial state, but not on the details of the quantum state |n . Guided by this result, we will now make the approximation that in the limit of Doppler cooling, γ > ν j , and for sufficiently large density of states, D(E) ≫ 1, the oscillations of the Franck-Condon coefficients with the vibrational numbers play a negligible role, whereas the average properties determine the cooling dynamics. We formulate this assumption by first introducing an average Franck-Condon coupling between shells of energy E n and E k , and then writing thereby neglecting the dependence of the LHS on the details of the state |n . We will discuss this assumption and possible alternatives in the next section. It will be one of the important results later to determine an explicit expression for Q (k) (E n , E k ) in Eq. (26). Furthermore, in (27) all ions of the crystal are assumed to be driven. In the next section we will discuss the case in which only a set is driven.
Applying (27) to Eq. (23), in the limit in which all ions are driven, we obtain In order to obtain the second line of (28) from the second line of (23) we have made a further approximation: First we write and then we set δ jl → 1/N . This corresponds to the assumption that, in the regime we are considering, the mechanical effect of the process of absorption + emission of a photon by one ion is equivalent to absorption by one ion and emission by another ion, weighted by the probability for the two ions to be the same. By defining the population densities P (E, t) of the energy shells through such that ∞ 0 dEP (E, t) = 1, and by using D(E) = g(E)∆E, we finally arrive at a rate equation as a function of the motional energy: where we have replaced the sums over the energy by integrals, valid when the average energy E(t) ≫ ∆E. The two parts of the rate equation (31) show how, after the various steps of averaging, the total population of a shell with energy E changes in time: Population is lost by excitation to shells with energy E 1 and subsequent emission into any shell, and population originally at E 2 is excited to a shell at E 1 and then scattered into a shell at energy E.
It should be noted that the restriction to a one-dimensional crystal enters into the number of modes and into the geometry of the laser beams, as well as into the pattern of the emitted radiation. However, the application of the "Random Phase Approximation" (22) and of assumption (27), which are at the basis of our derivation, are in no way restricted by the dimensionality of the problem. Actually, an extension to three dimensions would endorse the two approximations, since the number of modes increases, and with it the density of states in the spectrum of the motional energies. Therefore, an equation for the total motional energy of the crystal of the form of (31) can be derived in three dimensions using the same considerations done here for the one-dimensional case.
B. Discussion

Ergodic assumption
Since the distribution of the total mechanical energy over the system is sufficient for describing the cooling process, assumption (27), leading to Eq. (31), simplifies the description of the laser cooling dynamics of the crystal significantly by reducing dramatically the dimensionality of the problem. Equation (31) could have been also obtained by assuming that the population of all states in the same energy shell is equal, thus defining p(E) as the population of a state at energy E, such that (c.f. Eq. (30)) Equation (32) implies (27) and is an even stronger assumption. It corresponds to assuming that the system behaves ergodically, i.e. that the populations of the states in an energy shell equalize faster than the average quantities of the system evolve in time. It leads to the same rate equation (31), but it is a more natural assumption in the last stages of the evolution, when the system tends asymptotically to the steady state. In general, (31) describes the dynamics of cooling to the extent that assumption (27) is valid, i.e. when the average coupling between an eigenstate of the mechanical energy and the states of an energy shell is a smooth function of their respective energy. This is true when the average coupling of each state of one energy shell to the states of another energy shell is of the same order of magnitude. Outside the Lamb-Dicke regime, given the oscillatory behaviour of the Franck-Condon coefficients, this holds if each state couples appreciably to more than one state of the other energy shell, i.e. for D(E) ≫ 1 and γ sufficiently large. In the Lamb-Dicke regime analogous considerations can be applied, and they will be discussed in detail in the next section.
The ergodic regime can also be justified by the existence of a physical process whose main effect is to thermalize the states within one energy shell and which acts on a shorter time scale than the energy-changing processes. This assumption is at the basis of treatments in the kinetic theory of quantum gases [22], where the interatomic collisions lead to thermalization, and the gas can be considered in a thermal quasi-equilibrium distribution on the time scale in which it is cooled. The assumption of rapid thermalization is also central to an earlier study of laser cooling of Coulomb clusters [12]: There the effect of mode-mode coupling (anharmonicity) has been proposed as a possible agent, which does not explicitly appear in the equations but justifies condition (27). Direct evidence of this effect has been found in numerical studies in [17] for the case of exact degeneracy between the modes frequencies.

On the laser intensity distribution over the crystal
In deriving (31) we have assumed that all ions are uniformly driven. We discuss now the case in which only a subset {j} of ions in the crystal are illuminated. For simplicity, let us assume in Eq. (23) that Ω l = Ω for l ∈ {j}, and Ω l = 0 otherwise. Due to the geometry of the crystal, the coupling of each ion to a mode is a function of the mode and of the ion's position in the crystal [see Eq. (11)]. Some positions in the chain are more strongly coupled to a certain set of modes, and some positions are even decoupled from certain modes [17]. Then, the ergodic assumption is valid provided a suitable set of ions is driven, such that the coupling of each mode of collective motion to the radiation is of the same order of magnitude. In that case, the ergodic assumption, Eq. (27), can be expressed as: where M is the number of driven ions (M ≤ N ). Taking into account this more general case, the rate equation (31) has the form: Finally, note that the number of driven ions, M , appears in (31) and (35) as an overall scaling factor, making the optical pumping rate of the crystal M times that of a single ion. In fact, it has been shown earlier that below saturation the ions behave as independent scatterers, and their contributions to cooling simply add up [10,17].

IV. THE SEMICLASSICAL LIMIT AND THE LAMB-DICKE LIMIT
In this section, starting from Eq. (31) or (35) we derive the limit where the motion can be treated classically (semiclassical limit), and obtain a Fokker-Planck-equation for the energy which is analytically solvable. We compare our result with the well-known treatments of [9] and [12] and find full agreement with their theoretical predictions. In the second part of this section, starting from (15) we consider the Lamb-Dicke limit, and discuss the conditions under which an equation for the motional energy of the type of (31) can be derived for describing cooling of N ions. For clarity in the derivation, we rewrite Eq. (31) as: Here, p(E) is the population of a state of energy E defined in (32), and L(E 1 − E) is the Lorentzian distribution, From this definition, using (27) and (24)(25), it can be verified that f The analogous relations follow for f (kz) E (E ′ ) when θ is replaced by θ 0 , where k z = k cos θ 0 . Equations (40) and (41) are equivalent to (24) and (25) in the approximation of a smooth energy scale.

A. Derivation of a Fokker-Planck equation for the energy
We first consider the limit in which the recoil frequency is much smaller than the linewidth, ω R ≪ γ. In this limit, while f We now expand Eq. (36) around E up to first order in the parameter ω R /γ and, using (39-41), we obtain: Here, and α = d(cos θ) cos 2 θN (cos θ) .
We rescale the time as and define C = 2 cos 2 θ 0 (cos 2 θ 0 + α) Note that that with definition (44), C < 0 for red detunings δ < 0 (ω L < ω 0 ). Eq. (43) can now be written as: Using the energy dependence of the smoothed density of states g(E) as given in (6), and (33), we get a Fokker-Planck Equation of the form: where Eq. (49) can be easily solved. In the following, we calculate the steady state solution and the time evolution of the system, and compare the result with the existing treatments evaluated in the limit of one ion (N = 1).

Steady State
The steady state distribution P 0 (E) satisfies the equation d dτ P 0 (E) = 0. Thus, it is solution of the differential equation: (where the integration constant has been set to 0 for a convergent solution) and has the form with normalization constant F . For red detunings (C < 0) and if the wave vector of the cooling laser has a component along the motional axis (cos θ 0 = 0) the integral of (51) over the energies converges. Then, the value of F is found from ∞ 0 dEP 0 (E) = 1, yielding: Using these results, the steady state energy has the form: Eq. (53) contains the dependence of the cooling limit on the angle θ 0 between the cooling laser beam and the direction of the motion. The final energy is minimal for cos θ 0 = 1, i.e. when the laser propagates parallel to the trap axis, and it diverges for cos θ 0 = 0, when the laser is orthogonal to the trap axis and there is thus no laser cooling.
The minimum of the final energy vs. detuning is reached for δ = −γ/2, as in the case of one ion (see for example [3]). Inserting into (53) N = 1, α = 1/3 (which corresponds to spatially isotropic spontaneous emission), and cos θ 0 = 1, we find the same result as Javanainen and Stenholm in their semiclassical expansion for one ion [9]. Hence, the final energy for N ions is N times the steady state energy achieved by Doppler cooling of one ion. This general result has also been found earlier in special cases, such as a Coulomb cluster in the Lamb-Dicke regime [12], and a two-ion crystal treated in Ref. [23] by extending the method of [9].

Time Evolution
The steady state solution suggests to use the following ansatz for solving Eq. (49): where U (t) is a positive function of time and F (t) is a normalization factor, such that at any instant t the relation ∞ 0 dEF (t)P (E, t) = 1 is fulfilled, i.e.
This ansatz corresponds to assuming that the distribution P (E, t) is always thermal with average energy Substituting (54) into Eq. (49) and using (46) we get a differential equation for U : Since this equation must be fulfilled for all values of E, we find: which has the solution: where N U 0 = E(0) is the initial energy. For t → ∞ we recover from Eq. (59) the steady state solution (51). The rate at which the steady state is reached is: where we have used (44). Thus the cooling rate increases linearly with the number of driven ions M . It is largest if all ions are driven, as expected.
Since the results derived so far agree precisely with those found earlier in specific cases, the general procedure which lead us from the quantum mechanical equations to the rate equation for the energy can be considered the common basis which underlies and unifies these earlier treatments. Furthermore, we have shown that the final energy of an N -ion crystal is N times that of a single motional mode, while the cooling rate is M times that for the one-dimensional motion of a single ion, M being the number of ions driven by the laser.
A final remark should be made on the assumptions leading to (49). In deriving the Fokker-Planck equation we have assumed ergodicity and Eq. (42), i.e. that p(E) varies negligibly on the recoil energy scale. The latter condition corresponds to the second-order expansion inh of Ref. [9]. In that work the derivation of a Fokker-Planck equation was based on the limit of overdamped oscillation, γ ≫ ν, in order to adiabatically eliminate the excited state from the equations. Our derivation, however, does not necessarily imply this limit. Only in the case of N = 1 ion, which is not the main focus of our study, we must have ∆E ≫hν in order to fulfill the condition of a large density of states, D(E) ≫ 1, and thus for this special case the overdamped oscillator limit is a requirement for the validity of (49).

B. Lamb-Dicke regime
When the atomic motion is well localized with respect to the laser wavelength (Lamb-Dicke regime), the Franck-Condon coefficients in (16) can be approximated by their first order expansion in the parameter ( k · x) 2 . For a single ion excited below saturation, the dynamics are described by a rate equation for the populations of the states with vibrational number n, which can be analytically solved [3,8]. In this form, since n is proportional to the mechanical energy of the ion, the equation for cooling is an equation for the energy. For many ions, a set of rate equations can also be derived which describe cooling of each mode and which are decoupled, since a simultaneous change of the energy of two or more modes by scattering of a photon is of higher order in the Lamb-Dicke parameter (and thus takes place on a longer time scale) [10,17].
Yet we show in this section that under some conditions we can derive a single equation for the total energy of many ions. Given the analytical simplicity of the model, this example is also instructive in order to see the limits of validity of the ergodic equation.
Let us consider Eqs. (15) in the Lamb-Dicke regime [24]. Here the reduction of this set of equations into the rate equations (23) is possible without assumption (22), since the terms (19), (21) are negligible either because they are rapidly oscillating or because their coupling to the populations is of higher order in the Lamb-Dicke parameter [10]. In first order in the parameter ω R /ν β , with β = 1, ..., N , the modulus square of the Franck-Condon coefficient has the form: Thus, given the initial state of motion |n = |n 1 , ..., n N , scattering of a photon involves three possible transitions: (i) the so-called carrier transition, of zero order in the Lamb-Dicke parameter, with final state |k = |n ; (ii) the red-sideband transitions, where |k = |n 1 , ..., n β − 1, ..., n N , and (iii) the blue-sideband transitions, where |k = |n 1 , ..., n β + 1, ..., n N . The cases (ii) and (iii) are of second order in the parameter η β . When we substitute (61) into Eq. (23), we obtain: n + 1 β |ρ|n + 1 β ω R ν β (n β + 1) cos 2 θ 0 L(−ν β ) + αL(0) where have used j η β j 2 = ω R /ν β , and the vector 1 β is defined in the N -dimensional Hilbert space through (1 β ) α = δ α,β for α = 1, ..., N . Eq. (62) is linear in the vibrational numbers n β of the single modes, which are weighted by the factors ω R L(ν β )/ν β . Only when these weights are of approximately equal magnitude, a hypothesis similar to (27), which allows summation over states of equal energy, can be applied. This is true when ν β ≪ γ for all modes β because in this limit the Lorentzian atomic resonance curve varies very slowly over all red (blue) sidebands of the modes. Then, by expanding the terms in (62) at the second order in the parameter ν β /γ, we obtain a Fokker-Planck equation of the same form as (49). In that sense, the condition ν β ≪ γ can be considered a requirement for deriving an ergodic equation in the Lamb-Dicke regime.
Alternatively, in [12] a single equation for the total motional energy in the Lamb-Dicke regime has been justified by assuming that all modes thermalize on a faster time-scale and thus imposing Application of (63) to (62) yields an equation for the total mechanical energy which can be solved analytically. In our treatment above, in the case where we find a Fokker-Planck equation, condition (63) is also fulfilled, but there it is a consequence of the average coupling among energy shells due to light scattering, rather than an extra assumption.

C. Discussion
The Fokker-Planck equation (49) has been derived in two different ways: in the semiclassical case by starting from the ergodic equation (31), and in the Lamb-Dicke regime by starting from equation (23). Both derivations rely on the limit γ ≫ ν, ω R , while in the Lamb-Dicke regime there is the additional constraint that ω R ≪ ν.
Let us now summarize the solutions of (49). The form of the solution is the same as the one obtained for cooling of one ion, whereby here the number N of ions appears in the steady-state energy (53) as a scaling factor, and the cooling rate (60) scales as the ratio M/N where M is the number of driven ions. This result supports the conclusion drawn in [12], that the cluster is cooled like a single ion. It is instructive, at this regard, to compare two particular cases which exhibit analogies: The axial motion of two ions in a trap with collective modes at frequencies ν 1 , ν 2 , and the motion of a single ion in a two-dimensional harmonic oscillator potential with the same frequencies, ν 1 and ν 2 , on the two axes. Neglecting for the moment the different masses, the mechanical Hamiltonians of the two systems are equivalent. The laser-ion interaction terms are also equivalent, provided that one of the two ions in the chain is driven, and that the Lamb-Dicke parameter for each mode of the chain is the same as the Lamb-Dicke parameter for each axis of the trap of the single ion. This latter condition can be fulfilled by choosing, in the single-ion case, the proper angle of propagation of the laser in the two-dimensional plane.
The main difference between the two cases arises in the spontaneous emission: in the two-dimensional one-ion case, the photon is scattered at a random angle in the plane and the ratio between the projections of the recoil energy on each axis can vary. In contrast, in the two-ion case, the scattering angle is the same for both modes, hence the average energy transfer to the two modes has a fixed ratio in any scattering event. This difference appears as a the geometrical factor in the dynamical equations, as well as in the expression for the final energy (53) [15]. From that consideration we expect that the generalization of our treatment to a 3-dimensional Coulomb crystal will change the one-dimensional results only by numerical factors representing the different geometry of the problem.

V. EVALUATION OF THE SEMICLASSICAL ERGODIC EQUATION
So far we have assumed a grid of energy ∆E, which defines the coarse-graining on which the relevant system properties for cooling are defined. The grid ∆E has been chosen according to the conditionhγ ≫ ∆E, and an equation for the energy has been derived for the situation that the number of states at energy E is D(E) ≫ 1. Thus, ∆E represents a limited resolution which allows us to average over many quantum states. We now discuss the limit ∆E ≫hν N . This case is expected to correspond to the classical limit, since the resolution is such that all details of the quantum spectrum are averaged out. We will derive the classical limit starting from the average Franck-Condon coefficients Q (k) (E, E ′ ) in (26), from which we obtain an explicit form for the ergodic equation (31).
Let us consider the sum in Eq. (26). Using the property of the trace, we can write: The summation over the states belonging to the shell of energy E can be re-expressed as: and Eq. (64) aquires the form: Using the Fourier transform of the δ-function, we get: such that Eq. (66) can be rewritten as: where we have defined: By using coherent states |α β to calculate the trace, using (4), (10), and the properties of coherent states [21], Eq. (69) takes the form: In Eq. (68), the integrand is the product of A j (τ, τ ′ ) and the window function sin(x)/x of widthh/∆E. According to Eq. (70), A j (τ, τ ′ ) depends on τ, τ ′ only through e iν β τ and e iν β τ ′ , and thus τ is scaled by the mode frequencies.
A. Derivation of the classical limit We now consider the limit E, ∆E ≫hν β for all β. This corresponds to the situation n β ≫ 1, e.g. before the cooling limit is reached, and in general to the regime γ ≫ ν β . In this limit we can expand the exponentials in (70) in the parameterhν/∆E. Up to first order, Eq. (68) is equivalent to: This result can be interpreted physically, considering that τ is the Fourier conjugate of the energy E/h and has the dimension of a time. Thus, Eq. (69) can be seen as the overlap between the state |α and the state |α ′ , which corresponds to a displacement of |α , followed by free evolution for a time τ , then by a displacement back, and finally by free evolution for a time τ ′ . Integral (68) sums over all intervals of times τ , τ ′ , for all trajectories in the neighbourhood of the energy shells E at τ and E ′ at τ ′ , weighted by the time window of resolutionh/∆E. For τ = τ ′ = 0 the overlap is maximum, since the state is unchanged. Although the overlap shows a certain periodicity, e.g. for N = 1 and τ, τ ′ multiples of 2π/ν, A(τ, τ ) = A(0, 0), for ∆E ≫hν these recurrence times fall out of the window function interval. Therefore the only appreciable contribution to the integral comes from τ, τ ′ ∼ 0. In other words, the classical limit corresponds to large energy uncertainty, such that only the short time evolution of the wave packet contributes appreciably to the integral, since only in this limit the trajectories are phase coherent. Eq. (71) can now be rewritten as an integral in the classical phase space using the definition of coherent states whereq β ,p β are real numbers. Integrating over τ, τ ′ yields: Integrating in phase space [see Appendix C], and using the relation f E (E ′ ) = C(E, E ′ )/(D(E ′ )∆E) we obtain: The function f E (E ′ ) is real, and thus well-defined, on the interval of energies [E +hω R − √ 4hω R E, E +hω R + √ 4hω R E], and it is normalized with respect to E ′ . In Fig. 3 we plot Eq. (74) for N = 1, 10, 100. For N = 1 it has the wellknown form of the classical momentum distribution of a harmonic oscillator at a given energy E: In fact in this case E ′ = E +hω R +hkp/m, i.e. the probability for the oscillator to have final energy E ′ after the scattering is equal to the probability that the oscillator, at energy E, has momentum p when the scattering event occurs. For higher N the distribution peaks more and more narrowly around E +hω R . Finally, note that, from the definition of f E (E ′ ) in (74) and from (38): This function is symmetric in E and E ′ as in the quantum mechanical case (see (26)). In the limit in which this result holds, we eventually find the explicit form of the ergodic rate equation by inserting Eq. (75) into Eq. (31). An explicit form of the rate equation can also be evaluated in the limit of N ≫ 1 ions [25].

VI. CONCLUSIONS
We have studied Doppler cooling of a Coulomb crystal: Starting from the full Master Equation, in the low saturation regime and by consecutive steps of averaging we have derived a rate equation for the total energy of the crystal, reducing dramatically the number of degrees of freedom and thus simplifying considerably the complexity of the problem. The equation is defined on a coarse-grained energy scale with grid ∆E such that ∆E ≪hγ, where γ is the linewidth of the electronic transition resonant with laser light. Its derivation is based on the quasi-continuum characteristic of the spectrum of motional energies E on ∆E, and on the assumption that the coupling between states belonging to different energy shells is a smooth function of E.
Starting from the general form of the equation, we have studied the semiclassical limit and the Lamb-Dicke limit, and in both cases we have found a Fokker-Planck equation describing cooling of a Coulomb crystal of N ions. The general solution agrees with the results of previous treatments, which were developed in perturbation theory for these limiting cases [3,8,9,12]. As observed in [12], the dynamics of cooling of a Coulomb crystal can be scaled to the one of a single ion.
In the semiclassical limit we have derived the explicit form of the rate equation for the mechanical energy, by calculating from the full quantum mechanical expression the classical probability of scattering between motional states at different energies. An explicit form of the equation can also be derived in the quantum limit for the case of N ≫ 1 ions. This derivation will be presented in future work [25].
Our results, with marginal changes, can be applied to a Coulomb crystal in three dimensions. In that case, the dimensionality enters into the number of modes, which for a crystal of N ions is 3N , and into the spatial distribution of the scattered photons, thus affecting the coefficients of the energy rate equation. From the comparison between the cooling problem for a one-dimensional crystal and for one ion in more dimensions, we expect the final energy in 3D to differ from the one-dimensional result only by a geometrical factor.
To conclude, we would like to remark on the generality of the treatment. Other incoherent processes in physical systems can be studied in an analogous way, provided that the rate determining the dynamics of interest can be singled out, and that on the corresponding energy scale the spectrum of energy levels is characterized by a quasi-continuum.

VII. ACKNOWLEDGEMENTS
The authors are deeply grateful to H. Walther and P. Zoller, who have motivated and stimulated this work in its various stages. Discussions with P. Zoller have been highly appreciated. G.M. thanks J.I. Cirac, B.-G. Englert, H. Hoffman, W. Schleich, and S. Stenholm for many stimulating discussions and helpful comments. This work has been partly supported by the European Commission (TMR networks ERB-FMRX-CT96-0077 and ERB-FMRX-CT96-0087) and by the Austrian Science Fund FWF (SFB15).

APPENDIX A: MOMENTS DERIVATION
One ion. Let us consider a one-dimensional harmonic oscillator of frequency ν and number state |n . Considering | n| exp(ikx)|k | 2 as a distribution over the final states |k of energy E k , given the initial state |n with energy E n , the first and the second moments of the distribution are: Using the property a † a|k = k|k and the closure relation of the states {|n }, we can contract the sum over k in (A1), (A2), using: k m | n|e iη(a † +a) |k | 2 = n|e −iη(a † +a) a † a m e iη(a † +a) |n .
Expression (A3) can be further simplified by using the commutation properties of the bosonic operators, and we obtain: where we have used η 2 = ω R /ν.
N ions. Let us now take an N -ion chain. Using the definitions of Section II, the first moment of the distribution is: where the subscript j refers to the driven ion. Analogously, the second moment has the form: Averaging (A7) over the ions of the chain we find: We define: I(E,p 1 ) = dq 1 ...dq N dp 2 ...dp N h N δ(E − N β=1p 2 β 2m − V (q 1 , ..,q N )), and move to the coordinates q β where V is a diagonal quadratic form. By introducing the set of rescaled variables Q β = mν 2 β /2q β , P β = 1/2mp β , integral (B2) is the measure of the surface of a unitary hypersphere in 2N − 1 dimensions. Integrating, we obtain: Substituting now this expression into (B1) and integrating overp 1 we get: where It can be easily verified that (B5) satisfies the relations (38) and (39-41).