Normal modes of vibrations of an ion crystal

WARNING! This posts uses heavy 3D animations so it might be running slowly (or not at all) on “old” devices

Continuing the series of posts about trapped ions in Paul traps let's expand more on eigenfrequencies and eigenmodes of vibration of an ion crystal. Eigenfrequencies (also known as natural frequences) characterize how a complex system with multiple degrees of freedom oscillates naturally. For a system of $N$ particles, such as an ion crystal, the motion can be decomposed into two components: the center-of-mass motion and the internal motion. The internal motion encompasses vibrations (oscillatory motion), rotations, and deformations, depending on the system's physical properties. In ion crystals specifically, the internal motion is primarily characterized by vibrational eigenmodes, which describe how the ions collectively oscillate relative to their center of mass.

Loading 3D ion crystal simulation…

The potential energy of an ion crystal in a Paul trap can be expressed as:

V_{\mathrm{cry}} = \sum\limits_{i=1}^N \frac{1}{2} m \mathbf{\Omega}^2 \cdot \mathbf{r}_i^2 + \sum\limits_{\substack{j, k=1\\j \neq k}}^N \frac{Z^2 e^2}{8\pi\epsilon_0} \frac{1}{\lvert \mathbf{r}_{j} - \mathbf{r}_{k}\rvert}.

(1)

The first term of Eq. 1 represents the confining potential created by the Paul trap itself. Each ion of mass $m$ experiences a harmonic trapping potential characterized by the trap frequencies $\mathbf{\Omega} = \left[\omega_x, \omega_y, \omega_z\right]$ along the three spatial directions. The position of each ion is given by its vector $\mathbf{r} = \left[x, y, z\right]$ .

The second term of Eq. 1 represents the Coulomb interaction between ion pairs, where ions with charge $Ze$ repel each other with a force proportional to their inverse separation distance $|\mathbf{r}_j - \mathbf{r}_k|$ .

Quantum computing applications requires reaching the limit where the potential generated by the Paul trap and Coulomb repulsion can be quantized. This necessitates cooling the ions to millikelvin temperatures, achievable through laser cooling without requiring cryogenic apparatus. However, cooling the experimental setup to a few kelvin is often advantageous, as it enhances stability and ion trap lifetimes. While this added cooling increases the complexity of the experimental setup, it is not strictly necessary for performing quantum operations.

At low enough temperatures we can assume that the ion positions can be expressed as:

\mathbf{r}_i\left(t\right) \approx \mathbf{r}_i^{\left(0\right)} + \mathbf{q}_i\left(t\right),

(2)

where $\mathbf{r}_i^{\left(0\right)}$ is the equilibrium position of the $i$ -th ion and $\mathbf{q}_i$ a small displacement. We can calculate the equlibrium positions of an ion crystal by imposing that the force on all ions at their equilibrium positions is zero:

\frac{\partial V_{\mathrm{cry}}}{\partial \mathbf{r}_i}\biggr\lvert_{\mathbf{r}_i = \mathbf{r}_i^{\left(0\right)}} = 0,

(3)

For $N= 2,3$ the equilibrium positions can be calculated analytically.

For non-expert, let's consider a simplified 1-dimensional model. In this example we take two ions in a harmonic potential generated by a Paul trap. You can read more about Paul traps here. The potential generated by a Paul trap can be written as:

V_{\mathrm{PT}} = m \omega^2 x^2

(1.a)

This generates a force on the ions $F_{\mathrm{PT}} = - \nabla V_{\mathrm{PT}}$ that is proportional to the ion charge ( $Z = 1$ for single charged ions) and applied external electric potential $\mathbb{E}$ .

F_{\mathrm{PT}} = - m \omega^2 x = - Z e \mathbb{E} x \quad \mathrm{with} \quad \omega = \sqrt{\frac{e \mathbb{E}}{m}}

(2.a)

Let's consider two ions in such a potential. We can find their equilibrium potential by calculating the position at which the Coulomb repulsion between ions and trapping force equilibrate. Let's define $d$ as the distance between the two ions at their minima $x^{(0)}_1$ for the first ion and $x^{(0)}_2$ for the second. As the ions are identical in mass and charge and the system is invariant under $1 \rightarrow 2$ we can write $x^{(0)}_1 = d/2$ and $x^{(0)}_2 = -d/2$ . The rest condition for the ions is such that $F_{\mathrm{PT}} = F_{\mathrm{Coulomb}}$ , for one ion this can be written as:

m \omega^2 \frac{d}{2} = \frac{Z e^2}{4 \pi \epsilon_0 d^2}

(3.a)

from which we can obtain the equlibrium position $x^{(0)}_1$ as:

x^{(0)}_1 = d/2 \quad \mathrm{so: } \quad x^{(0)}_1 = \left(\frac{Z e^2}{16 \pi \epsilon_0 m \omega^2}\right)^{1/3}

(3.a)

It is convenient to define a length scale $l$ as:

l = \left(\frac{Z e^2}{4 \pi \epsilon_0 m \omega^2}\right)^{1/3},

(4.a)

such that we can define the dimensionless equlibrium positions as $u_m = x^{(0)}_m/l$ . For the two ions case we just described we have $u_1 = (1/2)^{2/3}$ and $u_2 = -(1/2)^{2/3}$ . In the case of $N=3$ the equilibrium positions are $u_1 = (5/4)^{2/3}$ and $u_2 = 0$ and $u_2 = -(5/4)^{2/3}$

In case of more ions we have to use numerical methods ¹. We can extract the equilibrium position from the simulation of a crystal of $N$ ions at the top of this post. You can change the shape of the crystal varying the trapping potentials $\omega_x$ , $\omega_y$ and $\omega_z$ or increase the number of ions $N$ .

The interplay between the external trapping potential and the Coulomb repulsion results in stable Coulomb crystals. The dimensionality of the crystal depends on the relative strength of the trapping potential along the different axes

To find the equilibrium positions, a cooling force $C_{\mathrm{cool}}$ is introduced into the simulation. This parameter controls the rate at which the ion crystal cools. A higher cooling force results in faster cooling allowing a faster convergence of the equilibrium positions. However, higher cooling rates may cause defects in the final crystal if the system becomes trapped in a local minimum of the potential energy.

Notice that, as the equilibrium positions are calculated at each time step of the simulation it might take a few seconds to converge to the exact value.

Ions' equilibrium positions

Ion N.	$x^{\left(0\right)}/l_z$	$y^{\left(0\right)}/l_z$	$z^{\left(0\right)}/l_z$

The equilibrium positions $x^{\left(0\right)}$ , $y^{\left(0\right)}$ and $z^{\left(0\right)}$ are given in units of the length scale:

l_z = \left(\frac{Z e^2}{4 \pi \epsilon_0 m \omega_z^2}\right)^{1/3}

(4)

For small displacements $\mathbf{q}$ from the equilibrium positions $\mathbf{r}^0$ we can expand total potential as:

V \approx V_0 + \sum\limits_{j=1}^N \mathbf{q}_j\frac{\partial V_{\mathrm{tot}}}{\partial \mathbf{r}_j}\biggr\lvert_{\mathbf{r}_j=\mathbf{r}^0_j} + \frac{1}{2} \sum\limits_{j,i=1}^N \mathbf{q}_j \mathbf{q}_i \frac{\partial V_{\mathrm{tot}}}{\partial \mathbf{r}_j\partial\mathbf{r}_i}\biggr\lvert_{\substack{\mathbf{r}_j=\mathbf{r}^0_j \\ \mathbf{r}_i=\mathbf{r}^0_i}} + \mathcal{O}(\mathbf{q}^3).

(5)

Note that, we can ignore the first term $V_0$ as is nothing more than a global energy offset. Futhermore, the equilibrium condtion of Eq. 3 imposes that also the second term has to be zero. We can then write the Lagrangian of the system at equilibrium as:

\mathcal{L} = \frac{m}{2}\sum\limits_{j=1}^N \dot{\mathbf{q}}_j^2 - \frac{1}{2} \sum\limits_{j,i=1}^N \mathbf{q}_j \mathbf{q}_i \mathcal{A}_{i,j},

(6)

with

\mathcal{A}_{i,j} = \frac{\partial V_{\mathrm{tot}}}{\partial \mathbf{r}_j\partial\mathbf{r}_i}\biggr\lvert_{\substack{\mathbf{r}_j=\mathbf{r}^0_j \\ \mathbf{r}_i=\mathbf{r}^0_i}},

(7)

the Hessian matrix calculated at the equilibrium positions.

As we know the ions potential and their equilibrium position we can calculate the Hessian matrix numerically. As this is not a trivial calculation we won't be updating it live as we did for the equilibrium positions, one has to put some extra effort and click the following button to generate the Hessian matrix.

Please note that, while you can update the matrix by clicking the button multiple times, the matrix acquires physical sense only when the ions reached their equilibrium positions (eg. their positions stopped fluctuating in the first table).

Note that while the first term in Eq. 7 is already diagonal, the Hessian matrix is not. We can then find the energy eigenvalues of the system by diagonalizing it. Also, since the matrix $\mathcal{A}_{i,j}$ is real, symmetric, and non-negative definite we are guaranteed that its eigenvalues are non-negative.

Let $A \in \mathbb{R}^{3N \times 3N}$ be the Hessian at equilibrium with eigenpairs $(\nu_m^2, \mathbf{b}^{(m)})$ for $m = 1, \ldots, 3N$ , where the eigenvectors are orthonormal:

A\mathbf{b}^{(m)} = \nu_m^2\,\mathbf{b}^{(m)}, \qquad \mathbf{b}^{(m)} \cdot \mathbf{b}^{(n)} = \delta_{mn}.

(8)

Collecting the eigenvectors into the orthogonal matrix,

B = \big[\,\mathbf{b}^{(1)}\ \mathbf{b}^{(2)}\ \cdots\ \mathbf{b}^{(3N)}\,\big], \qquad B^{\top} A B = \operatorname{diag}(\nu_1^2, \ldots, \nu_{3N}^2), \qquad B^{\top} B = I,

(9)

we can then define the normal coordinates $\mathcal{Q} = (\mathcal{Q}_1, \ldots, \mathcal{Q}_{3N})^{\top}$ by the orthogonal change of variables

\mathbf{q} = B\,\mathcal{Q}, \qquad \dot{\mathbf{q}} = B\,\dot{\mathcal{Q}}.

(10)

Let's calculate the Hessian for the two ions systems. We start with Eq. 1 the total potential of the system is

V_{\mathrm{2ions}} = \frac{1}{2} m \omega^2 x_1^2 + \frac{1}{2} m \omega^2 x_2^2 + \frac{Z^2 e^2}{8\pi\epsilon_0} \frac{1}{\lvert x_{1} - x_{2}\rvert} + \frac{Z^2 e^2}{8\pi\epsilon_0} \frac{1}{\lvert x_{2} - x_{1}\rvert}.

(1.b)

Before calculating the Hessian matrix let's compute a few derivatives first

\frac{\partial V_{\mathrm{2ions}}}{\partial x_1} = m\omega^2 x_1 - 2\frac{Z^2 e^2}{8\pi\epsilon_0} \frac{\left(x_1 - x_2\right)}{\lvert x_1 - x_2\rvert^3} \\[10pt] \frac{\partial V_{\mathrm{2ions}}}{\partial x_2} = m\omega^2 x_2 + 2\frac{Z^2 e^2}{8\pi\epsilon_0} \frac{\left(x_1 - x_2\right)}{\lvert x_1 - x_2\rvert^3}\\[10pt] \frac{\partial V_{\mathrm{2ions}}}{\partial x_1\partial x_2} = \frac{\partial V_{\mathrm{2ions}}}{\partial x_2\partial x_1} = - 4\frac{Z^2 e^2}{8\pi\epsilon_0} \frac{\left(x_1 - x_2\right)^2}{\lvert x_1 - x_2\rvert^5}\\[10pt] \frac{\partial V_{\mathrm{2ions}}}{\partial x_2\partial x_2} = \frac{\partial V_{\mathrm{2ions}}}{\partial x_1\partial x_1} = m\omega^2 + 4\frac{Z^2 e^2}{8\pi\epsilon_0} \frac{\left(x_1 - x_2\right)^2}{\lvert x_1 - x_2\rvert^5}\\[10pt]

(2.b)

From the previous section we know that $x_1 - x_2 = d$ with $d = \left(Z e^2 / \left(2 \pi \epsilon_0 m \omega^2\right)\right)^{1/3}$ . We can then compute the Hessian matrix at the equilibrium positions

\begin{bmatrix} \frac{\partial V_{\mathrm{2ions}}}{\partial x_1\partial x_1} & \frac{\partial V_{\mathrm{2ions}}}{\partial x_1\partial x_2}\\[10pt] \frac{\partial V_{\mathrm{2ions}}}{\partial x_2\partial x_1} & \frac{\partial V_{\mathrm{2ions}}}{\partial x_2\partial x_2} \end{bmatrix} = \begin{bmatrix} 2 m \omega^2 & - m \omega^2\\[10pt] - m \omega^2 & 2 m \omega^2 \end{bmatrix} = m\omega^2\begin{bmatrix} 2 & - 1\\[10pt] - 1 & 2 \end{bmatrix}

(3.b)

This Hessian matrix can be easily diagonalized obtaining the following eigenvectors:

\mathbf{b}_1 = \left[1, 1\right] \\[10pt] \mathbf{b}_1 = \left[-1, 1\right].

(3.b)

with corresponding eigenvalues $\nu_1^2 = 1$ and $\nu_2^2 = 3$ . The Lagrangian then reduces to:

\mathcal{L} = \frac{m}{2}\left( \dot{\mathcal{Q}}_1 + \nu_1^2 \mathcal{Q}_1 + \dot{\mathcal{Q}}_2 + \nu_2^2 \mathcal{Q}_2 \right),

(4.b)

which is identical to the Lagrangian of two decoupled harmonic oscillators (modes of vibration), one with frequency $\nu_1 = 1$ and $\nu_2 = \sqrt{3}$ .

The first mode $\mathcal{Q}_1$ represents all the ions oscillating in unison, moving back and forth as if rigidly connected (in phase). This mode is known as the center-of-mass mode.

The second mode $\mathcal{Q}_1$ involves the ions oscillating in opposite directions (in antiphase) and is referred to as the breathing mode.

Now we can diagonalise it and obtain its eigenspectra:

Eigenvalue Chart

Once diagonalized the Hessian matrix $\mathcal{A}_{i,j}$ , and changed our coordinates $\mathbf{b} \rightarrow \mathcal{Q}$ , the Lagrangian (Eq. 6) of the system reduces to:

\mathcal{L} = \frac{m}{2}\sum\limits_{m=1}^{3N} \left( \dot{\mathcal{Q}}_m^{2} - \nu_m^2 \mathcal{Q}_m^{2}\right),

(11)

which is analogous to the Lagrangian of $3N$ decoupled harmonic oscillators, each with a frequency $\nu_m$ .

If the particles have different masses $m_i$ , the kinetic term is

T = \tfrac{1}{2}\,\dot{\mathbf{q}}^{\top} M \dot{\mathbf{q}}, \qquad M = \operatorname{diag}(m_1, m_2, \ldots, m_{3N}),

(1.c)

with the Lagrangian,

L = \tfrac{1}{2}\,\dot{\mathbf{q}}^{\top} M \dot{\mathbf{q}} \;-\; \tfrac{1}{2}\,\mathbf{q}^{\top} A \mathbf{q}.

(2.c)

Let's define the mass-weighted coordinates,

\tilde{\mathbf{q}} = M^{1/2}\mathbf{q}, \qquad \dot{\tilde{\mathbf{q}}} = M^{1/2}\dot{\mathbf{q}}.

(3.c)

then,

L = \tfrac{1}{2}\,\dot{\tilde{\mathbf{q}}}^{\top}\dot{\tilde{\mathbf{q}}} \;-\; \tfrac{1}{2}\,\tilde{\mathbf{q}}^{\top}\big(M^{-1/2} A M^{-1/2}\big)\tilde{\mathbf{q}}.

(4.c)

Since $\tilde{A} = M^{-1/2} A M^{-1/2}$ is symmetric, there exists an orthogonal $\tilde{B}$ with

\tilde{B}^{\top} \tilde{A} \tilde{B} = \operatorname{diag}(\nu_1^2, \ldots, \nu_{3N}^2), \qquad \tilde{B}=\left[\tilde{\mathbf{b}}^{(1)}, \tilde{\mathbf{b}}^{(2)}, \cdots \right].

(5.c)

In the original coordinate system ( $\mathbf{q}$ ) the eigenvectors $\mathbf{b}^{(m)}$ are given by $\mathbf{b}^{(m)}=M^{-1/2}\tilde{\mathbf{b}}^{(m)}$ which satisfy the general eigenvalue problem $A\mathbf{b}^{(m)}=\nu_m^2 M \mathbf{b}^{(m)}$ and are orthonormal under the mass-weighted inner product $M\,\mathbf{b}^{(m)}\cdot \mathbf{b}^{(n)} = \delta_{mn}$ .

Defining $\tilde{\mathbf{q}} = \tilde{B}\,\mathcal{Q}$ gives,

L = \tfrac{1}{2}\sum_{m=1}^{3N}\left( \dot{\mathcal{Q}}_m^{\,2} - \nu_m^{2}\,\mathcal{Q}_m^{\,2} \right),

(6.c)

and the relation to physical displacements is,

\mathbf{q} = M^{-1/2}\,\tilde{B}\,\mathcal{Q}.

(7.c)

Moving back to our original coordinate system the physical displacement is,

\mathbf{q}(t) = \sum_{m=1}^{3N} \mathbf{b}^{(m)}\,\mathcal{Q}_m(t).

(12)

In this basis each harmonic oscillator is associated to a common mode of vibration of the ions. For each eigenvalue we can easily visualize its associated common mode of vibration:

Loading eigenmodes simulation…

Conclusions and Acknowledgements

This post aims to clarify how an ion crystal composed of $N$ ions can be modeled as $3N$ decoupled harmonic oscillators—one for each spatial dimension.

By expanding the potential around the equilibrium configuration and calculating its Hessian matrix, we derive the frequencies and vibrational modes of the system. This forms the foundation for generating controlled interactions between ions.

Consider two ions of the crystal, ion $i$ and $j$ . When exciting a particular vibrational mode the relative distance between them changes, as shown in the last animation. As the ion distances change, the Coulomb interaction modifies the system's potential energy.

By selectively exciting the correct modes at the appropriate frequencies using a laser beam, it becomes possible to engineer effective spin-spin interactions, a crucial component for trapped-ion quantum simulators.

In a future post, I’ll go into more detail about how to choose the "correct" modes and frequencies.

Finally, I’d like to thank Liam for providing a helpful code for computing Hessian matrices and Rene, whose Mathematica notebook and work inspired this post.

Footnotes

Quantum dynamics of cold trapped ions with application to quantum computation ↩

Tags

Ions' equilibrium positions

Eigenvalue Chart

Conclusions and Acknowledgements

Footnotes