Normal modes, rotational inertia, and thermal fluctuations of trapped ion crystals

Dubin, Daniel H. E.

doi:10.1063/5.0021732

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 ≫ 1$ ⁠. 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,²⁴ 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 - V$ ⁠, where T and V are kinetic and potential energies, respectively.²⁵ 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²⁸ 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 = (q_{1}, \dots, q_{N})$ and associated canonical momenta $p = (p_{1}, \dots, p_{N})$ whose Hamiltonian has the quadratic form

H (z, t) = \frac{1}{2} z \cdot H \cdot z + f (t) \cdot z,

(1)

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 $H_{i j} = H_{j i}$ ⁠.

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

\begin{matrix} \dot{z} = [z, H] \\ = [z, z] \cdot \frac{\partial H}{\partial z} \\ = J \cdot H \cdot z + J \cdot f \\ = D \cdot z + J \cdot f, \end{matrix}

(2)

where we introduce the dynamical matrix $D = J \cdot H$ as well as the fundamental symplectic matrix²⁸ $J \equiv [z, z]$ ⁠. The fundamental symplectic matrix is given, in block form, by

J = (\begin{matrix} 0, & 1 \\ - 1, & 0 \end{matrix}),

(3)

where 1 and 0 are the unit and zero tensors, respectively. The matrix is antisymmetric and expresses the basic Poisson bracket relationships $[q_{i}, q_{j}] = [p_{i}, p_{j}] = 0$ and $[q_{i}, p_{j}] = - [p_{j}, q_{i}] = δ_{i j}$ ⁠, where $δ_{i j}$ 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

z (t) = exp (- i ω t) u_{ω},

(4)

for some (possibly complex) frequency ω and some (time-independent, possibly complex) vector $u_{ω}$ ⁠. Substituting Eq. (4) into Eq. (2) and assuming $f = 0$ then yields an eigenvalue problem for ω and $u_{ω}$ ⁠,

- i ω u_{ω} = D \cdot u_{ω} .

(5)

Thanks to the Hamiltonian nature of the linear dynamical equations, the eigenfrequencies ω and eigenvectors $u_{ω}$ of Eq. (5) have the following three properties:

The eigenvectors $u_{ω}$ form an orthogonal set with respect to a generalized inner product defined for any complex vectors a and b as $(a, b) \equiv a^{*} \cdot H \cdot b$ ⁠:
$(u_{ω}, u_{\bar{ω}}) = 0$ ⁠, provided that $ω \neq {\bar{ω}}^{*}$ ⁠.
A given eigenvalue ω is real, provided that the corresponding eigenvector satisfies $(u_{ω}, u_{ω}) \neq 0$ ⁠.
For each eigenmode $(ω, u_{ω})$ for which $ω \neq 0$ ⁠, there is a second eigenmode $(- ω^{*}, u_{- ω^{*}})$ for which $u_{- ω^{*}} = u_{ω}^{*}$ ⁠. Thus, for real ω the $ω \neq 0$ eigenmodes come in $\pm ω$ 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

- i (- ω^{*}) u_{ω}^{*} = D \cdot u_{ω}^{*},

(6)

showing that $u_{ω}^{*}$ is also an eigenvector of D with eigenfrequency $- ω^{*}$ ⁠, which completes the proof of property 3.

Properties 1 and 2 follow from the fact that $i D$ 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.³⁴

The Hermitian property of the matrix $i D$ is defined by the relationship $(a, i D \cdot b) = {(b, i D \cdot a)}^{*}$ ⁠, 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^{*} \cdot H \cdot D \cdot b = - b \cdot H \cdot D \cdot a^{*} .

(7)

(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 \equiv H \cdot D = H \cdot J \cdot H$ that appears in the above expression. This matrix is antisymmetric: $L_{j i} = - L_{i j}$ ⁠. This follows from the symmetry and antisymmetry, respectively, of the matrices H and J:

L_{j i} = H_{j k} J_{k l} H_{l i} = H_{k j} (- J_{l k}) H_{i l} = - H_{i l} J_{l k} H_{k j} = - L_{i j} .

(8)

The antisymmetry of L proves Eq. (7), which in turn proves that $i D$ 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_{ω}, u_{ω}) \neq 0$ 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 (⁠ $\bar{ω} = ω$ 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 2N eigenvectors $u_{ω}$ then form an orthogonal set in the 2N 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:

z (t) = \sum_{ω} a_{ω} (t) u_{ω},

(9)

where the complex amplitude $a_{ω} (t)$ associated with each eigenvector $u_{ω}$ can be found by taking an inner product of both sides of Eq. (9), applying property 1 (orthogonality) of the eigenvectors:

a_{ω} (t) = \frac{(u_{ω}, z (t))}{(u_{ω}, u_{ω})} .

(10)

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

z (t) = \sum_{ω > 0} (a_{ω} (t) u_{ω} + a_{- ω} (t) u_{ω}^{*}) .

(11)

The real nature of the vector z then implies that $a_{- ω} (t) = a_{ω}^{*} (t)$ and we can then write Eq. (11) as

z (t) = \sum_{ω > 0} a_{ω} (t) u_{ω} + c . c .,

(12)

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

A differential equation for the time evolution of the complex mode amplitude $a_{ω} (t)$ follows by substitution of Eq. (9) into Eq. (2). Taking an inner product with respect to one of the eigenvectors $u_{ω}$ then yields

{\dot{a}}_{ω} (t) = - i ω a_{ω} (t) + f_{ω} (t),

(13)

where $f_{ω} (t) = (u_{ω}, J \cdot f (t)) / (u_{ω}, u_{ω})$ ⁠. The forcing coefficient $f_{ω}$ can also be written as

f_{ω} (t) = \frac{u_{ω}^{*} \cdot H \cdot J \cdot f (t))}{(u_{ω}, u_{ω})} = \frac{- i ω u_{ω}^{*} \cdot f (t)}{(u_{ω}, u_{ω})},

(14)

where we used $D = J \cdot H$ 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_{ω} (t)$ by substituting Eq. (9) into Eq. (1):

\begin{matrix} H = \frac{1}{2} \sum_{ω} \sum_{\bar{ω}} a_{ω} (t) a_{\bar{ω}} (t) u_{\bar{ω}} \cdot H \cdot u_{ω} + \sum_{ω} a_{ω} (t) f \cdot u_{ω} \\ = \frac{1}{2} \sum_{ω} \sum_{\bar{ω}} a_{ω} (t) a_{- \bar{ω}} (t) u_{\bar{ω}}^{*} \cdot H \cdot u_{ω} + \sum_{ω} a_{ω} (t) f \cdot u_{ω} \\ = \frac{1}{2} \sum_{ω} a_{ω} (t) a_{ω}^{*} (t) u_{ω}^{*} \cdot H \cdot u_{ω} + \sum_{ω} a_{ω} (t) f \cdot u_{ω} \\ = \sum_{ω > 0} H_{ω}, \end{matrix}

(15)

where in the second line we replaced the dummy index $\bar{ω}$ with $- \bar{ω}$ 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_{ω}$ in each eigenmode, given by

H_{ω} = | a_{ω} (t) |^{2} (u_{ω}, u_{ω}) + 2 f (t) \cdot R e (a_{ω} u_{ω}),

(16)

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_{ω}, u_{ω}) > 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_{ω}$ ⁠, 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_{ω}, a_{\bar{ω}}^{*}]$ ⁠. This bracket can be evaluated using Eq. (10) and the symmetry of the Hamiltonian matrix:

\begin{matrix} [a_{ω}, a_{\bar{ω}}^{*}] = \frac{u_{ω}^{*} \cdot H \cdot [z, z] \cdot H \cdot u_{\bar{ω}}}{(u_{ω}, u_{ω}) (u_{\bar{ω}}, u_{\bar{ω}})} \\ = \frac{u_{ω}^{*} \cdot H \cdot J \cdot H \cdot u_{\bar{ω}}}{(u_{ω}, u_{ω}) (u_{\bar{ω}}, u_{\bar{ω}})} \\ = \frac{u_{ω}^{*} \cdot H \cdot D \cdot u_{\bar{ω}}}{(u_{ω}, u_{ω}) (u_{\bar{ω}}, u_{\bar{ω}})} \\ = \frac{- i \bar{ω} u_{ω}^{*} \cdot H \cdot u_{\bar{ω}}}{(u_{ω}, u_{ω}) (u_{\bar{ω}}, u_{\bar{ω}})} \end{matrix}

(17)

= \frac{- i ω}{(u_{ω}, u_{ω})} δ_{ω \bar{ω}},

(18)

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

[a_{ω}, a_{\bar{ω}}] = 0,

(19)

for ω and $\bar{ω}$ greater than zero. In this case, both eigenvectors in Eq. (17) are starred. However, recall that $u_{\bar{ω}}^{*} = u_{- \bar{ω}}$ ⁠, 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_{ω}$ ⁠, we find it useful to impose the condition on these amplitudes that for $ω > 0$ and $\bar{ω} > 0$ ⁠,

[a_{ω}, a_{\bar{ω}}^{*}] = - i δ_{ω \bar{ω}},

(20)

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

(u_{ω}, u_{ω}) = ω .

(21)

Since both sides of this equation are positive, a normalization constant for $u_{ω}$ 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_{ω}, P_{ω})$ ⁠, defined by

a_{ω} = \frac{1}{\sqrt{2}} (Q_{ω} + i P_{ω}) .

(22)

In order to show that these are canonical pairs, invert Eq. (22) (and its complex conjugate) to give $Q_{ω} = 2^{- 1 / 2} (a_{ω} + a_{ω}^{*})$ and $P_{ω} = - i 2^{- 1 / 2} (a_{ω} - a_{ω}^{*})$ ⁠. Then,

\begin{matrix} [Q_{ω}, P_{\bar{ω}}] = - \frac{i}{2} ([a_{ω}, a_{\bar{ω}}] - {[a_{ω}, a_{\bar{ω}}]}^{*} - [a_{ω}, a_{\bar{ω}}^{*}] + [a_{ω}^{*}, a_{\bar{ω}}]) \\ = - \frac{i}{2} (0 - 0 + i δ_{ω \bar{ω}} + i δ_{ω \bar{ω}}) \\ = δ_{ω \bar{ω}}, \end{matrix}

(23)

and similarly, $[P_{ω}, P_{\bar{ω}}] = 0 = [Q_{ω}, Q_{\bar{ω}}]$ ⁠. We now see the point of Eq. (20): this choice determines that $[Q_{ω}, P_{ω}] = 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

H = \sum_{ω > 0} H_{ω},

(24)

H_{ω} = \frac{ω}{2} (Q_{ω}^{2} + P_{ω}^{2}) + f_{1 ω} Q_{ω} + f_{2 ω} P_{ω},

(25)

where $f_{1 ω} (t) = \sqrt{2} f (t) \cdot R e (u_{ω})$ and $f_{2 ω} (t) = - \sqrt{2} f (t) \cdot I m (u_{ω})$ ⁠. Note that for our choice of canonical pairs, the $Q_{ω}$ and $P_{ω}$ variables have the same dimensions of $\sqrt{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

\begin{matrix} {\dot{Q}}_{ω} = \frac{\partial H}{\partial P_{ω}} = ω P_{ω} + f_{2 ω}, \\ {\dot{P}}_{ω} = - \frac{\partial H}{\partial Q_{ω}} = - ω Q_{ω} - f_{1 ω}, \end{matrix}

(26)

which are seen to be equivalent to Eqs. (13) and (14) after application of Eqs. (22) and Eq. (21).

When $f_{1 ω} = f_{2 ω} = 0$ ⁠, the equations for the ω mode are unforced and the oscillator energy $H_{ω}$ 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

z = S \cdot Z,

(27)

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

S = \sqrt{2} (Re U, - Im U),

(28)

and the $2 N \times N$ matrix U has columns consisting of the N eigenvectors $u_{ω}$ with $ω > 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)$ ⁠,

⟨ F ⟩ = \frac{\int d z F (z) exp (- H (z) / T)}{\int d z exp (- H (z) / T)},

(29)

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_{ω}, P_{ω})$ can simplify such calculations. Since the transformations are canonical, the Jacobian of the transformation is unity and the averaging integrals become

⟨ F ⟩ = \frac{\int F (z) \prod_{ω > 0} d Q_{ω} d P_{ω} exp (- H_{ω} / T)}{\int \prod_{ω > 0} d Q_{ω} d P_{ω} exp (- H_{ω} / T)} .

(30)

For instance, it is easy to show that

⟨ Q_{ω} Q_{\bar{ω}} ⟩ = ⟨ P_{ω} P_{\bar{ω}} ⟩ = \frac{T}{ω} δ_{ω \bar{ω}},

(31)

⟨ Q_{ω} P_{\bar{ω}} ⟩ = 0,

(32)

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

⟨ H_{ω} ⟩ = T .

(33)

Similarly, the classical partition function $Z_{c} \equiv h^{- N} \int d z exp [- H (z) / T]$ (where h is Planck's constant) evaluates to

Z_{c} = \prod_{ω > 0} \frac{2 π T}{h ω},

(34)

the same form as the well-known classical limit for oscillators that do not have the position-momentum coupling considered in this paper.³⁵

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 $u_{0}$ ⁠. This zero-frequency eigenvector satisfies

D \cdot u_{0} = 0.

(35)

Thus, $u_{0}$ 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 $u_{0}$ 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 \cdot H$ ⁠:

\begin{matrix} J \cdot D \cdot u_{0} = J \cdot J \cdot H \cdot u_{0} \\ = - H \cdot u_{0} = 0, \end{matrix}

(36)

where, in the last step, we used the identity $J \cdot J = - 1$ ⁠.

Now, there are only $2 N - 1$ independent eigenvectors of the dynamical matrix D, the $2 (N - 1)$ eigenvectors with non-zero frequencies, and the single zero frequency eigenvector. Since the phase space has dimension 2N, the $2 N - 1$ 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 $2 N - 1$ eigenvectors. We will refer to this vector as ${\bar{u}}_{0}$ ⁠. 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₀. This constant is

P_{0} = u_{0} \cdot J \cdot z .

(37)

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

\begin{matrix} {\dot{P}}_{0} = u_{0} \cdot J \cdot \dot{z} \\ = u_{0} \cdot J \cdot (D \cdot z + J \cdot f) \\ = u_{0} \cdot J \cdot J \cdot (H \cdot z + f) \\ = - u_{0} \cdot H \cdot z - u_{0} \cdot f \\ = 0 \cdot z + 0 = 0, \end{matrix}

(38)

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

We can now find the vector ${\bar{u}}_{0}$ using the constancy of P₀ in the linear dynamics. Each oscillatory eigenmode with $ω \neq 0$ must satisfy $P_{0} = constant$ ⁠, and since P₀ is linear in $z$ the only possibility is $P_{0} = 0$ ⁠. Therefore, the eigenmodes all satisfy

0 = u_{0} \cdot J \cdot u_{ω} .

(39)

Note that this is also true for the ω = 0 eigenvector $u_{0}$ ⁠, since $0 = u_{0} \cdot J \cdot u_{0}$ by the antisymmetry of J. We will construct a vector ${\bar{u}}_{0}$ that is orthogonal to all of the eigenmodes by solving the equation

u_{0} \cdot J \cdot u_{ω} = {\bar{u}}_{0} \cdot H \cdot u_{ω},

(40)

for all eigenvectors $u_{ω}$ ⁠. By Eq. (39), such a vector will satisfy $({\bar{u}}_{0}, u_{ω}) = 0$ for all ω including ω = 0. A necessary and sufficient condition for solution of Eq. (40) is that ${\bar{u}}_{0}$ satisfies

H \cdot {\bar{u}}_{0} = u_{0} \cdot J .

(41)

However, since H has a vector $u_{0}$ in its nullspace, ${\bar{u}}_{0}$ 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 $u_{0} \cdot J \cdot u_{0} = 0$ ⁠, and this implies that only $2 N - 1$ of the equations in Eq. (41) are linearly independent: the system satisfies $u_{0} \cdot H \cdot {\bar{u}}_{0} = 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 $u_{0}$ ⁠, obtaining $2 N - 1$ independent equations. A unique particular solution for ${\bar{u}}_{0}, {\bar{u}}_{p}$ ⁠, can then obtained by specifying an extra condition on the solution that, for example, ${\bar{u}}_{0} \cdot u_{0} = 0$ ⁠. The vector ${\bar{u}}_{p}$ is real, since all coefficients appearing in the equations are real. The general solution is this particular solution added to the nullspace eigenvector:

{\bar{u}}_{0} = {\bar{u}}_{p} + C u_{0},

(42)

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$ ⁠:

z = \sum_{ω \neq 0} a_{ω} u_{ω} + a_{0} u_{0} + {\bar{a}}_{0} {\bar{u}}_{0} .

(43)

The vectors form a complete orthogonal set so that the coefficients $a_{ω}$ and ${\bar{a}}_{0}$ can be obtained by projection:

a_{ω} = \frac{(u_{ω}, z)}{(u_{ω}, u_{ω})},

(44)

{\bar{a}}_{0} = \frac{({\bar{u}}_{0}, z)}{({\bar{u}}_{0}, {\bar{u}}_{0})} .

(45)

However, a₀ cannot be found using the standard projection method because $u_{0}$ is orthogonal to itself: according to Eq. (36), $u_{0} \cdot H \cdot u_{0} = 0$ ⁠. Instead, a₀ can be determined using the properties of the fundamental symplectic matrix. Acting on both sides of Eq. (43) with ${\bar{u}}_{0} \cdot J$ ⁠, we obtain

{\bar{u}}_{0} \cdot J \cdot z = \sum_{ω \neq 0} a_{ω} {\bar{u}}_{0} \cdot J \cdot u_{ω} + a_{0} {\bar{u}}_{0} \cdot J \cdot u_{0} + {\bar{a}}_{0} {\bar{u}}_{0} \cdot J \cdot {\bar{u}}_{0} .

(46)

However, ${\bar{u}}_{0} \cdot J \cdot {\bar{u}}_{0} = 0$ due to the antisymmetry of $J$ ⁠, and also ${\bar{u}}_{0} \cdot J \cdot u_{ω} = 0$ for $ω \neq 0$ ⁠. This follows because, for $ω \neq 0$ ⁠,

{\bar{u}}_{0} \cdot J \cdot u_{ω} = \frac{{\bar{u}}_{0} \cdot J \cdot D \cdot u_{ω}}{- i ω} = \frac{{\bar{u}}_{0} \cdot J \cdot J \cdot H \cdot u_{ω}}{- i ω} = \frac{{\bar{u}}_{0} \cdot H \cdot u_{ω}}{i ω} = 0.

(47)

We are therefore left with

{\bar{u}}_{0} \cdot J \cdot z = a_{0} {\bar{u}}_{0} \cdot J \cdot u_{0} = - a_{0} {\bar{u}}_{0} \cdot H \cdot {\bar{u}}_{0},

(48)

where in the second form we employed Eq. (41). Thus, we obtain for a₀

a_{0} = - \frac{{\bar{u}}_{0} \cdot J \cdot z}{({\bar{u}}_{0}, {\bar{u}}_{0})} .

(49)

Returning to Eq. (43), the complex coefficients $a_{ω}$ still come in $\pm ω$ pairs satisfying $a_{- ω} = a_{ω}^{*}$ for $ω \neq 0$ ⁠. The time-dependence of these coefficients is still given by Eq. (13). The time-dependence of a₀ and ${\bar{a}}_{0}$ follows in the same way, by substitution of Eq. (43) into the equation of motion, Eq. (2). The result, after projecting out all the $ω \neq 0$ eigenvectors, is

\begin{matrix} {\dot{a}}_{0} u_{0} + {\bar{a} ̇}_{0} {\bar{u}}_{0} = {\bar{a}}_{0} D \cdot {\bar{u}}_{0} + Δ f \\ = {\bar{a}}_{0} u_{0} + Δ f, \end{matrix}

(50)

where we employed $D \cdot u_{0} = 0$ and where $Δ f = J \cdot f - \sum_{ω \neq 0} f_{ω} u_{ω}$ is the projection of $J \cdot f$ into the $(u_{0}, {\bar{u}}_{0})$ subspace. In the second line, we used

D \cdot {\bar{u}}_{0} = u_{0},

(51)

which follows from Eq. (41) by applying $J$ to both sides, and using $u_{0} \cdot J = - J \cdot u_{0}$ and $J \cdot J = - 1$ ⁠. Taking an inner product of Eq. (50) with respect to ${\bar{u}}_{0}$ then implies that

\begin{matrix} {\bar{a} ̇}_{0} = \frac{({\bar{u}}_{0}, Δ f)}{({\bar{u}}_{0}, {\bar{u}}_{0})} \\ = - \frac{f \cdot u_{0}}{({\bar{u}}_{0}, {\bar{u}}_{0})}, \end{matrix}

(52)

where in the second line we rewrote the numerator using $({\bar{u}}_{0}, Δ f) = ({\bar{u}}_{0}, J \cdot f) = (J \cdot f, {\bar{u}}_{0}) = - f \cdot J \cdot H \cdot {\bar{u}}_{0} = - f \cdot D \cdot {\bar{u}}_{0} = - f \cdot u_{0}$ ⁠, and in the last step we used Eq. (51).

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

P_{0} = {\bar{a}}_{0} {\bar{u}}_{0} \cdot H \cdot {\bar{u}}_{0},

(53)

where we employed Eqs. (39) and (41). Taking a time derivative and using Eq. (52) yield ${\dot{P}}_{0} = - f \cdot u_{0}$ ⁠.

Finally, the dynamics of a₀ follows by acting on both sides of Eq. (50) with ${\bar{u}}_{0} \cdot J$ ⁠, yielding

{\dot{a}}_{0} = {\bar{a}}_{0} + \frac{{\bar{u}}_{0} \cdot J \cdot Δ f}{{\bar{u}}_{0} \cdot J \cdot u_{0}},

(54)

Substituting for $Δ f$ ⁠, and using $J \cdot J = - 1$ ⁠, Eq. (47), and Eq. (41), we are left with

{\dot{a}}_{0} = {\bar{a}}_{0} + \frac{{\bar{u}}_{0} \cdot f}{({\bar{u}}_{0}, {\bar{u}}_{0})} .

(55)

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

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:

H = \sum_{ω > 0} H_{ω} + \frac{1}{2} {\bar{a}}_{0}^{2} ({\bar{u}}_{0}, {\bar{u}}_{0}) + a_{0} f \cdot u_{0} + {\bar{a}}_{0} f \cdot {\bar{u}}_{0},

(56)

where $H_{ω}$ is still given by Eq. (16). Poisson brackets involving a₀ and ${\bar{a}}_{0}$ follow from Eqs. (49) and (45):

\begin{matrix} [a_{0}, {\bar{a}}_{0}] = - \frac{{\bar{u}}_{0} \cdot J \cdot [z, z] \cdot H \cdot {\bar{u}}_{0}}{{({\bar{u}}_{0}, {\bar{u}}_{0})}^{2}} \\ = - \frac{{\bar{u}}_{0} \cdot J \cdot J \cdot H \cdot {\bar{u}}_{0}}{{({\bar{u}}_{0}, {\bar{u}}_{0})}^{2}} = \frac{{\bar{u}}_{0} \cdot H \cdot {\bar{u}}_{0}}{{({\bar{u}}_{0}, {\bar{u}}_{0})}^{2}} \\ = \frac{1}{({\bar{u}}_{0}, {\bar{u}}_{0})} . \end{matrix}

(57)

This implies that a₀ and the momentum $P_{0} = {\bar{a}}_{0} ({\bar{u}}_{0}, {\bar{u}}_{0})$ form a canonical pair [see Eq. (53)]. Similarly, one can also show that $[a_{0}, a_{ω}] = [P_{0}, a_{ω}] = 0$ for all $ω \neq 0$ ⁠.

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

H = \sum_{ω > 0} {\frac{ω}{2} (Q_{ω}^{2} + P_{ω}^{2}) + f_{1 ω} Q_{ω} + f_{2 ω} P_{ω}} + \frac{P_{0}^{2}}{2 ({\bar{u}}_{0}, {\bar{u}}_{0})} + a_{0} f \cdot u_{0} + P_{0} \frac{f \cdot {\bar{u}}_{0}}{({\bar{u}}_{0}, {\bar{u}}_{0})} .

(58)

The dynamical equations for $Q_{ω}$ and $P_{ω}$ remain unchanged from Eq. (26). The equations of motion for a₀ and P₀ are

{\dot{a}}_{0} = \frac{\partial H}{\partial P_{0}} = \frac{P_{0}}{({\bar{u}}_{0}, {\bar{u}}_{0})} + \frac{f \cdot {\bar{u}}_{0}}{({\bar{u}}_{0}, {\bar{u}}_{0})},

(59)

{\dot{P}}_{0} = - \frac{\partial H}{\partial a_{0}} = - f \cdot u_{0},

(60)

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

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 $u_{01}$ and $u_{02}$ in the nullspace of both the dynamical matrix D and the Hamiltonian matrix H. These eigenvectors produce two constants of the motion, $P_{01}$ and $P_{02}$ ⁠, given by

P_{01} = u_{01} \cdot J \cdot z,

(61)

P_{02} = u_{02} \cdot J \cdot z .

(62)

Two possibilities must be separately considered: (i) $u_{01} \cdot J \cdot u_{02} = 0$ and (ii) $u_{01} \cdot J \cdot u_{02} \equiv J_{12} \neq 0$ ⁠. 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 $[P_{01}, P_{02}]$ ⁠, using Eqs. (61) and (62):

\begin{matrix} [P_{01}, P_{02}] = - u_{01} \cdot J \cdot [z, z] \cdot J \cdot u_{02} \\ = - u_{01} \cdot J \cdot J \cdot J \cdot u_{02} \\ = u_{01} \cdot J \cdot u_{02} = J_{12} . \end{matrix}

(63)

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 ${\bar{u}}_{01}$ and ${\bar{u}}_{02}$ are required in order to describe a general phase-space vector z according to

z = \sum_{ω \neq 0} a_{ω} u_{ω} + a_{01} u_{01} + a_{02} u_{02} + {\bar{a}}_{01} {\bar{u}}_{01} + {\bar{a}}_{02} {\bar{u}}_{02} .

(64)

The need for the extra vectors ${\bar{u}}_{01}$ and ${\bar{u}}_{02}$ can be seen by acting on the equation with $u_{01} \cdot J$ ⁠:

\begin{matrix} u_{01} \cdot J \cdot z = \sum_{ω \neq 0} a_{ω} u_{01} \cdot J \cdot u_{ω} + a_{01} u_{01} \cdot J \cdot u_{01} \\ + a_{02} u_{01} \cdot J \cdot u_{02} + {\bar{a}}_{01} u_{01} \cdot J \cdot {\bar{u}}_{01} + {\bar{a}}_{02} u_{01} \cdot J \cdot {\bar{u}}_{02} . \end{matrix}

(65)

Conservation of $P_{01}$ implies that $u_{01} \cdot J \cdot u_{ω} = 0$ ⁠, while $u_{01} \cdot J \cdot u_{01} = 0$ by symmetry and $u_{01} \cdot J \cdot u_{02} = 0$ by assumption. Repeating the procedure with $u_{02} \cdot J$ ⁠, we have

u_{01} \cdot J \cdot z = {\bar{a}}_{01} u_{01} \cdot J \cdot {\bar{u}}_{01} + {\bar{a}}_{02} u_{01} \cdot J \cdot {\bar{u}}_{02},

(66)

u_{02} \cdot J \cdot z = {\bar{a}}_{02} u_{02} \cdot J \cdot {\bar{u}}_{01} + {\bar{a}}_{02} u_{02} \cdot J \cdot {\bar{u}}_{02} .

(67)

To satisfy these equations for a general vector z, we require nonzero values for both ${\bar{a}}_{01}$ and ${\bar{a}}_{02}$ ⁠, proving that both vectors ${\bar{u}}_{01}$ and ${\bar{u}}_{02}$ are required for a complete set.

Let us now consider the solution of these equations for ${\bar{a}}_{01}$ and ${\bar{a}}_{02}$ ⁠, along with the determination of $a_{01}$ and $a_{02}$ ⁠. Recall that these latter two coefficients can be found by applying ${\bar{u}}_{01} \cdot J$ and ${\bar{u}}_{02} \cdot J$ to Eq. (64). Noting that ${\bar{u}}_{01} \cdot J \cdot u_{ω} = {\bar{u}}_{02} \cdot J \cdot u_{ω} = 0$ [see Eq. (47)], and that ${\bar{u}}_{01} \cdot J \cdot {\bar{u}}_{01} = {\bar{u}}_{02} \cdot J \cdot {\bar{u}}_{02} = 0$ by symmetry, we are left with

{\bar{u}}_{01} \cdot J \cdot z = a_{01} {\bar{u}}_{01} \cdot J \cdot u_{01} + a_{02} {\bar{u}}_{01} \cdot J \cdot u_{02} + {\bar{a}}_{02} {\bar{u}}_{01} \cdot J \cdot {\bar{u}}_{02},

(68)

{\bar{u}}_{02} \cdot J \cdot z = a_{01} {\bar{u}}_{02} \cdot J \cdot u_{01} + a_{02} {\bar{u}}_{02} \cdot J \cdot u_{02} - {\bar{a}}_{01} {\bar{u}}_{01} \cdot J \cdot {\bar{u}}_{02} .

(69)

We can simplify the solution of Eqs. (66)–(69) by recalling that the vectors ${\bar{u}}_{01}$ and ${\bar{u}}_{02}$ are constructed to be orthogonal to all of the eigenvectors by solution of the (underdetermined) equations

H \cdot {\bar{u}}_{01} = u_{01} \cdot J,

(70)

H \cdot {\bar{u}}_{02} = u_{02} \cdot J .

(71)

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

{\bar{u}}_{01} \cdot H \cdot z = {\bar{a}}_{01} {\bar{u}}_{01} \cdot H \cdot {\bar{u}}_{01} + {\bar{a}}_{02} {\bar{u}}_{01} \cdot H \cdot {\bar{u}}_{02},

(72)

{\bar{u}}_{02} \cdot H \cdot z = {\bar{a}}_{02} {\bar{u}}_{02} \cdot H \cdot {\bar{u}}_{01} + {\bar{a}}_{02} {\bar{u}}_{02} \cdot H \cdot {\bar{u}}_{02},

(73)

{\bar{u}}_{01} \cdot J \cdot z = - a_{01} {\bar{u}}_{01} \cdot H \cdot {\bar{u}}_{01} - a_{02} {\bar{u}}_{01} \cdot H \cdot {\bar{u}}_{02} + {\bar{a}}_{02} {\bar{u}}_{01} \cdot J \cdot {\bar{u}}_{02},

(74)

{\bar{u}}_{02} \cdot J \cdot z = - a_{01} {\bar{u}}_{02} \cdot H \cdot {\bar{u}}_{01} - a_{02} {\bar{u}}_{02} \cdot H \cdot {\bar{u}}_{02} - {\bar{a}}_{01} {\bar{u}}_{01} \cdot J \cdot {\bar{u}}_{02} .

(75)

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 ${\bar{u}}_{p 1}$ and ${\bar{u}}_{p 2}$ of these equations:

{\bar{u}}_{01} = {\bar{u}}_{p 1} + α_{1} u_{01} + α_{2} u_{02},

(76)

{\bar{u}}_{02} = {\bar{u}}_{p 2} + β_{1} u_{01} + β_{2} u_{02} .

(77)

It is useful to choose the constants $α_{1}, α_{2}, β_{1}, β_{2}$ so that ${\bar{u}}_{01} \cdot J \cdot {\bar{u}}_{02} = 0$ ⁠. This can be accomplished, for example, by taking $α_{1} = β_{1} = β_{2} = 0$ and choosing α₂ such that ${\bar{u}}_{p 1} \cdot J \cdot {\bar{u}}_{02} + α_{2} u_{02} \cdot J \cdot {\bar{u}}_{02} = 0$ ⁠. Using Eq. (71), the solution for α₂ is

α_{2} = - \frac{{\bar{u}}_{p 1} \cdot J \cdot {\bar{u}}_{02}}{{\bar{u}}_{02} \cdot H \cdot {\bar{u}}_{02}} .

(78)

The condition ${\bar{u}}_{01} \cdot J \cdot {\bar{u}}_{02} = 0$ implies that we can drop the ${\bar{a}}_{01}$ and ${\bar{a}}_{02}$ terms in Eqs. (74) and (75), so they become

{\bar{u}}_{01} \cdot J \cdot z = - a_{01} {\bar{u}}_{01} \cdot H \cdot {\bar{u}}_{01} - a_{02} {\bar{u}}_{01} \cdot H \cdot {\bar{u}}_{02},

(79)

{\bar{u}}_{02} \cdot J \cdot z = - a_{01} {\bar{u}}_{02} \cdot H \cdot {\bar{u}}_{01} - a_{02} {\bar{u}}_{02} \cdot H \cdot {\bar{u}}_{02} .

(80)

Equations (79) and (80) can be solved for $a_{01}$ and $a_{02}$ ⁠, while Eqs. (72) and (73) can be solved for ${\bar{a}}_{01}$ and ${\bar{a}}_{02}$ ⁠. The solutions are

a_{01} = \frac{h_{12} j_{2 z} - h_{22} j_{1 z}}{h_{11} h_{22} - h_{12}^{2}}, a_{02} = \frac{h_{12} j_{1 z} - h_{11} j_{2 z}}{h_{11} h_{22} - h_{12}^{2}},

(81)

{\bar{a}}_{01} = \frac{h_{22} h_{1 z} - h_{12} h_{2 z}}{h_{11} h_{22} - h_{12}^{2}}, {\bar{a}}_{02} = \frac{h_{11} h_{2 z} - h_{12} h_{1 z}}{h_{11} h_{22} - h_{12}^{2}},

(82)

where $h_{i j} = {\bar{u}}_{0 i} \cdot H \cdot {\bar{u}}_{0 j}, h_{i z} = {\bar{u}}_{0 i} \cdot H \cdot z$ ⁠, and $j_{i z} = {\bar{u}}_{0 i} \cdot J \cdot z$ ⁠. We assume throughout that the Hamiltonian satisfies $h_{11} h_{22} - h_{12}^{2} \neq 0$ ⁠. (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.)

The system energy can be evaluated by substituting Eq. (64) into Eq. (1), yielding

\begin{matrix} H = \sum_{ω > 0} H_{ω} + \frac{1}{2} h_{11} {\bar{a}}_{01}^{2} + \frac{1}{2} h_{22} {\bar{a}}_{02}^{2} + h_{12} {\bar{a}}_{01} {\bar{a}}_{02} \\ + a_{01} f \cdot u_{01} + a_{02} f \cdot u_{02} + {\bar{a}}_{01} f \cdot {\bar{u}}_{01} + {\bar{a}}_{02} f \cdot {\bar{u}}_{02} . \end{matrix}

(83)

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 $[a_{01}, a_{02}] = [{\bar{a}}_{01}, {\bar{a}}_{02}] = 0$ and the nontrivial brackets

\begin{matrix} [a_{01}, {\bar{a}}_{01}] = \frac{h_{22}}{h_{11} h_{22} - h_{12}^{2}}, \\ [a_{02}, {\bar{a}}_{02}] = \frac{h_{11}}{h_{11} h_{22} - h_{12}^{2}}, \\ [a_{01}, {\bar{a}}_{02}] = - \frac{h_{12}}{h_{11} h_{22} - h_{12}^{2}} . \end{matrix}

(84)

Canonical pairs can be found by noting that the constants of the motion $P_{01}$ and $P_{02}$ are related to ${\bar{a}}_{01}$ and ${\bar{a}}_{02}$ via

P_{01} = h_{11} {\bar{a}}_{01} + h_{12} {\bar{a}}_{02},

(85)

P_{02} = h_{22} {\bar{a}}_{02} + h_{12} {\bar{a}}_{01},

(86)

where we substituted Eq. (64) into Eqs. (61) and (62) and applied Eqs. (70) and (71).

Using these variables along with Eq. (84), it is then an exercise to show that $[a_{01}, P_{01}] = [a_{02}, P_{02}] = 1$ while $[a_{02}, P_{01}] = [a_{01}, P_{02}] = [P_{01}, P_{02}] = 0$ ⁠. Thus, the canonical pairs are $(a_{01}, P_{01})$ and $(a_{02}, P_{02})$ ⁠. We now need only invert Eqs. (85) and (86),

{\bar{a}}_{01} = \frac{h_{22} P_{01} - h_{12} P_{02}}{h_{11} h_{22} - h_{12}^{2}},

(87)

{\bar{a}}_{02} = \frac{h_{11} P_{02} - h_{12} P_{01}}{h_{11} h_{22} - h_{12}^{2}},

(88)

and employ these results in the Hamiltonian, which becomes

\begin{matrix} H = \sum_{ω > 0} H_{ω} + \frac{1}{2} \frac{h_{22} P_{01}^{2} + 2 h_{12} P_{01} P_{02} + h_{11} P_{02}^{2}}{h_{11} h_{22} - h_{12}^{2}} \\ + \frac{h_{22} P_{01} - h_{12} P_{02}}{h_{11} h_{22} - h_{12}^{2}} f \cdot {\bar{u}}_{01} + \frac{h_{11} P_{02} - h_{12} P_{01}}{h_{11} h_{22} - h_{12}^{2}} f \cdot {\bar{u}}_{02} + a_{01} f \cdot u_{01} + a_{02} f \cdot u_{02} . \end{matrix}

(89)

The equations of motion for $P_{01}$ and $P_{02}$ are then

{\dot{P}}_{01} = - \frac{\partial H}{\partial a_{01}} = - f \cdot u_{01},

(90)

{\dot{P}}_{02} = - \frac{\partial H}{\partial a_{02}} = - f \cdot u_{02},

(91)

and the equations of motion for the amplitudes $a_{01}$ and $a_{02}$ are

{\dot{a}}_{01} = \frac{\partial H}{\partial P_{01}} = \frac{h_{22} P_{01} + h_{12} P_{02}}{h_{11} h_{22} - h_{12}^{2}} + \frac{h_{22} f \cdot {\bar{u}}_{01} - h_{12} f \cdot {\bar{u}}_{02}}{h_{11} h_{22} - h_{12}^{2}},

(92)

{\dot{a}}_{02} = \frac{\partial H}{\partial P_{02}} = \frac{h_{11} P_{02} + h_{12} P_{01}}{h_{11} h_{22} - h_{12}^{2}} + \frac{h_{11} f \cdot {\bar{u}}_{02} - h_{12} f \cdot {\bar{u}}_{01}}{h_{11} h_{22} - h_{12}^{2}} .

(93)

When the forcing f is zero, the Hamiltonian is independent of $a_{01}$ and $a_{02}$ ⁠, hence the canonical momenta $P_{01}$ and $P_{02}$ are constant, as expected, and $a_{01}$ and $a_{02}$ both have uniform rates of change.

3. Constants of the motion not in involution

When the two zero-frequency modes satisfy $J_{12} = u_{01} \cdot J \cdot u_{02} \neq 0$ ⁠, the constants of the motion $P_{01}$ and $P_{02}$ 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

z = \sum_{ω \neq 0} a_{ω} u_{ω} + a_{01} u_{01} + a_{02} u_{02} .

(94)

No extra vectors ${\bar{u}}_{01}$ or ${\bar{u}}_{02}$ 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 $u_{01}$ with Eq. (72):

0 = u_{01} \cdot H \cdot {\bar{u}}_{02} = u_{02} \cdot J \cdot u_{01} = - J_{12} \neq 0,

(95)

a contradiction.

The amplitude coefficients a₀₁ and a₀₂ in Eq. (94) can now be obtained by acting with $u_{02} \cdot J$ and $u_{01} \cdot J$ ⁠, respectively, yielding

a_{01} = - \frac{u_{02} \cdot J \cdot z}{J_{12}} = - \frac{P_{02}}{J_{12}},

(96)

a_{02} = \frac{u_{01} \cdot J \cdot z}{J_{12}} = \frac{P_{01}}{J_{12}},

(97)

where we used $u_{01} \cdot J \cdot u_{ω} = u_{02} \cdot J \cdot u_{ω} = 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_{ω}$ are still obtained with the usual inner product, see Eq. (10).

The Hamiltonian is found by applying Eq. (94) to Eq. (1), yielding

\begin{matrix} H = \sum_{ω > 0} H_{ω} + a_{01} f \cdot u_{01} + a_{02} f \cdot u_{02} \\ = \sum_{ω > 0} H_{ω} - \frac{P_{02}}{J_{12}} f \cdot u_{01} + \frac{P_{01}}{J_{12}} f \cdot u_{02}, \end{matrix}

(98)

where in the second line we employed Eqs. (96) and (97). The equations of motion for $P_{01}$ and $P_{02}$ then follow from Eq. (63):

{\dot{P}}_{01} = [P_{01}, H] = J_{12} \frac{\partial H}{\partial P_{02}} = - f \cdot u_{01},

(99)

{\dot{P}}_{02} = [P_{02}, H] = - J_{12} \frac{\partial H}{\partial P_{01}} = - f \cdot u_{02},

(100)

which are the same as when P₀₁ and P₀₂ 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₀₁ and a₀₂ are time-independent [see Eqs. (96) and (97)]. This behavior differs from the previous case, where a₀₁ and a₀₂ 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 $ω = Ω \equiv Ω_{r} + i γ$ ⁠, where $Ω_{r} > 0$ and $γ > 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 $Ω^{*} = Ω_{r} - i γ$ 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, $Ω_{r} \pm i γ$ ⁠.

Note that the above mode with frequency $Ω^{*}$ has a negative growth rate $- γ$ and is independent of the mode required by property 3, with frequency $- Ω^{*}$ ⁠. The latter mode has a positive growth rate and is needed (when $Ω_{r} \neq 0$ ⁠) 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 $- Ω$ and a negative growth rate is also required by property 3 and is the complex conjugate of the mode with frequency $Ω^{*}$ ⁠, in order to produce a real solution for z. (If $Ω_{r} = 0$ ⁠, only two of these four modes are required, as the others are redundant. In what follows, we assume that $Ω_{r} > 0$ ⁠. The $Ω_{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

z (t) = \sum_{ω > 0} a_{ω} (t) u_{ω} + a_{Ω} (t) u_{Ω} + a_{Ω^{*}} (t) u_{Ω^{*}} + c . c .,

(101)

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_{Ω}, u_{Ω}) = 0$ ⁠. It is therefore not possible to determine $a_{Ω}$ in terms of $z$ using the standard projection formula, Eq. (10). Fortunately, however, the complex mode with frequency $Ω^{*}$ can be used to determine $a_{Ω}$ 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

a_{Ω} (t) = \frac{(u_{Ω^{*}}, z (t))}{(u_{Ω^{*}}, u_{Ω})},

(102)

with an analogous expression for $a_{Ω^{*}}$ ⁠. 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:

{\dot{a}}_{Ω} (t) = - i Ω a_{Ω} (t) + f_{Ω} (t),

(103)

where $f_{Ω} (t) = (u_{Ω^{*}}, J \cdot f) / (u_{Ω^{*}}, u_{Ω})$ ⁠. This forcing coefficient can also be written as

f_{Ω} (t) = - i Ω \frac{u_{Ω^{*}}^{*} \cdot f}{(u_{Ω^{*}}, u_{Ω})},

(104)

using the same algebraic steps which led to Eq. (14). When $f_{Ω} = 0$ ⁠, the differential equation (103) is unforced and the solution grows exponentially with time at the growth rate γ, $a_{Ω} = A_{Ω} exp (- i Ω t) \propto exp (γ t)$ ⁠, where $A_{Ω}$ is an integration constant determined by initial conditions. The analogous equation for the mode with frequency $Ω^{*}$ implies a decaying mode amplitude $a_{Ω^{*}} = A_{Ω^{*}} exp (- i Ω^{*} t) \propto exp (- γ 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:

H = \sum_{ω > 0} H_{ω} + H_{u},

(105)

where H_u, the unstable mode contribution to the energy, is

H_{u} = 2 R e {a_{Ω} a_{Ω^{*}}^{*} (u_{Ω^{*}}, u_{Ω})} + 2 R e {a_{Ω} f \cdot u_{Ω} + a_{Ω^{*}} f \cdot u_{Ω^{*}}},

(106)

and where the stable mode contribution $H_{ω}$ 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_{Ω} (t)$ and $a_{Ω^{*}} (t)$ have differing time dependences $exp (- i Ω t)$ and $exp (- i Ω^{*} t)$ ⁠, respectively (the former growing and the latter decaying), the combination $a_{Ω} a_{Ω^{*}}^{*}$ 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_{Ω}, a_{ω}^{*}]$ for $R e ω > 0$ ⁠. Using Eq. (102) and the same series of steps as led to Eq. (18), we obtain

[a_{Ω}, a_{ω}^{*}] = {\begin{array}{l} \frac{- i Ω}{(u_{Ω^{*}}, u_{Ω})}, & ω = Ω^{*}, \\ 0 & otherwise, \end{array}

and similarly, $[a_{Ω}, a_{ω}] = 0$ for all eigenfrequencies ω with $R e ω > 0$ ⁠. For canonical coordinates, we therefore choose normalization $u_{Ω^{*}}$ and $u_{Ω}$ such that

(u_{Ω^{*}}, u_{Ω}) = Ω .

(107)

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_{Ω}, P_{Ω})$ and $(Q_{Ω^{*}}, P_{Ω^{*}})$ via the linear transformation

a_{Ω} = \frac{Q_{Ω} + i P_{Ω^{*}}}{\sqrt{2}},

(108)

a_{Ω^{*}} = \frac{Q_{Ω^{*}} + i P_{Ω}}{\sqrt{2}} .

(109)

With these choices, one can easily show that $[Q_{Ω}, P_{Ω}] = [Q_{Ω^{*}}, P_{Ω^{*}}] = 1$ and $[Q_{Ω}, Q_{ω}] = [P_{Ω}, P_{ω}] = 0$ ⁠. When written in terms of these variables, the unstable mode Hamiltonian is

\begin{matrix} H_{u} = Ω_{r} (Q_{Ω} Q_{Ω^{*}} + P_{Ω} P_{Ω^{*}}) + γ (Q_{Ω} P_{Ω} - Q_{Ω^{*}} P_{Ω^{*}}) \\ + \sqrt{2} f \cdot (Q_{Ω} R e u_{Ω} - P_{Ω^{*}} I m u_{Ω} + Q_{Ω^{*}} R e u_{Ω^{*}} - P_{Ω} I m u_{Ω^{*}}), \end{matrix}

(110)

where we have also employed Eq. (107) and have taken $Ω = Ω_{r} + i γ$ ⁠. Hamilton's equations for the complex-frequency modes then yield

\begin{matrix} {\dot{Q}}_{Ω} = \frac{\partial H_{u}}{\partial P_{Ω}} = Ω_{r} P_{Ω^{*}} + γ Q_{Ω} - \sqrt{2} f \cdot I m u_{Ω^{*}}, \\ {\dot{P}}_{Ω^{*}} = - \frac{\partial H_{u}}{\partial Q_{Ω^{*}}} = - Ω_{r} Q_{Ω} + γ P_{Ω^{*}} - \sqrt{2} f \cdot R e u_{Ω^{*}}, \\ {\dot{Q}}_{Ω^{*}} = \frac{\partial H_{u}}{\partial P_{Ω^{*}}} = Ω_{r} P_{Ω} - γ Q_{Ω^{*}} - \sqrt{2} f \cdot I m u_{Ω}, \\ {\dot{P}}_{Ω} = - \frac{\partial H_{u}}{\partial Q_{Ω}} = - Ω_{r} Q_{Ω^{*}} - γ P_{Ω} - \sqrt{2} f \cdot R e u_{Ω} . \end{matrix}

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 $Ω \to Ω^{*}$ in this equation.

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

z (t) = \sum_{ω > 0} a_{ω} (t) u_{ω} + c c . + a_{i γ} (t) u_{i γ} + a_{- i γ} (t) u_{- i γ},

(111)

with $(u_{i γ}, u_{i γ}) = (u_{- i γ}, u_{- i γ}) = 0$ according to property 2, but $(u_{- i γ}, u_{i γ}) \neq 0$ according to property 1. Following through with the rest of the algebra, we find that $a_{i γ (t)}$ and $a_{- i γ} (t)$ are real; that the normalization condition for canonical coordinates is $(u_{- i γ}, u_{i γ}) = γ$ ⁠; that the canonical coordinates can be chosen as $Q_{γ} = a_{i γ}, P_{γ} = a_{- i γ}$ ⁠; and that for this choice, the unstable mode Hamiltonian is

H_{u} = γ Q_{γ} P_{γ} + f \cdot (Q_{γ} u_{i γ} + P_{γ} u_{- i γ}) .

(112)

This Hamiltonian leads to the equations of motion

\begin{matrix} {\dot{Q}}_{γ} = \frac{\partial H_{u}}{\partial P_{γ}} = γ Q_{γ} + f \cdot u_{- i γ}, \\ {\dot{P}}_{γ} = - \frac{\partial H_{u}}{\partial Q_{γ}} = - γ P_{γ} - f \cdot u_{i γ}, \end{matrix}

which agree with Eqs. (103) and (104) when $Ω = \pm i γ$ ⁠. 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 = - B \hat{z}$ ⁠, with B > 0, and an electrostatic trap potential $ϕ_{0} (r, z)$ that is confining in the z direction for positive charges. In some experiments, this potential is nearly a pure quadrupole, $ϕ_{0} (r, z) = (1 / 2) Q (z^{2} - r^{2} / 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 $q_{i} > 0$ ⁠, and position r_i, $i = 1, \dots, N$ ⁠. (For negative charges in the trap, remove the −sign from B and add a −sign to $ϕ_{0}$ so that B and Q remain positive, and treat $q_{i} > 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 $ω_{r} > 0$ (i.e., the rotation is in the positive $ϕ$ direction). In a frame rotating with the charges, the system Hamiltonian is

H = \sum_{i} \frac{{(p_{i} - m_{i} Ω_{i} A (r_{i}))}^{2}}{2 m_{i}} + Φ (r_{1}, \dots, r_{N}) .

(113)

Here, the canonical momentum for particle i is $p_{i} = m_{i} {\dot{r}}_{i} + m_{i} Ω_{i} A (r_{i})$ where $Ω_{i} = q_{i} B / (m_{i} c) - 2 ω_{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 $\nabla \times A = - \hat{z}$ ⁠, where $\hat{z}$ is the unit vector in the z direction. A useful gauge choice for $A$ is the cylindrically symmetric gauge $A = - (1 / 2) \hat{z} \times r$ ⁠. This choice of gauge makes

p_{i} = m_{i} {\dot{r}}_{i} - \frac{1}{2} m_{i} Ω_{i} \hat{z} \times r_{i} .

(114)

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

Φ = \sum_{i > j} ϕ_{i j} + \sum_{i} q_{i} ϕ_{i} (r_{i}, z_{i}),

(115)

where $ϕ_{i j} = q_{i} q_{j} / | r_{i} - r_{j} |$ is the electrostatic Coulomb potential between particles i and j (neglecting, for simplicity, image charge effects in the surrounding electrodes), and

q_{i} ϕ_{i} (r_{,} z) = q_{i} ϕ_{0} (r, z) + \frac{ω_{r}}{2} (\frac{q_{i} B}{c} - m_{i} ω_{r}) r^{2}

(116)

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 $q_{i} ϕ_{i} (r_{,} z) = (1 / 2) q_{i} Q (z^{2} + β_{i} r^{2})$ ⁠, where the trap parameter $β_{i} \equiv ω_{r} (B / c - m_{i} ω_{r} / q_{i}) / Q - 1 / 2$ ⁠.

For positive values of the trap parameters β_i and Q [or, more generally, for a potential $q_{i} ϕ_{i} (r, z)$ that increases from the trap center with both increasing r and z], a (neutrally) stable ion crystal equilibrium exists with ${\dot{r}}_{i} = 0$ and $r_{i} = R_{i}$ for equilibrium positions R_i satisfying $\partial Φ / \partial R_{i} = 0, i = 1, \dots, N$ ⁠. Equilibrium configurations can be obtained numerically by minimization of the system potential energy $Φ$ ⁠, although for small numbers of charges one can find equilibria analytically via symmetry arguments.³⁶ For $N ≫ 1$ ⁠, 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 $\pm z$ axis, each at a distance $d = {(q / 4 Q)}^{1 / 3}$ from the origin. When $β < 1$ ⁠, the ions are on opposite sides of the origin in the x–y plane, each a distance $d / β^{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₂ or N₂, 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.

FIG. 1.

View large Download slide

Equilibrium of two identical charges in a quadrupole trap.

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 $β = 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 $β < 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 $β > 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

FIG. 2.

View large Download slide

r–z positions in a spheroidal Coulomb crystal of N = 236 identical charges in a quadrupolar Penning trap with trap parameter $β = 3 / 4$ ⁠. Distances are in terms of the distance ${(q / Q)}^{1 / 3}$ ⁠.

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 $r_{i} = R_{i} + δ r_{i}$ and momenta $p_{i} = m_{i} Ω_{i} A (R_{i}) + δ p_{i}$ ⁠. The Taylor expansion to second order in the small displacements from equilibrium then results in a linearized Hamiltonian system, with Hamiltonian

H = \frac{1}{2} (\sum_{i} \frac{{(δ p_{i} + m_{i} Ω_{i} \hat{z} \times δ r_{i} / 2)}^{2}}{m_{i}} + \sum_{i j} δ r_{i} \cdot V_{i j} \cdot δ r_{j}),

(117)

where $V_{i j} = \partial^{2} Φ / \partial R_{i} \partial R_{j}$ ⁠. The form of $V_{i j}$ 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 = (δ r_{1}, \dots, δ r_{N})$ and $p = (δ p_{1}, \dots, δ p_{N})$ ⁠, and with the symmetric Hamiltonian matrix given in block form by

H = (\begin{matrix} V + C & Ω \\ Ω^{t r} & M^{- 1} \end{matrix}) .

(118)

Here, $V, C, Ω$ ⁠, and $M^{- 1}$ are $3 N \times 3 N$ matrices. The matrix $M^{- 1}$ is the inverse of the diagonal mass matrix M for the system, with the diagonal elements $M_{i i}^{- 1}$ given by the vector $(m_{1}^{- 1}, m_{1}^{- 1}, m_{1}^{- 1}, \dots, m_{N}^{- 1}, m_{N}^{- 1}, m_{N}^{- 1})$ ⁠. The matrix $Ω$ 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:

Ω = (\begin{matrix} Ω_{1} & 0 \\ ⋱ \\ 0 & Ω_{N} \end{matrix}) .

(119)

Each diagonal block is, in dyadic notation, given by

Ω_{i} = \frac{1}{2} Ω_{i} (\hat{x} \hat{y} - \hat{y} \hat{x}), i = 1, \dots, N .

(120)

The symmetric matrix $V = \partial^{2} Φ / \partial R \partial R$ is the potential energy matrix given in block form by

V = (\begin{matrix} V_{11} & V_{12} & \dots & V_{1 N} \\ V_{21} & V_{22} & \dots & V_{2 N} \\ \dots \\ V_{N 1} & V_{N 2} & \dots & V_{N N} \end{matrix}) .

(121)

The matrix C is a magnetic potential contribution whose elements are zero everywhere except along the diagonal. The vector of diagonal elements $C_{i i}$ is given by

\frac{1}{4} (m_{1} Ω_{1}^{2}, m_{1} Ω_{1}^{2}, 0, \dots, m_{N} Ω_{N}^{2}, m_{N} Ω_{N}^{2}, 0) .

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

D = (\begin{matrix} - Ω & M^{- 1} \\ - V - C & - Ω \end{matrix}),

(122)

where we used the antisymmetry of the Lorentz matrix to write $Ω^{t r} = - Ω$ ⁠. The eigenvalues and eigenvectors of D provide us with the normal modes of the system, as per Eq. (5). As discussed in Sec. II these modes diagonalize the system energy.

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 $ω_{z} \equiv \sqrt{q Q / 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 $ω_{+}$ and $ω_{-}$ ⁠, respectively, where

ω_{\pm} = \sqrt{Ω^{2} + 2 ω_{⊥}^{2} \pm Ω \sqrt{Ω^{2} + 4 ω_{⊥}^{2}}} / \sqrt{2},

(123)

and where $ω_{⊥} \equiv \sqrt{β} ω_{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 $u_{0 z} = (r_{0 z}, p_{0 z})$ where

\begin{matrix} r_{0 z} = (\hat{z} \times R_{1}, \dots, \hat{z} \times R_{N}), \\ p_{0 z} = - \frac{1}{2} (m_{1} Ω_{1} \hat{z} \times \hat{z} \times R_{1}, \dots, m_{N} Ω_{N} \hat{z} \times \hat{z} \times R_{N}) \\ = \frac{1}{2} (m_{1} Ω_{1} R_{1} \hat{r}, \dots, m_{N} Ω_{N} R_{N} \hat{r}), \end{matrix}

(124)

and where R_j is the cylindrical radius of equilibrium position R_j for the jth ion. As was discussed previously in more general terms, this eigenmode corresponds to a constant of the motion, the momentum $P_{0 z}$ given by Eq. (37):

\begin{matrix} P_{0 z} = u_{0 z} \cdot J \cdot z = r_{0 z} \cdot p - p_{0 z} \cdot r \\ = \sum_{i} (R_{i} \hat{ϕ} \cdot δ p_{i} - \frac{1}{2} m_{i} Ω_{i} R_{i} \hat{r} \cdot δ r_{i}) \\ = \sum_{i} (m_{i} R_{i}^{2} δ {\dot{ϕ}}_{i} - m_{i} R_{i} Ω_{i} δ r_{i}), \end{matrix}

(125)

where in the last step we substituted for $δ p_{i}$ using Eq. (114). The constant $P_{0 z}$ is, of course, the perturbed total canonical angular momentum associated with rotations about the z axis.

The corresponding vector ${\bar{u}}_{0 z}$ ⁠, required for the rotational inertia $({\bar{u}}_{0 z}, {\bar{u}}_{0 z})$ [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 ${\bar{u}}_{0 z}$ 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),

{\bar{u}}_{0 z} = (- \frac{2 p_{0 z}}{3 m ω_{⊥}^{2}}, (1 + \frac{Ω^{2}}{6 ω_{⊥}^{2}}) r_{0 z})

(126)

and the rotational inertia $({\bar{u}}_{0 z}, {\bar{u}}_{0 z})$ is then given by the expression

({\bar{u}}_{0 z}, {\bar{u}}_{0 z}) = m \sum_{i} R_{i}^{2} (1 + \frac{Ω^{2}}{3 ω_{⊥}^{2}}) .

(127)

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

\frac{V}{m ω_{z}^{2}} = (\begin{matrix} β + \frac{1}{4} (1 - 3 cos 2 θ) & 0 & \frac{3}{4} sin 2 θ & - \frac{1}{4} (1 - 3 cos 2 θ) & 0 & - \frac{3}{4} sin 2 θ \\ β - 1 / 2 & 0 & 0 & β - 1 / 2 & 0 \\ \frac{1}{4} (5 + 3 cos 2 θ) & - \frac{3}{4} sin 2 θ & 0 & - \frac{1}{4} (1 + 3 cos 2 θ) \\ β + \frac{1}{4} (1 - 3 cos 2 θ) & 0 & \frac{3}{4} sin 2 θ \\ β - 1 / 2 & 0 \\ \frac{1}{4} (5 + 3 cos 2 θ) \end{matrix}) .

(128)

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

TABLE I.

Eigenmodes for two identical charges in a quadrupole trap.

ω	$u_{ω}$	$(u_{ω}, u_{ω})$
	COM modes:
ω_z	$(0, 0, 1, 0, 0, 1, 0, 0, - i m ω, 0, 0, - i m ω)$	$4 m ω^{2}$
$ω_{\pm}$	$(1, \pm i, 0, 1, \pm i, 0, - i m ω \pm i \frac{m Ω}{2}, \pm m ω - \frac{m Ω}{2}, 0, - i m ω \pm i \frac{m Ω}{2}, \pm m ω - \frac{m Ω}{2}, 0)$	$4 m (ω^{2} + ω_{⊥}^{2})$
	Other modes for $β > 1$ ⁠:
$\sqrt{3} ω_{z}$	$(0, 0, 1, 0, 0, - 1, 0, 0, - i m ω, 0, 0, i m ω)$	$4 m ω^{2}$
$ω_{r \pm}$	$(1, \pm i, 0, - 1, \mp i, 0, - i m ω \pm i \frac{m Ω}{2}, \pm m ω - \frac{m Ω}{2}, 0, i m ω \mp i \frac{m Ω}{2}, \mp m ω + \frac{m Ω}{2}, 0)$	$4 m (ω^{2} + (β - 1) ω_{z}^{2})$
	Other modes for $β < 1$ ⁠:
$ω_{z} \sqrt{1 - β}$	$(0, 0, 1, 0, 0, - 1, 0, 0, - i m ω, 0, 0, i m ω)$	$4 m ω^{2}$
$\sqrt{Ω^{2} + 3 ω_{⊥}^{2}}$	$(1, \frac{i Ω}{ω}, 0, - 1, - i \frac{Ω}{ω}, 0, - i m ω + i \frac{m Ω^{2}}{2 ω}, \frac{m Ω}{2}, 0, i m ω - i \frac{m Ω^{2}}{2 ω}, - \frac{m Ω}{2}, 0)$	$4 m ω^{2}$
0	$(0, 1, 0, 0, - 1, 0, \frac{m Ω}{2}, 0, 0, - \frac{m Ω}{2}, 0, 0)$	0
	Other modes for β = 1:
$ω_{θ \pm}$	$(r_{θ}, - r_{θ}, p_{θ}, - p_{θ})$	$\frac{4 m}{3} ω^{2} \frac{Ω^{4} - ω^{2} (Ω^{2} - 3 ω_{z}^{2}) - 3 Ω^{2} ω_{z}^{2} cos 2 θ}{ω_{z}^{2} {sin}^{2} θ}$
	$r_{θ} = (- i ω, Ω, - 2 i ω \frac{ω^{2} - Ω^{2} - 3 ω_{z}^{2} {sin}^{2} θ}{3 ω_{z}^{2} sin 2 θ})$
	$p_{θ} = m (\frac{Ω^{2}}{2} - ω^{2}, - i \frac{ω Ω}{2}, - 2 ω^{2} \frac{ω^{2} - Ω^{2} - 3 ω_{z}^{2} {sin}^{2} θ}{3 ω_{z}^{2} sin 2 θ})$
0	$u_{0 z} = (0, 1, 0, 0, - 1, 0, \frac{m Ω}{2}, 0, 0, - \frac{m Ω}{2}, 0, 0)$	0
0	$u_{0 y} = (cos θ, 0, - sin θ, - cos θ, 0, sin θ, 0, - \frac{m Ω}{2} cos θ, 0, 0, \frac{m Ω}{2} cos θ, 0)$	0

We first consider $β > 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, $δ r_{1} = - δ r_{2}$ ⁠. One of these is an axial stretch mode only along the z axis, with frequency $\sqrt{3} ω_{z}$ ⁠. The other two modes consist of circular motion in x and y at the frequencies $ω_{r +}$ and $ω_{r -}$ where

ω_{r \pm} = \sqrt{Ω^{2} + 2 (β - 1) ω_{z}^{2} \pm Ω \sqrt{Ω^{2} + 4 (β - 1) ω_{z}^{2}}} / \sqrt{2} .

(129)

As $β \to 1, ω_{r -}$ 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)]:

H = \sum_{ω > 0} | a_{ω} |^{2} (u_{ω}, u_{ω}),

(130)

with the coefficients $(u_{ω}, u_{ω})$ given in Table I.

On the other hand, for $β < 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 / β^{1 / 3}$ of the ions, defining $P'_{0 z} = P_{0 z} / R = δ p_{1 y} - δ p_{2 y} - m Ω (δ x_{1} - δ x_{2})$ ⁠. In addition, we introduce the normalized vector $\bar{u}'_{0} = {\bar{u}}_{0} / R$ with ${\bar{u}}_{0}$ given by Eq. (126). The two other modes are a tilt mode consisting only of axial motion with $δ z_{1} = - δ z_{2}$ ⁠, at frequency $ω_{z} \sqrt{1 - β}$ ⁠, and an “upper-hybrid” mode consisting of elliptical motion of the charges in the x–y plane at frequency $\sqrt{Ω^{2} + 3 ω_{⊥}^{2}}$ ⁠. Due to the neutral mode, the energy takes the form

H = \sum_{ω > 0} | a_{ω} |^{2} (u_{ω}, u_{ω}) + \frac{1}{2} \frac{P_{0 z}^{' 2}}{(\bar{u}'_{0}, \bar{u}'_{0})},

(131)

where the scaled moment of inertia $(\bar{u}'_{0}, \bar{u}'_{0}) = 2 m (1 + Ω^{2} / 3 ω_{⊥}^{2})$ ⁠, 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 $ω_{θ +}$ and $ω_{θ -}$ depend on θ:

ω_{θ \pm} = \sqrt{\frac{Ω^{2} + 3 ω_{z}^{2} \pm \sqrt{{(Ω^{2} - 3 ω_{z}^{2})}^{2} + 12 Ω^{2} ω_{z}^{2} {sin}^{2} θ}}{2}} .

(132)

As $θ \to 0$ ⁠, these frequencies approach $ω_{r +} = Ω$ and the stretch mode frequency $\sqrt{3} ω_{z}$ ⁠. As $θ \to π / 2, ω_{θ -} \to 0$ ⁠, and $ω_{θ +} \to \sqrt{Ω^{2} + 3 ω_{z}^{2}}$ ⁠.

1. Nonexistence of the rotational inertia for spherical confinement

The two zero frequency modes produce two constants of the motion, the angular momentum $P'_{0 z}$ due to rotation of the cluster about the z axis (again dividing the radii $d sin θ$ ⁠) and the angular momentum $P'_{0 y}$ due to rotations about the y axis (also dividing out distance d from the axis):

P'_{0 y} = cos θ (δ p_{x 1} - δ p_{x 2} + \frac{1}{2} m Ω (δ y_{1} - δ y_{2})) - sin θ (δ p_{z 1} - δ p_{z 2}) .

(133)

These two constants are not in involution, $J'_{x y} = [P'_{0 y}, P'_{0 z}] = 2 m Ω cos θ$ ⁠. 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 $β < 1$ system when the perturbed axial angular momentum $P'_{0 z}$ 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 $β < 1$ ⁠, Hamiltonian (131) implies that the angle variable $δ ϕ = a_{0}$ changes linearly with time according to

δ \dot{ϕ} = \frac{P_{0 z}^{'}}{(\bar{u}'_{0}, \bar{u}'_{0})},

(134)

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 $(\bar{u}'_{0}, \bar{u}'_{0})$ which relates the scaled angular momentum change $P_{0 z}^{'}$ to the rotation frequency change $δ \dot{ϕ}$ ⁠.

However, for β = 1, neither $P'_{0 z}$ nor $P'_{0 y}$ appear in the Hamiltonian. The phase space configuration $z$ evolves according to

z = d δ θ u_{0 y} + R δ ϕ u_{0 z} + \sum_{ω \neq 0} a_{ω} u_{ω},

(135)

where $u_{0 z}$ and $u_{0 y}$ are given in Table I, and where

δ θ = - \frac{P'_{0 z}}{J'_{x y} d},

(136)

and

δ ϕ = \frac{P'_{0 y}}{J'_{x y} R},

(137)

[see Eqs. (94), (96), and (97)]. Equation (137) shows that the angle $δ ϕ$ does not evolve in time as it did for $β < 1$ [see Eq. (134)]. This is because the axial angular momentum perturbation $P'_{0 z}$ does not change the rotation frequency of the cluster. Indeed, the moment of inertia $({\bar{u}}_{0 z}, {\bar{u}}_{0 z})$ is undefined since the vector ${\bar{u}}_{0 z}$ does not exist, as discussed in Sec. II B 2.

The axial angular momentum perturbation $P'_{0 z}$ is instead accomplished by a rotation $δ θ$ of the cluster about the y axis, as exhibited in Eq. (136). This can occur because for $Ω \neq 0$ 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 $Ω = 20 ω_{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 $ω_{⊥} = \sqrt{β} ω_{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.

FIG. 3.

View large Download slide

Eigenfrequencies (in units of ω_z) for the spheroidal Coulomb crystal shown in Fig. 2, for vortex frequency Ω = 0, with frequencies ordered from the highest frequency to the lowest non-negative frequency. The mode number n (⁠ $n = 1, \dots, 3 \times 236 = 708$ ⁠) labels the position of a given mode in this ordering.

FIG. 4.

View large Download slide

Eigenfrequencies (in units of ω_z) for the spheroidal Coulomb crystal shown in Fig. 2, for vortex frequency $Ω = 20 ω_{z}$ ⁠, ordered from the highest frequency to the lowest non-negative frequency. Note the breaks in the frequency axis, separating modes into cyclotron, axial, and E × B branches. As shown in the previous figure, the mode number n (⁠ $n = 1, \dots, 3 \times 236 = 708$ ⁠) labels the position of a given mode in the ordering.

For $Ω = 20 ω_{z}$ ⁠, the mode frequencies condense into three groups: a group of N cyclotron modes with $ω > Ω$ ⁠; a group of N axial modes with $ω \approx 1$ ⁠, and a group of N E × B modes with low frequencies. Variation of Ω indicates that these latter mode frequencies scale with Ω as $1 / Ω$ ⁠. 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.⁴⁰ 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.¹⁰ 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).¹⁰

A single example is displayed in Fig. 5. Here, we consider a fluid displacement of the form $δ r_{fluid} = {(δ x, δ y, δ z)}_{fluid} \propto (x, - y, 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 $ϕ$ direction. For the case of a spheroidal plasma with $β = 3 / 4$ ⁠, the frequencies of these two modes is predicted to be $ω_{fluid} = \sqrt{Ω^{2} / 4 + 0.925253 ω_{z}^{2}} \pm Ω / 2$ ⁠. In Fig. 5, we determine the mode amplitude $a_{ω}$ for each mode by projecting the phase-space displacement $z$ onto each mode using Eq. (10), with the eigenmodes for $Ω = 20 ω_{z}$ ⁠, and with $z$ given by $δ r_{fluid}$ for each particle, along with the associated canonical momentum $- m Ω \hat{z} \times δ r_{fluid} / 2$ ⁠: $z = (x_{1}, - y_{1}, 0, \dots, x_{N}, - y_{N}, 0, - m Ω y_{1} / 2, - m Ω x_{1} / 2, 0, \dots,, - m Ω y_{N} / 2, - m Ω x_{N} / 2, 0)$ ⁠. This phase space configuration corresponds to an initially-stationary elliptical distortion of the plasma crystal.

FIG. 5.

View large Download slide

The magnitude of eigenmode amplitudes $a_{ω}$ plotted on a logarithmic scale vs mode number n, for an initial condition corresponding to an ellipsoidal distortion of the crystal, as described in the text.

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, $ω / ω_{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 $ω_{fluid} / ω_{z} = 20.0462, 0.0462$ for the cyclotron and E × B mode, respectively. The weaker of the three peaks at $ω = 1.3063 ω_{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 $ω_{fluid} = 1.3037 ω_{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 $⟨ δ z_{j} δ z_{k} ⟩$ for two particles j and k. The jth and kth axial displacements appearing here correspond to the α and β components of the phase space vector $z$ where $α = 3 j$ and $β = 3 k$ ⁠. 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),

δ z_{j} = \sum_{ω > 0} (Q_{ω} + i P_{ω}) u_{ω α} / \sqrt{2} + c . c . + a_{0} u_{0 α} + P_{0} \frac{{\bar{u}}_{0 z α}}{({\bar{u}}_{0 z}, {\bar{u}}_{0 z})},

(138)

and similarly for $δ z_{k}$ ⁠. Using Eqs. (31) and (32), the thermal averages then yield

⟨ δ z_{i} δ z_{j} ⟩ = T \sum_{ω > 0} \frac{u_{ω α} u_{ω β}^{*}}{(u_{ω}, u_{ω})} + c . c + \frac{T}{({\bar{u}}_{0 z}, {\bar{u}}_{0 z})} {\bar{u}}_{0 z α} {\bar{u}}_{0 z β},

(139)

where we used Eq. (21) to replace ω by $(u_{ω}, u_{ω})$ ⁠, so that the expression is valid for unnormalized eigenvectors. The last term arises from the average over the zero frequency mode momentum $⟨ P_{0}^{2} ⟩$ ⁠, which we evaluated using Hamiltonian (58). There is no contribution from a₀ because the neutrally stable eigenmode is a pure rotation about the z axis by angle $δ ϕ = a_{0}$ ⁠, with eigenvector given by Eq. (124). This eigenvector has no axial component (i.e., $u_{0 α} = 0$ ⁠).

1. Thermal averages for the two particle Coulomb cluster

We can evaluate $⟨ δ z_{j} δ z_{k} ⟩$ for the two particle system using the information in Table I. For $β > 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 $β > 1$ ⁠, since charges in equilibrium are aligned along the z axis. We thus obtain, for $β > 1$ ⁠,

⟨ δ z_{1}^{2} ⟩ = ⟨ δ z_{2}^{2} ⟩ = \frac{2 T}{4 m ω_{z}^{2}} + \frac{2 T}{12 m ω_{z}^{2}} = \frac{2 T}{3 m ω_{z}^{2}},

(140)

⟨ δ z_{1} δ z_{2} ⟩ = \frac{2 T}{4 m ω_{z}^{2}} - \frac{2 T}{12 m ω_{z}^{2}} = \frac{T}{3 m ω_{z}^{2}},

(141)

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 $β < 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 ${\bar{u}}_{0}$ when the charges are trapped in the x–y plane, see Eq. (126). Then, Eq. (139) yields

⟨ δ z_{1}^{2} ⟩ = ⟨ δ z_{2}^{2} ⟩ = \frac{2 T}{4 m ω_{z}^{2}} + \frac{2 T}{4 m ω_{z}^{2} (1 - β)} = \frac{T}{2 m ω_{z}^{2}} \frac{2 - β}{1 - β},

(142)

⟨ δ z_{1} δ z_{2} ⟩ = \frac{2 T}{4 m ω_{z}^{2}} - \frac{2 T}{4 m ω_{z}^{2} (1 - β)} = - \frac{T}{2 m ω_{z}^{2}} \frac{β}{1 - β},

(143)

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 $β \to 1$ ⁠, as the tilt mode becomes zero frequency, allowing large fluctuations in the axial displacements of the charges. Note that the equations imply that $⟨ {(δ z_{1} + δ z_{2})}^{2} ⟩ / 4 = T / (2 m ω_{z}^{2})$ ⁠, independent of β. This is the mean square fluctuation in the axial center of mass position, with the expected form $T / (N m ω_{z}^{2})$ for a particle of mass $N m$ 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.²⁴ 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, $⟨ δ x_{j} δ x_{k} ⟩$ ⁠. This average is finite for an N = 2 cluster with equilibrium positions in the x–z plane because zero-frequency rotations through $ϕ$ 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, $α = 3 j - 2$ and $β = 3 k - 2$ ⁠. For the two particle Coulomb cluster, and taking $β > 1$ ⁠, we obtain

⟨ δ x_{1}^{2} ⟩ = ⟨ δ x_{2}^{2} ⟩ = \frac{T}{2 m} (\frac{1}{ω_{-}^{2} + ω_{⊥}^{2}} + \frac{1}{ω_{+}^{2} + ω_{⊥}^{2}} + \frac{1}{ω_{r -}^{2} + (β - 1) ω_{z}^{2}} + \frac{1}{ω_{r +}^{2} + (β - 1) ω_{z}^{2}}),

(144)

⟨ δ x_{1} δ x_{2} ⟩ = \frac{T}{2 m} (\frac{1}{ω_{-}^{2} + ω_{⊥}^{2}} + \frac{1}{ω_{+}^{2} + ω_{⊥}^{2}} - \frac{1}{ω_{r -}^{2} + (β - 1) ω_{z}^{2}} - \frac{1}{ω_{r +}^{2} + (β - 1) ω_{z}^{2}}) .

(145)

Each term on the right hand side depends explicitly on Ω, but their sum does not. In fact, $1 / (ω_{-}^{2} + ω_{⊥}^{2}) + 1 / (ω_{+}^{2} + ω_{⊥}^{2}) = 1 / ω_{⊥}^{2}$ and $1 / (ω_{r -}^{2} + (β - 1) ω_{z}^{2}) + 1 / (ω_{r +}^{2} + (β - 1) ω_{z}^{2}) = 1 / ((β - 1) ω_{z}^{2})$ ⁠. When these formulas are applied to Eqs. (144) and (145), we obtain the Ω-independent result

⟨ δ x_{1}^{2} ⟩ = ⟨ δ x_{2}^{2} ⟩ = \frac{T}{2 m ω_{z}^{2}} (\frac{1}{β} + \frac{1}{β - 1}),

(146)

⟨ δ x_{1} δ x_{2} ⟩ = \frac{T}{2 m ω_{z}^{2}} (\frac{1}{β} - \frac{1}{β - 1}),

(147)

where we used the relationship $ω_{⊥}^{2} = β ω_{z}^{2}$ ⁠.

Re-evaluating the averages for $β < 1$ ⁠, we obtain

⟨ δ x_{1}^{2} ⟩ = ⟨ δ x_{2}^{2} ⟩ = \frac{T}{2 m} (\frac{1}{ω_{-}^{2} + ω_{⊥}^{2}} + \frac{1}{ω_{+}^{2} + ω_{⊥}^{2}} + \frac{1}{Ω^{2} + 3 ω_{⊥}^{2}} + \frac{{(Ω / 3 ω_{⊥}^{2})}^{2}}{1 + Ω^{2} / 3 ω_{⊥}^{2}}),

(148)

⟨ δ x_{1} δ x_{2} ⟩ = \frac{T}{2 m} (\frac{1}{ω_{-}^{2} + ω_{⊥}^{2}} + \frac{1}{ω_{+}^{2} + ω_{⊥}^{2}} - \frac{1}{Ω^{2} + 3 ω_{⊥}^{2}} - \frac{{(Ω / 3 ω_{⊥}^{2})}^{2}}{1 + Ω^{2} / 3 ω_{⊥}^{2}}) .

(149)

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 ${\bar{u}}_{0}$ and Eq. (127) for $({\bar{u}}_{0}, {\bar{u}}_{0})$ ⁠. Again, each term depends on Ω, but when summed the results are Ω-independent:

⟨ δ x_{1}^{2} ⟩ = ⟨ δ x_{2}^{2} ⟩ = \frac{2 T}{3 m ω_{⊥}^{2}},

(150)

⟨ δ x_{1} δ x_{2} ⟩ = \frac{T}{3 m ω_{⊥}^{2}} .

(151)

As another test of the Bohr–van Leeuwen theorem, we consider the average $⟨ {(δ r_{1} \cdot {\hat{R}}_{1})}^{2} ⟩ = ⟨ {(δ x_{1} sin θ + δ z_{1} cos θ)}^{2} ⟩$ for the N = 2 cluster at β = 1, oriented at angle θ with respect to the z-axis. Here, ${\hat{R}}_{1} = (sin θ, 0, cos θ)$ is the unit vector in the direction of $R_{1}$ ⁠, 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 $ω_{θ \pm}$ enter the average:

\begin{matrix} ⟨ {(δ r_{1} \cdot {\hat{R}}_{1})}^{2} ⟩ = ⟨ δ x_{1}^{2} ⟩ {sin}^{2} θ + ⟨ δ z_{1}^{2} ⟩ {cos}^{2} θ + 2 ⟨ δ x_{1} δ z_{1} ⟩ sin θ cos θ \\ = {sin}^{2} θ \frac{T}{2 m} \sum_{ω = ω_{\pm}} \frac{1}{ω^{2} + ω_{⊥}^{2}} + {cos}^{2} θ \frac{T}{2 m ω_{z}^{2}} \\ + 2 T \sum_{ω = ω_{θ_{\pm}}} \frac{1}{(u_{ω}, u_{ω})} (ω^{2} {sin}^{2} θ + Δ^{2} {cos}^{2} θ + 2 ω Δ sin θ cos θ), \end{matrix}

(152)

where $Δ = 2 ω \frac{ω^{2} - Ω^{2} - 3 ω_{z}^{2} {sin}^{2} θ}{3 ω_{z}^{2} sin 2 θ}$ ⁠. 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 / (2 m ω_{z}^{2})$ ⁠, the expected radial thermal fluctuation for a single particle of mass 2m in a spherically symmetric harmonic well of frequency ω_z. Surprisingly, perhaps, considering its complexity, the last term sums to $T / (6 m ω_{z}^{2})$ ⁠. Thus we obtain

⟨ {(δ r_{1} \cdot {\hat{R}}_{1})}^{2} ⟩ = \frac{2 T}{3 m ω_{z}^{2}},

(153)

which agrees with Eq. (150) when $θ = π / 2$ and β = 1, and with Eq. (140) when θ = 0.

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: $\sum_{j = 1}^{N} {⟨ δ z_{j}^{2} ⟩}_{ω}$ ⁠, where ${⟨ δ z_{j}^{2} ⟩}_{ω} = 2 T | u_{ω 3 j} |^{2} / (u_{ω}, u_{ω})$ is the mean square fluctuation in axial position z_j caused by mode ω, as per Eq. (139).

FIG. 6.

View large Download slide

The thermal average $\sum_{j = 1}^{N} {⟨ δ z_{j}^{2} ⟩}_{ω}$ for each eigenmode, in units of $T / (m ω_{z}^{2})$ ⁠, for the spherioidal crystal shown in Fig. 2, and for two vortex frequencies, Ω = 0 and $Ω = 20 ω_{z}$ ⁠. For $Ω = 20 ω_{z}$ ⁠, the contributions of cyclotron modes (⁠ $1 \leq n \leq 236$ ⁠), axial modes (⁠ $236 < n \leq 472$ ⁠), and E × B modes (n > 472) are substantially different in magnitude.

We use this average to evaluate each mode's contribution to $⟨ δ z^{2} ⟩ = N^{- 1} \sum_{j = 1}^{N} ⟨ z_{j}^{2} ⟩$ via

N ⟨ δ z^{2} ⟩ = \sum_{ω > 0} \sum_{j = 1}^{N} {⟨ δ z_{j}^{2} ⟩}_{ω} + \frac{T}{({\bar{u}}_{0 z}, {\bar{u}}_{0 z})} \sum_{j = 1}^{N} {\bar{u}}_{0 z 3 j}^{2};

(154)

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 $Ω = 20 ω_{z}$ ⁠, cyclotron modes make a negligible contribution to the average, as one might also expect; we obtain $\sum_{ω_{cyclotron}} \sum_{j = 1}^{N} {⟨ δ z_{j}^{2} ⟩}_{ω} = 1.096 \times 10^{- 6} T / (m ω_{z}^{2})$ ⁠. Axial modes make a larger contribution, contributing to the average an amount $\sum_{ω_{axial}} \sum_{j = 1}^{N} {⟨ δ z_{j}^{2} ⟩}_{ω} = 385.0969 T / (m ω_{z}^{2})$ ⁠. Surprisingly, however, the low-frequency E × B drift modes dominate by a factor of 10, contributing $\sum_{ω_{E \times B}} \sum_{j = 1}^{N} {⟨ δ z_{j}^{2} ⟩}_{ω} = 3575.0989 T / (m ω_{z}^{2})$ ⁠. One normally thinks of $E \times B$ 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 $⟨ δ z^{2} ⟩ = 3690.4421 / N = 16.7815$ in units of $T / m ω_{z}^{2}$ ⁠. 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 ${\bar{u}}_{0 z}$ vanish, i.e., ${\bar{u}}_{0 z 3 j} = 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 $⟨ δ z^{2} ⟩$ are small compared to λ_eff. This could open the door for trapped-ion quantum simulation and sensing work^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 $Ω_{i} = 0$ ⁠. In this case, the dynamical matrix Eq. (122) reduces to

D = (\begin{matrix} 0 & M^{- 1} \\ - V & 0 \end{matrix}),

(155)

and the eigenvalue problem Eq. (5) corresponds to the coupled equations $- i ω r = M^{- 1} \cdot p, - i ω p = - V \cdot r$ which can be combined into a reduced eigenvalue problem for $ω^{2}$ ⁠, $ω^{2} r = M^{- 1} \cdot V \cdot r$ ⁠. 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 \cdot M \cdot b$ 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 $Ω_{i} ≫ ω_{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 $ω - Ω_{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

H^{(0)} = (\begin{matrix} C & Ω \\ - Ω & M^{- 1} \end{matrix}),

(156)

and

D^{(0)} = J \cdot H^{(0)} = (\begin{matrix} - Ω & M^{- 1} \\ - C & - Ω \end{matrix}),

(157)

while the first-order parts are

H^{(1)} = (\begin{matrix} V & 0 \\ 0 & 0 \end{matrix}),

(158)

and

D^{(1)} = J \cdot H^{(1)} = (\begin{matrix} 0 & 0 \\ - V & 0 \end{matrix}) .

(159)

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 6N dimensional eigenvector $u_{ω}$ 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_{ω} = u_{j, α}^{(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

\begin{matrix} u_{j α}^{(0)} = (0, \dots, & r_{j α}, 0, \dots & p_{j α}, 0, \dots & 0) . \\ i = 1, \dots, & j, \dots, & N + j, \dots & 2 N . \end{matrix}

(160)

The vector $(r_{j α}, p_{j α})$ solves the six-dimensional single-particle eigenvalue problem corresponding to the particle j elements in $D^{(0)}$ ⁠:

- i ω (r_{j α}, p_{j α}) = (\begin{matrix} 0 & - Ω_{j} / 2 & 0 & 1 / m_{j} & 0 & 0 \\ Ω_{j} / 2 & 0 & 0 & 0 & 1 / m_{j} & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 / m_{j} \\ - m Ω_{j}^{2} / 4 & 0 & 0 & 0 & - Ω_{j} / 2 & 0 \\ 0 & - m Ω_{j}^{2} / 4 & 0 & Ω_{j} / 2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}) \cdot (r_{j α}, p_{j α}) .

(161)

This eigenvalue problem has five independent eigenvectors, three of which correspond to zero-frequency modes, and two of which are cyclotron modes with frequencies $\pm Ω_{j}$ ⁠. The positive and negative frequency cyclotron eigenvectors for particle j are

(r_{j +}, p_{j +}) = \frac{1}{2 \sqrt{2} m_{j} Ω_{j}} (2, 2 i, 0, - i m_{j} Ω_{j}, m_{j} Ω_{j}, 0), ω = Ω_{j},

(162)

(r_{j -}, p_{j -}) = {(r_{j +}, p_{j +})}^{*}, ω = - Ω_{j} .

(163)

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 $\hat{ϕ}$ direction).

The three independent zero-frequency eigenvectors for particle j are

(r_{j X}, p_{j X}) = (1, 0, 0, 0, - m_{j} Ω_{j} / 2, 0),

(164)

(r_{j Y}, p_{j Y}) = (0, 1, 0, m_{j} Ω_{j} / 2, 0, 0),

(165)

(r_{j Z}, p_{j Z}) = (0, 0, 1, 0, 0, 0) .

(166)

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 $P_{j X}, P_{j Y}, P_{j Z}$ ⁠. For example, $P_{j X} = (r_{j X}, p_{j X}) \cdot J \cdot (δ x_{j}, δ y_{j}, δ z_{j}, δ p_{x j}, δ p_{y j}, δ p_{z j})$ ⁠, and so on. For $Ω_{j} \neq 0$ ⁠, one can check that $P_{j X}$ and $P_{j Y}$ are not in involution,

[P_{j X}, P_{j Y}] = m_{j} Ω_{j},

(167)

but $[P_{j Y}, P_{j Z}] = [P_{j X}, P_{j Z}] = 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:

({\bar{r}}_{j Z}, {\bar{p}}_{j Z}) = (0, 0, 0, 0, 0, 0, 1),

(168)

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

\begin{matrix} {\bar{u}}_{j Z} = (0, \dots & , 0, {\bar{p}}_{j Z}, 0, \dots 0) . \\ i = 1, \dots & N + j, \dots 2 N . \end{matrix}

(169)

These vectors are not eigenvectors, but instead satisfy

D^{(0)} \cdot {\bar{u}}_{j, Z} = u_{j, Z}^{(0)} / m_{j} .

(170)

The vectors are orthogonal to all of the zeroth-order eigenvectors $u_{k α}^{(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

{(a, b)}^{(0)} = a^{*} \cdot H^{(0)} \cdot b .

(171)

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_{ω}$ ⁠, writing it as a superposition of the degenerate zeroth order modes with frequency Ω_a:

u_{ω} = \sum_{j = 1}^{N_{a}} c_{j ω} u_{j, +}^{(0)} + u^{(1)},

(172)

where the sum over j sums only over particles of species a, the coefficients $c_{j ω}$ 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

- i ω \sum_{j = 1}^{N_{a}} c_{j ω} u_{j, +}^{(0)} - i ω u^{(1)} = \sum_{j = 1}^{N_{a}} c_{j ω} D^{(0)} \cdot u_{j, +}^{(0)} + \sum_{j = 1}^{N_{a}} c_{j ω} D^{(1)} \cdot u_{j, +}^{(0)} + D^{(0)} \cdot u^{(1)} + D^{(1)} \cdot u^{(1)} .

(173)

We drop the last term in this equation because it is second-order, use the fact that $D^{(0)} \cdot u_{j, +}^{(0)} = - i Ω_{a} u_{j, +}^{(0)}$ ⁠, multiply through by i, and take a zeroth-order inner product with respect to one of the eigenmodes $u_{k, +}^{(0)}$ ⁠. The orthogonality and the Hermitian nature of the matrix $i D^{(0)}$ then annihilate several terms in the equation, leaving us with

ω c_{k ω} {(u_{k, +}^{(0)}, u_{k, +}^{(0)})}^{(0)} = Ω_{a} c_{k ω} {(u_{k, +}^{(0)}, u_{k, +}^{(0)})}^{(0)} + \sum_{j = 1}^{N_{a}} c_{j ω} {(u_{k, +}^{(0)}, i D^{(1)} \cdot u_{j, +}^{(0)})}^{(0)} .

(174)

The inner product ${(u_{k, +}^{(0)}, u_{k, +}^{(0)})}^{(0)}$ involves only particle k and, using Eqs. (156), (160), and (162), evaluates to

{(u_{k, +}^{(0)}, u_{k, +}^{(0)})}^{(0)} = \frac{1}{m_{a}}

(175)

[recall that all particles in Eq. (174) are of species a]. The inner product ${(u_{k, +}^{(0)}, i D^{(1)} \cdot u_{j, +}^{(0)})}^{(0)}$ involves only particles j and k in species a and evaluates to

\begin{matrix} {(u_{k, +}^{(0)}, i D^{(1)} \cdot u_{j, +}^{(0)})}^{(0)} = u_{k, +}^{(0) *} \cdot H^{(0)} \cdot i D^{(1)} \cdot u_{j, +}^{(0)} \\ = u_{k, +}^{(0) *} \cdot H^{(0)} \cdot J \cdot i H^{(1)} \cdot u_{j, +}^{(0)} \\ = - i Ω_{a} u_{k, +}^{(0) *} \cdot i H^{(1)} \cdot u_{j, +}^{(0)} \end{matrix}

(176)

\begin{matrix} = Ω_{a} r_{k +}^{*} \cdot V_{k j} \cdot r_{j +} \\ = \frac{1}{2 m_{a}^{2} Ω_{a}} (1, - i, 0) \cdot V_{k j} \cdot (1, i, 0) \\ = \frac{1}{2 m_{a}^{2} Ω_{a}} (V_{kjxx} + V_{kjyy}), \end{matrix}

(177)

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 $V_{kjxx}$ is the $\hat{x} \hat{x}$ component of the potential matrix $V_{k j} = \partial^{2} Φ / \partial R_{j} \partial R_{k}$ and similarly for $V_{kjyy}$ ⁠. Equation (174) can then be written in vector form as

(ω - Ω_{a}) c_{ω} = F \cdot c_{ω},

(178)

where the matrix F has components

F_{j k} = \frac{V_{jkxx} + V_{jkyy}}{2 m_{a} Ω_{a}} .

(179)

Thus, the frequency shift $ω - Ω_{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 $N_{a} \times N_{a}$ ⁠, yielding N_a eigenfrequencies that differ from Ω_a by a small shift, proportional to $1 / Ω_{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 \cdot V \cdot r > 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 $v_{k} (1, i, 0)$ for some complex amplitude v_k. Then, consider the quantity $r^{*} \cdot V \cdot r = 2 m_{a} Ω_{a} v^{*} \cdot F \cdot v$ where $v$ is the vector of complex coefficients v_k. However, $r^{*} \cdot V \cdot r > 0$ because, for any complex vector $r = a + i b$ (for real $a$ and $b$ ⁠), $r^{*} \cdot V \cdot r = a \cdot V \cdot a + b \cdot V \cdot b$ (since V is symmetric) and both terms on the right hand side are positive. Therefore, the quantity $v^{*} \cdot F \cdot v$ 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^{*} \cdot F \cdot v$ ⁠.

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 $Ω = 20 ω_{z}$ ⁠, the fractional error between the exact mode frequencies and the approximate frequencies is quite small.

FIG. 7.

View large Download slide

Fractional difference between cyclotron frequencies evaluated using Eq. (178) and the exact frequencies evaluated using Eq. (5) and displayed in Fig. 4, for the same parameters as in that figure.

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 / \sqrt{2}) \sum_{ω > 0} a_{ω} \sum_{j = 1}^{N_{a}} c_{j ω} u_{j +}^{(0)} + c . c .$ ⁠, where $a_{ω}$ is the amplitude of each mode. This definition, with the extra factor of $1 / \sqrt{2}$ ⁠, allows us to identify $| c_{k ω} |$ as m_a multiplied by the speed of particle k in a given mode, for unit amplitude $a_{ω} = 1$ ⁠, because the unit amplitude results in a cyclotron radius for this particle of $| c_{k ω} | / m_{a} Ω_{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:

\begin{matrix} H = \frac{1}{2} \sum_{ω, \bar{ω} > 0} a_{ω} a_{\bar{ω}}^{*} \sum_{j, k = 1}^{N_{a}} c_{k ω} c_{j \bar{ω}} (u_{k +}^{(0) *} \cdot H^{(0)} \cdot u_{j +}^{(0)} + u_{k +}^{(0) *} \cdot H^{(1)} \cdot u_{j +}^{(0)}) \\ = \frac{1}{2} \sum_{ω, \bar{ω} > 0} a_{ω} a_{\bar{ω}}^{*} \sum_{j, k = 1}^{N_{a}} c_{k ω} c_{j \bar{ω}} (\frac{1}{m_{a}} δ_{j k} + \frac{1}{m_{a} Ω_{a}} F_{k j}) \\ = \sum_{ω, \bar{ω} > 0} a_{ω} a_{\bar{ω}}^{*} c_{ω} \cdot c_{\bar{ω}} (\frac{1}{2 m_{a}} + \frac{1}{2 m_{a} Ω_{a}} (\bar{ω} - Ω_{a})) \end{matrix}

(180)

= \frac{1}{2 m_{a} Ω_{a}} \sum_{ω > 0} ω | a_{ω} |^{2} c_{ω} \cdot c_{ω} .

(181)

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_{ω} |^{2} c_{ω} \cdot c_{ω}$ ⁠, also summed over the modes, and divided by two times the particle mass. The second term in the parenthesis, proportional to $\bar{ω} - Ω_{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_{ω}$ in terms of the zero-frequency eigenvectors along with the extra vectors ${\bar{u}}_{j Z}$ ⁠:

u_{ω} = \sum_{j = 1}^{N} (X_{j} u_{j X}^{(0)} + Y_{j} u_{j Y}^{(0)} + Z_{j} u_{j Z}^{(0)} + P_{j} {\bar{u}}_{j Z}^{(0)}) + u^{(1)},

(182)

where X_j, Y_j, and Z_j are displacement amplitudes for particle j in the x, y, and z directions, respectively, P_j is the axial momentum of particle 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

u_{j X}^{(0)} \cdot J \cdot u^{(1)} = u_{j Y}^{(0)} \cdot J \cdot u^{(1)} = u_{j Z}^{(0)} \cdot J \cdot u^{(1)} = {\bar{u}}_{j Z}^{(0)} \cdot J \cdot u^{(1)} = 0.

(183)

We substitute Eq. (182) into Eq. (5) and again break $D$ into zeroth-order and first-order parts:

\begin{matrix} - i ω \sum_{j = 1}^{N} ((X_{j} u_{j X}^{(0)} + Y_{j} u_{j Y}^{(0)} + Z_{j} u_{j Z}^{(0)} + P_{j} {\bar{u}}_{j Z}^{(0)}) - i ω u^{(1)} \\ = D^{(0)} \cdot \sum_{j = 1}^{N} (X_{j} u_{j X}^{(0)} + Y_{j} u_{j Y}^{(0)} + Z_{j} u_{j Z}^{(0)} + P_{j} {\bar{u}}_{j Z}^{(0)}) \\ + D^{(1)} \cdot \sum_{j = 1}^{N} (X_{j} u_{j X}^{(0)} + Y_{j} u_{j Y}^{(0)} + Z_{j} u_{j Z}^{(0)} + P_{j} {\bar{u}}_{j Z}^{(0)}) \\ + D^{(0)} \cdot u^{(1)} + D^{(1)} \cdot u^{(1)} . \end{matrix}

(184)

We drop the last term on the right hand side, since it is second order. Also, $D^{(0)} \cdot u_{j α}^{(0)} = 0$ for all the zero-frequency eigenmodes, so Eq. (184) simplifies to

\begin{matrix} - i ω \sum_{j = 1}^{N} ((X_{j} u_{j X}^{(0)} + Y_{j} u_{j Y}^{(0)} + Z_{j} u_{j Z}^{(0)} + P_{j} {\bar{u}}_{j Z}^{(0)}) - i ω u^{(1)} \\ = \sum_{j = 1}^{N} P_{j} u_{j Z}^{(0)} / m_{j} + D^{(1)} \cdot \sum_{j = 1}^{N} (X_{j} u_{j X}^{(0)} + Y_{j} u_{j Y}^{(0)} + Z_{j} u_{j Z}^{(0)} + P_{j} {\bar{u}}_{j Z}^{(0)}) \\ + D^{(0)} \cdot u^{(1)}, \end{matrix}

(185)

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 $- {\bar{u}}_{k Z}^{(0)} \cdot J = u_{k Z}^{(0)}$ and using the orthogonality conditions (183) yields

- i ω Z_{k} = P_{k} / m_{k} + {\bar{u}}_{k Z}^{(0)} \cdot H^{(1)} \cdot \sum_{j = 1}^{N} (X_{j} u_{j X}^{(0)} + Y_{j} u_{j Y}^{(0)} + Z_{j} u_{j Z}^{(0)} + P_{j} {\bar{u}}_{j Z}^{(0)}),

(186)

where we used the identity $J \cdot D = J \cdot J \cdot H = - H$ ⁠. However, one can use Eqs. (158) and (169) to check that ${\bar{u}}_{k Z}^{(0)} \cdot H^{(1)} = 0$ ⁠, which annihilates the sum on the right hand side so Eq. (186) simplifies to

- i ω Z_{k} = P_{k} / m_{k},

(187)

the standard relationship between the axial position and momentum. Now act on Eq. (185) with $u_{k Z}^{(0)} \cdot J = {\bar{u}}_{k Z}^{(0)}$ ⁠, which results in

\begin{matrix} - i ω P_{k} = - u_{k Z}^{(0)} \cdot H^{(1)} \cdot \sum_{j = 1}^{N} (X_{j} u_{j X}^{(0)} + Y_{j} u_{j Y}^{(0)} + Z_{j} u_{j Z}^{(0)} + P_{j} {\bar{u}}_{j Z}^{(0)}) \\ = - \sum_{j = 1}^{N} (V_{kjzx} X_{j} + V_{kjzy} Y_{j} + V_{kjzz} Z_{j}), \end{matrix}

(188)

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 $u_{k Y}^{(0)} \cdot J$ ⁠. Here, we will use $u_{k Y}^{(0)} \cdot J \cdot u_{k X}^{(0)} = - m_{k} Ω_{k}$ ⁠, which follows from Eqs. (167) and (63) or can be checked directly using Eqs. (160), (164), and (165). This projection results in

\begin{matrix} i ω m_{k} Ω_{k} X_{k} = - u_{k Y}^{(0)} \cdot H^{(1)} \cdot \sum_{j = 1}^{N} (X_{j} u_{j X}^{(0)} + Y_{j} u_{j Y}^{(0)} + Z_{j} u_{j Z}^{(0)} + P_{j} {\bar{u}}_{j Z}^{(0)}) \\ = - \sum_{j = 1}^{N} (V_{kjyx} X_{j} + V_{kjyy} Y_{j} + V_{kjyz} Z_{j}) . \end{matrix}

(189)

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 $- i ω X_{k}$ in the −x direction. Finally, act on Eq. (185) with $u_{k X}^{(0)} \cdot J$ ⁠. This projection results in

- i ω m_{k} Ω_{k} Y_{k} = - \sum_{j = 1}^{N} (V_{kjxx} X_{j} + V_{kjxy} Y_{j} + V_{kjxz} Z_{j}) .

(190)

The right hand side is the force in the x direction on particle k. This force produces an E × B drift velocity $- i ω Y_{k}$ 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 $V_{kjxz} = V_{kjyz} = 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, $Z_{j} = Z_{j}^{slow} + Z_{j}^{fast}$ ⁠. The slow evolution, on the E × B timescale, is conditioned on the axial particle positions being in axial force balance:

0 = - \sum_{j = 1}^{N} (V_{kjzx} X_{j} + V_{kjzy} Y_{j} + V_{kjzz} Z_{j}^{slow}) .

(191)

We can regard this force balance condition as a set of coupled linear equations for the slow axial displacements, written in vector form as $V_{z x} \cdot X + V_{z y} \cdot Y + V_{z z} \cdot Z^{slow} = 0$ ⁠, where the tensor $V_{x z} = \partial Φ / \partial X \partial Z$ and similarly for the other terms. This equation can be solved by matrix inversion,

Z^{slow} = - V_{z z}^{- 1} \cdot (V_{z x} \cdot X + V_{z y} \cdot Y) .

(192)

The fast motion in z is then evaluated using Eqs. (187) and (188) after canceling the slow terms:

\begin{matrix} - i ω Z_{k}^{fast} = P_{k} / m_{k}, \\ - i ω P_{k} = - \sum_{j = 1}^{N} V_{kjzz} Z_{j}^{fast} . \end{matrix}

(193)

This is the reduced eigenvalue problem for the axial modes.¹⁰ 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 $ω^{2}$ ⁠, of the familiar type encountered in the textbook Lagrangian theory of coupled oscillators,

ω^{2} M_{1} \cdot Z_{ω}^{fast} = V_{z z} \cdot Z_{ω}^{fast},

(194)

where $Z_{ω}^{fast}$ is an eigenvector of axial displacements and the matrix $M_{1}$ is the diagonal mass matrix of dimension N with diagonal elements $m_{j}, j = 1, \dots, N$ ⁠. The N eigenvalues $ω^{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_{ω}^{fast}$ are real and orthogonal with respect to the inner product ${(a, b)}_{z} = a \cdot M_{1} \cdot b$ ⁠. 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.

FIG. 8.

View large Download slide

Fractional difference between axial frequencies ω evaluated using Eq. (194) and the exact frequencies ω_exact evaluated using Eq. (5) and displayed in Fig. 4, for the same parameters as in that figure.

The axial eigenvectors diagonalize the axial energy H_z where

H_{z} = \frac{1}{2} \sum_{j = 1}^{N} P_{j}^{2} / m_{j} + \frac{1}{2} \sum_{j, k = 1}^{N} Z_{j} V_{jkzz} Z_{k} .

(195)

Writing $Z_{j} = \sum_{ω} a_{ω} Z_{j}^{fast}$ ⁠, where $a_{ω}$ 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

H_{z} = \sum_{ω} a_{ω}^{2} ω^{2} {(Z^{fast}, Z^{fast})}_{z} .

(196)

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 2N dimensional transverse displacement eigenvector $R_{⊥ ω} = (X, Y)$ ⁠, the equations become

- i ω R_{⊥ ω} = D_{⊥} \cdot R_{⊥ ω},

(197)

where $D_{⊥} = G \cdot J \cdot V_{⊥}$ is the dynamical matrix for E × B drift modes, $G$ is a diagonal $2 N \times 2 N$ matrix with diagonal elements $({(m_{1} Ω_{1})}^{- 1}, \dots, {(m_{N} Ω_{N})}^{- 1}, {(m_{1} Ω_{1})}^{- 1}, \dots, {(m_{N} Ω_{N})}^{- 1})$ ⁠, $J$ is the $2 N \times 2 N$ fundamental symplectic matrix [see Eq. (3)], and $V_{⊥}$ is the following symmetric $2 N \times 2 N$ potential energy tensor,

V_{⊥} = (\begin{matrix} V_{x x} & V_{x y} \\ V_{y x} & V_{y y} \end{matrix}) - (\begin{matrix} V_{x z} \cdot V_{z z}^{- 1} \cdot V_{z x} & V_{x z} \cdot V_{z z}^{- 1} \cdot V_{z y} \\ V_{y z} \cdot V_{z z}^{- 1} \cdot V_{z x} & V_{y z} \cdot V_{z z}^{- 1} \cdot V_{z y} \end{matrix}) .

(198)

This tensor determines the potential energy V in an E × B displacement of the form $(X, Y, Z^{slow})$ through the expression

V = \frac{1}{2} (X, Y) \cdot V_{⊥} \cdot (X, Y) .

(199)

This can be proven by writing V as

V = \frac{1}{2} (X, Y, Z^{slow}) \cdot (\begin{matrix} V_{x x} & V_{x y} & V_{x z} \\ V_{y x} & V_{y y} & V_{y z} \\ V_{z x} & V_{z y} & V_{z z} \end{matrix}) \cdot (X, Y, Z^{slow}),

(200)

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_{⊥}$ 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 $i D_{⊥}$ is Hermitian with respect to the inner product ${(a, b)}_{⊥} = a^{*} \cdot V_{⊥} \cdot b$ for any vector a and b. This follows from the fact that the matrix $L = V_{⊥} \cdot G \cdot J \cdot V_{⊥}$ is antisymmetric, which in turn follows from the symmetry and antisymmetry, respectively, of $V_{⊥}$ and J, along with the fact that $G \cdot J = J \cdot G$ ⁠. 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_{⊥}$ are real, the nonzero frequency modes come in $\pm ω$ 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.

FIG. 9.

View large Download slide

Fractional difference between E × B frequencies ω evaluated using Eq. (197) and the exact frequencies ω_exact evaluated using Eq. (5) and displayed in Fig. 4, for the same parameters as in that figure.

The energy of an E × B mode can be written in terms of the eigenvectors $R_{⊥ ω}$ ⁠. Since the E × B modes form an orthogonal set, a general E × B displacement of the form $\sum_{ω > 0} a_{ω} (R_{⊥ ω}, Z^{slow}) + c . c .$ ⁠, where $a_{ω}$ is the amplitude of a given mode, produces a diagonalized potential energy V given by

V = \sum_{ω > 0} | a_{ω} |^{2} R_{⊥ ω}^{*} \cdot V_{⊥} \cdot R_{⊥ ω} .

(201)

The kinetic energy contribution to E × B modes is negligible, scaling as $1 / B^{2}$ ⁠. 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.⁴⁷

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.¹⁷ When cooling and heating are applied, modes with disparate frequencies can come to quite different equilibrium energies, depending on the details of the driving.⁴⁸ 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 $β = 3 / 4$ crystal shown in Fig. 2, when Ω = 0 (the Brillouin limit), the rotation rate is $ω_{r} = \sqrt{5 / 4} ω_{z} = 1.118 ω_{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 $Ω = 20 ω_{z}$ ⁠, ω_r is either $20.0623 ω_{z}$ or $0.0623 ω_{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.

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 $V_{j k} = \partial^{2} Φ / \partial R_{j} \partial R_{k}$ for charges j and k with equilibrium positions R_j and R_k, respectively. The system potential energy $Φ$ is given by Eq. (115), and substituting for this, we obtain

V_{j k} = δ_{j k} \frac{\partial^{2}}{\partial R_{j}^{2}} (ϕ_{j} (R_{j}, Z_{j}) + \sum_{\begin{matrix} l = 1 \\ l \neq j \end{matrix}}^{N} \frac{q_{j} q_{l}}{| R_{j} - R_{l} |}) + (1 - δ_{j k}) \frac{\partial^{2}}{\partial R_{j} \partial R_{k}} \frac{q_{j} q_{k}}{| R_{j} - R_{k} |} .

(A1)

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

\begin{matrix} V_{j k} = δ_{j k} (\hat{r} \hat{r} \frac{\partial^{2}}{\partial R_{j}^{2}} + \hat{z} \hat{z} \frac{\partial^{2}}{\partial Z_{j}^{2}} + (\hat{r} \hat{z} + \hat{z} \hat{r}) \frac{\partial^{2}}{\partial R_{j} \partial Z_{j}}) ϕ_{j} (R_{j}, Z_{j}) \\ + δ_{j k} \sum_{\overset{l = 1}{l \neq j}}^{N} q_{j} q_{l} (\frac{3 (R_{j} - R_{l}) (R_{j} - R_{l})}{| R_{j} - R_{l} |^{3}} - \frac{1}{| R_{j} - R_{l} |^{5}}) \\ - (1 - δ_{j k}) q_{j} q_{k} (\frac{3 (R_{j} - R_{k}) (R_{j} - R_{k})}{| R_{j} - R_{k} |^{3}} - \frac{1}{| R_{j} - R_{k} |^{5}}), \end{matrix}

(A2)

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 $δ v_{j}, j = 1, \dots, N$ ⁠, are

δ {\dot{r}}_{i} = δ v_{i},

(B1)

m_{i} δ {\dot{v}}_{i} = - \sum_{j} V_{i j} \cdot δ r_{j} - m_{i} Ω_{i} δ v_{i} \times \hat{z},

(B2)

where the 3 × 3 symmetric tensor $V_{i j} = \nabla_{i} \nabla_{j} Φ$ ⁠.

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 6N dimensional configuration vector $η = (δ r_{1}, \dots, δ r_{N}, δ v_{1}, \dots, δ v_{N})$ ⁠, as was done in the Hamiltonian method. Then, the linearized equations of motion can be written, in analogy to Eq. (2), as

\dot{η} = D' \cdot η,

(B3)

where the $6 N \times 6 N$ dynamical matrix $D'$ consists of four $3 N \times 3 N$ blocks,

D' = (\begin{matrix} 0 & 1 \\ - M^{- 1} \cdot V & - 2 Ω \end{matrix}),

(B4)

where the Lorentz tensor $Ω$ ⁠, 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

E = \frac{1}{2} η \cdot E \cdot η,

(B5)

where the symmetric energy matrix $E$ is

E = (\begin{matrix} V & 0 \\ 0 & M \end{matrix}) .

(B6)

Eigenmodes of the form $η (t) = exp (- i ω t) ψ_{ω}$ satisfy the eigenvalue problem

- i ω ψ_{ω} = D' \cdot ψ_{ω} .

(B7)

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 $i D'$ is Hermitian, with respect to a modified inner product defined by $(a, b)' = a^{*} \cdot E \cdot b$ 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 \cdot D'$ ⁠. One can show that this matrix is antisymmetric by direct calculation using Eqs. (B4) and (B6):

L = (\begin{matrix} 0 & V \\ - V & - 2 M \cdot Ω \end{matrix}) .

(B8)

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^{*} \cdot E \cdot i D' \cdot b = {[b^{*} \cdot E \cdot i D^{'} \cdot a]}^{*}$ ⁠, so $i D^{'}$ 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)'$ ⁠, 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 $η$ in terms of the complete set of orthogonal eigenvectors

η = \sum_{ω} a_{ω} ψ_{ω},

(B9)

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

E = \sum_{ω > 0} | a_{ω} |^{2} (ψ_{ω}, ψ_{ω})' .

(B10)

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 $({\bar{u}}_{0 z}, {\bar{u}}_{0 z})$ ⁠, we require a solution for the vector ${\bar{u}}_{0 z}$ of the equation

H \cdot {\bar{u}}_{0 z} = u_{0 z} \cdot J,

(C1)

where $u_{0 z}$ 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 ${\bar{u}}_{0 z} = (\bar{r}, \bar{p})$ ⁠, which when used in Eq. (C1) yields two coupled equations

(V + C) \cdot \bar{r} + Ω \cdot \bar{p} = - p_{0 z} = - \frac{1}{2} m Ω R_{⊥},

(C2)

Ω^{t r} \cdot \bar{r} + m^{- 1} \bar{p} = r_{0 z} = \hat{z} \times R,

(C3)

where $R_{⊥}$ 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 $Ω$ ⁠, and the definition of $C$ ⁠, we can write these equations as

V \cdot \bar{r} + \frac{1}{4} m Ω^{2} {\bar{r}}_{⊥} - \frac{1}{2} Ω \hat{z} \times \bar{p} = - \frac{1}{2} m Ω R_{⊥},

(C4)

\frac{1}{2} Ω \hat{z} \times \bar{r} + m^{- 1} \bar{p} = \hat{z} \times R .

(C5)

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

\hat{z} \times \bar{p} = - m R_{⊥} + \frac{1}{2} m Ω {\bar{r}}_{⊥} .

(C6)

Applying this result to Eq. (C4) implies

- V \cdot \bar{r} = m Ω R_{⊥} .

(C7)

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

\bar{p} = m \hat{z} \times R - \frac{1}{2} m Ω \hat{z} \times \bar{r} .

(C8)

When these results are used to calculate the rotational inertia ${\bar{u}}_{0 z} \cdot H \cdot {\bar{u}}_{0 z}$ ⁠, the result is

({\bar{u}}_{0 z}, {\bar{u}}_{0 z}) = m R_{⊥} \cdot R_{⊥} + \bar{r} \cdot V \cdot \bar{r} .

(C9)

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 $- V \cdot R_{⊥} = - 3 m ω_{⊥}^{2} R_{⊥}$ ⁠. Therefore, the solution of Eq. (C7) is

\bar{r} = - \frac{Ω}{3 ω_{⊥}^{2}} R_{⊥},

(C10)

and together with Eq. (C5) this implies

\bar{p} = m \hat{z} \times R (1 + \frac{Ω^{2}}{6 ω_{⊥}^{2}}) .

(C11)

These are the results quoted in Eq. (126). The rotational inertia, Eq. (127), follows from the substitution of Eq. (C10) into Eq. (C9), using Eq. (C7).

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 $ψ = (c, c^{*})$ rather than the phase space coordinates $z = (q, p)$ ⁠. In terms of these pairs, the Hamiltonian for a linearized system is

H = \frac{1}{2} ψ^{*} \cdot h \cdot ψ,

(D1)

where the matrix $h$ is Hermitian, $h^{†} = h$ ⁠, and has the symmetric block form

h = (\begin{matrix} A & B \\ B^{*} & A^{*} \end{matrix}),

(D2)

where the N × N matrix $A = A^{†}$ is Hermitian and the N × N matrix $B = B^{t r}$ is symmetric. The creation and annihilation pairs are related to $z$ via the linear transformation

ψ = T \cdot z,

(D3)

where

T = \frac{1}{\sqrt{2}} (\begin{matrix} 1 & i 1 \\ 1 & - i 1 \end{matrix}),

(D4)

or in the component form $c = (q + i p) / \sqrt{2}, c^{*} = (q - i p) / \sqrt{2}$ ⁠. These creation/annihilation pairs have Poisson bracket relationships that may be succinctly expressed by the equation

[ψ, ψ^{*}] = - i σ,

(D5)

where the matrix σ is defined as

σ = (\begin{matrix} 1 & 0 \\ 0 & - 1 \end{matrix}) .

(D6)

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

H = T^{†} \cdot h \cdot T .

(D7)

The standard Bogoliubov approach is to find a transformation to new creation–annihilation pairs $ϕ = (a, a^{*})$ ⁠,

ψ = s \cdot ϕ,

(D8)

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

H = \frac{1}{2} ϕ^{*} \cdot k \cdot ϕ,

(D9)

with the new Hamiltonian matrix $k$ given by

k = s^{†} \cdot h \cdot s = diagonal .

(D10)

The added requirement that the transformation be canonical requires $[ϕ, ϕ^{*}] = - i σ$ ⁠, which using Eq. (D8) implies that $s$ must satisfy the symplectic condition in the creation/annihilation representation,

s \cdot σ \cdot s^{†} = σ .

(D11)

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):

\dot{ψ} = [ψ, H] = [ψ, ψ^{*}] \cdot h \cdot ψ = - i σ \cdot h \cdot ψ = - i d \cdot ψ,

(D12)

where we introduce the quantum dynamical matrix $d \equiv σ \cdot h$ ⁠. This dynamical matrix is analogous to the matrix $D$ introduced in Sec. II. Now, consider eigenmodes of the form $ψ = exp (- i ω t) w_{ω}$ for a mode of frequency ω, where $w_{ω}$ is the associated vector. Equation (D12) implies that these vectors are eigenvectors of $d$ with eigenvalues ω:

d \cdot w_{ω} = ω w_{ω} .

(D13)

These eigenvectors are related to the eigenvectors $u_{ω}$ of the dynamical matrix $D$ by Eq. (D3), which implies that

w_{ω} = T \cdot u_{ω} .

(D14)

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

The eigenvectors $w_{ω}$ form an orthogonal set with respect to a generalized inner product defined for any complex vectors $a$ and $b$ as $(a, b) \equiv a^{*} \cdot h \cdot b$ ⁠:
$(w_{ω}, w_{\bar{ω}}) = 0$ provided that $ω \neq {\bar{ω}}^{*}$ ⁠.
A given eigenvalue ω is real, provided that the corresponding eigenvector satisfies $(w_{ω}, w_{ω}) \neq 0$ ⁠.
For each eigenmode $(ω, w_{ω})$ for which $ω \neq 0$ ⁠, there is a second eigenmode $(- ω^{*}, w_{- ω^{*}})$ for which $w_{- ω^{*}} = Λ \cdot w_{ω}^{*}$ ⁠, where the matrix $Λ$ is
$Λ = (\begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix}) .$
(D15)

Thus, for real ω, the $ω \neq 0$ eigenmodes come in $\pm ω$ 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:

(a, d \cdot b) = {(b, d \cdot a)}^{*} .

(D16)

This requires that the matrix $l \equiv h \cdot d$ is a Hermitian matrix: $l = l^{†}$ ⁠, which can be proven using the same set of steps as shown in Eq. (8):

l_{j i} = h_{j k} σ_{k l} h_{l i} = h_{k j}^{*} σ_{l k} h_{i l}^{*} = h_{i l}^{*} σ_{l k} h_{k j} = l_{i j}^{*} .

(D17)

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):

h = Λ \cdot h^{*} \cdot Λ .

(D18)

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

\begin{matrix} d = σ \cdot Λ \cdot σ \cdot σ \cdot h^{*} \cdot Λ \\ = - Λ \cdot d^{*} \cdot Λ, \end{matrix}

(D19)

where we used the identity $σ \cdot Λ \cdot σ = - 1$ ⁠. Substituting Eq. (D19) into Eq. (D13), acting on both sides with $Λ$ and using $Λ \cdot Λ = 1$ then yields

d^{*} \cdot Λ \cdot w_{ω} = - ω Λ \cdot w_{ω} .

(D20)

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_{ω}, w_{ω}) \neq 0$ for all eigenmodes, so that all eigenfrequencies are real and are also nonzero and non-degenerate. First, we write a general vector $ψ = (c, c^{*})$ in terms of a linear combination of the eigenvectors $w_{ω}$ ⁠:

ψ (t) = \sum_{ω} a_{ω} (t) w_{ω} .

(D21)

This is merely another way to express the Bogoliubov transformation Eq. (D8), taking the components of the vector $ϕ (t)$ to be $a_{ω} (t)$ and the transformation matrix $s$ to have columns given by the eigenvectors:

s = (w_{ω_{1}}, w_{ω_{2}}, \dots, w_{ω_{2 N}}) .

(D22)

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

H = \frac{1}{2} \sum_{ω} a_{ω} a_{ω}^{*} (w_{ω}, w_{ω}) .

(D23)

To make further progress, we order the eigenfrequencies such that frequencies ω₁ 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

s = (w_{ω_{1}}, \dots, w_{ω_{N}}, Λ \cdot w_{ω_{1}}^{*}, \dots, Λ \cdot w_{ω_{N}}^{*}),

(D24)

and that we can write Eq. (D21) as

ψ (t) = \sum_{ω > 0} (a_{ω} (t) w_{ω} + a_{- ω} (t) Λ \cdot w_{ω}^{*}) .

(D25)

In order for this equation to match Eq. (D8) with $ϕ = (a, a^{*})$ ⁠, this requires

a_{- ω} = a_{ω}^{*},

(D26)

which implies that $ϕ$ has the required form

ϕ = (a_{ω_{1}}, \dots, a_{ω_{N}}, a_{ω_{1}}^{*}, \dots, a_{ω_{N}}^{*}) .

(D27)

We can use this result to simplify the Hamiltonian, summing only over positive frequencies $ω_{1}, \dots, ω_{N}$ ⁠:

H = \frac{1}{2} \sum_{ω > 0} a_{ω} a_{ω}^{*} {(w_{ω}, w_{ω}) + (Λ \cdot w_{ω}^{*}, Λ \cdot w_{ω}^{*})} .

(D28)

This can be further simplified using the identity

(Λ \cdot w_{ω}^{*}, Λ \cdot w_{ω}^{*}) = (w_{ω}, w_{ω}),

(D29)

which yields the simplified diagonalized Hamiltonian

H = \sum_{ω > 0} a_{ω} a_{ω}^{*} (w_{ω}, w_{ω}) .

(D30)

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

(Λ \cdot w_{ω}^{*}, Λ \cdot w_{ω}^{*}) = w_{ω} \cdot Λ \cdot h \cdot Λ \cdot w_{ω}^{*} = w_{ω} \cdot h^{*} \cdot w_{ω}^{*} = w_{ω}^{*} \cdot h \cdot w_{ω},

(D31)

where we used $Λ^{t r} = Λ$ in the first step and $h^{†} = h$ in the last step.

Next, we ensure that the transformation to the new $ϕ$ variables is canonical by requiring that their Poisson brackets satisfy $[ϕ, ϕ^{*}] = - i σ$ ⁠. In the component form, this requires $[a_{ω}, a_{\bar{ω}}^{*}] = - i δ_{ω, \bar{ω}}$ and $[a_{ω}, a_{\bar{ω}}] = 0$ (for ω and $\bar{ω}$ 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

a_{ω} (t) = \frac{(w_{ω}, ψ (t))}{(w_{ω}, w_{ω})} .

(D32)

Applying this to $[a_{ω}, a_{\bar{ω}}^{*}]$ ⁠, we use Eq. (D5) to obtain

\begin{matrix} [a_{ω}, a_{\bar{ω}}^{*}] = \frac{w_{ω}^{*} \cdot h \cdot (- i σ) \cdot h \cdot w_{\bar{ω}}}{(w_{ω}, w_{ω}) (w_{\bar{ω}}, w_{\bar{ω}})} = \frac{w_{ω}^{*} \cdot h \cdot (- i d) \cdot w_{\bar{ω}}}{(w_{ω}, w_{ω}) (w_{\bar{ω}}, w_{\bar{ω}})} \\ = \frac{w_{ω}^{*} \cdot h \cdot (- i \bar{ω} w_{\bar{ω}})}{(w_{ω}, w_{ω}) (w_{\bar{ω}}, w_{\bar{ω}})} \\ = \frac{- i ω δ_{ω, \bar{ω}}}{(w_{ω}, w_{ω})}, \end{matrix}

(D33)

where in the first step we used $h^{†} = 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_{ω}, a_{\bar{ω}}] = 0$ when ω and $\bar{ω}$ are greater than zero. Thus, the transformation is canonical, provided that we normalize the eigenvectors so that

(w_{ω}, w_{ω}) = ω, ω > 0,

(D34)

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

H = \sum_{ω > 0} a_{ω} a_{ω}^{*} ω .

(D35)

We can more directly connect this approach to the standard Bogoliubov approach by applying the vector $w_{\bar{ω}}^{*} \cdot σ$ to Eq. (D13):

w_{\bar{ω}}^{*} \cdot σ \cdot d \cdot w_{ω} = ω w_{\bar{ω}}^{*} \cdot σ \cdot w_{ω} .

(D36)

Using $σ \cdot σ = 1$ and $d = σ \cdot h$ then allows the left hand side to be written as an inner product:

(w_{\bar{ω}}, w_{ω}) = ω w_{\bar{ω}}^{*} \cdot σ \cdot w_{ω} .

(D37)

The orthogonality of the eigenmodes then yields

(w_{ω}, w_{ω}) δ_{ω, \bar{ω}} = ω w_{\bar{ω}}^{*} \cdot σ \cdot w_{ω} .

(D38)

For $ω > 0$ ⁠, we can substitute for the inner product using Eq. (D34), while for $ω < 0$ we can employ Eq. (D29) to see that $(w_{ω}, w_{ω}) = - ω, ω < 0$ ⁠, which when used in Eq. (D38) yields

w_{\bar{ω}}^{*} \cdot σ \cdot w_{ω} = sign (ω) δ_{ω, \bar{ω}} .

(D39)

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

s^{†} \cdot σ \cdot s = σ .

(D40)

Taking the complex conjugate of this equation and using $σ^{t r} = σ$ 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 $ϕ = T \cdot Z$ in Eq. (D8), which implies

S = T^{- 1} \cdot s \cdot T .

(D41)

Since T is a unitary transformation, $T^{- 1} = T^{†}$ and Eqs. (D14) and (D24) imply that $T^{- 1} \cdot s = (U, U^{*})$ ⁠, which when applied to Eq. (D41) leads back to Eq. (28), $S = \sqrt{2} (Re U, - Im U)$ ⁠.

References

1.

J. W.

Britton

,

B. C.

Sawyer

,

A. C.

Keith

,

C.-C. J.

Wang

,

J. K.

Freericks

,

H.

Uys

,

M. J.

Biercuk

, and

J. J.

Bollinger

,

Nature

484

,

489

(

2012

).

https://doi.org/10.1038/nature10981

Google Scholar

Crossref

PubMed

2.

B. C.

Sawyer

,

J. W.

Britton

,

A. C.

Keith

,

C.-C. J.

Wang

,

J. K.

Freericks

,

H.

Uys

,

M. J.

Biercuk

, and

J. J.

Bollinger

,

Phys. Rev. Lett.

108

,

213003

(

2012

).

https://doi.org/10.1103/PhysRevLett.108.213003

Google Scholar

Crossref

PubMed

3.

B. C.

Sawyer

,

J. W.

Britton

, and

J. J.

Bollinger

,

Phys. Rev. A

89

,

033408

(

2014

).

https://doi.org/10.1103/PhysRevA.89.033408

Google Scholar

Crossref

4.

J. G.

Bohnet

,

B. C.

Sawyer

,

J. W.

Britton

,

M. L.

Wall

,

A. M.

Rey

,

M.

Foss-Feig

, and

J. J.

Bollinger

,

Science

352

,

1297

(

2016

).

https://doi.org/10.1126/science.aad9958

Google Scholar

Crossref

PubMed

5.

K. A.

Gilmore

,

J. G.

Bohnet

,

B. C.

Sawyer

,

J. W.

Britton

, and

J. J.

Bollinger

,

Phys. Rev. Lett.

118

,

263602

(

2017

).

https://doi.org/10.1103/PhysRevLett.118.263602

Google Scholar

Crossref

PubMed

6.

B. J.

McMahon

,

C.

Volin

,

W. G.

Rellergert

, and

B. C.

Sawyer

,

Phys. Rev. A

101

,

013408

(

2020

).

https://doi.org/10.1103/PhysRevA.101.013408

Google Scholar

Crossref

7.

S. L.

Gilbert

,

J. J.

Bollinger

, and

D. J.

Wineland

,

Phys. Rev. Lett.

60

,

2022

(

1988

).

https://doi.org/10.1103/PhysRevLett.60.2022

Google Scholar

Crossref

PubMed

8.

J. N.

Tan

,

J. J.

Bollinger

,

B.

Jelenkovic

, and

D. J.

Wineland

,

Phys. Rev. Lett.

75

,

4198

(

1995

).

https://doi.org/10.1103/PhysRevLett.75.4198

Google Scholar

Crossref

PubMed

9.

T. B.

Mitchell

,

J. J.

Bollinger

,

D. H. E.

Dubin

,

X.-P.

Huang

,

W. M.

Itano

, and

R. H.

Baughman

,

Science

282

,

1290

(

1998

).

https://doi.org/10.1126/science.282.5392.1290

Google Scholar

Crossref

PubMed

10.

D. H. E.

Dubin

and

J. P.

Schiffer

,

Phys. Rev. E

53

,

5249

(

1996

).

https://doi.org/10.1103/PhysRevE.53.5249

Google Scholar

Crossref

11.

H.

Fukuyama

,

Solid State Commun.

17

,

1323

(

1975

).

https://doi.org/10.1016/0038-1098(75)90696-1

Google Scholar

Crossref

12.

H.

Bonsall

and

A. A.

Maradudin

,

Phys. Rev. B

15

,

1959

(

1977

).

https://doi.org/10.1103/PhysRevB.15.1959

Google Scholar

Crossref

13.

T.

Nagai

and

H.

Fukuyama

,

J. Phys. Soc. Jpn.

51

,

3431

(

1982

).

https://doi.org/10.1143/JPSJ.51.3431

Google Scholar

Crossref

14.

K. I.

Golden

,

G.

Kalman

, and

P.

Wynns

,

Phys. Rev. B

48

,

8882

(

1993

).

https://doi.org/10.1103/PhysRevB.48.8882

Google Scholar

Crossref

15.

T.

Ott

,

H.

Kahlert

,

A.

Reynolds

, and

M.

Bonitz

,

Phys. Rev. Lett.

108

,

255002

(

2012

).

https://doi.org/10.1103/PhysRevLett.108.255002

Google Scholar

Crossref

PubMed

16.

N. A.

Usov

,

Y. B.

Grebenshchikov

, and

F. R.

Ulinich

,

Sov. Phys. JETP

51

,

148

(

1980

).

Google Scholar

17.

S.-J.

Chen

, Ph.D. thesis (

University of California

,

San Diego

,

1994

).

18.

D. A.

Baiko

,

Phys. Rev. E

80

,

046405

(

2009

).

https://doi.org/10.1103/PhysRevE.80.046405

Google Scholar

Crossref

19.

N. N.

Bogoliubov

,

J. Phys. USSR

11

,

23

(

1947

).

Google Scholar

20.

T.

Holstein

and

H.

Primakoff

,

Phys. Rev.

58

,

1098

(

1940

).

https://doi.org/10.1103/PhysRev.58.1098

Google Scholar

Crossref

21.

A.

Fetter

and

J.

Walecka

,

Quantum Theory of Many Particle Systems

(

Dover

,

New York

,

2003

).

Google Scholar

22.

C. C.

Wang

,

A. C.

Keith

, and

J. K.

Freericks

,

Phys. Rev. A

87

,

029901

(

2013

).

https://doi.org/10.1103/PhysRevA.87.029901

Google Scholar

Crossref

23.

D. H. E.

Dubin

,

Phys. Rev. E

53

,

5268

(

1996

).

https://doi.org/10.1103/PhysRevE.53.5268

Google Scholar

Crossref

24.

N.

Bohr

, Doctoral dissertation (

Neils Bohr Collected Works

, edited by

L.

Rosenfeld

and

J.

Rud Nielsen

,

Elsevier

,

1972

).

Google Scholar

25.

H.

Goldstein

,

Classical Mechanics

(

Addison-Wesley

,

Reading, MA

,

1950

), Chap. 10.

Google Scholar

26.

K. S.

Fine

,

A. C.

Cass

,

W. G.

Flynn

, and

C. F.

Driscoll

,

Phys. Rev. Lett.

75

,

3277

(

1995

).

https://doi.org/10.1103/PhysRevLett.75.3277

Google Scholar

Crossref

PubMed

27.

L. J.

Campbell

and

R. M.

Ziff

,

Phys. Rev. B

20

,

1886

(

1979

).

https://doi.org/10.1103/PhysRevB.20.1886

Google Scholar

Crossref

28.

H.

Goldstein

,

C.

Poole

, and

J.

Safko

,

Classical Mechanics

(

Addison-Wesley

,

Reading, MA

,

2002

), Chap. 9.4.

Google Scholar

29.

D. H. E.

Dubin

,

Phys. Rev. Lett.

121

,

015001

(

2018

).

https://doi.org/10.1103/PhysRevLett.121.015001

Google Scholar

Crossref

PubMed

30.

D. H. E.

Dubin

,

Phys. Plasmas

26

,

102111

(

2019

).

https://doi.org/10.1063/1.5116375

Google Scholar

Crossref

31.

M. G.

Krein

,

Transl. Am. Math. Soc.

120

(

2

),

139

(

1983

).

https://doi.org/10.1090/trans2/120

Google Scholar

Crossref

32.

V. I.

Arnold

,

Mathematical Methods of Classical Mechanics

(

Springer-Verlag

,

Newy York

,

1980

), p.

225

.

Google Scholar

33.

R.

Johnson

,

R.

Obaya

,

S.

Novo

,

C.

Nunez

, and

E. R.

Fabbri

,

Nonautonomous Linear Hamiltonian Systems: Oscillation, Spectral Theory, and Control

(

Springer

,

New-York

,

2016

).

Google Scholar

Crossref

34.

P. C.

Shields

,

Elementary Linear Algebra

(

Worth, New York

,

1973

).

Google Scholar

35.

R. K.

Pathria

,

Statistical Mechanics

(

Pergamon Press

,

Oxford

,

1972

), p.

74

.

Google Scholar

36.

R.

Rafac

,

J. P.

Schiffer

,

J. S.

Hangst

,

D. H. E.

Dubin

, and

D. J.

Wales

,

Proc. Natl. Acad. Sci. U. S. A.

88

,

483

(

1991

).

https://doi.org/10.1073/pnas.88.2.483

Google Scholar

Crossref

PubMed

37.

D. H. E.

Dubin

and

T. M.

O'Neil

,

Phys. Rev. Lett.

60

,

511

(

1988

).

https://doi.org/10.1103/PhysRevLett.60.511

Google Scholar

Crossref

PubMed

38.

A.

Rahman

and

J. P.

Schiffer

,

Phys. Rev. Lett.

57

,

1133

(

1986

).

https://doi.org/10.1103/PhysRevLett.57.1133

Google Scholar

Crossref

PubMed

39.

D. H. E.

Dubin

and

T. M.

O'Neil

,

Rev. Mod. Phys.

71

,

87

(

1999

).

https://doi.org/10.1103/RevModPhys.71.87

Google Scholar

Crossref

40.

D. H. E.

Dubin

,

Phys. Rev. Lett.

66

,

2076

(

1991

).

https://doi.org/10.1103/PhysRevLett.66.2076

Google Scholar

Crossref

PubMed

41.

D. J.

Wineland

,

J. C.

Bergquist

,

W. M.

Itano

,

J. J.

Bollinger

, and

C. H.

Manney

,

Phys. Rev. Lett.

59

,

2935

(

1987

).

https://doi.org/10.1103/PhysRevLett.59.2935

Google Scholar

Crossref

PubMed

42.

E. A.

Cornell

,

K. R.

Boyce

,

D.

Fygenson

, and

D. E.

Pritchard

,

Phys. Rev. A

45

,

3049

(

1992

).

https://doi.org/10.1103/PhysRevA.45.3049

Google Scholar

Crossref