Tuning nonthermal to thermal distributions in time-dependent Paul traps

# Tuning nonthermal to thermal distributions in time-dependent Paul traps

H. Landa Institut de Physique Théorique, Université Paris-Saclay, CEA, CNRS, 91191 Gif-sur-Yvette, France
###### Abstract

We study the probability distribution of an atomic ion being laser-cooled in a radio-frequency Paul trap, in terms of Hamiltonian action-angle coordinates appropriate for the time-dependent potential. Close to the effective minimum in many Paul traps, the periodically-driven potential is approximately quadratic in the coordinates. We find that despite the microscopic nonequilibrium forces, a stationary thermal-like exponential distribution in the action is reached after coarse-graining the angles, with a mean action corresponding to the Doppler cooling limit. However, residual electric fields can push the ion from the origin of the quadrupole potential and set it into a large-amplitude oscillation. Above a threshold amplitude of such “excess micromotion”, the action distribution for motion expanded about the driven oscillation broadens and becomes distinctly nonthermal. We find that by a proper choice of the laser detuning the distribution can be made exponential again, with a mean action close to the Doppler cooling limit. This is important for quantum information processing and other applications, and in particular the derived approach can be immediately applied to crystals of trapped ions in planar configurations, where the driven motion of ions is unavoidable.

## I Introduction

Perhaps the most common method of trapping charged particles in vacuum is by use of electrodynamic Paul traps paul1990electromagnetic . The electric potential generated in these traps is time-dependent, as the applied electrode voltages are periodically modulated at a radio-frequency (rf) rate. Close enough to the effective potential minimum in many Paul trap variants, the potential is well approximated by a quadrupole, quadratic in the coordinates. A periodically-driven quadrupole potential results in time-dependent, linear, Mathieu equations of motion. The motion of the ions is then characterized by quasi-periodic harmonic oscillations at the so-called “secular” frequencies, superimposed with the fast “micromotion” driven by the time-dependent potential. Since the equations of motion are integrable (being linear), the description of the motion can be simplified using classical action-angle coordinates. For a given initial condition, the trajectory in the phase-space of coordinates and momenta is restricted to rotations on an invariant manifold. This is a torus whose dimension is half that of the phase-space. The ion’s position on the torus at any time is given by the angles which are the generalized coordinates, each evolving independently and increasing linearly in time (with a fixed angular frequency). The measure within the torus (the bounded hyper-volume in the phase-space) is determined by the actions which in the absence of perturbations are conserved quantities of the motion, taking the role of generalized momenta.

An important tool employed to control the dynamics of an ion in the phase-space is laser cooling. Laser cooling is a well established method for removing entropy, e.g. from a trapped ion’s motion wineland1979laser ; javanainen1980 ; javanainen1980a ; javanainen1981laser ; gordon1980motion ; stenholm1986semiclassical ; cirac1994laser ; aura2002 ; leibfried2003 ; wesenberg2007 ; epstein2007 ; marciante2010 ; PhysRevA.96.012519 . A stochastic process by its very nature, laser cooling is based on the absorption of photons with a narrow momentum bandwidth determined by the laser, followed by a spontaneous emission of photons with randomly directed momenta.

In this work we use a recently developed semiclassical framework for studying laser cooling dynamics of ions driven by micromotion at arbitrary frequency and amplitude rfcooling . We focus on a single ion’s distribution in the final stage of the cooling. Since the Paul trap potential is time-dependent, the distribution function in the positions and velocities is non-stationary; it is periodic with the period of the trap cirac1994laser . A key idea underlying this treatment is that by using the action, which is a time-independent quantity, the coherent, driven component in the kinetic energy is readily separated from the stochastic component, which makes the action-angle coordinates a very useful choice. In terms of the action, it is easy to distinguish an effective, stationary thermal-like exponential distribution from other, nonthermal (and in general much more broad) distributions possible for the ion in the final stage of the cooling, arising for reasons that will be discussed in the following.

In Sec. II we present the setup that is studied in detail in this work, of an ion being laser-cooled in a region of a Paul trap where the potential is approximately a quadrupole potential. We present analytic and numerical results for the simplified case of one-dimensional (1D) motion, with the detailed derivations for three-dimensional (3D) motion left for the Appendix (see below). The main observation that we follow is that for typical laser cooling rates we can consider the motion as Hamiltonian over a large number of rotations on the invariant torus, and the laser can be modelled as acting on the action alone (with the angles averaged over). The effect of the laser is to permit the ion to drift and diffuse between the invariant tori of the Hamiltonian phase-space. Using a Fokker-Planck (FP) equation for the probability distribution of the ion in terms of the action allows us to account transparently for the micromotion.

The main results of this work are presented in Sec. III. For an ion within a time-independent harmonic trap, the probability distribution in the final stage of the cooling is known to be thermal javanainen1980 . We show that even with periodically-driven dynamics where the ion accelerates between the photon absorption and the spontaneous emission, a nearly identical, exponential distribution is obtained with a mean action corresponding to the Doppler limit. We extend our study to include stray electric fields which may push the ion away from the origin of the quadrupole potential, leading to “excess micromotion”. Treating excess micromotion in detail, we study how the distribution changes from thermal to very distinctly nonthermal as a function of the excess micromotion. Remarkably, we find that for any amplitude of excess micromotion, the distribution can be narrowed down into an exponential distribution close to that of the Doppler cooling limit, by increasing appropriately the laser detuning.

We conclude with a brief summary and an outlook for possible applications and generalizations in Sec. IV. In the Appendix, we lay down the details of the action-angle treatment of a coupled system of Mathieu oscillators in an arbitrary dimension. The exact action-angle coordinates (derived in App. A) can be used for obtaining the FP coefficients for 3D motion of a single ion at any amplitude within quadrupole traps. The derived expressions can be applied to the linearized oscillations of a crystal of trapped ions about their minimum positions rfions . Closed form expressions are derived for Doppler cooling in the linear limit (App. B). Some derivations of the 1D results presented within the main text of the paper are given in App. C.

## Ii Model

### ii.1 Motion in a quadrupole trap with excess micromotion

We start with the motion of an ion in a 1D time-dependent quadrupole potential, displaced by a constant electric field ,

 Ve=12(az−2qzcos2t)z2−Ezz. (1)

The units are nondimensional, obtained by rescaling the time by half the trap’s driving frequency, , and rescaling the coordinate by a natural lengthscale (discussed below), absorbing also the ion charge and mass into the nondimensional parameters, , , and in a standard way leibfried2003 ; rfions ; rfchaos ; rfcooling , detailed also in App. A. The solution of the inhomogeneous linear Mathieu equation of motion (e.o.m) derived from can be written as rfmodes

 z(t)=¯z(t)+u(t),¯z=∑nB2nei2nt, (2)

where and . The term , being -periodic (since the rf drive frequency is 2 in rescaled units) is known as excess micromotion as it can (typically) be minimized (but not easily eliminated), by using controlled electric fields (that make effectively small). Although it cannot be cooled away [in contrast to ], is completely coherent and not stochastic devoe1989role ; zigzagexperiment , and in this work we assume that it is invariant under time-reversal. Equation (2) defines a time-dependent canonical transformation to the new coordinate , with the conjugate momentum

 p=˙u, (3)

(with the mass equal to 1). The transformed Hamiltonian becomes

 H0(u,p,t)=12p2+VM.o.(u,t), (4)

with the nondimensional Mathieu oscillator potential

 VM.o.(u,t)=12(az−2qzcos2t)u2. (5)

We can now introduce the action-angle variables and , defined by a second (time-dependent) canonical transformation, using two functions of the phase space and time,

 I=Λ(u,p,t),θ=Θ(u,p,t). (6)

The action-angle variables constitute a very useful choice of variables, since is conserved during the Hamiltonian motion in the time-dependent potential, in contrast to the energy. The exact transformation functions are given in App. C.1. To the leading order in we can write

 u≈√2I/νzcosθ,p≈−√2Iνzsinθ+qzsin(2t)u, (7)

with

 θ(t)=νzt+ϕ, (8)

where is determined by the initial conditions, and the characteristic exponent of the Mathieu equation can be approximated by

 νz≈√az+q2z/2,az,q2z≪1, (9)

which determines the secular oscillation frequency in physical units, . To the same accuracy we have for the coefficients of of Eq. (2),

 B0≈Ez/νz,B2≈−B0qz/4, (10)

with the rest of the series truncated (App. C.1). We define the amplitude of excess micromotion to be the peak-to-peak oscillation due to , given to this order by

 Az=4B2=2qzEz/νz, (11)

which is linear in .

### ii.2 Laser cooling in the finite lifetime treatment

In this subsection we review for completeness the laser cooling approach employed in this work. We consider laser cooling within a two-level approximation for the ion and a semiclassical framework rfcooling . It is based on conservation of energy and momentum at each photon absorption event and each spontaneous emission occurring after a random delay due to the nonzero lifetime of the electronic excited level. Between the absorption and emission, we assume that the ion moves completely classically on an invariant torus of the Hamiltonian phase-space, and is decoupled from the electromagnetic field. The presented theory is valid in the limit of a low saturation of the transition, when the ion spends most of the time in its electronic ground-state, which is often chosen in practice for allowing to reach the lowest cooling limit.

The following parameters define the cooling; , the detuning of the laser from the electronic transition frequency (for an ion at rest), , the spontaneous emission rate (the linewidth of the transition), , the on-resonance Rabi-frequency (with its squared magnitude proportional to the laser intensity), , the laser wavenumber, and , the photon recoil momentum, with being Planck’s constant. We treat laser cooling by assuming a weak coupling to the electromagnetic field (low laser intensity), that is expressed using the saturation parameter proportional to the laser intensity,

 s=2(|ΩR|/Γ)2≪1. (12)

Then we assume that the probability of the ion to absorb a photon at a given point in phase space takes the form

 ρ(p,t)=s/21+(2Δeff/Γ)2, (13)

with the effective detuning due to the Doppler shift,

 Δeff=Δ−kvz(p,t), (14)

that depends on the real-space velocity which is a function of the canonical momentum (expanded about , and the time . For our 1D treatment the velocity is given [using Eq. (2)] by

 vz(t)=˙z(t)=˙¯z(t)+p(t). (15)

As derived in rfcooling , when the stochastic dynamics are slow in comparison with the Hamiltonian motion, integrating over the angle allows one to obtain an effective Fokker-Planck equation for the probability distribution

 ∂P(I,t)∂t=−∂S(I,t)∂I≡−∂∂I[ΠIP]+12∂2∂I2[ΠIIP], (16)

with the probability flux. Denoting with an overbar the torus average over any function of the phase space which is assumed to have some arbitrary period ;

 Ξ[I,θ,t+T]=Ξ[I,θ,t], (17)

we define

 ¯¯¯¯Ξ(I)≡1T∫T0dt12π∫Ξ[I,θ,t]dθ. (18)

The cooling coefficients entering Eq. (16) for the laser cooling setup considered here are

 ΠI(I)=Γ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ρ[pr∂Λ∂p+12p2r(∂2Λ∂p2+μ⟨∂2Λ∂p2⟩Γ)], (19) ΠII(I)=Γp2r¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ρ[(∂Λ∂p)2+μ⟨(∂Λ∂p)2⟩Γ], (20)

where comes from the second moment of the dipole radiation rfcooling and is often denoted by javanainen1981laser or leibfried2003 in the existing literature. All terms in the right hand side of the equations above are functions of the phase space time and point where the absorption occurred,

 Za≡{u,p,t}, (21)

that is averaged on a given torus by the definition in Eq. (18). The required derivatives of the action-angle transformation function can be found in App. A, and to the leading order in [Eq. (9)], coincide with the simple formulae for the harmonic oscillator,

 ∂Λ∂p≈−√2I/νzsinθ,∂2Λ∂p2≈1νz. (22)

We stress that the approximation in Eq. (22) gives accurate results for the excess micromotion (with a neglected contribution approximately equal to ). The emission phase space point is averaged over through the integration of the waiting time distribution . Given an absorption that occurred at the phase-space point , the mean value of any function of the phase space, at the time of emission, is given by

 ⟨Ξ(Za)⟩Γ≡∫∞0Γe−Γt′Ξ(Z(t+t′;Z(t)=Za))dt′. (23)

The time integral in Eq. (23) is to be performed along the trajectory, denoted with the notation of Eq. (21) as , as it starts at the phase space point and evolves according to the Hamiltonian motion at fixed . The zero lifetime limit can be obtained from the treatment above if in Eq. (23) it can be assumed rfcooling that . In addition, an adiabaticity condition is assumed that justifies the averaging of Eq. (18). A simple criterion is obtained by requiring a small relative change in action due to both drift and diffusion, during a cycle of the motion:

 ΠI/νz≪I,ΠII/νz≪I2. (24)

A steady-state of the FP equation with appropriate boundary conditions can be obtained by setting the left-hand-side to 0, and then noting that the reflecting boundary condition at the origin [] implies that this would be a zero-current state []. Integrating the resulting equation gives the time-independent distribution

 P(I)=N[ΠII(I)]−1exp{2∫IΠI(I′)ΠII(I′)dI′}, (25)

with the normalization factor. This steady-state solution is relevant if the exponential decays fast enough (as a function of ), implying physically that indeed the ion remains trapped for a long time in a bounded region of phase space.

For the numerical calculations presented in the following we take a ion. As a specific example, the Mathieu oscillator potential of Eq. (1) can be realized at the center of a five-wire trap rfchaos , in the direction perpendicular to the electrodes’ plane. The rescaling parameters taken here are (the rf electrodes’ width) and . With a voltage of on the rf electrodes and an axial trapping frequency MHz, the Mathieu parameters are , , , and . The dimensional laser parameters that we consider are

 ~k≈2π/313nm−1,~Γ≈2π×19MHz, (26)

and the spontaneous emission coefficient is

 μ=2/5. (27)

## Iii Thermal and nonthermal distributions in the final stage

### iii.1 The cooling limit without excess micromotion

In the absence of excess micromotion, the ion velocities can be assumed to be small near the origin of the quadrupole trap. This defines a specific limit of the cooling, that can in fact be analyzed analytically for the general 3D motion as in App. B. In this limit the Lorentzian can be linearized in the velocity, subject to the condition

 I≪Ilinear=(Γ2+4Δ28kΔ)212νz. (28)

In this approximation, we can write for the action drift and diffusion coefficients (with ), using the expressions given in App. B and using Eq. (59),

 ΠlI=γI+12iz,ΠlII=izI, (29)

with

 iz=prFr(1+μ)cz, (30)

where the mean radiation force and the rate of momentum damping having the well-known forms (javanainen1980 ; leibfried2003 ; rfcooling , respectively,

 Fr=prΓs/21+(2Δ/Γ)2,γ=4kprsΔ/Γ[1+(2Δ/Γ)2]2, (31)

generalizes a similar coefficient for the harmonic oscillator rfcooling , and is defined in Eq. (110), while to the leading order in , using Eqs. (9) and (22), we have

 cz≈ν−1z. (32)

Then for the final stage of the cooling reduces to an exponential, thermal equilibrium-like distribution in the conserved action, that is nearly identical to the well-known Doppler cooling limit that is obtained in limit of a vanishing excited level lifetime within a static harmonic potential (the heavy particle javanainen1980 ). The equilibrium is determined by the balance of momentum dissipation and diffusive heating. It is a zero-current distribution in the action, given by

 P(I)=λze−λzI, (33)

where , which gives the mean and standard deviation of the action, , can be obtained from Eq. (92) and reads here

 Ilimit≡λ−1z=ℏΓ8cz(1+μ)[Γ2|Δ|+2|Δ|Γ]. (34)

Although the Doppler cooling limit leibfried2003 is derived by assuming an ion in a time-independent potential and being stationary in space during a photon absorption-emission cycle, the final distribution derived here, for motion with micromotion, turns out to be thermal as well, and coincides with that of the harmonic oscillator of the same secular frequency (up to the small correction due to not being exactly ). The fact that the ion can emit the photon at a random time along the trajectory after the absorption, and hence with the ion’s a momentum (or kinetic energy) being very different than that at the absorption, does not lead on average to an increase in the action diffusion due to the scattering events. The analysis in App. B indicates that this is a consequence of the memoryless decay process together with the linearity of the oscillations, and the linearity of the Lorentzian (in the small velocity limit). Under such conditions, the phase-space averaging of such scattering events (with the accessible ion momenta), cancels out the effects of micromotion (to leading order). It serves to show the strength of the action-angle as the right choice of coordinates for analyzing the motion in the time-dependent potential.

### iii.2 The cooling limit with excess micromotion

However, with fixed and increasing, the validity of this statement soon breaks, and the distribution becomes both nonthermal and with higher values of mean and variance. Figure 1 shows the mean action in the steady state of cooling, as a function of the excess micromotion and the detuning. extends to in this study, which corresponds for the presented parameters to nm when rescaled by the natural length m, obtained for V/m.

We define the minimal mean action that is obtained in the upper left corner of the figure,

 Imin≡⟨I⟩(Az=0,Δ=−Γ/2), (35)

and using Eq. (34), we have to an accuracy of order , . For increased values of (with fixed), the value of grows slowly at first; remains for nm (V/m), and only then rises sharply, and the distribution becomes nonthermal (as will be discussed in the following). We show [App. C.2] that indeed the linear term in the expansion of as a function of (at any fixed ) vanishes [for the model of Eq. (1), with a static stray electric field]. In the region with

 Alinear(Δ)=[Γ2+4Δ2]/(8k|Δ|), (36)

the effect of the micromotion is small, and grows approximately inversely as a function of . This curve, , is indicated by a dotted line in Fig. 1(a).

Notably, for detunings , we find that first decreases with , reaching a minimal value before growing sharply. This is a result of the coefficient at order in the expansion of being negative in that region, as shown in App. C.2. In Sec. III.3 we present simple (though implicit) expressions allowing to calculate , plotted as a solid line in Fig. 1(a), along which a thermal distribution in the action can be obtained for any amplitude of excess micromotion. can be seen to be a nearly straight line (curving only close to ), and the final action depends only weakly on ; for the presented parameters (up to nm, corresponding to V/m).

Thus we find that can be used as a control parameter to counteract, to a high extent, the effect of excess micromotion. As shown in rfcooling , increasing is also efficient for the cooling of an ion from the high amplitude motion regime of approximately integrable motion in a surface-electrode trap, wherein a low-detuning laser (with ) may more easily heat the ion past the trap barrier. However, it should be noted that a laser beam with a large detuning could actually capture the ion in a large amplitude motion away from the trap centre rfcycles .

### iii.3 Optimal cooling with excess micromotion

A precise calculation of the cooling parameters and the action distribution requires using the approach presented in Sec. II.2. We can however obtain a simpler approximation allowing to calculate also the function , by expanding the Loretzian in [see App. C.2] and maintaining the periodic excess micromotion part,

 prΓρ≈Fe(t)+γe(t)p, (37)

with

 Fe(t)=prΓs/21+(2Δe(t)/Γ)2,γe(t)=4kprsΔe(t)/Γ[1+(2Δe(t)/Γ)2]2, (38)

and

 Δe(t)=Δ−k˙¯z(t). (39)

In the absence of micromotion, the above expansion reduces to the coefficients in Eq. (31). Let us assume that and can be averaged in time independently of the torus averaging and of the spontaneous decay waiting time. We define

 ¯Fe=1π∫π0Fe(t)dt,¯γe=1π∫π0γe(t)dt. (40)

Then the resulting distribution will take again an exponential form

 P(I)=λee−λeI,λe=−2¯γe/he (41)

where

 he=pr¯Fe(1+μ)cz. (42)

The mean final action has a minimum when takes a maximum, and this condition defines the curve of optimal detuning as a function of the excess micromotion amplitude, which can be found by integrating the coefficients in Eq. (40) using of Eq. (2), whose time derivative can be approximated as

 ˙¯z(t)=−∑n>04nB2nsin(2nt)≈−Azsin(2t), (43)

with defined in Eq. (11). This approximation of averaging and independently of the torus averaging, amounts for our parameters to a small underestimate () of the final action and of the corresponding optimal detuning curve , shown in Fig. 1. A detailed study of the parameter dependence of this approximation is beyond the scope of the current work.

## Iv Summary and Outlook

The main result of the current work has two aspects; a practical one and a more conceptual one. On the practical side, our main result is that excess micromotion of a large amplitude does not prevent reaching a low, thermal-like distribution for the fluctuations expanded about the coherent driven motion, achievable by simply tuning the cooling laser frequency. This is useful in particular when considered in the context of ion crystals discussed below. This concrete result also demonstrates a broad conclusion that can be drawn, which is that even though the ion is strongly driven by the trap potential in combination with the stochastic photon scattering processes, by choosing the right frame, the driven and the stochastic dynamics can be separated, with each becoming considerably simpler. The action-angle coordinates form the right coordinates because the fluctuations in action represent the relevant stochastic, noncoherent part of the motion, and moreover the angles can be safely coarse grained, leaving a clear picture of the underlying complex dynamics in terms of the actions.

An important foreseeable application of the presented results is to a crystal of ions. In a crystal configuration for which some ions are not positioned at a point where the rf potential vanishes (e.g. in a planar configuration in a linear Paul trap), those ions will perform driven periodic motion akin to excess micromotion, determined by the interplay of the Coulomb interaction and the periodic drive rfions . It is the dynamic equivalent of an ion’s equilibrium position in a static crystal, and cannot be removed. By tailoring the laser detuning, the final action distribution of these ions can be significantly reduced. A setup consisting of a few lasers can be treated by adding the drift and diffusion coefficients calculated separately for each laser. Spatially inhomogeneous laser profiles and laser parameters which are modulated in time can be transparently treated. In combination with noise heating rfcooling , a detailed characterization of the final distribution of laser-cooled ions and the dynamics leading to it can be obtained, and ideas for its manipulation can be explored chupeau2018engineered . The ensuing dynamics can be studied using the analytic tools of App. A, by expansion in small linearized deviations about the driven periodic motion. This way, cooling and heating dynamics corresponding to Gaussian white noise, and the stationary distribution of chains of ions in 1D and crystals in 2D and 3D configurations, can be studied PhysRevA.64.063407 ; morigi2003ion ; morigi2001twospecies ; fogarty2016optomechanical ; laupretre2018controlling ; kamsap2017experimental ; PhysRevLett.119.043001 ; 1367-2630-13-4-043019 ; PhysRevX.8.021028 ; Mitchell13111998 ; Drewsen_Long_Range_Order ; SchilerProteinsPRL ; tabor2012suitability ; mavadia2013control ; bohnet2016quantum .

Interesting extensions of the theory could include more general electronic level structure PhysRevA.96.012519 ; janacek2018effect , and applying the action-angle framework to setups where power-law distributions in energy (in an averaged sense) were predicted for collisions of ions with neutral atoms devoe2009power ; rouse2017superstatistical ; rouse2018energy ; meir2018direct . The interplay of micromotion, noise and laser cooling is of significant importance for applications in quantum information processing and the operation of quantum gates and entanglement operations with trapped ions PhysRevLett.81.3631 ; berkeland1998minimization ; PhysRevLett.101.260504 ; rfions ; zigzagexperiment ; Landa2014 ; PhysRevA.90.022332 ; wang2015quantum ; arnold2015prospects ; keller2015precise ; yan2016exploring ; mielenz2016arrays ; bruzewicz2016scalable ; keller2017optical ; keller2018controlling ; welzel2018spin ; delehaye2018single ; PhysRevA.97.062325 . In particular, the actions standing at the heart of the current work correspond exactly to the quantum mechanical phonons with a periodically-driven harmonic potential Glauber1992 , in terms of Floquet-Lyapunov modes rfions ; rfmodes . Using the exact time-dependent wavefunctions in quadrupole traps rfmodes ; mihalcea2017study , would allow to extend the theory to the quantum limit.

###### Acknowledgements.
H.L. thanks Ananyo Maitra, Dietrich Leibfried, Denis Ullmo, and Roni Geffen for fruitful discussions, and acknowledges support by IRS-IQUPS of Université Paris-Saclay, and by LabEx PALM under grant number ANR-10-LABX-0039-PALM.

## Appendix A A coupled system of Mathieu oscillators

We consider an ion that has been cooled to the center of a Paul trap of a general type. We assume that the potential can be approximated as a quadrupole (this is however not necessarily the case in multipole traps MartinaRing ). With and being the vector coordinate and velocity of an ion in dimensions, we first rescale the time by half the micromotion frequency, and the coordinates by a natural unit of length, , relevant for the trap at hand,

 →r→→r/w,t→Ωt/2,→v→→v/(wΩ/2). (44)

Using the rescaling in Eq. (44) allows us also to define a nondimensional momentum using the ion mass and by absorbing it and the ion charge into the parameters, we define a nondimensional potential energy resulting from an electrostatic voltage ,

 V→V/[mw2Ω2/(4e)], (45)

a nondimensional electric field,

 E→E/[mwΩ2/(4e)], (46)

and a nondimensional Planck constant scaled according to

 ℏ→ℏ/(mw2Ω/2). (47)

In general, the quadrupole potential may be composed of a sum of a few quadrupole potential terms whose origin does not coincide, or there may be “stray” electric fields that push the ion from the origin of the quadrupole potential. The resulting motion contains a component known as “excess micromotion”, since it can (typically) be minimized by using additional DC electric fields. In this case it is useful to describe the motion by using generalized coordinates with conjugate momenta , differing from the real space position and velocity (with ) by a time-dependent displacement rfions , that eliminates the terms linear in the coordinates from the potential energy,

 →r(t)=¯r0(t)+→u(t),→v(t)=˙¯r0(t)+→p(t). (48)

We assume that the nondimensional potential is time-reversal invariant and -periodic (which can include the particular case where it is time-independent), and can be expanded into a system of coupled Mathieu oscillators rfions , obtaining the nondimensional Hamiltonian

 H0(→u,→p,t)=12(→p)2+VM.o.(→u,t), (49)

with

 VM.o.(→u,t)=12→ut(A−2Qcos2t)→u (50)

where denotes the transpose, and , are matrices that describe the linearized DC and rf voltages tensors. By the Laplace equation, . The equations of motion derived from the linearized potential form a coupled system of parametric oscillators rfmodes ; shaikh2012stability . If and commute, then they can be diagonalized to give a system of decoupled Mathieu equations. In the opposite case, no such simple transformation exists, and the three spatial directions will be mixed by the micromotion.

The most general solution of this motion is a sum over decoupled linear oscillators

 →u=∑j√Ij(2Re∑n→Cj2nei2nteiθj), (51) →p=−∑j√Ij(2Im∑n→Cj2n(2n+νj)ei2nteiθj), (52)

where the summation extends over and the coefficients for give the micromotion modulation. In Eqs. (51)-(52) we have defined

 θj=νjt+ϕj, (53)

with the characteristic exponents of the Mathieu system, that are related to the secular frequencies of motion in the trap (, in physical units), by

 ωj=νjΩ/2, (54)

and and are (for now) arbitrary constants related to the initial conditions. We assume that the motion is stable, i.e. that are all real, and are between 0 and 1.

As detailed in rfmodes , a time-dependent (Floquet-Lyapunov) linear transformation can be used to transform the real space coordinates and momenta, to new (complexified) coordinates and the canonically conjugate momenta ,

 (55)

where and , are complex matrices constructed from the -periodic part of the column vector solutions given in Eqs. (51)-(52), i.e.

 U=(∑→Cj2nei2nt...),V=(i∑(2n+νj)→Cj2nei2nt...), (56)

and is the complex conjugate matrix, and we will use for the hermitian conjugate matrix, and also that is the transposed matrix. These matrices are to be rescaled by multiplication by a constant diagonal matrix,

 U→U(−2iVt(0)U(0))−1/2, (57)

and similarly, guaranteeing the normalization

 Vt(0)U(0)=12i, (58)

which also implies the identity (that will be used later),

 Ut(t)V∗(t)−Vt(t)U∗(t)=−i. (59)

Then the transformation in Eq. (55) is canonical and its inverse is

 (60)

The Hamiltonian then transforms according to

 H0(→ξ,→χ)=12∑jνj(ξjχj+χjξj). (61)

The time evolution of the new canonical coordinates is given by (noting that )

 ξj(t)=√Ijei(νjt+ϕj),χj(t)=√Ije−i(νjt+ϕj). (62)

We can now introduce naturally a second canonical transformation to the action-angle coordinates ,

 Ij=ξjχj,θj=ilnχj, (63)

which justifies the notation in Eq. (53), since

 ξj=√Ijeiθj,χj=√Ije−iθj, (64)

whence the Hamiltonian obtains the simple form

 H0(→I)=∑jνjIj. (65)

The relations in Eq. (63) define implicitly the action-angle transformation from the real-space coordinates,

 Ij=Λj(→u,→p,t), θj=Θj(→u,→p,t), (66)

with the transformation functions and depending explicitly on time (and being -periodic), and can be constructed explicitly using Eq. (63) and Eq. (60). In the following we will not need the functions and explicitly, however we will use the partial derivatives

 ∂Λj∂pα=iUαjξj−iU∗αjχj,∂2Λj∂pα∂pβ=UαjU∗βj+U∗αjUβj. (67)

We note here that the entire derivation above would remain completely valid if Eq. (50) is to be replaced by a more general time-reversal invariant and -periodic linear system, i.e. a system of coupled Hill equations McLachlan ; Yakubovich ; rfions ; rfmodes , with higher harmonics of the fundamental frequency . In addition, we note that setting , the solutions for , reduce to coupled harmonic oscillators with nondimensional frequencies . This is the pseudopotential approximation (within the harmonic approximation), with the motion described by the potential

 Vh.o.(→u)=12→utN→u, (68)

with a time-independent coupling matrix.

## Appendix B The linear limit of cooling for a coupled Mathieu system

As derived in rfcooling , a Fokker-Planck equation in actions (equal to the space dimension) can be written using the probability flux vector whose components are given by

 Sj(→I,t)≡ΠjP−12∑k∂∂Ik[ΠjkP], (69)

with the FP equation taking the form

 ∂P(→I,t)∂t=−∑j∂Sj(→I,t)∂Ij= −∑j∂∂Ij[ΠjP]+12∑j,k∂2∂Ij∂Ik[ΠjkP]. (70)

For a multidimensional torus, generalizing Eq. (21), a phase-space point is defined by

 Za≡{→u,→p,t}, (71)

and defining the torus average over any function of the phase space, where is assumed to have an arbitrary period ,

 Ξ(→I,→θ,t+T)=Ξ(→I,→θ,t), (72)

we have

 ¯¯¯¯Ξ(→I)≡1T∫T0dt1(2π)D∫Ξ(→I,→θ,t)dD→θ, (73)

 Πj(→I)/νj≪Ij,Πjk(→I)/√νjνk≪IjIk. (74)

The mean drift and diffusion rates are given by

 Πj(→I)=Γ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ρ(Za)⟨δIj⟩,Πjk(→I)=Γ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ρ(Za)⟨δIjδIk⟩, (75)

with

 ⟨δIj⟩=pr∑α^kα∂Λj(Za)∂pα+12p2r∑α,β[^kα^kβ∂2Λj(Za)∂pα∂pβ+μαβ⟨∂2Λj(Za)∂pα∂pβ⟩Γ], (76)

and

 ⟨δIjδIk⟩=p2r∑α,β[^kα^kβ∂Λj∂pα∂Λk∂pβ+μαβ⟨∂Λj∂pα∂Λk∂pβ⟩Γ], (77)

where all terms on the r.h.s above are functions of . The absorption probability for takes the form

 ρ(→p,t)=s/21+(2Δeff/Γ)2,Δeff=Δ−→k⋅→v(→p,t), (78)

with showing that can be a function of the canonical momentum and time [Eq. (48)].

Following the treatment of App. A we consider motion described by the coupled Mathieu system of Eq. (50), and in this section we assume for simplicity that excess micromotion is negligible, i.e. we set

 ¯r0(t)=0. (79)

Under these conditions, that the velocities are small and that the potential is quadrupole, we can derive the following expressions for the linear limit of the cooling, which turns out to be identical within both the zero lifetime limit and the finite lifetime assumptions of the derivation.

Having assumed small velocities in this limit, we set in the rescaled units. Then, linearizing the Lorentzian in velocity we can write the coefficients in a well-known form javanainen1980 ; leibfried2003 ; rfcooling ,

 prΓρ(→p)≈Fr+γ∑β^kβpβ, (80)

with

 Fr=prΓs/21+(2Δ/Γ)2,γ=4kprsΔ/Γ[1+(2Δ/Γ)2]2. (81)

Here gives a mean radiation force (for ), and is a rate of damping. The condition for the validity of this linearization is

 →k⋅→v≪[Γ2+4Δ2]/(8|Δ|). (82)

Using Eqs. (55), (63) and (67), the tori averages in Eq. (75) would then contain bilinear terms of the oscillators and . Let us assume first that the ion can be assumed stationary in space between the absorption and emission, which can be expressed by substituting

 Γe−Γt′≈δ(t′) (83)

in of Eq. (23). We will see below that the result without this restriction is in fact identical. Assuming that are nondegenerate 111For an analysis of laser-induced correlations between two nearly-degenerate harmonic oscillator modes, see javanainen1980 ., the only surviving combination will involve , due to Eq. (64). We hence find that these terms become constants, or linear in the actions. Defining the coefficients

 cjαβ=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(UαjU∗βj+U∗αjUβj), (84) djαβ=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(iUαjV∗βj−iU∗αjVβj), (85)

allows us to put the FP coefficients in the form

 Πlj(→I)=gjIj+hj/2,Πljk(→I)=~hjIjδj,k, (86)

with

 gj=γ∑α,β^kα^kβdjαβ,hj=~hj=∑α,βDlαβcjαβ, (87)

and

 Dlαβ=prFr(^kα^kβ+μαβ). (88)

We note that the fact that the same coefficient appears in both and after the averaging (where they result from averaging different terms), is a consequence of the linearity of the oscillations, and is important in the following. In addition, for an isotropic laser, i.e. if , then Eq. (59) implies that .

For a zero current state, with the probability flux vector defined in Eq. (69), we require

 ΠljP=12∂∂Ij(ΠljjP). (89)

Substituting the exponential ansatz

 P=(∏jλj)exp{−∑jλjIj}, (90)

the zero-current condition becomes

 gjIj+12hj=12[~hj−λj~hjIj], (91)

and here we see that in order for the constant term to cancel on both sides of the equation, we must have , i.e. the same coefficient in both the drift and the diffusion. The distribution is then solved by

 λj=−2gj/hj. (92)

This is a thermal-like, equilibrium distribution with mean action

 ⟨Ij⟩=1/λj, (93)

provided that . This requires and also that none of the normal modes decouples from the laser.

Relaxing the instantaneous decay assumption of Eq. (83), we now carry out the integration over the waiting time (exponential) distribution for decay from the excited level. The tori averages will result in terms which can be written in the form

 (94)

In the above integral, the order of integration does not matter, and hence the integration over can be performed last. In that case, the expression in Eq. (94) reduces, by the arguments leading to Eq. (84), to

 (95)

and similarly,

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨∂2Λj∂pα∂pβ⟩Γ=cjαβ. (96)

This change of integration order does not hold (a-priori) if higher order terms in the Lorentzian expansion of Eq. (80) have to be included, because they will appear outside of the integration (the term linear in the momentum integrates to 0 in any case).

Hence the final stage of the cooling is described by an equilibrium distribution in the action coordinates, even when the ion cannot be treated as frozen between the absorption and emission, and even if it is driven by the high frequency rf trap to a large kinetic energy. Sufficient conditions for the validity of this approximation are the linearization [Eq. (82)] and the nondegeneracy of the modes, in addition to the conditions of the derivation of the finite lifetime particle limit, presented in rfcooling .

## Appendix C Cooling dynamics in 1D

### c.1 Hamiltonian dynamics

Using the general derivation of App. A, the solution of the linear and homogeneous, Mathieu oscillator equation of motion derived from Eq. (4) is given by

 u=√I(U(t)eiθ+U(t)∗e−iθ), (97)

and

 p=√I(V(t)eiθ+V(t)∗e−iθ). (98)

The functions and are both -periodic functions,

 U(t)=∑nC2nei2nt,V(t)=i∑nC2n(2n+νz)ei2nt, (99)

and the normalization of the coefficients is given by

 V(0)U(0)=12i. (100)

The partial derivatives of the action angle transformation of Eq. (6) are given by

 ∂Λ∂p=√I(iU(t)eiθ−iU(t)∗e−iθ),∂2Λ∂2p=2|U(t)|2. (101)

For the particular solution of the inhomogenous equation derived from of Eq. (1), we can substitute of Eq. (2), obtaining the recursion relations

 [B2n(az−4n2)−qz(B2n−2+B2n+2)]=Ezδn,0. (102)

To get a continued fraction expansion we use the assumed time-reversal invariance of the solution, , and define , getting , so

 B0=Ez/(az−2qzc0), (103)
 c2n−2=qz/(az−4n2−qzc2n),n≥1. (104)

To the leading order [Eq. (9)] we have for the coefficients of ,

 B0≈Ez/νz,B2≈−B0qz/4,B4≈0. (105)

### c.2 Excess Micromotion

In the absence of micromotion, the final stage of the cooling in the Mathieu oscillator potential has been studied in rfcooling and here, above. We now generalize this treatment to include the excess micromotion (focusing on 1D motion).

We assume that the Lorentzian can be expanded in the velocity, subject to conditions of validity that will be elucidated later. We carry the expansion to third order in , and then we set in the rescaled units (with the mass set to 1), and assume in addition that remains small, obtaining

 prΓρ≈Fr+γvz+κv2≈Fr+γ˙¯z+κ˙¯z2+γp, (106)

with and defined in Eq. (31), and the nonlinear damping coefficient is

 κ=2k2prs(1−12(Δ/Γ)2)/Γ[1+(2Δ/Γ)2]3. (107)

Plugging Eq. (106) and Eq. (101) into Eqs. (19)-(20), using Eq. (59) and performing the averaging of Eq. (73) by exploiting the fact that the order of the integration of the exponential decay of Eq. (23) can be interchanged with the averaging in this limit of linearization [App. B], we get

 ΠI=γI+jz/2,ΠII=jzI, (108)

where the coefficient of excess micromotion is

 jz=pr(1+μ)(Frcz+γdz+κez), (109)

with

 cz=