The normal modes of a trapped ion crystal are derived using an approach based on the Hermitian properties of the system's dynamical matrix. This method is equivalent to the standard Bogoliubov method, but for classical systems, it is arguably simpler and more general in that canonical coordinates are not necessary. The theory is developed for stable, unstable, and neutrally stable systems. The method is then applied to ion crystals in a Penning trap. Reduced eigenvalue problems for the case of large applied magnetic fields are developed, for which the spectrum breaks into E × B drift modes, axial modes, and cyclotron modes. Thermal fluctuation levels in these modes are analyzed and shown to be consistent with the Bohr–van-Leeuwen theorem, provided that neutrally stable modes associated with crystal rotations are included in the analysis. An expression for the rotational inertia of the crystal is derived, and a magnetic contribution to this inertia, which dominates in large magnetic fields, is described. An unusual limit is discovered for the special case of spherically symmetric confinement, in which the rotational inertia does not exist and changes in angular momentum leave the rotation frequency unaffected.

## I. INTRODUCTION

This paper examines the normal modes of an ion crystal confined in the static electric and magnetic fields of a Penning trap. Such ion crystals, consisting of anywhere from a few to thousands of ions, are employed in a variety of applications ranging from fundamental studies of quantum entanglement, quantum simulation, and frequency standards^{1–6} to studies of the properties of strongly coupled plasmas.^{7–10} Linear normal modes of oscillation in these crystals are used as a diagnostic and manipulation tool in some of these studies, and a detailed understanding of the modes is essential to the success of this work.

In this paper, we lay out a general theory for the normal modes. In previous work, normal modes of a periodic correlated Coulomb lattice in a uniform magnetic field were found using Fourier methods,^{11–15} taking advantage of the periodicity of the system, and in Refs. 16–18 by using a Bogoliubov transformation.^{19–21} In Ref. 22, the normal modes for a nonuniform 2D planar ion crystal in Penning trap geometry were found, using a method that employed a preliminary diagonalization of the potential matrix and that required the solution of a quadratic eigenvalue problem, but that can nevertheless be connected to a non-canonical version of the approach employed in this paper (see Appendix B). Here, we consider the general case of a nonuniform 3-dimensional magnetized ion crystal, from which the previous 2D results can be obtained as a limit. Some of the modes in 3D nonuniform crystals have been previously described,^{10,23} but a general theory for all of the modes has not been previously published to our knowledge.

In the first half of the paper (Sec. II), we develop a theory for the normal modes of a general linearized classical Hamiltonian system. The theory differs from the standard Bogoliubov approach mentioned above, focusing on the Hermitian properties of linearized Hamiltonians. There are some advantages to our approach: canonical coordinates are not required and extra transformations to/from the creation/annihilation representation of the dynamics (necessary in the Bogoliubov method) are avoided. We include an analysis of neutrally stable modes, since such modes often occur in trapped ion crystals, associated with rotations in symmetric external trap fields. This requires a discussion of the differing forms of the diagonalized Hamiltonian when constants of the motion are in involution (i.e., their Poisson bracket vanishes) or are not in involution. For completeness, we also consider the modes of an unstable Hamiltonian system. We find that canonical coordinates remain mixed in pairs of exponentially growing and decaying mode amplitudes in such a way that the system energy remains conserved as mode amplitudes grow and decay.

In the paper's second half (Sec. III), we apply this theory to determine the modes of an ion crystal and consider two examples in detail, a Coulomb cluster consisting of two charges and a 3D Coulomb crystal with $N\u226b1$. For the former system, all of the normal modes can be evaluated analytically. In the latter case, the modes are evaluated numerically. Averages over thermal fluctuations are discussed. The fluctuation energy associated with vibrations as well as rotations is analyzed. Expressions for thermal fluctuation amplitudes are shown to be consistent with the Bohr–van-Leewen theorem,^{24} provided that any contributions from zero-frequency rotational modes are included in the averages.

In relation to the rotational modes, expressions for rotational inertia of an ion crystal are developed, including a novel magnetic addition to the rotational inertia. This magnetic addition, arising from the vector-potential portion of the angular momentum, is a dominant contribution in many experiments that employ large magnetic fields. A surprising magnetic effect is discussed for the case of a spherically symmetric trap potential, in which the rotational inertia ceases to exist. Under these conditions, variations in angular momentum leave the rotation frequency unaffected but instead change the crystal's orientation with respect to the magnetic field. This phenomenon is connected to the occurrence of constants of the motion (components of the angular momentum) that are not in involution.

In the large magnetic field limit of interest in many of the experiments, reduced eigenvalue problems are developed that separately describe cyclotron, axial, and *E *×* B* drift eigenmodes. As far as we know, the reduced eigenvalue problems for cyclotron and *E *×* B* modes have not be written down previously for a general nonuniform crystal structure, although these limits have been considered for periodic lattices.^{16,17} We also describe the coupling between axial and *E *×* B* modes that occurs in three-dimensional crystals. These reduced eigenvalue problems provide some intuition as to the form of the eigenmodes and eigenfrequencies in large magnetic fields.

## II. NORMAL MODES OF A LINEARIZED HAMILTONIAN SYSTEM

The determination of the normal modes of a system of coupled linear oscillators is a venerable problem in mathematical physics, with applications in a variety of scientific fields. The standard solution is a textbook problem in classical mechanics, in which normal modes are found by solving for the eigenvectors and eigenvalues of combined kinetic and potential energy matrices that arise in a linearized Lagrangian of the form $L=T\u2212V$, where *T* and *V* are kinetic and potential energies, respectively.^{25} The coordinate transformation relating particle displacements to a sum over the eigenmode amplitudes (the “normal coordinates” for the oscillators) is a point transformation, and as such can be easily handled in the context of Lagrangian mechanics. However, for more general Lagrangians in which coordinates and velocities are mixed (for instance, through velocity-dependent potentials that arise in the application of a magnetic field), point transformations are no longer adequate. Problems of this nature arise, for example, in the modes of an ion crystal in a Penning trap, the study of vibrational modes of molecules in applied magnetic fields, and in the normal modes of crystals of point vortices in 2D Euler flow.^{26,27} A transformation mixing both positions and momenta is now required in order to diagonalize the Hamiltonian.

A general approach to the solution of this problem is the Bogoliubov transformation,^{19–21} a linear transformation in which the symplectic condition^{28} for canonical transformations is imposed, and conditions required to diagonalize the Hamiltonian are then determined. The standard Bogoliubov transformation was developed with quantum problems in mind, and formulated in terms of creation and annihilation operators which are related to the position and momentum by a linear transformation. This approach is reviewed in Appendix D.

In this paper, we use a method that is formulated specifically for the classical problem. To our surprise, we have not come across a discussion of this method in any previous publication, so we lay out the theory in some detail (however, the literature on Hamiltonian methods is vast and it is possible, even likely, that this approach is not novel). This method focuses on the Hermitian properties of the matrix operators that appear in linearized Hamiltonian systems, rather than on the symplectic condition. Normal modes are obtained as eigenvectors and eigenfrequencies of the dynamical matrix **D** that determines the linear equations of motion. We derive a Hermitian property of this matrix (with respect to an inner product involving the system Hamiltonian) that is then used to show that the eigenvectors form a complete orthogonal set under certain conditions involving the system stability with respect to small perturbations. Diagonalization of the system energy is easily accomplished using these eigenvectors, without imposing the symplectic condition on the linear transformation.

The added requirement that the transformation to normal mode coordinates be canonical is then satisfied through a specific normalization condition on the eigenvectors. Once this condition is imposed, the transformation is equivalent to the Bogoliubov transform. The equivalence to the Bogoliubov method is discussed in detail in Appendix D of the paper.

Although our diagonalization method is equivalent to the Bogoliubov method, for classical systems, it is arguably more straightforward in both its derivation and its application, as extra transformations from the creation/annihilation representation to the position/momentum representation are not required. Furthermore, it is possible to apply our method to linearized non-canonical Hamiltonian systems, since our method does not require canonical coordinates. We describe an example of this approach to diagonalizing a non-canonical system in Appendix B: the ion crystal in a position-velocity space representation. We have also used a continuum version of this method in calculations involving the normal modes of a non-canonical conservative system of fluid equations.^{29,30}

A linearized Hamiltonian system is a dynamical system with *N* coordinates $q=(q1,\u2026,qN)$ and associated canonical momenta $p=(p1,\u2026,pN)$ whose Hamiltonian has the quadratic form

where $z=(q,p)$ is the system phase-space configuration vector, $f(t)$ is a “forcing” vector, and **H** is the Hamiltonian matrix, a matrix of coefficients independent of **z** and which, for our purposes, is also assumed to be time-independent. (There is also considerable interest in the time-dependent problem, particularly in the area of linear control theory and in generalizations of parametric resonance to multidimensional systems.^{31–33}) Suitable choices of the off diagonal coefficients in **H** allow the matrix to be of symmetric form, satisfying $Hij=Hji$.

The linear equations of motion that arise from this Hamiltonian are, in vector form, given in terms of the Poisson bracket $[\xb7,\xb7]$ as

where we introduce the dynamical matrix $D=J\xb7H$ as well as the fundamental symplectic matrix^{28} $J\u2261[z,z]$. The fundamental symplectic matrix is given, in block form, by

where **1** and **0** are the unit and zero tensors, respectively. The matrix is antisymmetric and expresses the basic Poisson bracket relationships $[qi,qj]=[pi,pj]=0$ and $[qi,pj]=\u2212[pj,qi]=\delta ij$, where $\delta ij$ is the Kronecker delta.

We will consider the normal modes of this linearized Hamiltonian system. The normal modes are unforced (i.e., $f=0$) solutions of Eq. (2) that are assumed to have a time dependence of the form

for some (possibly complex) frequency *ω* and some (time-independent, possibly complex) vector $u\omega $. Substituting Eq. (4) into Eq. (2) and assuming $f=0$ then yields an eigenvalue problem for *ω* and $u\omega $,

Thanks to the Hamiltonian nature of the linear dynamical equations, the eigenfrequencies *ω* and eigenvectors $u\omega $ of Eq. (5) have the following three properties:

The eigenvectors $u\omega $ form an orthogonal set with respect to a generalized inner product defined for any complex vectors

**a**and**b**as $(a,b)\u2261a*\xb7H\xb7b$:$(u\omega ,u\omega \xaf)=0$, provided that $\omega \u2260\omega \xaf*$.

A given eigenvalue

*ω*is real, provided that the corresponding eigenvector satisfies $(u\omega ,u\omega )\u22600$.For each eigenmode $(\omega ,u\omega )$ for which $\omega \u22600$, there is a second eigenmode $(\u2212\omega *,u\u2212\omega *)$ for which $u\u2212\omega *=u\omega *$. Thus, for real

*ω*the $\omega \u22600$ eigenmodes come in $\xb1\omega $ pairs.

It is straightforward to prove these properties. Property 3 arises from the fact that the dynamical matrix **D** has real coefficients. Taking the complex conjugate of Eq. (5) then yields

showing that $u\omega *$ is also an eigenvector of **D** with eigenfrequency $\u2212\omega *$, which completes the proof of property 3.

Properties 1 and 2 follow from the fact that $iD$ is a Hermitian (self-adjoint) matrix with respect to the above-defined inner product, which we prove below. It is well-known that the eigenvalues and eigenvectors of a Hermitian matrix satisfy properties 1 and 2. These properties of Hermitian matrices are often referred to as the spectral theorem, and a proof may be found in many linear algebra textbooks.^{34}

The Hermitian property of the matrix $iD$ is defined by the relationship $(a,iD\xb7b)=(b,iD\xb7a)*$, which must be satisfied for all vectors **a** and **b**. Dividing out the factor of *i* and using the definition of the inner product, this Hermitian relationship can be expressed as

(A matrix **D** that satisfies this equation is sometimes referred to as “anti-Hermitian” with respect to **H**, due to the negative sign in the equation.) Consider the matrix $L\u2261H\xb7D=H\xb7J\xb7H$ that appears in the above expression. This matrix is antisymmetric: $Lji=\u2212Lij$. This follows from the symmetry and antisymmetry, respectively, of the matrices **H** and **J**:

The antisymmetry of **L** proves Eq. (7), which in turn proves that $iD$ is Hermitian.

### A. Stable system

We can use properties 1–3 in order to analyze the evolution of the solution $z(t)$ to the dynamical equations. First consider the simplest case, of a *stable* system for which $(u\omega ,u\omega )\u22600$ for all modes, so that all eigenfrequencies *ω* are real (property 2). Assume also (for simplicity) that there are no *ω* = 0 modes, and that all mode frequencies are different. Cases which have one or more neutrally stable (*ω* = 0) modes introduce certain technical issues that are addressed in Sec. II B, and examples of such neutral modes will be considered in Sec. III. Degeneracies ($\omega \xaf=\omega $ for two or more separate eigenvectors) can be handled easily by orthogonalizing degenerate eigenvectors within the subspace created by these vectors; see, for example, Ref. 25.

Under these assumptions, the set of 2*N* eigenvectors $u\omega $ then form an orthogonal set in the 2*N* dimensional vector space for the phase-space vector **z**, spanning the vector space and thus forming a complete set. We can therefore construct a representation of the vector $z(t)$ in terms of the eigenvectors:

where the complex amplitude $a\omega (t)$ associated with each eigenvector $u\omega $ can be found by taking an inner product of both sides of Eq. (9), applying property 1 (orthogonality) of the eigenvectors:

Since eigenmodes come in pairs (property 3), Eq. (9) can also be written as

The real nature of the vector **z** then implies that $a\u2212\omega (t)=a\omega *(t)$ and we can then write Eq. (11) as

where $c.c.$ stands for complex conjugate.

A differential equation for the time evolution of the complex mode amplitude $a\omega (t)$ follows by substitution of Eq. (9) into Eq. (2). Taking an inner product with respect to one of the eigenvectors $u\omega $ then yields

where $f\omega (t)=(u\omega ,J\xb7f(t))/(u\omega ,u\omega )$. The forcing coefficient $f\omega $ can also be written as

where we used $D=J\xb7H$ and the complex conjugate of Eq. (5), along with the symmetry and antisymmetry, respectively, of **H** and **J**.

The system energy can be written in terms of the mode amplitudes $a\omega (t)$ by substituting Eq. (9) into Eq. (1):

where in the second line we replaced the dummy index $\omega \xaf$ with $\u2212\omega \xaf$ using property 3, in the third line we used property 1 (orthogonality) of the modes, and in the last line we identified the energy $H\omega $ in each eigenmode, given by

using property 3 to sum only over positive eigenfrequencies.

Note also that in many, if not all, applications, an unforced stable oscillator has positive energy compared to the equilibrium $z=0$, in which case Eq. (16) implies that $(u\omega ,u\omega )>0$.

Equations (9), (13), (15), and (16) provide a complete description of the dynamics of a stable linearized Hamiltonian system. We have already diagonalized the system energy in Eq. (15), without consideration of canonical variables in the linear transformation from **z** to mode amplitude variables $a\omega $, and we have found equations for the evolution of each mode amplitude.

However, Eq. (16) can also be thought of as a Hamiltonian for mode *ω*, provided that we introduce the proper canonical variables. These variables can be constructed using the following argument. Consider the Poisson bracket $[a\omega ,a\omega \xaf*]$. This bracket can be evaluated using Eq. (10) and the symmetry of the Hamiltonian matrix:

where in the fourth line we used Eq. (5) and in the last line we used orthogonality of the eigenvectors. Similarly, one can show that

for *ω* and $\omega \xaf$ greater than zero. In this case, both eigenvectors in Eq. (17) are starred. However, recall that $u\omega \xaf*=u\u2212\omega \xaf$, the eigenvector for a negative frequency mode. This mode is orthogonal to all positive frequency modes by property 1, proving Eq. (19).

Now, to define canonical variables based on the complex amplitudes $a\omega $, we find it useful to impose the condition on these amplitudes that for $\omega >0$ and $\omega \xaf>0$,

(The reason for this condition will become clear in a moment.) According to Eq. (18), we therefore choose normalizations of the eigenvectors such that

Since both sides of this equation are positive, a normalization constant for $u\omega $ can be found to satisfy this equation. Note that only the magnitude of this constant is determined. The phase of the constant can be chosen arbitrarily, which allows a certain degree of latitude in the canonical transformation. As an aside, note also that Eqs. (19) and (20) are analogous to the commutator relationships required for the creation and annihilation operators in the Bogoliubov method; see Appendix D.

We can now introduce *N* real-valued canonical pairs $(Q\omega ,P\omega )$, defined by

In order to show that these are canonical pairs, invert Eq. (22) (and its complex conjugate) to give $Q\omega =2\u22121/2(a\omega +a\omega *)$ and $P\omega =\u2212i2\u22121/2(a\omega \u2212a\omega *)$. Then,

and similarly, $[P\omega ,P\omega \xaf]=0=[Q\omega ,Q\omega \xaf]$. We now see the point of Eq. (20): this choice determines that $[Q\omega ,P\omega ]=1$; a different choice would lead to a value other than 1 on the right hand side of this Poisson bracket relationship.

Applying Eq. (22) to Eqs. (15) and (16) and using the normalization condition Eq. (21) yields the diagonalized system Hamiltonian

where $f1\omega (t)=2f(t)\xb7Re(u\omega )$ and $f2\omega (t)=\u22122f(t)\xb7Im(u\omega )$. Note that for our choice of canonical pairs, the $Q\omega $ and $P\omega $ variables have the same dimensions of $energy/frequency$, but other choices are of course possible via a secondary canonical transformation.

Hamilton's equations of motion applied to Eq. (24) then yield

When $f1\omega =f2\omega =0$, the equations for the *ω* mode are unforced and the oscillator energy $H\omega $ is a constant of the motion.

Finally, we note that Eqs. (11) and (22) imply that the linear transformation from phase-space variables $z=(q,p)$ to new variables $Z=(Q,P)$ can be written as a matrix equation

where the $2N\xd72N$ symplectic transformation matrix $S$ is given by

and the $2N\xd7N$ matrix **U** has columns consisting of the *N* eigenvectors $u\omega $ with $\omega >0$, normalized as per Eq. (21). A similar result holds for the Bogoliubov transformation, although there the phase space variables **z** and **Z** are replaced by creation/annihilation pairs (see Appendix D), and extra linear transformations between these pairs and the phase space variables must be performed to obtain Eq. (28).

#### 1. Thermal averages

It is well-known that the diagonalized Hamiltonian simplifies many calculations involving the energy. For example, consider the thermal average of a phase space function $F(z)$,

where *T* is the temperature. In what follows, we assume that the forcing coefficient $f$ is zero.

In many cases, transformation to the canonical variables $(Q\omega ,P\omega )$ can simplify such calculations. Since the transformations are canonical, the Jacobian of the transformation is unity and the averaging integrals become

For instance, it is easy to show that

which implies that the mean energy in each normal mode is the temperature,

Similarly, the classical partition function $Zc\u2261h\u2212N\u222bdz\u2009exp\u2009[\u2212H(z)/T]$ (where *h* is Planck's constant) evaluates to

the same form as the well-known classical limit for oscillators that do not have the position-momentum coupling considered in this paper.^{35}

### B. Neutrally stable system

We now return to Eq. (9) and consider a modification to it that is necessary when there is a neutrally stable (*ω* = 0) eigenmode of Eq. (5). This modification affects the form of the Hamiltonian and involves a new constant of the motion associated with the neutral mode. Let us assume that there is only one such mode, whose eigenvector we label $u0$. This zero-frequency eigenvector satisfies

Thus, $u0$ is in the nullspace of **D**. By assumption, it is the only vector in the nullspace. This eigenvector is real, since **D** is real.

The eigenvector $u0$ is also in the nullspace of the Hamiltonian matrix **H**. This follows by applying to both sides of Eq. (35) the fundamental symplectic matrix **J** and using $D=J\xb7H$:

where, in the last step, we used the identity $J\xb7J=\u22121$.

Now, there are only $2N\u22121$ independent eigenvectors of the dynamical matrix **D**, the $2(N\u22121)$ eigenvectors with non-zero frequencies, and the single zero frequency eigenvector. Since the phase space has dimension 2*N*, the $2N\u22121$ eigenvectors no longer form a complete set that can be used to represent general phase-space vectors $z$. However, we require such a representation in order to fully diagonalize the Hamiltonian. We therefore need one more vector that is orthogonal to the $2N\u22121$ eigenvectors. We will refer to this vector as $u\xaf0$. It is not an eigenvector of $D$, but can instead be obtained by consideration of the constants of the motion. When there is a zero frequency eigenvector, Eq. (36) implies that the Hamiltonian has a symmetry that produces a new constant of the motion, the momentum *P*_{0}. This constant is

The time derivative of *P*_{0} can be shown to equal zero using Eq. (2), if we also assume that $f\xb7u0=0$:

where, in the last step, we used Eq. (36) and the symmetry of the Hamiltonian matrix.

We can now find the vector $u\xaf0$ using the constancy of *P*_{0} in the linear dynamics. Each oscillatory eigenmode with $\omega \u22600$ must satisfy $P0=constant$, and since *P*_{0} is linear in $z$ the only possibility is $P0=0$. Therefore, the eigenmodes all satisfy

Note that this is also true for the *ω* = 0 eigenvector $u0$, since $0=u0\xb7J\xb7u0$ by the antisymmetry of **J**. We will construct a vector $u\xaf0$ that is orthogonal to all of the eigenmodes by solving the equation

for all eigenvectors $u\omega $. By Eq. (39), such a vector will satisfy $(u\xaf0,u\omega )=0$ for all *ω* including *ω* = 0. A necessary and sufficient condition for solution of Eq. (40) is that $u\xaf0$ satisfies

However, since **H** has a vector $u0$ in its nullspace, $u\xaf0$ cannot be obtained via a standard matrix inversion solution to Eq. (41) because the inverse of **H** does not exist. In fact, the solution to Eq. (41) is underdetermined. The right-hand-side is perpendicular (in the usual dot-product sense) to the nullspace of **H** since $u0\xb7J\xb7u0=0$, and this implies that only $2N\u22121$ of the equations in Eq. (41) are linearly independent: the system satisfies $u0\xb7H\xb7u\xaf0=0$. The solution to such a problem can be obtained in a number of ways. For example, one can project Eq. (41) onto the subspace that is perpendicular to $u0$, obtaining $2N\u22121$ independent equations. A unique particular solution for $u\xaf0,\u2009u\xafp$, can then obtained by specifying an extra condition on the solution that, for example, $u\xaf0\xb7u0=0$. The vector $u\xafp$ is real, since all coefficients appearing in the equations are real. The general solution is this particular solution added to the nullspace eigenvector:

for any value of the constant *C*. The value of *C* can be chosen arbitrarily without affecting any of our subsequent results.

We can now use this extended system of vectors to represent a general phase space vector $z$:

The vectors form a complete orthogonal set so that the coefficients $a\omega $ and $a\xaf0$ can be obtained by projection:

However, *a*_{0} cannot be found using the standard projection method because $u0$ is orthogonal to itself: according to Eq. (36), $u0\xb7H\xb7u0=0$. Instead, *a*_{0} can be determined using the properties of the fundamental symplectic matrix. Acting on both sides of Eq. (43) with $u\xaf0\xb7J$, we obtain

However, $u\xaf0\xb7J\xb7u\xaf0=0$ due to the antisymmetry of $J$, and also $u\xaf0\xb7J\xb7u\omega =0$ for $\omega \u22600$. This follows because, for $\omega \u22600$,

We are therefore left with

where in the second form we employed Eq. (41). Thus, we obtain for *a*_{0}

Returning to Eq. (43), the complex coefficients $a\omega $ still come in $\xb1\omega $ pairs satisfying $a\u2212\omega =a\omega *$ for $\omega \u22600$. The time-dependence of these coefficients is still given by Eq. (13). The time-dependence of *a*_{0} and $a\xaf0$ follows in the same way, by substitution of Eq. (43) into the equation of motion, Eq. (2). The result, after projecting out all the $\omega \u22600$ eigenvectors, is

where we employed $D\xb7u0=0$ and where $\Delta f=J\xb7f\u2212\u2211\omega \u22600f\omega u\omega $ is the projection of $J\xb7f$ into the $(u0,u\xaf0)$ subspace. In the second line, we used

which follows from Eq. (41) by applying $J$ to both sides, and using $u0\xb7J=\u2212J\xb7u0$ and $J\xb7J=\u22121$. Taking an inner product of Eq. (50) with respect to $u\xaf0$ then implies that

where in the second line we rewrote the numerator using $(u\xaf0,\Delta f)=(u\xaf0,J\xb7f)=(J\xb7f,u\xaf0)=\u2212f\xb7J\xb7H\xb7u\xaf0=\u2212f\xb7D\xb7u\xaf0=\u2212f\xb7u0$, and in the last step we used Eq. (51).

When $f\xb7u0=0$ Eq. (52) is an expression of the conservation of the momentum *P*_{0} in the dynamics. This can be seen by substituting for **z** from Eq. (43) into Eq. (37), yielding

where we employed Eqs. (39) and (41). Taking a time derivative and using Eq. (52) yield $P\u03070=\u2212f\xb7u0$.

Finally, the dynamics of *a*_{0} follows by acting on both sides of Eq. (50) with $u\xaf0\xb7J$, yielding

Thus, for $f=0$, *a*_{0} increases linearly in time with a rate given by $a\xaf0$.

We now return to the question of the diagonalization of the system energy and the proper choice of canonical pairs. Substituting **z** from Eq. (43) into the system energy Eq. (2), the same argument which led to Eq. (15) now results in some new terms involving the zero-frequency modes:

where $H\omega $ is still given by Eq. (16). Poisson brackets involving *a*_{0} and $a\xaf0$ follow from Eqs. (49) and (45):

This implies that *a*_{0} and the momentum $P0=a\xaf0(u\xaf0,u\xaf0)$ form a canonical pair [see Eq. (53)]. Similarly, one can also show that $[a0,a\omega ]=[P0,a\omega ]=0$ for all $\omega \u22600$.

Finally, we write the diagonalized Hamiltonian in terms of the canonical pairs as

The dynamical equations for $Q\omega $ and $P\omega $ remain unchanged from Eq. (26). The equations of motion for *a*_{0} and *P*_{0} are

which agree with Eqs. (55) and (52). In these equations, the inner product $(u\xaf0,u\xaf0)$ can be interpreted as the inertia associated with the momentum *P*_{0}.

This completes the derivation of the diagonalized Hamiltonian and the canonical coordinates for a neutrally stable linearized Hamiltonian system with a single zero frequency mode.

#### 1. Two neutral modes

Before we move on, a few remarks must be made regarding the case when there is more than one zero frequency mode. This special case actually arises more often than one might suspect; examples include systems with spherical symmetry, systems with translational or cylindrical symmetry that are also undergoing a second order structural phase transition, or systems with translational symmetry in more than one dimension. Let us consider the case where there are two neutrally stable modes; the case of more than two can be understood by an extension of this example.

Now, there are two independent eigenvectors $u01$ and $u02$ in the nullspace of both the dynamical matrix **D** and the Hamiltonian matrix **H**. These eigenvectors produce two constants of the motion, $P01$ and $P02$, given by

Two possibilities must be separately considered: (i) $u01\xb7J\xb7u02=0$ and (ii) $u01\xb7J\xb7u02\u2261J12\u22600$. In case (i), the constants of the motion are in involution (their Poisson bracket vanishes), while in case (ii) they are not. This can be seen by evaluating the Poisson bracket $[P01,P02]$, using Eqs. (61) and (62):

#### 2. Constants of the motion in involution

When the constants are in involution, the eigenvectors by themselves are not a complete set and two vectors $u\xaf01$ and $u\xaf02$ are required in order to describe a general phase-space vector **z** according to

The need for the extra vectors $u\xaf01$ and $u\xaf02$ can be seen by acting on the equation with $u01\xb7J$:

Conservation of $P01$ implies that $u01\xb7J\xb7u\omega =0$, while $u01\xb7J\xb7u01=0$ by symmetry and $u01\xb7J\xb7u02=0$ by assumption. Repeating the procedure with $u02\xb7J$, we have

To satisfy these equations for a general vector **z**, we require nonzero values for both $a\xaf01$ and $a\xaf02$, proving that both vectors $u\xaf01$ and $u\xaf02$ are required for a complete set.

Let us now consider the solution of these equations for $a\xaf01$ and $a\xaf02$, along with the determination of $a01$ and $a02$. Recall that these latter two coefficients can be found by applying $u\xaf01\xb7J$ and $u\xaf02\xb7J$ to Eq. (64). Noting that $u\xaf01\xb7J\xb7u\omega =u\xaf02\xb7J\xb7u\omega =0$ [see Eq. (47)], and that $u\xaf01\xb7J\xb7u\xaf01=u\xaf02\xb7J\xb7u\xaf02=0$ by symmetry, we are left with

We can simplify the solution of Eqs. (66)–(69) by recalling that the vectors $u\xaf01$ and $u\xaf02$ are constructed to be orthogonal to all of the eigenvectors by solution of the (underdetermined) equations

We first use these two equations to replace Eqs. (66)–(69) by the equivalent equations,

Also, we note that since Eqs. (70) and (71) are underdetermined, any linear combination of the nullspace vectors of **H** can be added to particular solutions $u\xafp1$ and $u\xafp2$ of these equations:

It is useful to choose the constants $\alpha 1,\alpha 2,\beta 1,\beta 2$ so that $u\xaf01\xb7J\xb7u\xaf02=0$. This can be accomplished, for example, by taking $\alpha 1=\beta 1=\beta 2=0$ and choosing *α*_{2} such that $u\xafp1\xb7J\xb7u\xaf02+\alpha 2u02\xb7J\xb7u\xaf02=0$. Using Eq. (71), the solution for *α*_{2} is

The condition $u\xaf01\xb7J\xb7u\xaf02=0$ implies that we can drop the $a\xaf01$ and $a\xaf02$ terms in Eqs. (74) and (75), so they become

Equations (79) and (80) can be solved for $a01$ and $a02$, while Eqs. (72) and (73) can be solved for $a\xaf01$ and $a\xaf02$. The solutions are

where $hij=u\xaf0i\xb7H\xb7u\xaf0j,\u2009hiz=u\xaf0i\xb7H\xb7z$, and $jiz=u\xaf0i\xb7J\xb7z$. We assume throughout that the Hamiltonian satisfies $h11h22\u2212h122\u22600$. (This requirement has an origin similar to the requirement that the determinant of the inertia tensor of a rigid body must be nonzero. It is a requirement on any physical Hamiltonian system.)

Canonical variables must be found in order to use this expression as a Hamiltonian. The Poisson brackets of the zero frequency amplitudes can be found using Eqs. (81) and (82). We obtain $[a01,a02]=[a\xaf01,a\xaf02]=0$ and the nontrivial brackets

Canonical pairs can be found by noting that the constants of the motion $P01$ and $P02$ are related to $a\xaf01$ and $a\xaf02$ via

Using these variables along with Eq. (84), it is then an exercise to show that $[a01,P01]=[a02,P02]=1$ while $[a02,P01]=[a01,P02]=[P01,P02]=0$. Thus, the canonical pairs are $(a01,P01)$ and $(a02,P02)$. We now need only invert Eqs. (85) and (86),

and employ these results in the Hamiltonian, which becomes

The equations of motion for $P01$ and $P02$ are then

and the equations of motion for the amplitudes $a01$ and $a02$ are

When the forcing **f** is zero, the Hamiltonian is independent of $a01$ and $a02$, hence the canonical momenta $P01$ and $P02$ are constant, as expected, and $a01$ and $a02$ both have uniform rates of change.

#### 3. Constants of the motion not in involution

When the two zero-frequency modes satisfy $J12=u01\xb7J\xb7u02\u22600$, the constants of the motion $P01$ and $P02$ are not in involution [Eq. (63)]. This case is easier to deal with than the previous case of constants in involution. The Poisson bracket relationships and the Hamiltonian now take a different form. Now the eigenvectors of **D** by themselves form a complete set for any phase space vector **z**, allowing us to write

No extra vectors $u\xaf01$ or $u\xaf02$ are needed. If such vectors were needed, they would satisfy Eqs. (71) and (72), but these equations no longer have solutions. This can be seen by, for example, taking a dot product of $u01$ with Eq. (72):

a contradiction.

The amplitude coefficients *a*_{01} and *a*_{02} in Eq. (94) can now be obtained by acting with $u02\xb7J$ and $u01\xb7J$, respectively, yielding

where we used $u01\xb7J\xb7u\omega =u02\xb7J\xb7u\omega =0$ [see Eq. (47)] and where the second forms in terms of the constants of the motion follow from Eqs. (61) and (62). The coefficients $a\omega $ are still obtained with the usual inner product, see Eq. (10).

where in the second line we employed Eqs. (96) and (97). The equations of motion for $P01$ and $P02$ then follow from Eq. (63):

which are the same as when *P*_{01} and *P*_{02} are in involution [see Eqs. (90) and (91)]. However, now these two variables form a canonical set. This implies that when $f=0$, the amplitudes *a*_{01} and *a*_{02} are time-independent [see Eqs. (96) and (97)]. This behavior differs from the previous case, where *a*_{01} and *a*_{02} continued to evolve at a fixed rate when $f=0$ [Eqs. (92) and (93)].

### C. Unstable System

We now consider the normal modes in an unstable Hamiltonian system. Unstable conservative systems are of importance in several contexts, such as in the study of ideal fluid and plasma instabilities. Here, we will consider the case of a system with an unstable mode with complex frequency $\omega =\Omega \u2261\Omega r+i\gamma $, where $\Omega r>0$ and $\gamma >0$ are the real frequency and growth rate, respectively, with both assumed to be greater than zero. In addition to this complex mode, a second mode with complex frequency $\Omega *=\Omega r\u2212i\gamma $ must also occur. This follows because the dynamical matrix in the eigenmode problem, Eq. (5), has only real coefficients, which implies that, in order to solve the characteristic polynomial in *ω*, all complex mode frequencies must come in pairs, $\Omega r\xb1i\gamma $.

Note that the above mode with frequency $\Omega *$ has a negative growth rate $\u2212\gamma $ and is independent of the mode required by property 3, with frequency $\u2212\Omega *$. The latter mode has a positive growth rate and is needed (when $\Omega r\u22600$) in order to construct a real solution for the phase space configuration $z$, in analogy to the argument accompanying Eq. (12). A fourth mode with frequency $\u2212\Omega $ and a negative growth rate is also required by property 3 and is the complex conjugate of the mode with frequency $\Omega *$, in order to produce a real solution for **z**. (If $\Omega r=0$, only two of these four modes are required, as the others are redundant. In what follows, we assume that $\Omega r>0$. The $\Omega r=0$ case will be briefly discussed at the end of the subsection.)

The modes of the unstable system form an orthogonal set according to property 1, and we use the eigenvectors associated with the modes to describe the phase space vector **z** via

just as was done for a stable system. However, some differences become apparent. According to property 2, a mode with complex eigenfrequency Ω must be orthogonal to itself: $(u\Omega ,u\Omega )=0$. It is therefore not possible to determine $a\Omega $ in terms of $z$ using the standard projection formula, Eq. (10). Fortunately, however, the complex mode with frequency $\Omega *$ can be used to determine $a\Omega $ via projection. According to property 1, this mode is orthogonal to all other eigenmodes (as well as itself), *except* for the mode with complex frequency Ω. Therefore, for modes with complex frequencies Eq. (10) are replaced by

with an analogous expression for $a\Omega *$. The dynamics of a mode with complex frequency then follows by substitution of Eq. (101) into the equation of motion Eq. (2), followed by projection, just as for the stable modes:

where $f\Omega (t)=(u\Omega *,J\xb7f)/(u\Omega *,u\Omega )$. This forcing coefficient can also be written as

using the same algebraic steps which led to Eq. (14). When $f\Omega =0$, the differential equation (103) is unforced and the solution grows exponentially with time at the growth rate *γ*, $a\Omega =A\Omega \u2009exp\u2009(\u2212i\Omega t)\u221d\u2009exp\u2009(\gamma t)$, where $A\Omega $ is an integration constant determined by initial conditions. The analogous equation for the mode with frequency $\Omega *$ implies a decaying mode amplitude $a\Omega *=A\Omega *\u2009exp\u2009(\u2212i\Omega *t)\u221d\u2009exp\u2009(\u2212\gamma t)$.

The system energy can be found in terms of the mode amplitudes by substitution of Eq. (101) into Eq. (1), just as for a stable system. However, when the orthogonality of the eigenmodes is applied, the energy is no longer perfectly diagonalized:

where *H _{u}*, the unstable mode contribution to the energy, is

and where the stable mode contribution $H\omega $ is unchanged, given by Eq. (16). The non-diagonal form of *H _{u}* is required by energy conservation. For an unforced system, one can see that although $a\Omega (t)$ and $a\Omega *(t)$ have differing time dependences $exp\u2009(\u2212i\Omega t)$ and $exp\u2009(\u2212i\Omega *t)$, respectively (the former growing and the latter decaying), the combination $a\Omega a\Omega **$ is time-independent, as required for an energy-conserving system.

Just as for a stable system, the energy can be formulated as a Hamilitonian when the proper canonical coordinates are introduced. First, we consider the Poisson bracket $[a\Omega ,a\omega *]$ for $Re\omega >0$. Using Eq. (102) and the same series of steps as led to Eq. (18), we obtain

and similarly, $[a\Omega ,a\omega ]=0$ for all eigenfrequencies *ω* with $Re\omega >0$. For canonical coordinates, we therefore choose normalization $u\Omega *$ and $u\Omega $ such that

It is possible for the inner product in Eq. (107) to evaluate to the complex frequency Ω because the vectors appearing in the inner product are different, with different normalization coefficients. We then introduce real-valued canonical pairs $(Q\Omega ,P\Omega )$ and $(Q\Omega *,P\Omega *)$ via the linear transformation

With these choices, one can easily show that $[Q\Omega ,P\Omega ]=[Q\Omega *,P\Omega *]=1$ and $[Q\Omega ,Q\omega ]=[P\Omega ,P\omega ]=0$. When written in terms of these variables, the unstable mode Hamiltonian is

where we have also employed Eq. (107) and have taken $\Omega =\Omega r+i\gamma $. Hamilton's equations for the complex-frequency modes then yield

The first two equations represent the dynamics of the unstable mode and may be seen to agree with Eqs. (103) and (104), once one applies Eqs. (107)–(109). Similarly, the last two equations determine the dynamics of the exponentially decaying mode and may be seen to agree with Eq. (103) upon replacing $\Omega \u2192\Omega *$ in this equation.

Finally, we briefly mention a few salient points regarding the special case $\Omega r=0$. In this case, the eigenvalue problem for the unstable mode with frequency $\Omega =i\gamma $ can be written as $\gamma ui\gamma =D\xb7ui\gamma $. Since *γ* and $D$ are real, the unstable eigenvector $ui\gamma $ is also real. The corresponding mode with frequency $\Omega *=\u2212i\gamma $ also has a real eigenvector $u\u2212i\gamma $. Since $\Omega =\u2212\Omega *$, these two modes are already paired according to property 3, and there are no other associated complex eigenmodes. Thus, Eq. (101) becomes

with $(ui\gamma ,ui\gamma )=(u\u2212i\gamma ,u\u2212i\gamma )=0$ according to property 2, but $(u\u2212i\gamma ,ui\gamma )\u22600$ according to property 1. Following through with the rest of the algebra, we find that $ai\gamma (t)$ and $a\u2212i\gamma (t)$ are real; that the normalization condition for canonical coordinates is $(u\u2212i\gamma ,ui\gamma )=\gamma $; that the canonical coordinates can be chosen as $Q\gamma =ai\gamma ,P\gamma =a\u2212i\gamma $; and that for this choice, the unstable mode Hamiltonian is

This Hamiltonian leads to the equations of motion

which agree with Eqs. (103) and (104) when $\Omega =\xb1i\gamma $. The first equation represents the dynamics of the unstable growing mode, with the second equation corresponding to the exponentially decaying mode.

This completes our discussion of the modes of an unstable linearized Hamiltonian system.

## III. NORMAL MODES OF AN ION CRYSTAL

As an example of the Hamiltonian approach outlined in Sec. II, consider the dynamics of *N* positive charges confined in the fields of a Penning trap: a uniform magnetic field $B=\u2212Bz\u0302$, with *B *>* *0, and an electrostatic trap potential $\varphi 0(r,z)$ that is confining in the *z* direction for positive charges. In some experiments, this potential is nearly a pure quadrupole, $\varphi 0(r,z)=(1/2)Q(z2\u2212r2/2)$, where *Q *>* *0, and this form will be used in our examples. However, this quadrupolar form is not necessary in the general theory described below. Each particle has mass *m _{i}*, charge $qi>0$, and position

**r**

_{i}, $i=1,\u2026,N$. (For negative charges in the trap, remove the −sign from

**B**and add a −sign to $\varphi 0$ so that

*B*and

*Q*remain positive, and treat $qi>0$ as the magnitude of each charge. This preserves the signs for all the subsequent coefficients and formulas used in this section.)

The charges are assumed to rotate about the *z* axis with some mean rotation frequency $\omega r>0$ (i.e., the rotation is in the positive $\varphi $ direction). In a frame rotating with the charges, the system Hamiltonian is

Here, the canonical momentum for particle *i* is $pi=mir\u02d9i+mi\Omega iA(ri)$ where $\Omega i=qiB/(mic)\u22122\omega r$ is the “vortex frequency” for particle *i* (the cyclotron frequency shifted by Coriolis effects) and $A(r)$ is the scaled magnetic vector potential, defined so that $\u2207\xd7A=\u2212z\u0302$, where $z\u0302$ is the unit vector in the z direction. A useful gauge choice for $A$ is the cylindrically symmetric gauge $A=\u2212(1/2)z\u0302\xd7r$. This choice of gauge makes

The function $\Phi $ is the total electrostatic potential energy of the system (as seen in the rotating frame) given by

where $\varphi ij=qiqj/|ri\u2212rj|$ is the electrostatic Coulomb potential between particles *i* and *j* (neglecting, for simplicity, image charge effects in the surrounding electrodes), and

is the effective external potential energy for charge *i* as seen in the rotating frame, including both the force from rotation through the magnetic field and centrifugal force. For a quadrupole trap potential, this can be written as $qi\varphi i(r,z)=(1/2)qiQ(z2+\beta ir2)$, where the trap parameter $\beta i\u2261\omega r(B/c\u2212mi\omega r/qi)/Q\u22121/2$.

For positive values of the trap parameters *β _{i}* and

*Q*[or, more generally, for a potential $qi\varphi i(r,z)$ that increases from the trap center with both increasing

*r*and

*z*], a (neutrally) stable ion crystal equilibrium exists with $r\u02d9i=0$ and $ri=Ri$ for equilibrium positions

**R**

_{i}satisfying $\u2202\Phi /\u2202Ri=0,i=1,\u2026,N$. Equilibrium configurations can be obtained numerically by minimization of the system potential energy $\Phi $, although for small numbers of charges one can find equilibria analytically via symmetry arguments.

^{36}For $N\u226b1$, there are typically many such equilibria corresponding to different crystalline configurations with slightly different energies and arrangements of the charges. These crystal equilibria have been discussed in some detail in several previous publications.

^{9,37–39}We will consider two examples in detail. Figure 1 displays the simplest nontrivial crystal consisting of two identical charges in a quadrupole trap, a 2-ion Coulomb cluster.

^{36,41,42}When the trap parameter

*β*is greater than one, the equilibrium has the charges on the $\xb1z$ axis, each at a distance $d=(q/4Q)1/3$ from the origin. When $\beta <1$, the ions are on opposite sides of the origin in the

*x*–

*y*plane, each a distance $d/\beta 1/3$ from the origin, and for

*β*= 1 the charges can be at any angle

*θ*with respect to the

*z*axis. This Coulomb cluster can be thought of as a classical version of a symmetric molecule such as

*H*

_{2}or

*N*

_{2}, in which the electrons are replaced by a neutralizing background charge (the “plum-pudding” model of J. J. Thompson). Later in this section, we will analytically evaluate the normal modes for this system, including the effect of the magnetic field.

An equilibrium configuration with larger *N* is displayed in Fig. 2, which shows the *r*–*z* positions of a crystal (local minimum energy state) consisting of *N *=* *236 identical charges in a quadrupole trap with trap parameter $\beta =3/4$. The charges tend to arrange themselves in spheroidal shells^{37,38} with an average density that is determined by the rotation rate and the external trap fields. For $\beta <1$, the system tends to form an oblate spheroid, which for sufficiently small *β* collapses into the *z *=* *0 plane. This particular regime of a single-plane plasma crystal is currently of interest as a useful system for the purposes of quantum simulation.^{1,2,43} For $\beta >1$, the system forms a prolate spheroid and for sufficiently large *β* the system forms a one-dimensional Coulomb string of charges distributed along the *z* axis.^{44,45}

These equilibria are all neutrally stable with respect to rotations about the *z* axis. When *β* = 1, the spherical symmetry of the effective trap potential implies that rotations about *x* and *y* axes are also neutral modes.

### A. Perturbed Hamiltonian and the dynamical matrix for crystal modes

In what follows, we assume that particles are displaced only slightly from one such equilibrium, with positions $ri=Ri+\delta ri$ and momenta $pi=mi\Omega iA(Ri)+\delta pi$. The Taylor expansion to second order in the small displacements from equilibrium then results in a linearized Hamiltonian system, with Hamiltonian

where $Vij=\u22022\Phi /\u2202Ri\u2202Rj$. The form of $Vij$ is given in Appendix A.

This Hamiltonian can be put in the matrix form of Eq. (1) with $f=0$, with phase space vector $z=(r,p)$ where $r=(\delta r1,\u2026,\delta rN)$ and $p=(\delta p1,\u2026,\delta pN)$, and with the symmetric Hamiltonian matrix given in block form by

Here, $V,\u2009C,\u2009\Omega $, and $M\u22121$ are $3N\xd73N$ matrices. The matrix $M\u22121$ is the inverse of the diagonal mass matrix **M** for the system, with the diagonal elements $Mii\u22121$ given by the vector $(m1\u22121,m1\u22121,m1\u22121,\u2026,mN\u22121,mN\u22121,mN\u22121)$. The matrix $\Omega $ is the Lorentz matrix coupling positions and momenta in the Hamiltonian. For our symmetric choice of vector potential, the Lorentz matrix is antisymmetric and is zero everywhere except in 3 × 3 blocks along the diagonal:

Each diagonal block is, in dyadic notation, given by

The symmetric matrix $V=\u22022\Phi /\u2202R\u2202R$ is the potential energy matrix given in block form by

The matrix **C** is a magnetic potential contribution whose elements are zero everywhere except along the diagonal. The vector of diagonal elements $Cii$ is given by

The dynamical matrix $D=J\xb7H$ corresponding to this Hamiltonian matrix is

### B. Center of mass modes and rotational modes

It is well-known that a few of the eigenmodes have simple analytic descriptions. In a pure quadrupole trap with a single species, there are three “center of mass” (COM) modes that consist of a displacement of the entire crystal. The axial COM mode consists of an oscillation in the *z* direction and has frequency *ω _{z}* where $\omega z\u2261qQ/m$ is referred to in the literature as the single particle axial frequency; it is the frequency at which a single trapped particle oscillates in

*z*when displaced from the origin. The cyclotron and

*E*×

*B*COM modes consist of rotational motions of the center of mass on the

*x*–

*y*plane, with frequencies $\omega +$ and $\omega \u2212$, respectively, where

and where $\omega \u22a5\u2261\beta \omega z$ is the single particle transverse frequency in an unmagnetized trap. All three COM mode frequencies are independent of the number of charges in the trap.

In addition, in a cylindrically symmetric trap potential (quadrupolar or not), there is a zero frequency eigenmode which is a pure rigid rotation about the z axis, with eigenvector $u0z=(r0z,p0z)$ where

and where *R _{j}* is the cylindrical radius of equilibrium position

**R**

_{j}for the

*j*th ion. As was discussed previously in more general terms, this eigenmode corresponds to a constant of the motion, the momentum $P0z$ given by Eq. (37):

where in the last step we substituted for $\delta pi$ using Eq. (114). The constant $P0z$ is, of course, the perturbed total canonical angular momentum associated with rotations about the *z* axis.

The corresponding vector $u\xaf0z$, required for the rotational inertia $(u\xaf0z,u\xaf0z)$ [see Eq. (58)], is the solution of Eq. (41). In general, this equation requires a numerical solution but it might also be of interest to note that there is a case where $u\xaf0z$ can be evaluated analytically: when the equilibrium consists of identical charges trapped in a quadrupolar trap with trap parameter *β* chosen to be sufficiently small so that the crystal equilibrium is a planar crystal confined to the *z *=* *0 plane. In this case (see Appendix C),

and the rotational inertia $(u\xaf0z,u\xaf0z)$ is then given by the expression

The first term in the parenthesis gives the usual kinetic inertia associated with rigid rotation, while the second term (which can dominate in strong magnetic fields) arises from the electrostatic energy associated with compression/expansion of the crystal as the rotation rate changes.

### C. Modes of a two particle Coulomb cluster

A few other cases allow an analytic solution for all of the modes. One example is the *N *=* *2 Coulomb cluster shown in Fig. 1. Assuming that the particles are aligned in equilibrium in the *x*–*z* plane, with particle 1 above the *z *=* *0 plane and particle 2 below the plane, the potential matrix $V$ evaluates to

The frequencies and corresponding eigenvectors are provided in Table I. These eigenvectors are not normalized.

ω
. | $u\omega $ . | $(u\omega ,u\omega )$ . |
---|---|---|

COM modes: | ||

ω _{z} | $(0,0,1,0,0,1,0,0,\u2212im\omega ,0,0,\u2212im\omega )$ | $4m\omega 2$ |

$\omega \xb1$ | $(1,\xb1i,0,1,\xb1i,0,\u2212im\omega \xb1im\Omega 2,\xb1m\omega \u2212m\Omega 2,0,\u2212im\omega \xb1im\Omega 2,\xb1m\omega \u2212m\Omega 2,0)$ | $4m(\omega 2+\omega \u22a52)$ |

Other modes for $\beta >1$: | ||

$3\omega z$ | $(0,0,1,0,0,\u22121,0,0,\u2212im\omega ,0,0,im\omega )$ | $4m\omega 2$ |

$\omega r\xb1$ | $(1,\xb1i,0,\u22121,\u2213i,0,\u2212im\omega \xb1im\Omega 2,\xb1m\omega \u2212m\Omega 2,0,im\omega \u2213im\Omega 2,\u2213m\omega +m\Omega 2,0)$ | $4m(\omega 2+(\beta \u22121)\omega z2)$ |

Other modes for $\beta <1$: | ||

$\omega z1\u2212\beta $ | $(0,0,1,0,0,\u22121,0,0,\u2212im\omega ,0,0,im\omega )$ | $4m\omega 2$ |

$\Omega 2+3\omega \u22a52$ | $(1,i\Omega \omega ,0,\u22121,\u2212i\Omega \omega ,0,\u2212im\omega +im\Omega 22\omega ,m\Omega 2,0,im\omega \u2212im\Omega 22\omega ,\u2212m\Omega 2,0)$ | $4m\omega 2$ |

0 | $(0,1,0,0,\u22121,0,m\Omega 2,0,0,\u2212m\Omega 2,0,0)$ | 0 |

Other modes for β = 1: | ||

$\omega \theta \xb1$ | $(r\theta ,\u2212r\theta ,p\theta ,\u2212p\theta )$ | $4m3\omega 2\Omega 4\u2212\omega 2(\Omega 2\u22123\omega z2)\u22123\Omega 2\omega z2\u2009cos\u20092\theta \omega z2\u2009sin2\theta $ |

$r\theta =(\u2212i\omega ,\Omega ,\u22122i\omega \omega 2\u2212\Omega 2\u22123\omega z2\u2009sin2\theta 3\omega z2\u2009sin\u20092\theta )$ | ||

$p\theta =m(\Omega 22\u2212\omega 2,\u2212i\omega \Omega 2,\u22122\omega 2\u2009\omega 2\u2212\Omega 2\u22123\omega z2\u2009sin2\theta 3\omega z2\u2009sin\u20092\theta )$ | ||

0 | $u0z=(0,1,0,0,\u22121,0,m\Omega 2,0,0,\u2212m\Omega 2,0,0)$ | 0 |

0 | $u0y=(cos\u2009\theta ,0,\u2212sin\u2009\theta ,\u2212cos\u2009\theta ,0,\u2009sin\u2009\theta ,0,\u2212m\Omega 2cos\u2009\theta ,0,0,m\Omega 2cos\u2009\theta ,0)$ | 0 |

ω
. | $u\omega $ . | $(u\omega ,u\omega )$ . |
---|---|---|

COM modes: | ||

ω _{z} | $(0,0,1,0,0,1,0,0,\u2212im\omega ,0,0,\u2212im\omega )$ | $4m\omega 2$ |

$\omega \xb1$ | $(1,\xb1i,0,1,\xb1i,0,\u2212im\omega \xb1im\Omega 2,\xb1m\omega \u2212m\Omega 2,0,\u2212im\omega \xb1im\Omega 2,\xb1m\omega \u2212m\Omega 2,0)$ | $4m(\omega 2+\omega \u22a52)$ |

Other modes for $\beta >1$: | ||

$3\omega z$ | $(0,0,1,0,0,\u22121,0,0,\u2212im\omega ,0,0,im\omega )$ | $4m\omega 2$ |

$\omega r\xb1$ | $(1,\xb1i,0,\u22121,\u2213i,0,\u2212im\omega \xb1im\Omega 2,\xb1m\omega \u2212m\Omega 2,0,im\omega \u2213im\Omega 2,\u2213m\omega +m\Omega 2,0)$ | $4m(\omega 2+(\beta \u22121)\omega z2)$ |

Other modes for $\beta <1$: | ||

$\omega z1\u2212\beta $ | $(0,0,1,0,0,\u22121,0,0,\u2212im\omega ,0,0,im\omega )$ | $4m\omega 2$ |

$\Omega 2+3\omega \u22a52$ | $(1,i\Omega \omega ,0,\u22121,\u2212i\Omega \omega ,0,\u2212im\omega +im\Omega 22\omega ,m\Omega 2,0,im\omega \u2212im\Omega 22\omega ,\u2212m\Omega 2,0)$ | $4m\omega 2$ |

0 | $(0,1,0,0,\u22121,0,m\Omega 2,0,0,\u2212m\Omega 2,0,0)$ | 0 |

Other modes for β = 1: | ||

$\omega \theta \xb1$ | $(r\theta ,\u2212r\theta ,p\theta ,\u2212p\theta )$ | $4m3\omega 2\Omega 4\u2212\omega 2(\Omega 2\u22123\omega z2)\u22123\Omega 2\omega z2\u2009cos\u20092\theta \omega z2\u2009sin2\theta $ |

$r\theta =(\u2212i\omega ,\Omega ,\u22122i\omega \omega 2\u2212\Omega 2\u22123\omega z2\u2009sin2\theta 3\omega z2\u2009sin\u20092\theta )$ | ||

$p\theta =m(\Omega 22\u2212\omega 2,\u2212i\omega \Omega 2,\u22122\omega 2\u2009\omega 2\u2212\Omega 2\u22123\omega z2\u2009sin2\theta 3\omega z2\u2009sin\u20092\theta )$ | ||

0 | $u0z=(0,1,0,0,\u22121,0,m\Omega 2,0,0,\u2212m\Omega 2,0,0)$ | 0 |

0 | $u0y=(cos\u2009\theta ,0,\u2212sin\u2009\theta ,\u2212cos\u2009\theta ,0,\u2009sin\u2009\theta ,0,\u2212m\Omega 2cos\u2009\theta ,0,0,m\Omega 2cos\u2009\theta ,0)$ | 0 |

We first consider $\beta >1$, where the charges align along the *z* axis in equilibrium, with *θ* = 0. In this case, there is no zero frequency rotational mode. In addition to the three center of mass modes, there are three other modes in which the charges perform opposite motions, $\delta r1=\u2212\delta r2$. One of these is an axial stretch mode only along the *z* axis, with frequency $3\omega z$. The other two modes consist of circular motion in *x* and *y* at the frequencies $\omega r+$ and $\omega r\u2212$ where

As $\beta \u21921,\omega r\u2212$ approaches zero frequency and the eigenvector corresponds to a sum of rotations about the *x* and *y* axes, which are neutral modes in the spherically symmetric *β* = 1 limit.

The energy takes the standard form for a stable system [see Eq. (15)]:

with the coefficients $(u\omega ,u\omega )$ given in Table I.

On the other hand, for $\beta <1$, the equilibrium is in the *x*–*y* plane. Now, there is a zero frequency rotation about the *z* axis in addition to the three center of mass modes. The associated constant of the motion is the angular momentum given by Eq. (125). For the two ion system, it is convenient to divide out the radii $R=d/\beta 1/3$ of the ions, defining $P\u20320z=P0z/R=\delta p1y\u2212\delta p2y\u2212m\Omega (\delta x1\u2212\delta x2)$. In addition, we introduce the normalized vector $u\xaf\u20320=u\xaf0/R$ with $u\xaf0$ given by Eq. (126). The two other modes are a tilt mode consisting only of axial motion with $\delta z1=\u2212\delta z2$, at frequency $\omega z1\u2212\beta $, and an “upper-hybrid” mode consisting of elliptical motion of the charges in the *x*–*y* plane at frequency $\Omega 2+3\omega \u22a52$. Due to the neutral mode, the energy takes the form

where the scaled moment of inertia $(u\xaf\u20320,u\xaf\u20320)=2m(1+\Omega 2/3\omega \u22a52)$, see Eq. (127).

Finally, when *β* = 1, the system has spherical symmetry and there are now equilibria oriented at any angle *θ* with respect to the *z* axis. In addition to the center of mass modes, there are two zero frequency modes consisting of free rotations about the *y* and *z* axes, and two modes whose frequencies $\omega \theta +$ and $\omega \theta \u2212$ depend on *θ*:

As $\theta \u21920$, these frequencies approach $\omega r+=\Omega $ and the stretch mode frequency $3\omega z$. As $\theta \u2192\pi /2,\u2009\omega \theta \u2212\u21920$, and $\omega \theta +\u2192\Omega 2+3\omega z2$.

#### 1. Nonexistence of the rotational inertia for spherical confinement

The two zero frequency modes produce two constants of the motion, the angular momentum $P\u20320z$ due to rotation of the cluster about the *z* axis (again dividing the radii $d\u2009sin\u2009\theta $) and the angular momentum $P\u20320y$ due to rotations about the *y* axis (also dividing out distance *d* from the axis):

These two constants are not in involution, $J\u2032xy=[P\u20320y,P\u20320z]=2m\Omega \u2009cos\u2009\theta $. Therefore, the constants do not appear in the energy, so it takes the form given in Eq. (130).

It is instructive to compare the evolution of the *β* = 1 system to that of the $\beta <1$ system when the perturbed axial angular momentum $P\u20320z$ is nonzero. This comparison illustrates a physical difference between systems with constants of the motion in involution and those for which the constants are not in involution. For $\beta <1$, Hamiltonian (131) implies that the angle variable $\delta \varphi =a0$ changes linearly with time according to

because a change in angular momentum corresponds to a change in the rotation frequency of the cluster. There is a finite scaled moment of inertia $(u\xaf\u20320,u\xaf\u20320)$ which relates the scaled angular momentum change $P0z\u2032$ to the rotation frequency change $\delta \varphi \u0307$.

However, for *β* = 1, neither $P\u20320z$ nor $P\u20320y$ appear in the Hamiltonian. The phase space configuration $z$ evolves according to

where $u0z$ and $u0y$ are given in Table I, and where

and

[see Eqs. (94), (96), and (97)]. Equation (137) shows that the angle $\delta \varphi $ does not evolve in time as it did for $\beta <1$ [see Eq. (134)]. This is because the axial angular momentum perturbation $P\u20320z$ does not change the rotation frequency of the cluster. Indeed, the moment of inertia $(u\xaf0z,u\xaf0z)$ is undefined since the vector $u\xaf0z$ does not exist, as discussed in Sec. II B 2.

The axial angular momentum perturbation $P\u20320z$ is instead accomplished by a rotation $\delta \theta $ of the cluster about the *y* axis, as exhibited in Eq. (136). This can occur because for $\Omega \u22600$ the canonical angular momentum depends on the cylindrical radius *R* of the charges (through the vector potential term), which varies as *θ* varies—see Fig. 1(c).

The surprising nonexistence of the rotational inertia in magnetized spherically symmetric crystals also occurs for larger Coulomb clusters, provided that a cluster's canonical angular momentum depends on its orientation for spherical confinement. A few cases were discussed in Ref. 39 in the context of a study of the configurations of minimum energy as angular momentum is varied.

### D. Modes of an N = 236 particle crystal

Turning now to an example with more charges, in Figs. 3 and 4, we display the mode frequencies for the case of the oblate spheroidal Coulomb crystal with *N *=* *236 that was shown in Fig. 2. We choose Ω = 0 in Fig. 3 and $\Omega =20\omega z$ in Fig. 4. There is one zero frequency mode corresponding to a rotation about the *z* axis, with an associated constant of the motion corresponding to the perturbed angular momentum. There are no degeneracies in the spectrum, except for the two COM modes consisting of oscillations of the whole crystal in the *x*–*y* plane at frequency $\omega \u22a5=\beta \omega z$. This degeneracy is broken by the application of a magnetic field (i.e., finite Ω); see Eq. (123). The density of states for these modes is similar to that displayed in Fig. 10 of Ref. 10. There are several modes with quite low frequencies, corresponding to torsional motions of the crystal with weak restoring forces.

For $\Omega =20\omega z$, the mode frequencies condense into three groups: a group of *N* cyclotron modes with $\omega >\Omega $; a group of *N* axial modes with $\omega \u22481$, and a group of *N E *×* B* modes with low frequencies. Variation of Ω indicates that these latter mode frequencies scale with Ω as $1/\Omega $. These mode groupings have been identified in previous publications.^{10,17,22,40} The cyclotron modes consist primarily of cyclotron motion in the *x*–*y* plane. The axial modes consist primarily of axial oscillations, and the *E *×* B* modes are drifts of the particle guiding centers in the collective electric field of the other particles as well as the trap field. There is one COM mode in each group, and there is one neutral mode in the *E *×* B* group, associated with rotation about the *z* axis.

Apart from a few special cases such as the center of mass and neutral modes, the eigenmodes have numerical forms the details of which vary depending on the precise crystal structure. Nevertheless, some of these modes are close to the sorts of modes predicted in a cold fluid theory of the normal mode oscillations.^{40} In the cold fluid theory, normal modes for a plasma confined in a quadrupolar trap were worked out analytically in terms of Legendre functions. Frequency shifts to these modes due to interparticle correlations were also worked out and were found to decrease in magnitude as the mode wavelength increases.^{10} An initial condition consisting of displacements associated with a given cold fluid normal mode can be described as a superposition of exact crystal eigenmodes, with only a small number of these eigenmodes dominating, provided that the mode is of low order (with a wavelength large compared to the interparticle spacing).^{10}

A single example is displayed in Fig. 5. Here, we consider a fluid displacement of the form $\delta rfluid=(\delta x,\delta y,\delta z)fluid\u221d(x,\u2212y,0)$. Such a displacement creates an ellipsoidal distortion of the crystal that, in fluid theory, is associated with two modes: a cyclotron frequency mode and an *E *×* B* diocotron mode, in which this distortion propagates around the plasma in the $\varphi $ direction. For the case of a spheroidal plasma with $\beta =3/4$, the frequencies of these two modes is predicted to be $\omega fluid=\Omega 2/4+0.925253\omega z2\xb1\Omega /2$. In Fig. 5, we determine the mode amplitude $a\omega $ for each mode by projecting the phase-space displacement $z$ onto each mode using Eq. (10), with the eigenmodes for $\Omega =20\omega z$, and with $z$ given by $\delta rfluid$ for each particle, along with the associated canonical momentum $\u2212m\Omega z\u0302\xd7\delta rfluid/2$: $z=(x1,\u2212y1,0,\u2026,xN,\u2212yN,0,\u2212m\Omega y1/2,\u2212m\Omega x1/2,0,\u2026,,\u2212m\Omega yN/2,\u2212m\Omega xN/2,0)$. This phase space configuration corresponds to an initially-stationary elliptical distortion of the plasma crystal.

The resulting plot of the magnitude of the amplitude of each eigenmode displays three strong peaks near *n *=* *90, 250, 480. These frequency peaks dominate the dynamics, as the system evolves from this initial condition. The frequencies corresponding to each peak are, respectively, $\omega /\omega z=20.0464,1.3063,0.0461$, which are, respectively, in the cyclotron, axial, and *E *×* B* range (see Fig. 4). These frequencies may be compared to *ω _{fluid}* which evaluates to $\omega fluid/\omega z=20.0462,0.0462$ for the cyclotron and

*E*×

*B*mode, respectively. The weaker of the three peaks at $\omega =1.3063\omega z$ does not correspond to either of these fluid modes, but is instead close to the cylindrically-symmetric axial fluid mode [the (2,0) mode]

^{40,46}with fluid frequency $\omega fluid=1.3037\omega z$ for these conditions. In the fluid theory, this cylindrically symmetric mode has no overlap with the

*θ*-dependent fluid displacement used here. Evidently, a weak coupling to this mode occurs due to the finite number of charges in the crystal. There also appears to be strong coupling to some of the very low frequency

*E*×

*B*modes (the broad peak on the right side of the figure), which is also not predicted in fluid theory.

### E. Thermal averages and zero frequency modes

Thermal averages over various functions are of importance in many applications. For example, consider the average $\u27e8\delta zj\delta zk\u27e9$ for two particles *j* and *k*. The *j*th and *k*th axial displacements appearing here correspond to the *α* and *β* components of the phase space vector $z$ where $\alpha =3j$ and $\beta =3k$. Then, remembering that there is a zero-frequency mode in this system, the position *z _{j}* can be written in terms of the eigenmodes via Eqs. (43), (22), and (53),

where we used Eq. (21) to replace *ω* by $(u\omega ,u\omega )$, so that the expression is valid for unnormalized eigenvectors. The last term arises from the average over the zero frequency mode momentum $\u27e8P02\u27e9$, which we evaluated using Hamiltonian (58). There is no contribution from *a*_{0} because the neutrally stable eigenmode is a pure rotation about the *z* axis by angle $\delta \varphi =a0$, with eigenvector given by Eq. (124). This eigenvector has no axial component (i.e., $u0\alpha =0$).

#### 1. Thermal averages for the two particle Coulomb cluster

We can evaluate $\u27e8\delta zj\delta zk\u27e9$ for the two particle system using the information in Table I. For $\beta >1$, only the axial COM mode and the stretch mode contribute. There is no contribution from the zero frequency mode because it does not exist for $\beta >1$, since charges in equilibrium are aligned along the *z* axis. We thus obtain, for $\beta >1$,

where the first term on the right hand side of each equation is from the axial COM mode and the second term is from the stretch mode.

For $\beta <1$, only the axial COM and the tilt mode contribute, and again the rotational mode does not contribute because there is no axial component of $u\xaf0$ when the charges are trapped in the *x*–*y* plane, see Eq. (126). Then, Eq. (139) yields

where the first term is again from the axial COM mode and now the second term is from the tilt mode. These averages diverge when $\beta \u21921$, as the tilt mode becomes zero frequency, allowing large fluctuations in the axial displacements of the charges. Note that the equations imply that $\u27e8(\delta z1+\delta z2)2\u27e9/4=T/(2m\omega z2)$, independent of *β*. This is the mean square fluctuation in the axial center of mass position, with the expected form $T/(Nm\omega z2)$ for a particle of mass $Nm$ in a harmonic well.

Also note that none of these averages depend on the magnetic field strength (i.e., on Ω), as expected from the Bohr–van Leeuwen theorem.^{24} Of course, in this example, none of the modes contributing to the averages depended on Ω. A less trivial application of the theorem arises in the evaluation of a different average, $\u27e8\delta xj\delta xk\u27e9$. This average is finite for an *N *=* *2 cluster with equilibrium positions in the *x*–*z* plane because zero-frequency rotations through $\varphi $ are in the *y* direction and do not affect the average. The formula for this average is the same as Eq. (139), except that now, $\alpha =3j\u22122$ and $\beta =3k\u22122$. For the two particle Coulomb cluster, and taking $\beta >1$, we obtain

Each term on the right hand side depends explicitly on Ω, but their sum does not. In fact, $1/(\omega \u22122+\omega \u22a52)+1/(\omega +2+\omega \u22a52)=1/\omega \u22a52$ and $1/(\omega r\u22122+(\beta \u22121)\omega z2)+1/(\omega r+2+(\beta \u22121)\omega z2)=1/((\beta \u22121)\omega z2)$. When these formulas are applied to Eqs. (144) and (145), we obtain the Ω-independent result

where we used the relationship $\omega \u22a52=\beta \omega z2$.

Re-evaluating the averages for $\beta <1$, we obtain

The first two terms are from the COM modes, the third term is from the upper hybrid mode, and the last term is from the zero-frequency rotational mode, where we used Eq. (126) for $u\xaf0$ and Eq. (127) for $(u\xaf0,u\xaf0)$. Again, each term depends on Ω, but when summed the results are Ω-independent:

As another test of the Bohr–van Leeuwen theorem, we consider the average $\u27e8(\delta r1\xb7R\u03021)2\u27e9=\u27e8(\delta x1\u2009sin\u2009\theta +\delta z1\u2009cos\u2009\theta )2\u27e9$ for the *N *=* *2 cluster at *β* = 1, oriented at angle *θ* with respect to the *z*-axis. Here, $R\u03021=(sin\u2009\theta ,0,\u2009cos\u2009\theta )$ is the unit vector in the direction of $R1$, the equilibrium position of charge 1. This average is finite because zero-frequency rotations of the cluster about the *y* or *z* axes have no effect—see Table I. Only the COM modes and the modes at frequencies $\omega \theta \xb1$ enter the average:

where $\Delta =2\omega \omega 2\u2212\Omega 2\u22123\omega z2\u2009sin2\theta 3\omega z2\u2009sin\u20092\theta $. The first two terms on the right hand side of the second line arise from the three COM modes, and for *β* = 1 they sum to $T/(2m\omega z2)$, the expected radial thermal fluctuation for a single particle of mass 2*m* in a spherically symmetric harmonic well of frequency *ω _{z}*. Surprisingly, perhaps, considering its complexity, the last term sums to $T/(6m\omega z2)$. Thus we obtain

#### 2. Thermal averages for the N = 236 particle crystal

As an example of thermal averages in a larger *N* system, in Fig. 6 we numerically evaluate the following thermal average for the *N *=* *236 crystal: $\u2211j=1N\u27e8\delta zj2\u27e9\omega $, where $\u27e8\delta zj2\u27e9\omega =2T|u\omega \u20093j|2/(u\omega ,u\omega )$ is the mean square fluctuation in axial position *z _{j}* caused by mode

*ω*, as per Eq. (139).

We use this average to evaluate each mode's contribution to $\u27e8\delta z2\u27e9=N\u22121\u2211j=1N\u27e8zj2\u27e9$ via

see Eq. (139).

From the figure, one can see that for Ω = 0 the lowest frequency torsional modes dominate the thermal average, as one would expect (in the figure the mode number *n* is ordered from highest to lowest frequency as shown in Fig. 3). For the case of a large magnetic field $\Omega =20\omega z$, cyclotron modes make a negligible contribution to the average, as one might also expect; we obtain $\u2211\omega cyclotron\u2211j=1N\u27e8\delta zj2\u27e9\omega =1.096\xd710\u22126T/(m\omega z2)$. Axial modes make a larger contribution, contributing to the average an amount $\u2211\omega axial\u2211j=1N\u27e8\delta zj2\u27e9\omega =385.0969T/(m\omega z2)$. Surprisingly, however, the low-frequency *E *×* B* drift modes dominate by a factor of 10, contributing $\u2211\omega E\xd7B\u2211j=1N\u27e8\delta zj2\u27e9\omega =3575.0989T/(m\omega z2)$. One normally thinks of $E\xd7B$ drift modes as motions in the *x*–*y* plane, but there can also be substantial axial motion in these modes, as we will see; and low-frequency torsional *E *×* B* motions make a large contribution to the axial fluctuations. Finally, the zero-frequency rotational mode (the last term in Eq. 154) contributes 0.2463 to the right hand side of the equation, for a total mean square fluctuation of $\u27e8\delta z2\u27e9=3690.4421/N=16.7815$ in units of $T/m\omega z2$. This fluctuation is independent of magnetic field strength as expected from the Bohr–van-Leeuwen theorem; we have checked that precisely the same value can be computed by summing over the Ω = 0 mode contributions shown in the figure. Note that for Ω = 0, there is no contribution to this thermal average from the zero frequency rotational mode, since for Ω = 0, the axial components of the vector $u\xaf0z$ vanish, i.e., $u\xaf0z\u20093j=0$; see Appendix C.

Calculations like those displayed in Fig. 6 may impact the prospects of extending sensitive techniques^{2,3} developed for the spectroscopy and thermometry of normal modes with single-plane ion crystals to large three-dimensional ion crystals. The basic technique uses an oscillating spin-dependent force with a frequency *μ* to map motion parallel to the magnetic field and with the same frequency *μ* onto precession of an internal spin-degree of freedom of the trapped ion. The spin precession can then be readout with a high signal-to-noise ratio. Technically, the spin-dependent force is an optical dipole force created at the intersection of two laser beams and is characterized by a wavelength *λ _{eff}*. Accurate spectroscopy and thermometry require the so-called Lamb–Dicke confinement criteria, where the ion axial fluctuations $\u27e8\delta z2\u27e9$ are small compared to

*λ*. This could open the door for trapped-ion quantum simulation and sensing work

_{eff}^{1,4,5}with large three-dimensional ion crystals. The large axial excursions of the low-frequency

*E*×

*B*modes should improve the prospects for laser Doppler cooling of these modes relative to that possible for single-plane crystals.

### F. Limiting cases

We now consider several limiting cases for which the eigenmode problem for charges in a trap simplifies. We first consider the unmagnetized case $\Omega i=0$. In this case, the dynamical matrix Eq. (122) reduces to

and the eigenvalue problem Eq. (5) corresponds to the coupled equations $\u2212i\omega r=M\u22121\xb7p,\u2212i\omega p=\u2212V\xb7r$ which can be combined into a reduced eigenvalue problem for $\omega 2$, $\omega 2r=M\u22121\xb7V\xb7r$. This is the standard Hermitian eigenvalue problem for coupled unmagnetized oscillators presented in many physics textbooks and referred to in the introduction to Sec. II. The real eigenvectors $r$ are found to form an orthogonal set with respect to the reduced inner product $(a,b)=a\xb7M\xb7b$ for real vectors $a$ and $b$. The mode frequencies *ω* are found to be real provided that the crystal equilibrium is stable (or neutrally stable), just as in the more general problem discussed in Sec. II. Other features of the unmagnetized eigenmodes for an ion crystal in a trap have been examined in previous papers and we will not comment further on this special case.

We next turn to the case of a large magnetic field, such that $\Omega i\u226b\omega r$. As shown in Fig. 4, the normal modes now separate into three frequency groupings. There are *N E *×* B* drift modes with low frequencies that scale with magnetic field strength *B* as $1/B$; *N* intermediate frequency axial modes with frequencies that are independent of *B* (for large *B*), and *N* high frequency cyclotron modes with frequencies close to (but slightly larger than) Ω_{i}. If there are several species of charge in the trap with different values of the vortex frequency Ω_{i}, there are modes near each value. The number of modes is the number of particles with that vortex frequency. In general, the frequency difference $\omega \u2212\Omega i$ for these modes scales as $1/B$.

Many characteristics of these strongly magnetized modes have been considered in previous publications.^{10,13,16–18,22,23} Here, we present reduced eigenvalue problems for each mode type and consider a few special cases in more detail.

The eigenmodes in a large magnetic field can be evaluated using degenerate perturbation theory applied to Eq. (5). We break up the dynamical matrix $D$ and the Hamiltonian matrix $H$ into zeroth-order and first order parts. The zeroth-order parts are

and

while the first-order parts are

and

The zeroth-order matrices describe the dynamics of non-interacting charges in a magnetic field. Because the charges are non-interacting, the eigenvectors and eigenvalues can be worked out for each particle separately. The 6*N* dimensional eigenvector $u\omega $ for one of these modes consists of zeros in every element except for those elements corresponding to particle *j*. We label this particular eigenvector $u\omega =uj,\alpha (0)$ where the superscript (0) indicates that it is a zeroth order eigenvector, the subscripts *j* and *α* label the particle and the mode type, respectively. There are five different mode types as we will see in a moment. This eigenvector has the form

The vector $(rj\alpha ,pj\alpha )$ solves the six-dimensional single-particle eigenvalue problem corresponding to the particle *j* elements in $D(0)$:

This eigenvalue problem has five independent eigenvectors, three of which correspond to zero-frequency modes, and two of which are cyclotron modes with frequencies $\xb1\Omega j$. The positive and negative frequency cyclotron eigenvectors for particle *j* are

The negative frequency cyclotron eigenvector is the complex conjugate of the positive frequency eigenvector, as expected from property 3. Both eigenvectors are needed to describe the real phase space vector $z$ corresponding to cyclotron motion for particle *j*. Both modes correspond to rotation of the particle position and momentum vectors in the counterclockwise sense (the positive $\varphi \u0302$ direction).

The three independent zero-frequency eigenvectors for particle *j* are

These eigenvectors correspond to displacements in the *x*, *y*, and *z* directions, respectively, and are labeled as such. The nonzero momentum components arise because the canonical momentum depends on the position. There are corresponding constants of the motion $PjX,PjY,PjZ$. For example, $PjX=(rjX,pjX)\xb7J\xb7(\delta xj,\delta yj,\delta zj,\delta pxj,\delta pyj,\delta pzj)$, and so on. For $\Omega j\u22600$, one can check that $PjX$ and $PjY$ are not in involution,

but $[PjY,PjZ]=[PjX,PjZ]=0$.

These five eigenvectors are *not* sufficient to form a complete set to describe the motion of particle *j*; we require a sixth vector. This problem should be familiar as it was covered in Sec. II B in the discussion of neutrally stable systems. We require a vector that is orthogonal to the other five eigenvectors. It could be found by solution of the linear algebra problem given by Eq. (41), but we can identify the solution without any algebra:

corresponding to a constant velocity in the *z* direction. There is one of these vectors for each particle. We refer to the corresponding 6*N* dimensional vector as $u\xafj,Z$, with zeroes in all elements except for those corresponding to the *j*th particle:

These vectors are not eigenvectors, but instead satisfy

The vectors are orthogonal to all of the zeroth-order eigenvectors $uk\alpha (0)$, both with respect to a standard dot product as well as with respect to the inner product defined by the zeroth order Hamiltonian matrix as

#### 1. Reduced eigenvalue problem for cyclotron modes

We now employ the zeroth-order eigenvectors in order to construct a cyclotron-frequency eigenmode that includes the effect of particle interactions to the first order. Let us assume that there are *N _{a}* particles of species

*a*, all of which have identical vortex frequency Ω

_{a}, and we will assume there are no other particles with this vortex frequency. We will then use degenerate perturbation theory to solve for the perturbed eigenmode $u\omega $, writing it as a superposition of the degenerate zeroth order modes with frequency Ω

_{a}:

where the sum over *j* sums only over particles of species *a*, the coefficients $cj\omega $ are to be determined, and where $u(1)$ is a small correction to the eigenmode. By assumption, this correction is orthogonal to the zeroth-order eigenmodes. (It will turn out that we do not need to calculate $u(1)$.) We substitute this eigenvector into Eq. (5) and write $D=D(0)+D(1)$ to obtain

We drop the last term in this equation because it is second-order, use the fact that $D(0)\xb7uj,+(0)=\u2212i\Omega auj,+(0)$, multiply through by *i*, and take a zeroth-order inner product with respect to one of the eigenmodes $uk,+(0)$. The orthogonality and the Hermitian nature of the matrix $iD(0)$ then annihilate several terms in the equation, leaving us with

The inner product $(uk,+(0),uk,+(0))(0)$ involves only particle *k* and, using Eqs. (156), (160), and (162), evaluates to

[recall that all particles in Eq. (174) are of species *a*]. The inner product $(uk,+(0),iD(1)\xb7uj,+(0))(0)$ involves only particles *j* and *k* in species *a* and evaluates to

where we used the conjugate transpose of Eq. (5) in the third line, Eqs. (158) and (160) in the fourth line, and Eq. (162) in the fifth line, and where $Vkjxx$ is the $x\u0302x\u0302$ component of the potential matrix $Vkj=\u22022\Phi /\u2202Rj\u2202Rk$ and similarly for $Vkjyy$. Equation (174) can then be written in vector form as

where the matrix **F** has components

Thus, the frequency shift $\omega \u2212\Omega a$ is an eigenvalue of the matrix $F$. The matrix is real and symmetric, so its eigenvalues are real and the eigenvectors are also real, and they form a complete orthogonal set.

Equations (178) and (179) are the reduced eigenvalue problem for the cyclotron modes of species *a*. Note that only particles of this species enter the matrix $F$ so the matrix is $Na\xd7Na$, yielding *N _{a}* eigenfrequencies that differ from Ω

_{a}by a small shift, proportional to $1/\Omega a$.

For a stable system the shift is positive; the extra restoring force from the oscillator potentials tends to increase the mode frequencies. One can see this from the following argument: any displacement **r** of the charges from equilibrium must increase the potential energy: $r\xb7V\xb7r>0$ for any real nonzero vector **r**. Consider now a complex displacement $r$ created by a superposition of cyclotron eigenvectors for species *a* particles only. From Eq. (162), each particles position change is complex, of the form $vk(1,i,0)$ for some complex amplitude *v _{k}*. Then, consider the quantity $r*\xb7V\xb7r=2ma\Omega av*\xb7F\xb7v$ where $v$ is the vector of complex coefficients

*v*. However, $r*\xb7V\xb7r>0$ because, for any complex vector $r=a+ib$ (for real $a$ and $b$), $r*\xb7V\xb7r=a\xb7V\xb7a+b\xb7V\xb7b$ (since

_{k}**V**is symmetric) and both terms on the right hand side are positive. Therefore, the quantity $v*\xb7F\xb7v$ must also be greater than zero for any vector

**v**, which implies that the eigenvalues of

**F**must be positive. This follows by writing $v$ as a superposition of the real eigenvectors of

**F**and evaluating $v*\xb7F\xb7v$.

In Fig. 7, we compare the eigenfrequencies evaluated using Eq. (178) to the exact eigenfrequencies obtained using Eq. (5), as a test of the perturbation theory. For the large magnetic field $\Omega =20\omega z$, the fractional error between the exact mode frequencies and the approximate frequencies is quite small.

The energy in cyclotron modes can be determined in terms of the approximate cyclotron frequencies and eigenvectors. [Of course, the energy of any given mode is also determined exactly via Eq. (16).] Consider a phase space displacement due to species *a* cyclotron modes, $z=(1/2)\u2211\omega >0a\omega \u2211j=1Nacj\omega uj+(0)+c.c.$, where $a\omega $ is the amplitude of each mode. This definition, with the extra factor of $1/2$, allows us to identify $|ck\omega |$ as *m _{a}* multiplied by the speed of particle

*k*in a given mode, for unit amplitude $a\omega =1$, because the unit amplitude results in a cyclotron radius for this particle of $|ck\omega |/ma\Omega a$, according to Eqs. (162) and (163). Applying this phase space displacement to the energy, Eq. (1), breaking the Hamiltonian matrix into zeroth and first order parts, and using orthogonality of the cyclotron modes along with Eqs. (175), (176), and (178), the cyclotron energy diagonalizes:

The first term in the parenthesis in Eq. (180) gives the total kinetic energy associated with free particle cyclotron motion: a sum of the squares of the particle kinetic momenta in a given mode, $|a\omega |2c\omega \xb7c\omega $, also summed over the modes, and divided by two times the particle mass. The second term in the parenthesis, proportional to $\omega \xaf\u2212\Omega a$, is the relatively small positive correction to the kinetic energy due to interactions between the charges.

#### 2. Reduced eigenvalue problem for the axial and E × B modes

Let us now turn to the axial modes and the *E *×* B* modes. These modes can also be described using a reduced eigenvalue problem that stems from degenerate perturbation theory applied to Eq. (5). Now, however, we expand an eigenvector $u\omega $ in terms of the zero-frequency eigenvectors along with the extra vectors $u\xafjZ$:

where *X _{j}*,

*Y*, and

_{j}*Z*are displacement amplitudes for particle

_{j}*j*in the

*x*,

*y*, and

*z*directions, respectively,

*P*is the axial momentum of particle

_{j}*j*, and $u(1)$ is a small correction (which we will avoid having to evaluate in what follows). This correction is assumed to satisfy orthogonality conditions

We drop the last term on the right hand side, since it is second order. Also, $D(0)\xb7uj\alpha (0)=0$ for all the zero-frequency eigenmodes, so Eq. (184) simplifies to

where we also applied Eq. (170).

Now, recall from Sec. II B that we can project out zero-frequency modes using the fundamental symplectic matrix $J$ rather than the inner product, since such modes are orthogonal to themselves. Acting on Eq. (185) with $\u2212u\xafkZ(0)\xb7J=ukZ(0)$ and using the orthogonality conditions (183) yields

where we used the identity $J\xb7D=J\xb7J\xb7H=\u2212H$. However, one can use Eqs. (158) and (169) to check that $u\xafkZ(0)\xb7H(1)=0$, which annihilates the sum on the right hand side so Eq. (186) simplifies to

the standard relationship between the axial position and momentum. Now act on Eq. (185) with $ukZ(0)\xb7J=u\xafkZ(0)$, which results in

where the second line applied Eqs. (158), (160), and (164)–(169). The right hand side is the axial force on particle *k* caused by other particle displacements from equilibrium.

Next, act on Eq. (185) with $ukY(0)\xb7J$. Here, we will use $ukY(0)\xb7J\xb7ukX(0)=\u2212mk\Omega k$, which follows from Eqs. (167) and (63) or can be checked directly using Eqs. (160), (164), and (165). This projection results in

The right hand side is the force in the *y* direction on particle *k* caused by displacements of other particles. The force produces an *E *×* B* drift velocity $\u2212i\omega Xk$ in the −*x* direction. Finally, act on Eq. (185) with $ukX(0)\xb7J$. This projection results in

The right hand side is the force in the *x* direction on particle *k*. This force produces an *E *×* B* drift velocity $\u2212i\omega Yk$ in the *y* direction.

Equations (187)–(190) constitute a reduced eigenvalue problem for axial and *E *×* B* modes that has projected out the cyclotron modes. However, the problem still mixes the *E *×* B* modes with the axial plasma modes. One can see from Eqs. (189) and (190) that axial displacements *Z _{j}* are coupled to

*x*and

*y*drifts through the

*x*and

*y*forces such axial displacements can produce. Also, axial accelerations can be caused by

*X*and

*Y*displacements, as seen in Eq. (188).

Now, there are circumstances where this coupling vanishes. For example, when the crystal equilibrium is a single lattice plane in the *z *=* *0 plane, symmetries of this equilibrium imply that $Vkjxz=Vkjyz=0$. The coupling between axial and transverse motions also vanishes when the crystal is a one-dimensional line of charges along the *z* axis. However, for more general crystal equilibria the coupling is nonzero.

#### 3. Decoupling the axial and E × B modes

In order to further decouple the *E *×* B* modes from the axial modes, we must, in general, resort to an asymptotic two-timescale analysis based on the different frequencies of these modes. The *E *×* B* modes, with frequencies scaling as $1/B$, are low frequency compared to the axial modes provided that *B* is sufficiently large. This regime implies that we may write *Z _{j}* as a sum of a slowly evolving and a rapidly evolving contribution, $Zj=Zjslow+Zjfast$. The slow evolution, on the

*E*×

*B*timescale, is conditioned on the axial particle positions being in axial force balance:

We can regard this force balance condition as a set of coupled linear equations for the slow axial displacements, written in vector form as $Vzx\xb7X+Vzy\xb7Y+Vzz\xb7Zslow=0$, where the tensor $Vxz=\u2202\Phi /\u2202X\u2202Z$ and similarly for the other terms. This equation can be solved by matrix inversion,

This is the reduced eigenvalue problem for the axial modes.^{10} When the approximations used in its derivation are poor, the full eigenvalue problem is available for exact results, or the reduced eigenvalue problem consisting of Eqs. (187)–(190) could be employed.

Equation (193) can be combined into a standard generalized eigenvalue problem for $\omega 2$, of the familiar type encountered in the textbook Lagrangian theory of coupled oscillators,

where $Z\omega fast$ is an eigenvector of axial displacements and the matrix $M1$ is the diagonal mass matrix of dimension *N* with diagonal elements $mj,j=1,\u2026,N$. The *N* eigenvalues $\omega 2$ are the squares of the axial mode frequencies, and can be shown to be real and positive for a stable equilibrium, using standard arguments. The eigenvectors $Z\omega fast$ are real and orthogonal with respect to the inner product $(a,b)z=a\xb7M1\xb7b$. In Fig. 8, we compare the frequencies determined using Eq. (194) to those obtained from the exact analysis, as a test of the theory. Just as for the cyclotron modes, the errors are small when the magnetic field is large.

The axial eigenvectors diagonalize the axial energy *H _{z}* where

Writing $Zj=\u2211\omega a\omega Zjfast$, where $a\omega $ is the amplitude of mode *ω*, and applying Eqs. (193) and (194) to Eq. (195), the kinetic and potential terms contribute equally, yielding for the axial energy the expression

Turning to the *E *×* B* modes, we apply Eq. (192) for the slow axial motion to Eqs. (189) and (190). These equations can then be combined into a vector form. Defining a 2*N* dimensional transverse displacement eigenvector $R\u22a5\omega =(X,Y)$, the equations become

where $D\u22a5=G\xb7J\xb7V\u22a5$ is the dynamical matrix for *E *×* B* drift modes, $G$ is a diagonal $2N\xd72N$ matrix with diagonal elements $((m1\Omega 1)\u22121,\u2026,(mN\Omega N)\u22121,(m1\Omega 1)\u22121,\u2026,(mN\Omega N)\u22121)$, $J$ is the $2N\xd72N$ fundamental symplectic matrix [see Eq. (3)], and $V\u22a5$ is the following symmetric $2N\xd72N$ potential energy tensor,

This tensor determines the potential energy *V* in an *E *×* B* displacement of the form $(X,Y,Zslow)$ through the expression

This can be proven by writing *V* as

and applying Eq. (192) along with the symmetry of the matrix elements under interchange of the *x*, *y*, and *z* subscripts.

Equation (197) is the reduced eigenvalue problem for *E *×* B* drift modes. The frequencies all scale as $1/B$ since the dynamical matrix $D\u22a5$ is proportional to $1/B$ through its dependence on **G**. The dynamical matrix has properties that mirror the general Hamiltonian matrices discussed in Sec. II, which in turn determine properties of the eigenmodes. First, the matrix $iD\u22a5$ is Hermitian with respect to the inner product $(a,b)\u22a5=a*\xb7V\u22a5\xb7b$ for any vector **a** and **b**. This follows from the fact that the matrix $L=V\u22a5\xb7G\xb7J\xb7V\u22a5$ is antisymmetric, which in turn follows from the symmetry and antisymmetry, respectively, of $V\u22a5$ and **J**, along with the fact that $G\xb7J=J\xb7G$. The proof follows along the same path as in Eq. (8).

As discussed in relation to properties 1 and 2 in Sec. II, the Hermitian nature of the dynamical matrix implies that the eigenfrequencies are real, provided that the system is at least neutrally stable, and that non-degenerate complex eigenvectors are orthogonal with respect to the above inner product. Also, since the matrix elements of $D\u22a5$ are real, the nonzero frequency modes come in $\xb1\omega $ pairs, as per property 3.

In Fig. 9, we compare the *E *×* B* frequencies determined using Eq. (197) to those obtained from the exact analysis, as a test of the theory. Just as for the cyclotron and axial modes, the errors are small.

The energy of an *E *×* B* mode can be written in terms of the eigenvectors $R\u22a5\omega $. Since the *E *×* B* modes form an orthogonal set, a general *E *×* B* displacement of the form $\u2211\omega >0a\omega (R\u22a5\omega ,Zslow)+c.c.$, where $a\omega $ is the amplitude of a given mode, produces a diagonalized potential energy *V* given by

The kinetic energy contribution to *E *×* B* modes is negligible, scaling as $1/B2$. For a cylindrically symmetric trap potential, there is an additional contribution to the potential energy from a zero frequency mode. The analysis of this contribution follows the same procedure as was developed in Sec. II B. The extra energy from this mode is the potential energy change from radial compression of the plasma due to a change in the rotation rate.

## IV. DISCUSSION

In this paper, we have detailed a method of diagonalizing the Hamiltonian for a general linearized Hamiltonian system. The method relies on the Hermitian properties of the dynamical matrix $D$ and, for a stable system, requires only the evaluation of the eigenmodes of this matrix. For a neutrally stable system, we found that the form of the Hamiltonian depends on whether or not constants of the motion associated with neutral modes are in involution. The normal mode form of the Hamiltonian was also derived for an unstable system.

In applying this formalism to determine the normal modes of a magnetized two-ion Coulomb cluster, we found that for spherically symmetric confinement the rotational inertia of the cluster is undefined, and we related this surprising result to the fact that constants of the motion associated with rotations of this system are not in involution. In this case, a change in angular momentum produces a change in crystal orientation rather than a change in rotation frequency.

Thermal fluctuations were also considered, and in particular, fluctuations in axial position were evaluated in order to make contact with ongoing experiments which depend, in part, on keeping these fluctuations below the level set by the Lamb–Dicke confinement.^{1–5} A somewhat surprising result of our analysis is the dominant contribution of low-frequency *E *×* B* modes to these axial fluctuations in 3D Coulomb crystals.

This paper focused on linear modes, but nonlinear interactions between modes are also important. One simple example of this coupling is the modulation of axial mode frequencies as ions slowly shift positions due to low-frequency *E *×* B* modes. This nonlinear behavior will be the subject of future investigations.^{47}

The fact that modes separate into three disparate frequency groups in a strong magnetic field also has consequences for the energy equilibration in such a system.^{17} When cooling and heating are applied, modes with disparate frequencies can come to quite different equilibrium energies, depending on the details of the driving.^{48} In particular, laser cooling of *E *×* B* modes may not be as efficient as for axial and cyclotron modes because the velocities associated with these modes are low, so their energies may not be well-equilibrated with the other mode branches. In 3D crystals, this problem may be somewhat alleviated by the axial motion associated with some *E *×* B* modes, which can increase coupling between these modes and axial modes. These effects will be subjects of further study.

The general theory presented in Sec. II kept forcing terms **f** in the oscillator equations. This is because time-dependent forcing terms are often present in ion crystal experiments, although they were not considered in Sec. III. In some experiments, oscillating fields are applied in order to purposely interact with and excite certain modes, such as the axial center of mass mode. Forcing can also be caused by non-axisymmetric field errors in the external potential that are static in the laboratory frame and by non-axisymmetric forces from laser beams (also stationary in the lab frame). In the frame of the rotating crystal, these forces oscillate at the rotation frequency *ω _{r}*, which can resonate with normal modes. For example, for the $\beta =3/4$ crystal shown in Fig. 2, when Ω = 0 (the Brillouin limit), the rotation rate is $\omega r=5/4\omega z=1.118\omega z$, which is within the spectrum of mode frequencies (see Fig. 3). It is therefore possible for forcing that is static in the lab frame to resonantly excite normal modes in this crystal to large amplitude. On the other hand, when $\Omega =20\omega z$,

*ω*is either $20.0623\omega z$ or $0.0623\omega z$ in the fast and slow rotation branches, respectively, and both values lie just beyond the frequency spectrum limits shown in Fig. 4. Resonant interactions with modes will lead to nonlinear effects, plasma heating, and mode damping which are beyond the scope of this paper's linear analysis, and probably require a simulation approach. The effect of interaction of modes with resonant external forcing will be further examined in future work.

_{r}## ACKNOWLEDGMENTS

The author acknowledges useful discussions with Dr. John Bollinger, Dr. Matt Affolter, Professor Dan Arovas, and Professor Scott Parker. This work is supported by AFOSR Contract No. FA 9550–19-1–0999, DOE Grant No. DE-SC0018236, and NSF Grant No. PHY1805764.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: THE POTENTIAL MATRIX

In this appendix, we work out the form of the 3 × 3 potential matrix $Vjk=\u22022\Phi /\u2202Rj\u2202Rk$ for charges *j* and *k* with equilibrium positions **R**_{j} and **R**_{k}, respectively. The system potential energy $\Phi $ is given by Eq. (115), and substituting for this, we obtain

Evaluating the derivatives yields the following dyadic form for the matrix:

where **1** is the 3 × 3 unit matrix.

### APPENDIX B: DIAGONALIZING THE ENERGY USING NON-CANONICAL VARIABLES

In order to evaluate the normal modes of oscillation of a Coulomb crystal, we employed a Hamiltonian approach using canonical variables in Sec. III. An alternate approach instead uses particle velocities rather than canonical momenta, but employs most of the same techniques as in the Hamiltonian method. In this Appendix, we outline this non-canonical method. It is closer to the Lagrangian approach used in Ref. 22 and is arguably easier to apply, provided that canonical coordinates are not required.

The linearized equations of motion in the rotating frame of the crystal equilibrium, when written in terms of the perturbed velocities $\delta vj,j=1,\u2026,N$, are

where the 3 × 3 symmetric tensor $Vij=\u2207i\u2207j\Phi $.

In order to solve for the linear normal modes of oscillation of this system, we combine Eqs. (B1) and (B2) into a single vector equation for the 6*N* dimensional configuration vector $\eta =(\delta r1,\u2026,\delta rN,\delta v1,\u2026,\delta vN)$, as was done in the Hamiltonian method. Then, the linearized equations of motion can be written, in analogy to Eq. (2), as

where the $6N\xd76N$ dynamical matrix $D\u2032$ consists of four $3N\xd73N$ blocks,

where the Lorentz tensor $\Omega $, the mass tensor $M$, and the potential tensor $V$ are the same as in Eq. (118).

The system energy *E* is a conserved quantity, given by

where the symmetric energy matrix $E$ is

Eigenmodes of the form $\eta (t)=exp\u2009(\u2212i\omega t)\psi \omega $ satisfy the eigenvalue problem

This eigenvalue problem is essentially identical to the secondary eigenvalue problem used to solve the quadratic eigenvalue problem in Ref. 22.

This modified non-canonical dynamical matrix still has the property that $iD\u2032$ is Hermitian, with respect to a modified inner product defined by $(a,b)\u2032=a*\xb7E\xb7b$ for vectors $a$ and $b$. This can be proven in an analogous manner to the proof for the Hamiltonian problem in Eqs. (7) and (8). Consider the matrix $L=E\xb7D\u2032$. One can show that this matrix is antisymmetric by direct calculation using Eqs. (B4) and (B6):

The dot product of the diagonal matrix $M$ and the antisymmetric block-diagonal matrix Ω is clearly antisymmetric, and therefore $L$ is antisymmetric; also, it is real. This implies that $a*\xb7E\xb7iD\u2032\xb7b=[b*\xb7E\xb7iD\u2032\xb7a]*$, so $iD\u2032$ is Hermitian. Therefore, the eigenmodes satisfy properties 1, 2, and 3 of Sec. II. The eigenvectors form an orthogonal set with respect to the inner product $(a,b)\u2032$, and so on.

The energy of the system is diagonalized by the eigenmodes. For simplicity, we consider only the case where the system is stable with no zero-frequency modes. Then, we may write a general phase space configuration $\eta $ in terms of the complete set of orthogonal eigenvectors

where $a\omega $ is the complex amplitude of mode *ω*. As mentioned before, property 3 of Sec. II together with the real nature of *η* implies that $a\u2212\omega =a\omega *$. Applying this to the energy in Eq. (B5) and using orthogonality and property 3 then yields the diagonalized energy

### APPENDIX C: ROTATIONAL INERTIA FOR A COULOMB CRYSTAL

In this appendix, we evaluate the rotational inertia of a Coulomb crystal consisting of identical ions. We then specialize the result to a quadrupolar trap in which the ions are confined in the *x*–*y* plane.

In order to evaluate the rotational inertia $(u\xaf0z,u\xaf0z)$, we require a solution for the vector $u\xaf0z$ of the equation

where $u0z$ is given by Eq. (124) and the Hamiltonian matrix **H** is given by Eq. (118). Taking advantage of the block form of **H**, we write $u\xaf0z=(r\xaf,p\xaf)$, which when used in Eq. (C1) yields two coupled equations

where $R\u22a5$ is the projection of the equilibrium positions $R$ onto the *x*–*y* plane, and where we have also imposed the assumption of a single species plasma for simplicity, but have not yet assumed a single-plane structure to the equilibrium. Using Eqs. (119) and (120) for $\Omega $, and the definition of $C$, we can write these equations as

Taking a cross-product of Eq. (C5) then yields

Applying this result to Eq. (C4) implies

The left hand side is the electrostatic force due to a displacement $r\xaf$ of the charges. The equation requires that this force must be purely radial. Once a solution is obtained, then $p\xaf$ is determined by Eq. (C5)

When these results are used to calculate the rotational inertia $u\xaf0z\xb7H\xb7u\xaf0z$, the result is

The first term is the usual kinetic rotational inertia of a rigid body consisting of identical masses and the second term is the extra inertia associated with potential energy from compression of the crystal.

In general, Eq. (C7) must be solved numerically, but for a single plane equilibrium in a quadrupole trap, the solution is available analytically. In this case, it is well-known that a purely radial perturbation in the position of each charge produces a radial restoring force, as required (this occurs in the radial breathing mode).^{10,22} According to Appendix B in Ref. 10, the restoring force from a radial expansion is $\u2212V\xb7R\u22a5=\u22123m\omega \u22a52R\u22a5$. Therefore, the solution of Eq. (C7) is

and together with Eq. (C5) this implies

### APPENDIX D: THE BOSONIC BOGOLIUBOV METHOD REVISITED

In this appendix, we review the Bogoliubov method and show that it is equivalent to the classical Hermitian method used in Sec. II, while describing a version of the Bogoliubov method that more closely follows the Hermitian method.

The Bogoliubov method for diagonalizing a linearized Hamiltonian is couched in terms of “creation and annihilation” pairs $\psi =(c,c*)$ rather than the phase space coordinates $z=(q,p)$. In terms of these pairs, the Hamiltonian for a linearized system is

where the matrix $h$ is Hermitian, $h\u2020=h$, and has the symmetric block form

where the *N *×* N* matrix $A=A\u2020$ is Hermitian and the *N *×* N* matrix $B=Btr$ is symmetric. The creation and annihilation pairs are related to $z$ via the linear transformation

where

or in the component form $c=(q+ip)/2,c*=(q\u2212ip)/2$. These creation/annihilation pairs have Poisson bracket relationships that may be succinctly expressed by the equation

where the matrix *σ* is defined as

Equation (D5) is equivalent to, and follows from, the Poisson bracket relationships $[z,z]=J$. The relationship between the Hamiltonian matrix $H$ of Eq. (1) and the matrix $h$ is found by substitution of Eq. (D3) into Eq. (D1), yielding

The standard Bogoliubov approach is to find a transformation to new creation–annihilation pairs $\varphi =(a,a*)$,

for some matrix $s$ such that the Hamiltonian is diagonalized:

with the new Hamiltonian matrix $k$ given by

The added requirement that the transformation be canonical requires $[\varphi ,\varphi *]=\u2212i\sigma $, which using Eq. (D8) implies that $s$ must satisfy the symplectic condition in the creation/annihilation representation,

Equations (D10) and (D11) are the Bogoliubov equations for $s$, whose solution provides the canonical transformation that diagonalizes the Hamiltonian. In what follows, we solve for $s$ using an approach similar to that used in Sec. II.

Consider the equation of motion for *ψ* that follows from Hamiltonian (D1):

where we introduce the quantum dynamical matrix $d\u2261\sigma \xb7h$. This dynamical matrix is analogous to the matrix $D$ introduced in Sec. II. Now, consider eigenmodes of the form $\psi =exp\u2009(\u2212i\omega t)w\omega $ for a mode of frequency *ω*, where $w\omega $ is the associated vector. Equation (D12) implies that these vectors are eigenvectors of $d$ with eigenvalues *ω*:

These eigenvectors are related to the eigenvectors $u\omega $ of the dynamical matrix $D$ by Eq. (D3), which implies that

We can now prove three properties of these eigenmodes that are directly analogous to the three properties described in Sec. II:

The eigenvectors $w\omega $ form an orthogonal set with respect to a generalized inner product defined for any complex vectors $a$ and $b$ as $(a,b)\u2261a*\xb7h\xb7b$:

$(w\omega ,w\omega \xaf)=0$ provided that $\omega \u2260\omega \xaf*$.

A given eigenvalue

*ω*is real, provided that the corresponding eigenvector satisfies $(w\omega ,w\omega )\u22600$.- For each eigenmode $(\omega ,w\omega )$ for which $\omega \u22600$, there is a second eigenmode $(\u2212\omega *,w\u2212\omega *)$ for which $w\u2212\omega *=\Lambda \xb7w\omega *$, where the matrix $\Lambda $ is(D15)$\Lambda =(0110).$

Thus, for real *ω*, the $\omega \u22600$ eigenmodes come in $\xb1\omega $ pairs.

As described in Sec. II, the first two properties are a consequence of the spectral theorem for Hermitian matrices. Here, as mentioned before, the quantum dynamical matrix $d$ is Hermitian with respect to the above inner product:

This requires that the matrix $l\u2261h\xb7d$ is a Hermitian matrix: $l=l\u2020$, which can be proven using the same set of steps as shown in Eq. (8):

The third property takes a bit more work than in Sec. II. Here, we use the following property of the quantum Hamiltonian matrix $h$ that follows from its special form, Eq. (D2):

Acting on both sides of the equation with $\sigma $ and using the identity $\sigma \xb7\sigma =1$ gives

where we used the identity $\sigma \xb7\Lambda \xb7\sigma =\u22121$. Substituting Eq. (D19) into Eq. (D13), acting on both sides with $\Lambda $ and using $\Lambda \xb7\Lambda =1$ then yields

The complex conjugate of this equation proves property 3.

We can now diagonalize the Hamiltonian using these eigenvectors, proceeding as described in Sec. II. We will assume for simplicity that $(w\omega ,w\omega )\u22600$ for all eigenmodes, so that all eigenfrequencies are real and are also nonzero and non-degenerate. First, we write a general vector $\psi =(c,c*)$ in terms of a linear combination of the eigenvectors $w\omega $:

This is merely another way to express the Bogoliubov transformation Eq. (D8), taking the components of the vector $\varphi (t)$ to be $a\omega (t)$ and the transformation matrix $s$ to have columns given by the eigenvectors:

Substituting Eq. (D21) into Eq. (D1) and using orthogonality of the eigenvectors (property 1) lead immediately to the diagonal form

To make further progress, we order the eigenfrequencies such that frequencies *ω*_{1} to *ω _{N}* are greater than zero, and the next set of

*N*frequencies are their paired opposites as per property 3. This implies that $s$ takes the form

and that we can write Eq. (D21) as

In order for this equation to match Eq. (D8) with $\varphi =(a,a*)$, this requires

which implies that $\varphi $ has the required form

We can use this result to simplify the Hamiltonian, summing only over positive frequencies $\omega 1,\u2026,\omega N$:

This can be further simplified using the identity

which yields the simplified diagonalized Hamiltonian

The identity can be proven with the aid of the complex conjugate of Eq. (D18):

where we used $\Lambda tr=\Lambda $ in the first step and $h\u2020=h$ in the last step.

Next, we ensure that the transformation to the new $\varphi $ variables is canonical by requiring that their Poisson brackets satisfy $[\varphi ,\varphi *]=\u2212i\sigma $. In the component form, this requires $[a\omega ,a\omega \xaf*]=\u2212i\delta \omega ,\omega \xaf$ and $[a\omega ,a\omega \xaf]=0$ (for *ω* and $\omega \xaf$ both greater than zero). We satisfy these equations in the same way as described in Sec. II. Equation (D21) and orthogonality of the modes imply that

Applying this to $[a\omega ,a\omega \xaf*]$, we use Eq. (D5) to obtain

where in the first step we used $h\u2020=h$, in the third step we used Eq. (D13), and in the last step we used the orthogonality of the eigenmodes (property 1). A similar argument [see Eqs. (18) and (19)] shows that $[a\omega ,a\omega \xaf]=0$ when *ω* and $\omega \xaf$ are greater than zero. Thus, the transformation is canonical, provided that we normalize the eigenvectors so that

a result analogous to Eq. (21). Applying this to Eq. (D30) yields the diagonalized Hamiltonian in the canonical form,

We can more directly connect this approach to the standard Bogoliubov approach by applying the vector $w\omega \xaf*\xb7\sigma $ to Eq. (D13):

Using $\sigma \xb7\sigma =1$ and $d=\sigma \xb7h$ then allows the left hand side to be written as an inner product:

The orthogonality of the eigenmodes then yields

For $\omega >0$, we can substitute for the inner product using Eq. (D34), while for $\omega <0$ we can employ Eq. (D29) to see that $(w\omega ,w\omega )=\u2212\omega ,\u2009\u2009\omega <0$, which when used in Eq. (D38) yields

However, since the eigenvectors are columns of the transformation matrix $s$ [see Eq. (D24)], this equation is equivalent to

Taking the complex conjugate of this equation and using $\sigma tr=\sigma $ yields the symplectic condition, Eq. (D11). Thus, our transformation matrix $s$, given by Eq. (D24) along with eigenvector normalizations (D34), is of the required symplectic form, and also diagonalizes the Hamiltonian as required by the Bogoliubov equations.

Finally, we can rederive the symplectic transformation $S$ between $z$ and $Z$ [see Eq. (27)] by employing Eq. (D3) along with $\varphi =T\xb7Z$ in Eq. (D8), which implies