1. Introduction
Mathieu’s equation, which holds a prominent place in mathematical physics, has a fascinating historical origin dating back to the 19th century. Emile Mathieu conducted research on the vibrations of elliptical drums. His work was motivated by the broader scientific context of understanding the wave phenomena, particularly in the study of acoustics and vibrations. In 1868, Mathieu published his seminal work and his investigations are rooted in the Helmholtz equation [1, 2]
$\begin{eqnarray}{{\rm{\nabla }}}^{2}W+{k}^{2}W=0,\end{eqnarray}$
in which the solutions are defined as W(u, v, z) = f(u)g(v)φ(z). In the elliptic cylinder coordinates $x=\rho \cosh (u)\cos (v)$, $y=\rho \sinh (u)\sin (v)$, and z = z. In the above equation, the curve by u = constant gives the confocal ellipses and the curve by v = constant gives the orthogonal hyperbolas. If we define Y = f and t = u, or Y = g and t = v, then we will obtain the same ordinary differential equation, now known as the Mathieu equation, which emerged as a solution to describe the complex oscillatory patterns observed in elliptical hoops [3]. For this reason, the Mathieu and the modified Mathieu equation can emerge from any differential equations involving the Helmholtz equations expressed in the elliptic cylinder coordinates. It has been discussed in details by Paul from 1953 to 1958; see [4] and references therein. For a long time, the Mathieu equation has been one of the most important linear differential equations studied in mathematics and physics [5, 6], and can be written as the following second-order linear ordinary differential equation [4, 7] $\begin{eqnarray}\frac{{{\rm{d}}}^{2}Y}{{\rm{d}}{t}^{2}}+(a-2q\cos (2t))Y=0.\end{eqnarray}$
Here t is the time, Y is the vibration amplitude, and the coefficients a and q, known as the Mathieu parameters, completely determine the stability of the dynamics. Obviously, the differential operator itself is a periodic function, with period T = π, however the solution is not necessarily periodic, but instead, a quasi-periodic solution is allowed. This kind of quasi-periodicity is a typical feature of differential equations in dynamical systems [8].The above Mathieu equation is a special condition of the Hill equation, which can be written as
$\begin{eqnarray}\frac{{{\rm{d}}}^{2}Y}{{\rm{d}}{t}^{2}}+\left(a+2\displaystyle \sum _{j=1}^{\infty }{q}_{j}\cos (2jt)\right)Y=0,\end{eqnarray}$
where qj are constants. The solution is related to a determinant of an infinite matrix. When only q1 and q2 are involved and qj = 0 for all j > 2, it will be reduced to the Whittaker–Hill equation [9–11].In [1], a lot of important applications have been summarized, including elliptic drums, inverted pendulum, radio frequency quadrupole, modulated LC circuit, floating body, etc. In recent years, the Paul trap for charged particles and mirror trap for neutral particles has gained some major attention. In [12–15], applications of the Mathieu equation in wave propagation in pips, electromagnetic wave guides and oscillations of water in a lake have also been presented. In these applications, the physics in the Paul trap is the major concern in this work for its potential applications in quantum computation, following the scheme by Cirac and Zoller [16].
The Mathieu model provides a mathematical framework to understand the behavior of charged particles in various types of traps. Trapped ions are confined spatially by electric or magnetic fields, and their motion within these traps is described by the Mathieu equation [4, 17]. In this case, the dynamics are fully determined by several parameters: a, q and the damping rate γ (see below). Moreover, the trap geometry and field strength are also determined from the values of these parameters [4, 18, 19]. Motivated by this research, in this work, we extended the single particle Mathieu equation to the many-body Mathieu equation, and study its dynamics based on classical equations. This work is organized in the following way. The dynamics in the quantum regime need to be considered elsewhere. We first start by presenting the criterion for stability and instability based on Floquet theorem in section 2 . Then, in section 3 , we study the dynamics with two to four trapped ions, in which in the large damping limit, the dynamics can be solved analytically. In this case, we show how the damping rate influences the relative phase and the vibration amplitude around their equilibrium positions. In section 4 , we study the dynamics with N trapped ions with a finite damping rate. In section 5 , we examine the stability and the phase chart in the many-body Mathieu equation, which has been argued in the previous section. We show that the phase chart is independent of the many-body interaction. In section 6 , the fate of the phase chart and the stability is examined in the presence of anharmonic nonlinear potential, showing that the parametric resonance is absent in the presence of a cubic nonlinear coefficient. Finally, we discuss the experimental relevance of our results in section 7 . We conclude this work and discuss the unique feature of this many-body Mathieu model and its potential applications in section 8 .
2. Stability of the Mathieu equation and the phase chart
The Mathieu equation is the simplest equation to study the stability and instability of a system with a periodic driving force, which is known as parametric resonance. This behavior is totally different from the forced resonance in the linear differential equation. In the forced resonance, the resonance can only happen when the driving frequency is the same as the intrinsic vibrational frequency of the system, in which the vibrational amplitude increases linearly with time t. However, in the parametric resonance, the resonance can be found in a wide range of parameters, with a vibrational amplitude increasing exponentially with time t [5, 6, 20]. In the following, we discuss the phase chart of the Mathieu equation with and without damping [21].
We can use the Floquet theory to understand the stability of this model. Let us consider the following equation with periodic coefficients [22–24]
$\begin{eqnarray}\frac{{{\rm{d}}}^{2}Y}{{\rm{d}}{t}^{2}}+2\gamma \frac{{\rm{d}}Y}{{\rm{d}}t}+(a-2q\cos (2t))Y=0,\end{eqnarray}$
where γ is the damping rate [25, 26]. This equation can be defined as [27] $\begin{eqnarray}\left[\begin{array}{c}\dot{Y}\\ \ddot{Y}\\ \end{array}\right]=\left[\begin{array}{cc}0 & 1\\ -(a-2q\cos (2t)) & -2\gamma \\ \end{array}\right]\left[\begin{array}{c}Y\\ \dot{Y}\\ \end{array}\right]=A\left[\begin{array}{c}Y\\ \dot{Y}\\ \end{array}\right].\end{eqnarray}$
We will obtain two independent solutions, that can be denoted as ${({Y}_{1}(t),{\dot{Y}}_{1}(t))}^{{\rm{T}}}$, and ${({Y}_{2}(t),{\dot{Y}}_{2}(t))}^{{\rm{T}}}$, with the initial conditions $\begin{eqnarray}\left[\begin{array}{c}{Y}_{1}(0)\\ {\dot{Y}}_{1}(0)\end{array}\right]=\left[\begin{array}{c}1\\ 0\end{array}\right],\quad \left[\begin{array}{c}{Y}_{2}(0)\\ {\dot{Y}}_{2}(0)\end{array}\right]=\left[\begin{array}{c}0\\ 1\end{array}\right].\end{eqnarray}$
Furthermore, we can construct the solutions as $\begin{eqnarray}{\rm{\Psi }}=\left[\begin{array}{cc}{Y}_{1}(t) & {Y}_{2}(t)\\ {\dot{Y}}_{1}(t) & {\dot{Y}}_{2}(t)\\ \end{array}\right],\end{eqnarray}$
with $\dot{{\rm{\Psi }}}=A{\rm{\Psi }}$, which is based on the Floquet theory, its solution should be $\Psi$(t + T) = $\Psi$(t)C. One can check this solution $\dot{{\rm{\Psi }}}(t+T)=A{\rm{\Psi }}(t+T)=A{\rm{\Psi }}(t)C=\dot{{\rm{\Psi }}}(t)C$ [28]. Using the initial condition that $\Psi$(0) = 1, we can immediately obtain $\begin{eqnarray}C=\left[\begin{array}{cc}{Y}_{1}(T) & {Y}_{2}(T)\\ {\dot{Y}}_{1}(T) & {\dot{Y}}_{2}(T)\\ \end{array}\right].\end{eqnarray}$
In this way, we expect $\Psi$(nT) = $\Psi$(0)Cn. Thus, the system is stable when and only when the eigenvalues of C, which are assumed by β, satisfy ∣β∣ ≤ 1. We can obtain the eigenvalues of the C matrix using $\begin{eqnarray}{\beta }^{2}-{\rm{Tr}}(C)\beta +{\rm{\det }}(C)=0.\end{eqnarray}$
The above equation is used to determine the phase boundary. In the above solution, we have ${Y}_{i}\in {\mathbb{R}}$ and ${\dot{Y}}_{i}\in {\mathbb{R}}$, thus C is a real matrix; yet non-Hermite. (I) When the solution is real, we should have a phase boundary at $1\pm {\rm{Tr}}(C)+{\rm{\det }}(C)=0$; (II) When β is a complex, the two solutions are complex valued, and we expect ${\rm{\det }}(C)=1$. The second condition will corresponds to the condition without damping γ = 0 [29, 30].The phase boundary is fully determined by the two parameters ${\rm{Tr}}(C)$ and ${\rm{\det }}(C)$. We find that the phase boundaries are determined by the regime enclosed by the following condition
$\begin{eqnarray}\begin{array}{r}1\,\geqslant \,\,\rm{det}\,{(C)}^{2}\,\geqslant \,| 1\pm {\rm{Tr}}(C){| }^{2}.\end{array}\end{eqnarray}$
With this theory, we can determine the phase boundary as a function of a and q. In figure 1, we present the phase chart without damping [31, 32], and in figure 2, we present the results with damping [33, 34]. The roots of the equation without damping are given by $\begin{eqnarray}\beta =\frac{{\rm{Tr}}(C)\pm \sqrt{{\rm{Tr}}{(C)}^{2}-4}}{2},\end{eqnarray}$
and the stable phase is given by $| {\rm{Tr}}(C)| \leqslant 2$ [35, 36]. In this condition, β is a complex number, and we can prove exactly that ∣β∣ = 1, in this condition we can define $\beta =\exp ({\rm{i}}\theta )$.
Figure 1. Stability chart with stable and unstable regions (tongues). The shadowed regimes correspond to the condition ∣β∣ > 1 in equation ( |
Figure 2. Comparative analysis of stable and unstable regions of the damped and undamped Mathieu equations. |
We see that damping has the effect of reducing the unstable regions of the system, which is consistent with our intuition; see [22, 37] and [5]. Whenever there is a dissipation of energy, the system goes to rest, which naturally means that it becomes stable. However, in the following section, we will present a different role played by the dissipative term on the relation between the vibration amplitude and the equilibrium position in the strong damping limit.
3. Dimensionless and large damping limits
Now, we turn to the major results that will be presented in this work. Our aim is to study the dynamics of trapped ions in the presence of a strong Coulomb interaction. The Hamiltonian for N trapped ions can be written as
$\begin{eqnarray}H={H}_{0}+U,\end{eqnarray}$
where $\begin{eqnarray}{H}_{0}=\displaystyle \sum _{i=1}^{N}\frac{{p}_{i}^{2}}{2m}+\frac{1}{2}m{\omega }^{2}(a-2q\cos ({\rm{\Omega }}\tau )){q}_{i},\end{eqnarray}$
$\begin{eqnarray}U=\frac{{e}^{2}}{4\pi {\epsilon }_{0}}\displaystyle \sum _{i\lt j}\frac{1}{| {q}_{i}-{q}_{j}| }.\end{eqnarray}$
In the above model, m is the mass of the trapped ions, ε0 is the vacuum permittivity, Ω is the modulation frequency, and ω is the system’s frequency. The equation of motion of this model can then be written as $\begin{eqnarray}\begin{array}{rcl}m{\ddot{q}}_{i} & = & -2\alpha {\dot{q}}_{i}-m{\omega }^{2}(a-2q\cos ({\rm{\Omega }}\tau )){q}_{i}\\ & & -\frac{2{e}^{2}}{4\pi {\epsilon }_{0}}\displaystyle \sum _{j\ne i}\frac{1}{{({q}_{i}-{q}_{j})}^{2}}.\end{array}\end{eqnarray}$
Here, α is introduced for the damping effect.Now, we use the following transformation to make the equation dimensionless
$\begin{eqnarray}\begin{array}{rcl}\tau & = & \frac{2t}{{\rm{\Omega }}},\quad {q}_{i}={\omega }^{2}\sqrt{{ \mathcal Q }}{y}_{i},\quad \gamma =\frac{\alpha }{m\omega },\\ \omega & = & \frac{{\rm{\Omega }}}{2},\quad { \mathcal Q }=\frac{2{e}^{2}}{4\pi {\epsilon }_{0}m{\omega }^{3}}.\end{array}\end{eqnarray}$
Substituting the new variables and parameters into the original equation results in the following interacting Mathieu equation $\begin{eqnarray}{\ddot{y}}_{i}=-\hat{L}({y}_{i})-\displaystyle \sum _{j\ne i}\frac{1}{{({y}_{i}-{y}_{j})}^{2}}.\end{eqnarray}$
Hereafter, we have defined a differential operator as $\begin{eqnarray}\hat{L}(y)=2\gamma \dot{y}+(a-2q\cos (2t))y.\end{eqnarray}$
Now all the quantities became dimensionless, and their values are of the order of unity. In the last term, the minus sign represents the repulsive interaction between the charged ions. Obviously, when 2q > a and γ is very small, the trapping potential may become unbounded from below, which may lead to instability.3.1. The case of two ions N = 2
The dynamics of a system of two ions can be effectively analyzed through the decomposition of the system into two components: the motion of the center of mass and the motion of the relative position of the ions. This decomposition simplifies the analysis by reducing the problem from a complex two-body interaction to two independent single-body equations. To this end, we define 2R = y1 + y2, and r = y1 − y2. The system yields two equations
$\begin{eqnarray}\ddot{R}\pm \frac{\ddot{r}}{2}+\hat{L}\Space{0ex}{2.5ex}{0ex}(R\pm \frac{r}{2}\Space{0ex}{2.5ex}{0ex})-\frac{1}{{r}^{2}}=0.\end{eqnarray}$
These two equations lead to two independent equations as following $\begin{eqnarray}\ddot{R}+\hat{L}(R)=0,\quad \ddot{r}+\hat{L}(r)-\frac{2}{{r}^{2}}=0.\end{eqnarray}$
Here the first equation (the motion of the center of mass) corresponds to the damped Mathieu equation, and the second one (the motion of the relative position) corresponds to the same damped Mathieu equation but with an extra Coulomb interaction term. In the instability phase when r → ∞, the second equation will be reduced to the linear Mathieu equation. For this reason, these two equations have the same phase boundary for stability-instability transitions. However, in the stable phase, their dynamics will be totally different. We find that R = 0 is the solution of the center of mass, but r = constant is not the solution of the relative position. For this reason, even in the presence of damping, the motion of the relative position will always oscillate in some way. Moreover, due to the non-penetration effect of the trapped ions, the motion of r will be restricted to r < 0 or r > 0, depending on its initial value. The results for various damping rates are presented in figure 3, showing that in the stable phase when R → 0, r will always oscillate periodically around its equilibrium position with the same period of the Mathieu equation.
Figure 3. Evolution of a system of two ions for a fixed pair of a = 1, q = 0.3. (a)–(d) represent the evolution with γ = 0, 0.1, 0.2 and 0.7, respectively. Since r(0) = 1, during the evolution, r(t) > 0 from the non-penetration effect of the charged ions. (a) and (b) are the dynamics in the unstable phase; and (c) and (d) are the results in the stable phase due to the damping. In the large γ limit, the solution of r can be described by equation ( |
In the case of strong damping (see figure 3(d)), we may assume r = re + δr, where re is the equilibrium position, one may find that δr oscillates in the way of the driving dynamics. When γ → ∞, we expect $a{r}_{e}-2/{r}_{e}^{2}=0$, yielding re = (2/a)1/3 = (2)1/3 for a = 1 used in figure 3. Thus, we assume that the dynamics are the following
$\begin{eqnarray}r={\left(\frac{2}{a}\right)}^{1/3}+A\sin (2t+\theta )+\cdots \,{,}\end{eqnarray}$
where θ is a phase introduced by the damping effect. Inserting this solution into the above equation, and to the leading term, we get $\begin{eqnarray}4{\left(\frac{a}{2}\right)}^{1/3}A\gamma -2q\cos \left(\theta \right)=0,\end{eqnarray}$
$\begin{eqnarray}\frac{3a-4}{{\left(2/a\right)}^{1/3}}A+2q\sin \left(\theta \right)=0.\end{eqnarray}$
These two terms correspond to the coefficients of $\cos (2t)$ and $\sin (2t)$, which should be equal to zero. The higher-order terms such as $\cos (4t)$ and $\sin (4t)$ are neglected. The solutions of A and θ can be written as $\begin{eqnarray}A=\frac{q}{2\gamma {(a/2)}^{1/3}}\cos (\theta ),\end{eqnarray}$
and $\begin{eqnarray}\frac{(3a-4)q}{2\gamma }\cos (\theta )+2q\sin (\theta )=0.\end{eqnarray}$
In the large γ limit, we expect θ ∼ 0, or $\sin (\theta )=(4-3a)/4\gamma $. In this limit, we expect the solution to be $\begin{eqnarray}r={\left(\frac{2}{a}\right)}^{1/3}+\frac{q}{{2}^{2/3}{a}^{1/3}\gamma }\sin (2t).\end{eqnarray}$
This solution can be understood by assuming r = (2/a)1/3 + y, where for small y and ∣(2/a)1/3∣ ≫ ∣y∣, we have 26 ). For this reason, the effect of damping influences the vibration amplitude, instead of driving the ions to rest. This picture applies to the conditions of more trapped ions.
$\begin{eqnarray}-\frac{2q\cos (2t)}{{(a/2)}^{1/3}}+(ay+\ddot{y}+2\gamma \dot{y}-2qy\cos (2t))=0.\end{eqnarray}$
For the four terms between the brackets, the term $2\gamma \dot{y}$ dominates and we have $\begin{eqnarray}-\frac{2q\cos (2t)}{{(a/2)}^{1/3}}+2\gamma \dot{y}=0,\end{eqnarray}$
yielding the following solution $\begin{eqnarray}y=\frac{q{r}_{e}}{2\gamma }\sin (2t)=\frac{q}{2{(a/2)}^{1/3}\gamma }\sin (2t),\end{eqnarray}$
in equation (The above results can also be understood from the balance between the energy Edriving introduced by the driving force and the energy Edamped dissipated by the damping force [38]. We have 29 ))
$\begin{eqnarray}{E}_{\,\rm{driving}\,}={\int }_{0}^{\pi }-2q\cos (2t)x\dot{x}{\rm{d}}t=-2\pi Aq{r}_{e},\end{eqnarray}$
$\begin{eqnarray}{E}_{\,\rm{damped}\,}={\int }_{0}^{\pi }2\gamma {\dot{x}}^{2}{\rm{d}}t=4\pi {A}^{2}\gamma .\end{eqnarray}$
A direct calculation using Edriving + Edamped = 0 will yield the above relation, with (see equation ( $\begin{eqnarray}A=\frac{q{r}_{e}}{2\gamma }.\end{eqnarray}$
This method has been used in [5] in the determination of the vibration amplitude of the damped oscillator, and the same idea can be used to understand the vibration amplitude in our model. As a result, when Edriving + Edamped > 0 the vibration amplitude will increase due to the absorption of the energy from the driving field; otherwise, it will decrease. This picture can be used to understand the quasi-periodic motion of the trapped ions, as discussed below. Thus, the quasi-periodic solution can also be found in nonlinear driven equations.3.2. The case of three ions N = 3
The same method can be applied to three trapped ions, by defining r1 = y1 − R, r2 = y3 − R, and 3R = y1 + y2 + y3. This results in a system of three coupled equations
$\begin{eqnarray}\begin{array}{l}\ddot{R}-{\ddot{r}}_{1}-{\ddot{r}}_{2}+\hat{L}(R-{r}_{1}-{r}_{2})\\ \quad -\frac{1}{{(-{r}_{1}-2{r}_{2})}^{2}}+\frac{1}{{(2{r}_{1}+{r}_{2})}^{2}}=0,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\ddot{R}+{\ddot{r}}_{1}+\hat{L}(R+{r}_{1})\\ \quad -\frac{1}{{({r}_{1}-{r}_{2})}^{2}}-\frac{1}{{(2{r}_{1}+{r}_{2})}^{2}}=0,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\ddot{R}+{\ddot{r}}_{2}+\hat{L}(R+{r}_{2})\\ \quad +\frac{1}{{(-{r}_{1}-2{r}_{2})}^{2}}+\frac{1}{{({r}_{1}+{r}_{2})}^{2}}=0.\end{array}\end{eqnarray}$
With this definition, we can obtain the equation of the center of mass as $\ddot{R}+\hat{L}(R)=0$, which is exactly the same as the driving Mathieu equation. This conclusion can be extended to arbitrary N ≥ 2 ions. The equations of motion for r1 and r2 can be written as $\begin{eqnarray}{\ddot{r}}_{1}+\hat{L}({r}_{1})=\frac{1}{{({r}_{1}-{r}_{2})}^{2}}+\frac{1}{{(2{r}_{1}+{r}_{2})}^{2}},\end{eqnarray}$
$\begin{eqnarray}{\ddot{r}}_{2}+\hat{L}({r}_{2})=-\frac{1}{{({r}_{1}-{r}_{2})}^{2}}-\frac{1}{{({r}_{1}+2{r}_{2})}^{2}}.\end{eqnarray}$
We first start by defining the equilibrium position equations for a strong damping r1 = re + δr, and r2 = − re − δr. In the large γ limit, we neglect δr, and would obtain
$\begin{eqnarray}a{r}_{e}=\frac{1}{4{r}_{e}^{2}}+\frac{1}{{r}_{e}^{2}}\to {r}_{e}=\frac{{5}^{1/3}}{{2}^{2/3}{a}^{1/3}}.\end{eqnarray}$
Therefore, the solutions would be $\begin{eqnarray}\begin{array}{r}{r}_{1}=-{r}_{2}=\frac{{5}^{1/3}}{{2}^{2/3}{a}^{1/3}}+A\sin (2t+\theta ).\end{array}\end{eqnarray}$
Using the perturbation theory, and assuming A to be small, we obtain from its equation of motion that $\begin{eqnarray}4A\gamma -2q\frac{{5}^{1/3}}{{2}^{2/3}{a}^{1/3}}\cos (\theta )=0,\end{eqnarray}$
$\begin{eqnarray}\frac{(5a-4)q\cos (\theta )}{2\gamma }+2q\sin (\theta )=0.\end{eqnarray}$
The solutions of A and $\sin (\theta )$ can be written as 32 ). This solution has been verified using numerical data. In figure 4, we present the results for three trapped ions, showing that in the large γ limit, indeed, when R → 0, r1 and r2 will vibrate around their equilibrium positions, which can be well described by the above analytical solution 43 ).
$\begin{eqnarray}\begin{array}{r}A=\frac{{5}^{1/3}q\cos (\theta )}{2{a}^{1/3}{2}^{2/3}\gamma },\quad \sin (\theta )=\frac{(4-5a)}{4\gamma }\cos (\theta ).\end{array}\end{eqnarray}$
In the large γ limit, θ ∼ 0, expecting the solution to be $\begin{eqnarray}r=\frac{{5}^{1/3}}{{2}^{2/3}{a}^{1/3}}+\frac{{5}^{1/3}q}{2{a}^{1/3}{2}^{2/3}\gamma }\sin (2t).\end{eqnarray}$
Obviously, the vibration amplitude in this solution also satisfies the energy balance condition of equation ( $\begin{eqnarray}\begin{array}{r}2\gamma \dot{y}-\frac{{5}^{1/3}2q}{{2}^{2/3}{a}^{1/3}}\cos (2t)=0,\end{array}\end{eqnarray}$
which yields $y=\frac{{5}^{1/3}q}{{2}^{2/3}2{a}^{1/3}\gamma }\sin (2t)$; see equation (
Figure 4. Evolution of a system of 3 ions for a fixed pair of a = 1, q = 0.3. (a)–(d) represent the evolution with γ = 0, 0.1, 0.2 and 0.7, respectively. |
3.3. The case of four ions N = 4
This method can be extended to any arbitrary number of trapped ions, and here we will only present results for N = 4. In a similar way, for i = 1, 4 and j = 2, 3 we can define 32 ). From these results, we expect that in the condition of N ≥ 4, we should have ${r}_{e}^{i}\propto 1/{a}^{1/3}$ and Ai ∝ q/a1/3γ, up to some universal constants to be determined numerically.
$\begin{eqnarray}\begin{array}{r}{y}_{i}=\pm {r}_{e}^{1}\pm {A}_{1}\sin (2t),\quad {y}_{j}=\pm {r}_{e}^{2}\pm {A}_{2}\sin (2t).\end{array}\end{eqnarray}$
This result is based on numerical simulation, showing that in the large γ limit, trapped ions 1 and 4 (and 2 and 3) can have almost opposite dynamics. From the equation of motion and assuming Ai to be small in comparison to ${r}_{e}^{i}$, we obtain $\begin{eqnarray}{r}_{e}^{1}=\frac{1.4368}{{a}^{1/3}},\quad {r}_{e}^{2}=\frac{0.454379}{{a}^{1/3}},\end{eqnarray}$
which can not be solved analytically. We also find using numerical fitting that $\begin{eqnarray}{A}_{1}=1.4368\frac{q}{2\gamma {a}^{1/3}},\quad {A}_{2}=0.454379\frac{q}{2\gamma {a}^{1/3}},\end{eqnarray}$
which can be understood using equation (So far, we focused on studying systems in the strong damping limit, and what we noticed is that it have barely any influence on the total energy of the ions. This means that despite the presence of damping from the environment or the laser field, the total energy of the system remains relatively conserved. The damping does not extract a significant large energy from the system, that is why it remains stable over time; see discussion in the next section. As a result, in this regime, the strong damping only influences the oscillation amplitudes of the ions around their equilibrium positions. In the following, we will only focus on the dynamics with finite damping, which will yield some new interesting features.
3.4. Anharmonic effect
The above analysis—based on the center of mass and the relative coordinates—depends strongly on the linearity of the equation of motion. In the presence of anharmonic effects in the trapping potential, the calculation is much more complicated. It has some profound impact on the dynamics of the system. Firstly, the method of the center of mass and the relative positions can not be applied anymore; and secondly, the phase chart and stability can be changed even in the nonlinear Mathieu equation; and thirdly, the criterion for stability based on the Floquet theorem will no longer be valid. For these reasons, the quasi-periodic solution may not be observed.
To this end, let us consider the following model with N = 3 with strong damping. The equation of motion can be written as 38 ) when η = 0, and will be solved using numerical methods. Using perturbation theory again, we will find the leading terms of the equations (with A1 = − A3) 39 ), with a finite relative phase. These solutions are the same as equation (32 ). The equilibrium position re can be numerically determined using equation (49 ). Obviously, with the increasing of η (assuming η > 0), re will decrease monotonically. In the large η limit, we should have re ∼ (5/4η)1/5 ≃ 1/η1/5.
$\begin{eqnarray}{\ddot{y}}_{i}+\hat{L}({y}_{i})+\eta {y}_{i}^{3}-\displaystyle \sum _{j\ne i}\frac{1}{{({y}_{i}-{y}_{j})}^{2}}=0,\quad i=1-3.\end{eqnarray}$
Let us assume that ${y}_{1,3}=\pm {r}_{e}\pm {A}_{1,3}\sin (2t)$, and y2 = 0, where by symmetry the equilibrium position of the first and third ions should be at ± re and the equilibrium position of the second ion should be at zero. This solution is supported by numerical evidences in figure 5. In the case of strong damping γ, ∣A1,3∣ ≪ ∣re∣, thus we have $\begin{eqnarray}a{r}_{e}+\eta {r}_{e}^{3}=\frac{5}{4{r}_{e}^{2}}.\end{eqnarray}$
This equation is the same as equation ( $\begin{eqnarray}\begin{array}{r}4\gamma {A}_{1}-2q{r}_{e}=0,\end{array}\end{eqnarray}$
for the coefficient of $\cos (2t)$, allowing us to obtain the expressions of $\begin{eqnarray}\begin{array}{r}{A}_{1}=-\frac{q{r}_{e}}{2\gamma },\quad {A}_{3}=\frac{q{r}_{e}}{2\gamma }.\end{array}\end{eqnarray}$
The coefficient of $\sin (2t)$ is neglected. A much better approximation of the solutions is using equation (
Figure 5. Evolution of a system of three ions for a = 1, q = 0.3, η = 0.01. (a)-(d) represent the evolution with γ = 0, 0.1, 0.5 and 1, respectively. |
Figure 6. The dynamics evolution of a system of two ions with asymmetric initial conditions, with a = 5, q = 1.7 and without damping. (a) and (b) depict a detailed view of the ions’ positions and velocities over a duration of t = 100. (c) and (d) represent the evolution of the force and the energy over a period of time of t = 100. Meanwhile, (e), (f), (g) and (h) depict the positions, the velocities, the forces and the energy of the ions over a shorter period of time t = [20 − 30], respectively. In (g), the force diverges when the distance between the ions approach zero simultaneously from the synchronization dynamics. |
We see that in the condition of anharmonic effects in the trapping potential, we still have the relation in equation (32 ) for the balance between the driving force and the damping force averaged in one full period. Therefore, the anharmonic effects, η, in the trapping potential only influences the vibration amplitude in a way of equation (49 ). A simple and straightforward argument will show that this conclusion to be correct even with many trapped ions.
4. Exchange of velocities and Newton’s cradle
4.1. Without damping
Then, we ask what will happen in the presence of finite damping? Generally speaking, the dynamics and vibration amplitudes depend on the competition of the driving forces, which may introduce some energy to the system, and the damping force, which can extract energy from the system. In the large γ limit, these two forces are balanced, leading to periodic vibrations. In the long time limit, we find that the dynamics of the trapped ions are independent of their initial positions and velocities (see equation (45 )). In the weak damping limit, the physics will be totally changed. For example, in the model without damping, this system should support a quasi-periodic solution, instead of a periodic solution, and large vibrations of the amplitude of the ions will be possible, in which the distance between the ions may become very small, leading to large forces from the repulsive Coulomb interaction. In this way, the dynamics will exhibit some new features, which resembles to the dynamics of Newton’s cradle; see figure 7. This motivates our following investigations.
Figure 7. Newton’s cradle and the sudden exchange of velocities between distant balls during collisions when all of them have the same mass and geometry. |
We will focus on the following quantities: the positions yi, the velocities ${v}_{i}={\dot{y}}_{i}$, the forces fi = − ∇iU, and the energy E = T + U. Here, E is the total energy of the system defined by the kinetic energy, the trapping energy and the Coulomb energy. When two ions approach each other, the force may diverge, leading to dramatic changes of velocities—exchange of velocities—[39], which is the same as the dynamics in Newton’s cradle. Furthermore, when the total energy E increases, the system absorbs energy from the driving force; otherwise, the driving force will extract energy from the system.
The results for yi, vi, fi and E for two trapped ions are presented in figure 6 without damping. We focus on the stable regime. The dynamics of y1 and y2 exhibit some interesting features not presented in the strong damping limit. Firstly, it presents some kind of micro-motion with period π from the driving field; and secondly, the system shows some global structure with much larger period, not necessary to the commensurate with the period of the driving field. In this way, the system displays some quasi-periodic dynamics. The velocity also exhibits similar dynamics in figure 6(b). We find that the vibration amplitude and the velocity are greatly enhanced, indicating that the system continuously absorbing and releasing energy from the driving field. In figure 6(f), we plot the detailed dynamics of the velocity, showing that when two ions are approaching each other (see figure 6(e)), there is a sudden exchange of velocities, which is the same as that in Newton’s cradle (figure 7). Our model may be the smallest system for the realization of Newton’s cradle. At this point, the force will become divergent or significantly large; see figures 6(c) and (g). In figure 6(d), we calculate the total energy of this model, showing that the total energy has the same quasi-periodic structure for the motion and the velocity. These results represent some unique features of the interacting Mathieu equation with the many-body Coulomb force.
Next, we are in a position to understand what will happen in the presence of many trapped ions. In Newton’s cradle, the exchange of velocities happens between two colliding balls; and the higher number of balls involved, the more complicated the dynamics become. To this end, we present the results for three trapped ions in figure 8 and for five trapped ions in figure 9 without damping. In these two figures, we always find the sudden exchange of velocities when all the ions are colliding with each other. The most striking feature is that in the presence of many trapped ions (N ≥ 3), the motion of the ions will be synchronized in such a way that they will collide into each other almost at the same time. We have verified that similar features can be found for much larger systems (for example, N ≥ 10), and these synchronization dynamics are independent of their initial positions and initial velocities. As a result, for three trapped ions, the exchange of velocities happens between y1 and y3 (assuming y1 < y2 < y3); and for five trapped ions, the exchange of velocities happen between y1 and y5, and y2 and y4. The synchronization dynamics are a unique feature in our model.
Figure 8. The dynamics of a system of three ions with asymmetric initial conditions with a = 5, q = 1.7 and without damping. The meanings of the sub-figures are the same as that in figure 6. |
Figure 9. The dynamics of a system of five ions with asymmetric initial conditions with a = 5, q = 1.7 and without damping. The meanings of the sub-figures are the same as that in figure 6. |
These results should not be mixed up with the melting of ions, in which the structure is totally destroyed [39]. In this case, the distance between the ions will increase, leading to the destruction of the crystal structures. In our model, we indeed found the dramatic increase of the distance between the ions, however, these ions are still trapped in the potential, and their distances will be decreased after some proper time, leading to coherent quasi-periodic dynamics. In this way, the structure of the ions will not be destroyed.
4.2. Effect of damping
With the above dynamics, we need to understand the role of damping in the many-body Mathieu equation, in which we expect a transition from the undamped physics discussed in the above subsection to the strong damping physics discussed in the large γ limit. These physics are intriguing, not only because of the effect of damping on the exchange of velocities, but also from the transition from a quasi-periodic solution to a periodic vibrational solution. Here, the strong damping and weak damping limits are only used to define the two limits with and without quasi-periodic solutions, yet in a true system, this is actually a smooth transition process.
By focusing on the dynamics of a system of five ions, we can much further observe the ion-ion interactions and highlight the unique features of damping. In the stable regime, the results for yi, vi, fi and E for five trapped ions with a weak damping and a strong one are presented in figures 10 and 11, respectively. In the case of weak damping, the system behaves almost exactly the same as without damping. This is seen in figures 10(a) and (e), the ions oscillate in a quasi-periodic way where they approach each other, the Coulomb force becomes strong in figures 10(c) and (g), leading to a sudden exchange of velocities, as observed before. In this condition, each ion will mirror the furthest ion in the system, thus we observe that ion y1 exchanges its velocity with ion y5, and y2 exchanges its velocity with ion y4, leaving ion y3 to be almost unchanged. This is a unique feature from Newton’s cradle, in which the kinetic energy is transferred from the left furthest ball to the right furthest ball through collisions, and everything in between stays stable. This feature is not changed significantly by the weak damping.
Figure 10. The dynamics of a system of five ions with asymmetric initial conditions with a = 5, q = 1.7, and weak damping γ = 0.05. The meanings of the sub-figures are the same as that in figure 6. |
Figure 11. The dynamics evolution of a system of five ions with asymmetric initial conditions with a = 5, q = 1.7, and a strong damping γ = 0.4. The meanings of the sub-figures are the same as that in figure 6. |
However, with the increasing of the damping rate, say γ = 0.4, as shown in figure 11, we notice some total different dynamics. We find that the previous structure of the motion of ions is destroyed; see figures 11(a) and (b). Ions will no longer oscillate quasi-periodically, because the damping will introduce some stability to the system, and each ion now tends to oscillate around its equilibrium position. This will cause the amplitude of each ion to be small in the large γ limit. As a result, we will no longer have any exchange of velocities, as demonstrated in figures 11(b) and (c). Now the dynamics of the system become more stable, and in this case, the distance between the trapped ions is large enough for the Coulomb force to be insignificant. The exchange of velocities will be observed only when the distance between the ions is approaching zero. In the strong γ limit, their dynamics will be reduced to that discussed in section 3 , showing that the anharmonic effect term η only influences the equilibrium position re, with relation between vibrational amplitude and equilibrium position given by equation (32 ). In this work, the validity of this relation is found in many different conditions.
5. Many-body physics stability chart
Finally, we are in a position to understand the role of the Coulomb interaction on the stability and the phase chart in this model using the following argument. Let us assume that the system is unstable. In this case, yi approaches infinity, as observed in the single particle Mathieu equation. When the distance between the ions is sufficiently large, the Coulomb interaction is negligible, thus the instability of this model should be fully determined by the phase chart of the single particle Mathieu model. Thus, we have sufficient reason to believe that the phase chart in this interacting Mathieu model should be the same as the single particle Mathieu equation. This is shown in figure 12, however, the dynamics can still be fundamentally different. In figure 12, we focus on the motion and the velocity in the unstable region for two groups of ions. The first group is composed of five ions in figures 12(a) and (c), and the second group is composed of a single ion in figures 12(b) and (d). Both systems will oscillate exponentially, but have different magnitudes. The difference between the amplitude of the motion of the first ion and the fifth ion is around 102, however, the difference between their velocities is around the order of 103. The vibration magnitude of the position and the velocity are much smaller in a single particle model. Despite of this difference, both systems display similar patterns of oscillations. When a group of ions are sufficiently far from each other, the Coulomb forces decrease rapidly to zero, resulting in no influence on their motion. In such cases, the ions act like independent particles, which gives rise to a collective behavior similar to that of a single particle.
Figure 12. Evolution of the system over a long period of time for a fixed pair of a = 15, q = 9.148 in the unstable region, without damping (γ = 0). (a) and (c) are the results for five ions, and (b) and (d) are the results for one ion, respectively. |
Next, we study the motion and the velocity of five ions and of one ion in the stable region, which are presented in figure 13, with a = 15 and qc ≃ 9.1473. Similar features have also been observed. We find that these two conditions have the same parameters and they have the same stability. However, in the many-body model, the vibration amplitude is greatly enhanced by about two orders of magnitude, as compared with the single particle one. This result also means that in the many-body model the larger N is, the larger the vibration amplitude will be. For this reason, the trapping of N particles (where N is large) without damping will become a significant issue. However, in the presence of damping, the enhanced vibration amplitude will be greatly suppressed.
Figure 13. Evolution of the system over a long period of time for a fixed pair of a = 15, q = 9.146 in the stable region without damping (γ = 0). The sub-figures depict the same information as figure 12. |
This observed phenomenon is fascinating in physics, because of the fact that the collective interactions of multiple particles gives rise to behaviors similar to those of individual particles. This parallelism highlights not only the fundamental principles regarding the behavior of trapped ion systems, but also the interconnections of their dynamics. We hope these physics can be observed in experiments [39]. In the case of strong damping, the single particle will cease to oscillate, yet the many-body one will oscillate around its equilibrium position, with position and amplitude satisfying equation (32 ).
6. Mathieu equation with cubic anharmonic term
The Mathieu equation with a cubic nonlinear term, in the presence of many trapped ions can be written as [40, 41]
$\begin{eqnarray}{\ddot{y}}_{i}+\hat{L}({y}_{i})+\eta {y}_{i}^{3}-\displaystyle \sum _{j\ne i}\frac{1}{{({y}_{i}-{y}_{j})}^{2}}=0,\end{eqnarray}$
where η is the anharmonic coefficient. The cubic nonlinear coefficient introduces a nonlinearity to the equation. Therefore, the response of the system will no longer be directly proportional to the applied force, which naturally will make the system exhibit much more complicated behaviors. The cubic nonlinear coefficient in the Mathieu equation for trapped ions typically originates from the anharmonicity in the trapping potential that is responsible of the interaction between the ion’s motion and the trapping field [42, 43]. Numerical calculations show that in the presence of a weak η, the vibration amplitude is $\begin{eqnarray}A\propto \frac{1}{\sqrt{\eta }},\end{eqnarray}$
thus for all a and q the system is always stable. This relation can be determined by the minimal of a/2A2 + η/4A4, with a < 0, yielding A2 = − 2a/η. Hence the cubic nonlinear term can fundamentally change the fate of the phase chart discussed before. In the following, we are interested in the effect of this nonlinear term on the vibration amplitude with strong damping.In this section, we will explore the effect of the nonlinear term on the collisions of ions. Assume A to be the vibration amplitude, then we require 52 ), with η = 0.01, we present the dynamics of five trapped ions in figure 14. While the pattern of the ions looks organized, the pattern of the velocities is not. As compared with figure 9, we find that the quasi-periodic dynamics no longer exist. Moreover, the exchange of velocities between the first and the last ion will not be observed; instead, the collision of ions will happen at different times (see figures 14(c) and (d)), leading to an exchange of velocities between adjacent ions. For this reason, the synchronization dynamics are not observed anymore. This is a typical feature in a nonlinear trap [41], in which the ions in different positions may feel some slightly different trapping frequency, therefore their dynamics and period will be different. We find that when η ∼ 0.001, the synchronization and the associated exchange of velocities can still be observed. However, in the much stronger nonlinear trap, the dynamics will become much more complicated. Furthermore, in the case of strong damping, the ions will vibrate around their equilibrium positions, leading to periodic dynamics, as observed in section 3 . In this case, the positions of ions may appear to have different vibration amplitudes and different effective periods. Thus, the synchronization comes from the parametric driving field and the harmonic potential, in which all ions tend to have the same behavior, including the same relative phase, which maybe measured in the future experiments.
$\begin{eqnarray}\eta {A}^{2}\sim a\sim { \mathcal O }(1),\end{eqnarray}$
in which condition the nonlinear term is important. By solving equation (
Figure 14. Dynamics evolution of a system of five ions with asymmetric initial conditions in the stable region, with a = 5, q = 1.7, γ = 0.0 and η = 0.01. |
The micro-motion of trapped ions, introduced by the driving field, could have an influence on the dynamics of the trapped ions. Under ideal conditions, quasi-periodic dynamics would result in continuous energy absorption and release from the driving field, facilitating the exchange of velocities. However, in realistic setups, micro-motion could disrupt the precise timing of the collisions of ions, particularly in many-body systems with strong interactions. This latter may also lead to collisions with background gas particles, transferring enough energy to the ions to induce a transition to a melted state, known as an ion cloud [39]. In a Paul trap, ions without cooling mechanisms can gain enough energy to escape, and when multiple ions are stored, excess micro-motion can increase their secular motion (rf heating). However, the micro-motion’s impact is negligible when $| {a}_{i}| \gg {q}_{i}^{2}$, with i ∈ {x, y}. In high-resolution spectroscopy, micro-motion primarily causes a second-order Doppler shift and a Stark shift, both of which can affect the accuracy of quantum information experiments and precision measurements [44]. Additionally, environmental noise, could impact the precision and coherence of ion dynamics, potentially disrupting synchronization, velocity exchange, and damping the amplitude of oscillations [45].
7. Experiments and possible signatures
When working with numerical simulations, it is necessary to use dimensionless variables to simplify the equations and reducing the complexity of the problem. However, to interpret these numerical results in real life, it is crucial to convert these dimensionless variables back into dimensional terms, and estimate their values. Therefore, we will have to estimate the values of a, q and γ possibly to be used in experiments.
In a Paul trap, there is typically both dc and RF voltages applied to the electrode. Therefore, the total voltage should be $V(t)={V}_{\rm{dc}\,}+{V}_{\,\rm{RF}}\cos ({\rm{\Omega }}t)$. Created by this total voltage is the electric potential, which in one dimension reads as [46–50]
$\begin{eqnarray}\begin{array}{r}{\rm{\Phi }}(y,t)=\frac{{V}_{\rm{dc}\,}}{2{R}_{\rm{dc}}^{2}}{y}^{2}+\frac{{V}_{\rm{RF}}}{2{R}_{\,\rm{dc}}^{2}}{y}^{2}\cos ({\rm{\Omega }}t),\end{array}\end{eqnarray}$
then, the electric field can be defined as E = − ∇Φ, which will allow us to find the equation of motion of the ions using F = eE, yielding [4, 51–55]. $\begin{eqnarray}\begin{array}{r}m\ddot{y}=-2\alpha \dot{y}-\frac{e}{{R}_{\,\rm{dc}\,}^{2}}({V}_{\rm{dc}\,}+{V}_{\,\rm{RF}}\cos ({\rm{\Omega }}\tau ))y.\end{array}\end{eqnarray}$
Therefore, the expression of the parameters in the dimensional system [17, 53, 56–59] $\begin{eqnarray}\begin{array}{rcl}a & = & \frac{4e{V}_{\rm{dc}\,}}{m{{\rm{\Omega }}}^{2}{R}_{\,\rm{dc}}^{2}},\quad | q| =\frac{2e{V}_{\rm{RF}\,}}{m{{\rm{\Omega }}}^{2}{R}_{\,\rm{dc}}^{2}},\\ \gamma & = & \frac{\alpha }{m\omega },\quad t=\frac{{\rm{\Omega }}\tau }{2},\quad {\rm{\Omega }}=2\omega .\end{array}\end{eqnarray}$
Here, e is the charge of the ions, Vdc is the applied dc voltage, VRF is the applied radio-frequency voltage, m is the mass of the ions, ω is the system frequency, Ω is the driving force frequency and Rdc is a parameter that defines the geometry of the trap and represents the distance from the trap center to the electrode.In an experimental setup for the trapped ions, the system is bound to experience some heating effect due to multiple reasons, like collisions between the ions, exchange of energy with the environment, and fluctuations of the applied electric fields. Therefore, the trapped ions need to be cooled down, in order to maintain their stability. Many techniques can be used for this. In [39, 51, 60, 61], two different mechanisms for the damping force were discussed, the first one is the drag force due to residual rarefied gas, and the second one is the interaction with the blue-detuned laser beams. In both cases, the damping rate γ can be tuned in a wide range and in general the damping from the laser field dominates due to the Doppler effect. Especially, this damping force is essential for the cooling of trapped ions [62]. That is why a damping term was introduced in equation (56 ) [63]. Furthermore, a and q can be estimated knowing the values of experimental variables Vdc, VRF, e, m and Rdc. Using the above relations, we get 32 ), the strong damping regime can be achieved when ∣A∣ ≪ re, yielding 58 ) by changing of VRF. Much more details can be found in [46]. Obviously, these values can be tuned in a wide range in experiments, thus the phase chart can be studied using this platform. From equation (59 ), we require ∣γ∣ ≫ 1, assuming $| q| \sim { \mathcal O }(1)$, and from α ∼ ω, we expect that the cooling of the trapped ions can be reached in the time scale of ${ \mathcal O }(1/\omega )$ in the strong damping limit; see figures 5 and 11.
$\begin{eqnarray}\frac{a}{| q| }=\frac{{V}_{\rm{dc}\,}}{{V}_{\,\rm{RF}}},\end{eqnarray}$
which can be tuned in a wide range of experiments. Furthermore, using equation ( $\begin{eqnarray}\left|\frac{q}{2\gamma }\right|\ll 1.\end{eqnarray}$
Thus, we expect the physics discussed in this work to be researched using trapped ions. For example, in an experiment with trapped 171Yb+ ions [46, 64] with Ω = 2π × 2 MHz, Rdc = 5 mm, RRF = 1.5 mm, hc = 1.24 eV · μm, and mc2 = 171 × 1860 × 0.511 MeV, we can obtain the following dimensionless value $\begin{eqnarray}a=3.5886.\end{eqnarray}$
The similar magnitude of q can be obtained using equation (The interacting Mathieu equation for trapped ions has broader implications for systems such as quantum computation, where synchronization and velocity exchange dynamics are relevant. These phenomena are important to those observed in many-body quantum systems, such as those with interacting qubits [65], and could be crucial for enhancing qubit coherence and designing more robust quantum gates and algorithms [66]. Additionally, similar synchronization dynamics may be explored in other physical systems, such as coupled mechanical oscillators, potentially advancing applications like pattern recognition in engineering [67].
8. Conclusions
To summarize, we present a complete analysis of the interacting Mathieu equation with the Coulomb interaction, which can be realized using trapped ions. We investigate the physics in various conditions. In the strong damping limit, the dynamics can be solved using perturbation theory, showing that the damping coefficient can play the role of reducing the vibration amplitude about their equilibrium positions. In the weak damping limit, we observe the synchronization dynamics and the associated exchange of velocities between distant trapped ions, which resembles to the dynamics of Newton’s cradle. In this case, the phase chart for the stability and instability dynamics is independent of the Coulomb interaction, thus, it is the same as the single particle Mathieu equation. Finally, we study the effect of the nonlinear anharmonic term in the trapping potential, in which the system will always be stable for all the parameters. In this case, the desynchronization prohibits the exchange of velocities between distant ions. We estimated that these results are within the reach of the current experiments with trapped ions and expect these interesting dynamics to be observed.
During the preparation of this work, we were constantly asking the fundamental question of what are the most unique features in this many-body classical model? The generalization of the exact solvable single particle Mathieu equation to the many-body models and even to the quantum realms are quite obviously important and necessary [17, 68], which should exhibit some interesting physics. In most of the cases, the quasi-periodic solutions are the major concern in these investigations. From these investigations with the Coulomb interaction in this work, it was quite possible to provide some concrete answers to our previous question, which are, synchronization dynamics with coherent phases, collision-induced exchange of velocity, desynchronization dynamics from the anharmonic effect, and periodic oscillations in the strong damping limits. These results are not limited to small amplitude vibrations. Some of these results may be beyond the scope of the Mathieu equation (thus without the Floquet theorem), indicating the necessary investigation of the interacting Mathieu equation in future. These results can be found for other types of inter-particle interactions (including Van der Waals interactions) [46], in which when r → ∞, V(r) = 0. It is important that some of these predictions can be immediately verified in experiments.
The model considered in this work can be regarded as the smallest Newton’s cradle that has ever existed. The results found might contribute to the study of the dynamics of trapped ions, which can be generalized to higher dimensions [69–72], due to the current research of trapped ions from the perspective of quantum computations and quantum simulations. The interacting Mathieu equation is probably much more important due to the wide range study of this model in dynamical systems in terms of the Floquet theorem in mathematics. It is hoped that the experimental progress for the realization of this model can stimulate the investigation of this interacting Mathieu model from a pure mathematical perspective beyond the quasi-periodic solutions [73, 74]. Meanwhile, this model is also interesting in the presence of noise at finite temperature with both dissipation and fluctuation, which is connected by the Einstein relation in the Brownian motion between the temperature, damping rate and diffusion constant, which is always presented in experiments when the temperature is not sufficiently low [45, 75, 76].