Modes of Oscillation in Radiofrequency Paul Traps

We examine the time-dependent dynamics of ion crystals in radiofrequency traps. The problem of stable trapping of general three-dimensional crystals is considered and the validity of the pseudopotential approximation is discussed. We derive analytically the micromotion amplitude of the ions, rigorously proving well-known experimental observations. We use a method of infinite determinants to find the modes which diagonalize the linearized time-dependent dynamical problem. This allows obtaining explicitly the ('Floquet-Lyapunov') transformation to coordinates of decoupled linear oscillators. We demonstrate the utility of the method by analyzing the modes of a small `peculiar' crystal in a linear Paul trap. The calculations can be readily generalized to multispecies ion crystals in general multipole traps, and time-dependent quantum wavefunctions of ion oscillations in such traps can be obtained.


I. INTRODUCTION
The trapping of charged particles by electromagnetic fields is an essential tool for many investigations within physics [1,2] and chemistry [3,4], as well as for studies of bio-molecules [5,6]. In the Paul trap, charged particles are confined by oscillating multipole radiofrequency (rf) fields [3]. Different types of Paul traps have been proposed and constructed over the years, and an increasing number of experimental and theoretical work is dedicated to the improvement of these devices (see, e.g. [7]). With the advent of laser-cooling techniques [8], a lot of effort is focused on the study of crystalline forms of trapped ions . Single trapped ions have been cooled to the quantum ground state of motion [38,39], and crystals of a few ions have been cooled to the ground state in at least a few of the motional modes [40]. Long-range order has been observed with large structures of thousands of ions in Penning [21,22] and Paul traps [37,41].
In this contribution we analyze the dynamics of crystals of charged particles confined in rf traps. In Sec. II we review some previous results on the trapping of ions in rf traps. In Sec. III we discuss the linear stability of ion crystals in general, and consider limits of validity of the pseudopotential (time-independent) approximation, in relation to symmetries of the trapped crystal. We derive the micromotion amplitude of ions in a general periodic solution in a Paul trap, and relate our findings to recent experimental results. We then expand the motion of the trapped ions in coordinates of small displacements about the periodic solution, which is the dynamic equivalent of a minimum of the trapping potential. This expansion is completely general and can be applied to any periodic trapping field. In Sec. IV we solve the coupled motion of the ions in the time-dependent potential. We find the modes which diagonalize the dynamical problem and obtain explicitly a time-dependent transformation to coordinates in which the motion is that of decoupled linear oscillators. Using this expansion the exact time-dependent wavefunctions of ions in rf traps can be obtained. We present a numerical study of a small crystal in a linear Paul trap in Sec. V, and conclude in Sec. VI with comments on further possible applications and research.

II. OVERVIEW OF PAUL TRAPPING
The first crystallization of a cloud of charged aluminum microparticles has been reported in [42]. Wigner crystals of 2 to approximately 100 trapped ions in a Paul trap were reported in [9,10] and further investigated in [43] both experimentally and numerically. The simulations included the rf trapping, Coulomb interaction, laser cooling and random noise. Depending on trap parameters, the ions have been found to equilibrate either as an apparently chaotic cloud or in an ordered structure. The latter is defined as the 'crystal' solution when it is a simple limit cycle, with the ions oscillating at the rf frequency about well-defined average points. The transition between the two phases has been investigated, it was shown that both phases can coexist and hysteresis in the transition has been observed [43].
The motion of two ions in the Paul trap has been investigated in detail in various publications [44][45][46][47][48][49][50][51][52][53]. In addition to the aforementioned phases, frequency-locked periodic attractors (where the nonlinearity pulls the motional frequencies into integral fractions of the external rf frequency) were found in numerical simulations and experiments. These solutions are different from the crystal in that the ions move in extended (closed) orbits in the trap, whose period is a given multiple of the rf period. However many of these frequency locked solutions are unstable, especially those of a large period, and perturbations such as those coming from the nonlinearity of the laser-cooling mechanism, tend to destroy them.
The reason that such crystals were termed 'peculiar', is that for a two-ion crystal with no angular momentum about the axial x-axis, there are only two possibilities in the harmonic pseudopotential approximation (excluding the case of degenerate secular frequencies). The crystal can either align with the x-axis or lie in the y-z plane. In the latter case we may assume that the crystal is aligned with the y-axis, and in both cases the z coordinate can be eliminated. Therefore a crystal which is at angle to both reduced axes x and y, i.e. one which forms an angle with both the axial axis and the plane of symmetry, is peculiar.
In [56,57], an improved pseudopotential approximation was derived for two ions in the hyperbolic Paul trap (and in [59] for the linear Paul trap). This nonlinear pseudopotential can be used to derive approximately the areas of stability in parameter space of the axial and radial two-ion crystal, and the orientation of the peculiar crystal. However, even the improved nonlinear pseudopotential cannot reproduce the areas of no stable crystal.
As mentioned above, the crystal solution in all of these works, is taken to be a periodic solution of the equations of motion with a period equal to that of the rf potential. This is a limit cycle of the equations, the simplest of the frequency locked solutions mentioned above, and it is an attractive one when cooling is present and it is stable. Its stability is analyzed in [58] by looking at the Poincaré map associated with this solution. The Poincaré map is a mapping of the phase space into itself, in which every inital point {x (0) , p (0)} is evolved according to the equations of motion and mapped to {x (T ) , p (T )} after one rf period T . The crystal is by definition a fixed point of the Poincaré map, and its linear stability is determined by the linearized dynamics in its vicinity, specifically by the eigenvalues of the so-called monodromy matrix. The linearization takes into account the leading (harmonic) coupling between the ions, expanded about the periodic solution. Further analysis (both linear and nonlinear) of the periodic orbit of two ions in a Paul trap appears in [60][61][62][63].

A. The Trapping Potential
We start with the potential energy of identical ions in the most general single, nonsegmented quadrupole rf trap. By choosing a specific frequencyω, which is characteristic of the secular frequencies in the problem, and measuring time and distances in units of 1/ω and d = (e 2 /mω 2 ) 1/3 , respectively (with m and e the ion's mass and charge), we write the nondimensional potential as where R i = {x i , y i , z i } is the vector coordinate of ion i, the trapping terms are given by Λ α = Ω 2 4 (a α − 2q α cos Ωt), α ∈ {x, y, z}, with a α and q α being the nondimensional Mathieu parameters of the respective coordinates [64], and Ω being the rf frequency (in units ofω).
Regarding the trapping parameters, the Laplace equation implies that Two commonly used types of traps are the hyperbolic [65] and the linear [64] Paul traps.
Taking the x-axis as the axial direction, the hyperbolic trap can be described by setting a y = a z ≡ a, a x = −2a and q y = q z ≡ q, q x = −2q, while in the linear trap we have a y = a z ≡ a, a x = −2a and q y = −q z ≡ q, q x = 0 (so a must be negative to obtain stable trapping). In the next subsection we will relate to symmetries of these two trapping geometries, but for now we keep the discussion completely general. The trapping potential V trap when considered alone, gives rise to decoupled Mathieu equations in each ion coordinate [66]. The corresponding characteristic exponents are β α (a α , q α ) and the secular frequencies are then ω α ≡ β α Ω 2 . In these units, V trap and V Coulomb are both of order unity for a crystal, when the trapping balances the Coulomb repulsion. These units allow naturally the introduction of the parameter Ifω is the (dimensional) secular frequency along some specific axis α, we have ω α = 1 and ǫ = β 2 α . Using the familiar lowest order approximation [64], β α ≈ a α + q 2 α /2, we see that the limit ǫ → 0 is equivalent a → 0, q 2 → 0. This parameter will be used in the following.
When the Mathieu parameters belong to a single-particle stability-zone, the potential of eq. (1) allows for bounded motion and large ion clouds can be trapped for extraordinary long times without cooling. This is given by the limit R i → ∞ in eq. (1), which means that ions are not crystallized, rather their motion is that of decoupled trapped particles. Despite the large amplitude motion, at low enough density the ions hardly interact.
However, stable single-particle parameters do not guarantee that a stable crystal solution exists, even with two ions, as discussed extensively in Sec. II, which has also been recently investigated with crystals consisting of upto several hundreds of ions [67,68]. We recall that a stable crystal solution is considered as a periodic solution with the same period of the rf, which is linearly stable under small perturbations. The existence of a stable crystal solution is a property of the fully nonlinear and time-dependent problem. There may also exist stable multi-periodic solutions (with a period which is some large, arbitrary multiple of the rf period), the frequency-locked solutions described in Sec. II.
In experiments, even at the presence of cooling, crystals may exist only metastably, changing completely after a given time or when parts of the crystal are moving with respect to the rest [37]. In this case, if the timescale of the change of the crystal is long enough, the analysis in the following subsections, using modes expanded about the (quasi-) periodic solution may still be useful, even though at least one mode would be unstable.

B. Linearization using the Pseudopotential Modes
In this subsection we linearize the crystal motion using the pseudopotential modes. Therefore this treatment assumes that a crystal solution exists, and is close to the crystal of the pseudopotential. As known from the example of two-ions in a hyperbolic Paul trap, this assumption may not hold, even for arbitrarily small values of a and q (where the crystal is 'peculiar' [56][57][58]). The linearized secular modes of a one-dimensional Coulomb chain were treated in detail in [69][70][71].
Let us define the squared ratio of secular radial and axial frequencies γ y = ω 2 y /ω 2 x , γ z = ω 2 z /ω 2 x and γ x = ω 2 x /ω 2 x = 1, thus choosingω as the axial trapping frequency, so we can rewrite the potential as the sum We expand about the minimum-configuration locations, R 0 i , obtained from the secular part of V in eq. (4), V pseudo + V Coulomb , and change to the normal modes Θ j by setting where D j i,α is the matrix of normal mode vectors, with rows indexed by the N ion numbers i and 3 directions α, and columns by the 3N normal mode numbers j. Writing the potential in terms of the normal modes, V = V harmonic + V rf + V nonlinear , we keep only the first two terms to get and the linearized equation of motion (e.o.m) derived from eq. (6) is Before treating the coupled system described by eq. (7), we first note that when the rf trapping potential is symmetric with respect to the axes of symmetry of the crystal, the equations of motion are in fact diagonal. For the following few paragraphs, let us therefore limit the discussion to the hyperbolic and the linear Paul traps, whose parameters were described following eq. (2). Given these two trapping geometries, the current expansion will be diagonal if the crystal is two-dimensional and lies in the plane of cylindrical symmtery of a hyperbolic trap, or the crystal is one-dimensional and aligned with one of the trap axes, in either the hyperbolic or the linear trap. In that case, we have R 0 i,α = 0 for the coordinates transverse to the crystal, and also the normal modes decouple in these directions from the modes tangential to the crystal (i.e. D j i,α is divided into blocks in the index α). We then can use two relations which hold in the plane of the crystal or along its axis, where hereafter,α runs on the symmetry directions , and b denotes the breathing mode. The identity on the left in eq. (8) is the completeness of the normal modes along the symmetry directions. To get the second identity we use the fact that in a harmonic trap, the breathing mode vector is exactly proportional to the minimumconfiguration locations R 0 i,α of the pseudopotential (see, e.g. App. B of [15]), so that and therefore the vector R 0 i,α is orthogonal to all other normal mode vectors. By using eq. (8) we can replace the r.h.s in eq. (7) with the simple expression We find that the equations of motion for modes in the crystal plane or along its axis are, after multiplying by ǫ = 4/Ω 2 of eq. (3) and rescaling t → Ωt/2, where γ, a and q are defined along the symmetry axes. Similar equations hold in the directions transverse to the crystal, without the inhomogenous r.h.s and with the corresponding γ, a and q.
Eq. (10) shows explicitly how the isotropy of the potential along the axes of symmtery of the crystal allows to decouple the equations of the modes. This decoupling puts the conditions for linear stability of the crystal in terms of diagonal Mathieu equations for each of the modes. In addition, with this symmetry the only mode with an inhomogenous r.h.s is the breathing mode, i.e. it is the only mode on which the rf potential acts as a driving force. This driven motion is in fact the π-periodic motion of the crystal as whole, about the static minimum-configuration locations of the pseudopotential, R 0 i . Returning to the case of a general trap, in order to put eq. (7) in a simple form we multiply it by ǫ = 4/Ω 2 of eq. (3) and rescale t → Ωt/2. Defining the two matrices and the vectors we rewrite eq. (7) in vector notation as We solve the homogeneous l.h.s of eq. (13) in Sec. IV. Since the pseudopotential modes are expanded aboout static configuration points, eq. (13) is an inhomogenous equation with a π-periodic r.h.s, which has a unique π-periodic solution (except for possibly a region of measure-zero in a α , q α space). This periodic solution of driven motion of the normal modes, corresponds to the exact π-periodic solution which defines the crystal in the rf trapping potential. Details of the solution of the inhomogeneous equation can be found in [72].
Eq. (13) shows that in general, the rf couples the pseudopotential normal modes and also acts as a driving force. Under this coupling, the true modes of oscillation of the system may in general have different frequencies (and even lose stability), and different oscillation directions than the pseudopotential modes upon which this expansion is based. Indeed, the linearization starting from the pseudopotential approximation may not be adequate in the general case. We investigate this point further in the following two subsections.

C. The Periodic Crystal Solution
We now abandon the pseudopotential approximation and turn to studying the timedependent potential directly. In this subsection we derive analytically the micromotion amplitude of the ions in typical crystals in Paul traps. The e.o.m for the ion coordinates, derived from eq. (1) after rescaling by t → Ωt/2, is Eq. (14) has π-periodic coefficients and is time-reversal invariant. We assume the existence of a crystal in the sense of Sec. II, i.e. a π-periodic and time-reversal invariant solution, which obtains the genreal form In this form, the average ion location is B 0,i (which is different from R π i (t = 0)). We wish now to see what can be said about R π i in typical Paul trapping experiments. We first define the dynamic matrix and write using eq. (15) its Taylor and Fourier expansion around B 0,i − B 0,j as Substituting the solution ansatz eq. (15) into eq. (14) we get in eq. (18) by its leading order, time-independent term from eq. (17), G 0,ij = G ij B 0,i , and require that the above relation holds for every t. We get, by using B 2n = B −2n and neglecting B 4n ≈ 0 (which is of the next order in q α ), the equation coming from the n = 0 term, and the coupled equation from the coefficient of e 2it , This sytem of equations can be seen as a linear homogenous system in the 2N variables B 0,i,α , B 2,i,α , i = 1, ..., N, for fixed α, since the equations are diagonal along the three different axes. Of course, this is not really a linear system since the matrix G 0,ij depends on B 0,i , but since we are not trying to actually solve the system, this will not matter. Let us fix the index α to some axis (supressing it in the following), and define the N-component vectors u 0,i = B 0,i,α and u 2,i = B 2,i,α . Then the above system can be written in block-matrix form where we have neglected a − ǫG 0 as compared with −4 in the lower-right block of eq. (21). Since we assume that the system has a solution (which approximates the π-periodic crystal), the above matrix must be singular. We can expand its determinant, Taking u 0 to be the vector from the kernel of (a − ǫG 0 + q 2 /2) which obeys (a − ǫG 0 ) u 0 = − (q 2 /2) u 0 , we find that the solution of eq. (21) is u 2 = − (q/4) u 0 . We therefore have obtained that in a general quadrupole trap, in the π-periodic crystal solution of eq. (15), every ion coordinate obeys (at least to leading order in a α /4, q α /4, ǫ α /4), i.e. that the micromotion amplitude in each coordinate is q α /2 of the respective average position, and at π phase with respect to the rf drive. In the hyperbolic trap the corresponding motion has been imaged as early as in [42]. In simulations of a generic trap (with different q and a parameters for the three axes), eq. (23) seems to hold accurately to within a few percent, for q up to ∼ 0.7, which is consistent with a deviation of order (q α /4) 2 .
The relation in eq. (23) loses its accuracy when either B 0,i,α ≪ 1 or q α ≪ 1. In the former case, for an ion near the zero of one of the rf axes, the corresponding micromotion amplitude in fact seems in the cases we have checked, to be lower (this is similar to a single trapped ion at the origin of a Paul trap, for which the unique π-periodic solution is the trivial one).
For the linear Paul trap case of q α ≪ 1 along the axial axis, the first-order expression in eq. (23) loses its meaning. We must therefore add to eq. (20) the second term in eq. (17), G 2,ij , which rotates at the frequency e 2it . Setting q = 0, eq. (21) is replaced with We now examine what can be deduced about G 2 . The off-diagonal elements are in general In particular, using eq. (23) with q y = −q z = q and q x = 0, we get that for the linear Paul trap Since the diagonal elements are the negative of the row sum, as in eq. (16), G 2,ii = − m =i G 2,im , we immediately find that for a crystal which is invariant (up to a permutation of the ions) under y ↔ ±z, G 2,ii = 0 + O (q 2 /4 2 ). In fact, for the linear Paul trap, the e.o.m, eq. (14) is invariant under y ↔ −z, t → t + π/2, so given one crystal solution, there is also a solution related by this transformation. Depending on the number of ions and trapping parameters, both solutions may actually be the same crystal. As the crystal size grows, by the same symmetry arguments, G 2,ii → 0 + O (q 2 /4 2 ).
The second row of eq. (24) gives B 2,i,x = −ǫ/4 G 2,ij B 0,j,x , and we argue that by the above symmtery arguments, for a typical symmetric or large crystal in a linear Paul trap, when the crystal configuration is also symmetric under x ↔ −x, the first order terms in the above summation will cancel, and (using ǫ ≈ q/2) If the symmetry alluded to above does not hold, eq. (27) should be replaced with an expression which is one order less, namely B 2,i,x = O ǫ 4 · q 4 B 0,i,x . Eq. (27) explains why in the linear Paul trap, there is essentially no micromotion excitation along the axial direction despite the strong Coulomb interaction. In addition, since q y = −q z in the linear trap, the oscillation described by eq. (23) is exactly the (2,2) mode of cold-fluid theory [14] which has been observed both in simulations [20] and recently discussed in connection with experiments [29]. In fact, in simulations we have performed (Sec. V), eq. (23) seems to hold extremely accurately radially (to within half a percent), and eq. (27) gives an accurate estimate for the axial micromotion amplitude (also matching e.g., [20,73]).

D. Linearization about the Periodic Crystal, and the Pseudopotential Limit
In this subsection we expand the modes of oscillation of the ions about the periodic crystal solution, eq. (15). The following expansion is applicable to general crystal configurations, does not rely on the pseudopotential modes or on a small parameter, and can be readily generalized to settings not considered in this paper, e.g. segmented traps and higher-order multipole traps. We will discuss its relation to the expansion of Sec. III B, which used the pseudopotential modes.
Expanding the potential of eq. (1) about the time-dependent solution R π i (t) , in the coordinates of small oscillations r i,α = R i,α − R π i,α (t), and defining the π-periodic matrix we find the linearized e.o.m This equation is a linearly coupled system for the ion coordinates with π-periodic coefficients. We see that the interaction term in eq. (29) is of the same order (in ǫ, or equivalently a and q 2 ) as the diagonal term. Expanding the matrix K στ ij (t) in a Fourier series in the form we first define the two matrices We can now switch notation to a simpler one with dynamical variables u m whose single index m stands for both indices {i, σ} of r i,σ , with the corresponding indices replaced in eq. (31), and rewrite eq. (29) in vector notation, with only the two leading harmonics of the Fourier series, as¨ u + [A − 2Q cos 2t] u = 0.
This equation describes linealrized coupled perturbations about the π-periodic solution which defines the crystal in the time-dependent potential. In Sec. IV we will solve this coupled system and find its decoupled modes of oscillation. Eq. (32) is seen to be identical in form to the l.h.s of eq. (13), which is given in terms of the pseudopotential modes. However, the current expansion is based on the exact force acting between the ions during their motion along the periodic trajectory of the full potential. Higher harmonics from the Fourier series of this force can be added, and in Sec. V we will in fact include the next harmonic (cos 4t), to obtain very accurate results for the modes. When the nonlinear motion is not exactly π-periodic, some of the linearized modes of the above expansion may not be oscillatory (as detailed in Sec. IV A), which describes (locally) the aperiodic motion of the crystal. If this change of the crystal is slow, the linearization of eq. (29) can still be useful.
Considering the limit of taking the Mathieu parameters to zero, this does not guarantee that the crystal will be (linearly) stable. If there are unstable modes, they may remain unstable even as a, q 2 , ǫ → 0. As for approaching the pseudopotential crystal in this limit, if the average ion locations B 0,i of the periodic solution tend to the pseudopotential minimum locations, then the limit of the pseudopotential modes is regained on the l.h.s of eq. (29). Characterizing the conditions for this to occur requires further investigations, however for simple cases the results of Sec. III C may provide some hints.
In particular, if the rf potential is isotropic with respect to the configuration (as discussed in Sec. III B for specific cases), and we neglect K 2 in eq. (31), the matrix Q becomes proportional to the identity. Then we can diagonalize eq. (32) by an orthogonal transformation (which is exactly the transformation to the normal modes) and obtain eq. (10). There will now not be an inhomogeneous r.h.s, since the periodic crystal motion is already accounted for. At the presence of driven breathing oscillations, the ions feel stronger repulsive nonlinearity when they oscillate radially inwards, so we may expect the real crystal to be a bit more expanded (with the ion distances somewhat larger) as compared with the pseudopotential approximation.
We note that for an axial chain of ions in the linear Paul trap, indeed K 2 ≈ 0 even if there is a small rf leak in axial direction, or when the ions do not lie exactly along the rf null axis, and eq. (10) is almost exact. For the axial modes, the correction to the pseudopotential modes will be small. For the radial modes, the correction will typically be small, provided that the rf frequency is large compared with the radial frequencies. However, for calculation of small effects which might be of importance in high-accuracy experiments such as for quantum information processing [64], the full time-dependent solution must be used.
As a final note, a generalization of the isotropy of the rf potential with respect to the crystal configuration, would be the existence of a constant orthogonal transformation which diagonalizes eq. (32), and thus decouples the modes in the non-isotropic case. Such a transformation would simultaneously diagonalize the matrices A and Q (which must commute), and does not exist in the general case.

A. The Floquet Problem
Eq. (32) is a linear differential equation with periodic coefficients and therefore amenable to treatment using Floquet theory, which we briefly review here. For more details see [72] and references within. For the Newtonian problem with f degrees of freedom (f = 3N for N ions in three-dimensions), the corresponding Floquet problem is stated in terms of coordinates in 2f -dimensional phase space by the definitions where 1 f is the f -dimensional identity matrix. The e.o.m is written in standard form aṡ In the following, an f -dimensional vector u will be denoted by a lower case Latin letter with an arrow. f -dimensional matrices will be denoted by capital Latin letters (Q). A 2f -dimensional vector φ will be denoted by a lower case Greek letter. Capital Greek letters (unitalicized) will denote 2f -dimensional matrices (Π, B).
A fundamental matrix solution to eq. (34) has 2f linearly independent column solutions and obeys the matrix equationΦ (t) = Π (t) Φ (t). A fundamental matrix solution which equals the identity matrix at t = 0, i.e. obeys Φ (0) = 1 2f , is known as the matrizant of eq. (34) (and is obviously unique). The matrizant can always be written in the form Φ (t) = Γ (t) e Bt Γ −1 (0) where Γ (t + T ) = Γ (t) is periodic with the period T of Π (t), and the constant matrix B is diagonal, with entries known as the characteristic exponents of the Floquet problem, The time-dependent linear coordinate change, known as the 'Floquet-Lyapunov' transformation, defined by transforms the e.o.m eq. (34) into the time-independent diagonal forṁ whose solutions are the Floquet modes, We note that in the general case, Γ mixes the coordinates u and their derivatives, and in addition is not unitary, although certain highly symmetric cases, e.g., as discussed in Secs. III B and III D, when the Mathieu equations decouple completely, are an exception.

B. Solution Using an Expansion in Infinite Continued Matrix Inversions
We expand the solutions of the e.o.m (32) by using an analytical expansion which utilizes infinite-continued matrix inversions [74], see [72] for details and further references. This expansion allows to obtain the frequencies and the coefficients of the solution vectors, to arbitrary precision. The solution is given in a form which is immediately suitable for obtaining the Floquet-Lyapunov transformation, as will be shown in the next subsection.
We seek for solutions of eq. (32), in the form of a sum of two linearly independent complex columns vectors, where b and c are complex constants determined by the initial conditions. Stable modes will be described by β taking a real nonintegral value (we exclude the case of integral β).
Following the discussion in §4.70 of [66], for trapping parameters in the first stability zone of the Mathieu equation, β can be chosen in the range 0 < β < 1 for all stable modes. By defining R 2n = A − (2n + β) 2 we can write infinite recursion relations for C 2n , and obtain two independent expansions in infinite continued matrix inversions and Multiplying eq. (41) by Q and defining we find that all characteristic exponents β are zeros of the determinant of Y 2,β (which is a function of β). If there are degenerate β's they will appear as degenerate zeros of this determinant. The vector C 0 for each β is an eigenvector of Y 2,β with eigenvalue 0. Since A and Q are symmetric, Y 2,β is symmetric as well, and so its kernel will be of dimension equal to the algebraic multiplicity of the β root. The vector C 2 can be obtained by an application of T 2,β Q to C 0 , for n = −1 we use C −2 = [T −2,β ] −1 Q C 0 , and so on for the other vectors. We note that the different vectors C 2n,β are not orthogonal in general, and the vectors at every order in n mix different coordinates. The general term of the expansion vanishes. Either A or Q may be singular and the expansion can still be applied in general. Even if both are singular, the expansion is valid if there are no integral values of β, a case which we do not tackle as noted above. Excluding perhaps isolated values of β (and atypically in the a−q parameter-space), all matrices which are inverted in the above expressions will be invertible, and while employing the algorithm in practice, the invertibility of the matrices is of course easily verified at each step. In Sec. V we use in fact a generalization of the above expansion [72] which includes also the next Fourier harmonic (cos 4t, ommited from eq. (30)).

C. The Floquet-Lyapunov Transformation For Stable Modes
We now further assume that all Floquet modes are stable, i.e. that the 2f linearlyindependent solutions of eq. (34) are oscillatory and thus come in complex conjugate pairs. This simplifies many expressions. We therefore take B of eq. (35) in the block form where β j are positive. We define the f -dimensional matrix U whose columns are constructed from the series of f -dimensional vectors C 2n,β j obtained from the recursion relations for the solutions of eq. (39), i.e.
where in the above expression and for the rest of this section, the summation is over n ∈ Z.
We similarly define the f -dimensional matrix V composed from column-vectors as The matrices U and V can be chosen to obey the normalization condition which is a rescaling imposed by multiplication with a (diagonal) matrix, such that and V accordingly. As shown in [72], the Floquet-Lyapunov transformation and its inverse can then be obtained in closed form, and are given block-wise by (where U * denotes the complex conjugate of the matrix U), This transformation is a canonical transformation from the Hamiltonian coordinates u and their conjugate momenta p =˙ u, to the variables given by where ξ is the new f -dimensional vector of coordinates and −iζ are the new conjugate momenta (we here break the notation a little). Using the realness of u and p it is easy to verify that ξ = ζ * . The time dependence of these modes is (with only a 1% DC asymmetry in the radial plane), such that the center-of-mass frequencies are degenerate to within 1%. The Mathieu parameters are q y = 0.41 and a x = 0.05766. In the pseudopotential approximation, the corresponding minimum-configuration is a (nearlyregular) octahedron, with two ions sitting on each axis of the trap. The rf crystal on the other hand, may well deserve to be called 'peculiar' (in the sense described in Sec. II), since it is not oriented with the axes. There are two ions lying along the y-axis, along which the confinement is strongest, but the two other pairs of ions sit along lines which are rotated in x−z plane. The large amplitude oscillation at the rf-period, and the absence of micromotion along the axial direction, confirm very accurately the results presented in Sec. III C.
Numerically, the periodic crystal solution can be obtained by starting from a simulation of the full e.o.m's with a friction (cooling) term and slowly turning it off (the adiabatic shutting-down of the damping term is important). The crystal is then followed for a period and the periodic solution and force matrix can be Fourier expanded to obtain the matrices A and Q. Since the crystal is 'peculiar', it has no corresponding pseudopotential limit (normal modes of regular polyhderons were investigated in [75], and see also references in [76]). For a very accurate description of the modes in the Paul trap, we add to eq. (32) the next Fourier harmonic to get¨ u + [A − 2Q 2 cos 2t − 2Q 4 cos 4t] u = 0, and expand it in the method of continued matrix inversions, using the formulas given in App. B of [72], which generalize the formulas presented in Sec. IV. This modification is required in order to obtain accurately the low frequency modes of nearly-degenerate configurations, as in a peculiar crystal.

VI. CONCLUDING COMMENTS
In this contribution we have investigated the dynamics of ion crystals in rf traps. We repeat the main results presented herein. In eq. (7) we show that in general, the rf couples the pseudopotential normal modes of the crystal, and also acts as a driving force. For a crystal configuration which has the same symmetry as the trapping potential (e.g. a chain of ions or a planar crystal in a trap with cylindrical symmetry), we find that the equations of motion for the modes become decoupled Mathieu equations as in eq. (10), for which the stability analysis is trivial (but may be different from the pseudopotential approximation).
In eqs. (23) and (27) we derive results which have been observed in experiments and numerical simulations, namely that in a general ion crystal the micromotion amplitude in each coordinate is q α /2 of the respective average position, and that in the linear Paul trap, axial micromotion is negligibly small.
When the crystal solution differs from the pseudopotential limit (as in the peculiar crystals discussed above), or when a very accurate expansion of the modes is desired, the derivation of eq. (32) takes into account the full rf crystal solution and ion interactions along the periodic trajectory. Sec. IV describes briefly how to solve for the decoupled modes of oscillations of the ions, allowing obtaining explicitly the mode frequencies and solution vectors.
We have focused on single-species crystals in quadrupole traps, specifically the linear and hyperbolic Paul traps. However, generalization to other cases is easy. Crystals in multipole traps [77][78][79] can be treated by expanding the motion around a suitable periodic solution R π i as in Sec. III C. Keeping only the leading terms will lead to the equations treated above. Segmented traps and trap arrays [80][81][82][83] can be handled simply by changing the Mathieu parameters felt by each ion (assuming that the rf drive has identical frequency). Crystals of nonidentical ions and other types of driving can be treated similarly, and various transformations [84] can be used to handle more general linear system similar to eq. (32), such as systems with first-order derivatives (e.g. linear damping, gyroscopic forces or magnetic fields [69,85]).
The framework presented here for calculating the classical normal modes and their frequencies can find immediate application in most studies involving trapped ion Coulomb crystals, including studies of how ion Coulomb crystals differ from uniformly charged liquid models [14], and possibly even crystalline beams [86,87]. The analysis is based on linearization of the nonlinear solution about a periodic solution, and the nonlinear correction terms can be written as a series expansion in the Floquet modes about the crystal solution. In a quadrupole trap the nonlinear terms would come only from the Coulomb interaction, and in a multipole trap, there will be terms coming from the nonlinear trapping potential. Such an expansion may serve as a starting point for study of nonlinear collective phenomena, e.g. the important phenomenon of rf heating (e.g. [20,73,88,89] and many more).
In addition, an exact quantum description of the modes, and wavefunctions in the configuration space of the ions, are presented in [72], utilizing the Floquet-Lyapunov transformation described in Sec. IV. The nontrivial time-dependent wavefunctions could, e.g., eventually become important for understanding trapped ion chemistry at ultracold temperatures [90].