Solutions of the Mathieu-Hill Equation for the Trapped Ion Harmonic Oscillator. Qualitative Discussion

1. Introduction

Parametric excitation of second order non-linear systems [3,4] has been the subject of intense investigations by means of different approaches [5,6], among which one has to also mention the method of normal forms. Ref. [7] uses such method to investigate the dynamical stability of a non-linear system characterized by a nonlinear Mathieu-Duffing equation of motion [8,9], in case of parametric excitation [10]. It was demonstrated that the intrinsic non-linear nature of the system induces a subharmonic region that was not reported for systems characterized by linear MH equations [9,11]. Normal modes analysis is also used in [12] to describe dynamical stability for trapped ion systems.

Investigations of the stability of physical systems that undergo periodic parametric driving [13], such as ions confined by oscillating electric fields (Paul traps) [14], is a subject of vivid scientific interest. The behaviour of these systems can be better explained based on an approach that employs the pseudopotential approximation, which explains well trapping of charged particles within a quadratic potential and resonances that arise out of parametric excitation [12]. In the pseudopotential approximation (for large values of the radiofrequency), both in the classical and quantum cases, the time periodic Hamilton function that characterizes an ion confined within an electrodynamic trap is replaced by an autonomous (time independent). Hamilton function. Furthermore, linear ion traps (LIT) that operate at two RF values represent a versatile tool to simultaneously confine two ion species of interest [15]. Parametric excitation is also employed to carry out mass-selective removal of ions out of an electrodynamic trap [16,17,18].

It is well established that a single trapped ion can be employed as a paradigm to test a choice of essential physical models realized as time-dependent quantum harmonic oscillators (QHO) [13,19,20,21,22,23,24,25,26,27,28,29,30]. A numerical modelling for a segmented Paul trap that exhibits four blade electrodes which generate the trapping field and two biasing rods, where the latter are employed to compensate micromotion [31], is performed in Ref. [32]. Such an approach delivers enhanced optical access for both fluorescence spectroscopy and individual ion addressing, a mandatory requirement to achieve ion crystals or perform quantum engineering and manipulation [33] of trapped ion quantum states [34,35]. Applications span atomic clocks [36,37,38], which are excellent tools to search for physics Beyond the Standard Model (BSM) [39,40,41], or novel Quantum Technologies (QT) based on ultracold ions [42,43]. Multipole linear Paul traps (LPTs) are also versatile tools to investigate cold ion-atom collisions and bring new experimental evidence of the phenomena involved [44,45,46,47].

Ref. [4] investigates Hamiltonian dynamics [48] of a single QHO [13,49,50] in presence of dissipation and parametric driving, while it establishes that a time-dependent parametric frequency incessantly drives the system out of its dynamical equilibrium state. In addition, fine tuning of the system parameters enables one to control the contention between dissipation and parametric driving, in a similar way in which Ref. [51] investigates the competition between multipole anharmonicities of the trap electric potential and a periodic kicking term (laser field), thus driving the system from a regular motion regime towards chaotic dynamics [52]. To achieve control of a QHO, a key issue lies in characterizing and suppressing the inherent noise [53,54]. The intrinsic noise spectrum of a trapped ion can be characterized and experimental demonstrations of QHO control are recently reported in literature [27].

This paper is intended as a follow-up of the review paper recently published in Photonics [2]. Section 2 starts from the analytical model introduced in [1]. The paper investigates the issue of the solutions of the MH equation that explains ion dynamics within an electrodynamic ion trap (EIT). We show that such equation can be expressed as the equation of a (quantum) HO. Then, we suggest two classes of complex solutions in order to solve the MH equation whilst we separate the micromotion from the secular motion of the ion. It is well established that the solution of the MH equation can be expressed as a Hill series [55,56,57,58]. We demonstrate such solution can be cast as a Fourier series, whereas the MH equation of motion turns into the equation of a (quantum) HO. We introduce two methods to solve the MH equation with complex solutions. The coefficients in the Fourier series are determined based on the initial conditions. We infer a system of linear equations that enables one to determine the constants that we introduce in the Fourier series solution. In case of the second method suggested, for a particular solution of the MH equation, we derive a recursive relationship between the coefficients that we introduce. By using this method we reproduce the MH equation which determines the Floquet exponent as a function of these coefficients, which also depend on the values of the characteristic value (eigenvalue) a and q parameter in the MH equation, which are real-valued. Both approaches are original. Section 3 applies the results previously derived in Section 2 and supplies an analytical model to characterize ion dynamics as bound or unbound. By considering the residual interaction, we supply the solutions of the MH equation and the corresponding coefficients along with the value of $a_0$ as functions of the Floquet exponent, which represents another original result. Finally, we demonstrate the MH equation can be regarded as a HO equation.

In Section A.1 we establish that a trapped ion can be assimilated with a kicked-damped HO model. We supply the solution of the damped HO in three cases, depending on the value of the discriminant for the equation of motion we investigate while we also perform a qualitative discussion. The value of the discriminant shows how many roots the equation exhibits. In Section B we analyse solutions of the MH equation for an ion confined within an electrodynamic (or radiofrequency - RF) trap, considered as a parametric oscillator (PO). Based on previous results [2], we review the Hill method employed to determine the Floquet coefficient. We discuss the solutions of the MH equation for a trapped ion treated as a HO and emphasize they can be bounded or unbounded.

The results in the paper apply to electrodynamic traps employed in ultrahigh-resolution spectroscopy, mass spectrometry (MS) and in the domain of quantum technologies (QT) based on ultracold trapped ions [59,60,61,62], with an emphasize on quantum sensors and quantum metrology.

2. Solutions of the Mathieu-Hill Equation

An iconic paper that reports on an elaborate analysis and discussion of the transition of few ion crystals into ion clouds is Ref. [1]. The paper demonstrates how the crystalline phase prevails until the Mathieu stability limit is attained, emphasizing that the corresponding transition can not be regarded as an order → chaos transition which occurs when a control parameter reaches a critical value. In addition, in the vicinity of this limit the system exhibits sensitivity to perturbations, which enables the experimenter to investigate crystal melting substantially forward the instability limit. The influence of the ion micromotion on ion crystal stability is also discussed in [1], where numerical modelling is used to characterize dynamical stability and discriminate between four possible regimes. Representations of large Coulomb crystals composed of trapped ions are reported in [44], where a Molecular Dynamics (MD) numerical modeling is employed to account for cold ion-atom elastic collisions that occur between Coulomb crystals and very light virtual atoms. A review of strongly coupled Coulomb systems consisting of ultracold trapped ions is performed in [63], while 2D crystals in zig-zag configurations and structural transitions are discussed in [64,65,66]. Recent approaches report direct observation of micromotion for a many body Coulomb crystal consisting of ultracold trapped ions, while novel techniques to measure the micromotion amplitude are implemented [67]. A numerical modelling is employed to determine phase transitions in small-ion Coulomb crystals confined in an electrodynamic ion trap (EIT) [68], where the phase transitions points are introduced as extremes of the interpolated functions employed. The novel conceptualization of fractal quasi-Coulomb crystals is introduced in [69], in case of surface electrodynamic traps (SET) characterized by a Cantor Dust electrode configuration. New experimental evidence of orientational melting that arises in a 2D crystal consisting of up to 15 ions is reported in [70].

An analytical and numerical modelling focused on dynamical stability in case of trapped ion systems is performed in [71], by extending the model introduced in [1,72]. Numerical modelling illustrates that the associated dynamics is either quasiperiodic or periodic, depending on the initial conditions. The evolution in time for a system of two coupled oscillators levitated in an RF trap is described using an analytical model that depends on the chosen control parameters, while ion dynamics is shown to be integrable only for discrete values of the ratio between the axial and the radial frequencies of the ion secular motion. A qualitative discussion of system dynamical stability is then performed by employing the Morse theory. The results are further applied to many body strongly coupled Coulomb systems (trapped ions), locally investigated in the vicinity of equilibrium configurations that designate ordered structures. These equilibrium configurations exhibit a large interest for Coulomb ion crystals (atomic clocks) or to achieve quantum information processing (QIP).

For example, the Kibble-Zurek mechanism (KZM) generally describes non-equilibrium dynamics along with topological defect occurrence for physical systems that undergo second order phase transitions. A generalized KZM that explains defect emergence in trapped ion systems is discussed in [73,74], where the analytical models proposed enable investigation of KZ physics in case of inhomogeneous systems.

We consider the MH equation for an ion confined in an EIT [1]

$$\ddot{x}+f\left(\tau \right)x=0,f\left(\tau \right)=a+2qcos\left(2\tau \right).$$

(1)

Eq. 1 is a linear differential equation with a periodic coefficient, as it is assumed that f is a periodic function $f\left(\tau +T_0\right)=f\left(\tau \right)$, with the particular case $T_0=\pi$. Because we refer to a Paul (RF) trap, $\tau =\Omega t/2$ is the dimensionless time [1,59], where $\Omega$ is the radiofrequency of the oscillating trapping voltage supplied between the electrodes. Moreover, the function $f:\mathbb{R}\to \mathbb{R}$ is real and continuous, of period $T_0>0$. Furthermore, we express the f function as a Fourier series [2,55,75,76]

$$f\left(\tau \right)=A_0+\sum _{k>0}{A}_{k}cosk\tau +\sum _{k>0}{B}_{k}sink\tau ,$$

(2)

where $A_0$ is a constant, whilst the $A_k$ and $B_k$ coefficients are considered as real

$$\sum _{k>0}\left(A_k^2+B_k^2\right)<\infty ,A_k,B_k\in \mathbb{R}.$$

(3)

The Floquet theory [58] explains that the general solution (which is a complex function) can be cast as [1]

$$x_1=Q\left(\tau \right)\Phi \left(\tau \right),$$

(4)

where the function

$$Q\left(\tau \right)=e^{i\mu \tau },\mu \in \mathbb{C},$$

(5)

is associated to the slow (or secular) ion motion, whilst $\Phi \left(\tau \right)=\Phi \left(\tau +T_0\right)$ represents a complex, periodical and twice differentiable function, that characterizes the micromotion [2,77]. In case of a MH equation, ion dynamics is stable as long as the Floquet exponent is a complex number $\mu \in \mathbb{C}$ [78,79]. Hence, the solution is written as

$$x=e^{i\mu \tau }\Phi \left(\tau \right)+e^{-i\mu \tau }{\Phi }^{*}\left(\tau \right)\in \mathbb{R},$$

(6)

which depends on the coefficients $A_k$ and $B_k$, with $\left(k>0\right)$ and $A_0$. ${\Phi }^{*}\left(\tau \right)$ denotes the complex conjugate of the previously introduced $\Phi \left(\tau \right)$ function. One can also expand $\Phi \left(\tau \right)$ as a Fourier series

$$\Phi \left(\tau \right)=C_0+\sum _{k>0}{C}_{k}cosk\tau +\sum _{k>0}{D}_{k}sink\tau .$$

(7)

If

$$\sum _{k>0}\left(A_k^2+B_k^2\right)\ll A_0,$$

then eq. (1) can be expressed as the equation of a (quantum) harmonic oscillator (HO) [28,29,30,80,81]

$$\ddot{x}+A_{0}x=0,$$

(8)

where $A_0$ is a known constant and the solution can be expressed as

$$x=Mcos{\omega }_1\tau +Nsin{\omega }_1\tau ,{\omega }_1=\sqrt{A_0},$$

(9)

with $M,N$ constants that are determined from the initial conditions. Further on, we illustrate two different methods to solve the MH equation:

• One chooses a solution expressed as $x=x_2+x_2^{*}$, where $x_2^{*}$ is the complex conjugate of $x_2$ so that

  $$x_2=e^{i\mu \tau }\left(C_0+C_{1}cos\tau +D_{1}sin\tau \right).$$

  (10)
  A double differentiation with respect to $\tau$ yields:

  $${\dot{x}}_2=e^{i\mu \tau }\left[i\mu C_0+cos\tau \left(D_1+i\mu C_1\right)+sin\tau \left(i\mu D_1-C_1\right)\right]$$

  (11)

  $$\begin{array}{c}\hfill {\ddot{x}}_2=e^{i\mu \tau }\left\{-{\mu }^2{C}_0+cos\tau \left[2i\mu D_1-C_1\left(1+{\mu }^2\right)\right]-sin\tau \left[D_1\left(1+{\mu }^2\right)+2i\mu C_1\right]\right\}.\end{array}$$

  (12)
  We now revert to the eq. (1) and introduce the expression of $f\left(\tau \right)$ supplied by eq. (2)

  $${\ddot{x}}_2+\left(A_0+A_{1}cos\tau +B_{1}sin\tau \right)x_2=0.$$

  (13)
  Eq. (13) is also verified for the complex conjugate $x_2^{*}$. Then, one uses eqs. (11) and (12) which are introduced in eq. (13) to further derive

  $$\begin{array}{c}-{\mu }^2{C}_0+cos\tau \left[2i\mu D_1-C_1\left(1+{\mu }^2\right)\right]-sin\tau \left[D_1\left(1+{\mu }^2\right)+2i\mu C_1\right]\hfill \\ \hfill +\left(A_0+A_{1}cos\tau +B_{1}sin\tau \right)\left(C_0+C_{1}cos\tau +D_{1}sin\tau \right)=0.\end{array}$$

  (14)
  Eq. (14) can also be cast as

  $$\begin{array}{c}-{\mu }^2{C}_0+A_0{C}_0+A_1{C}_1{cos}^2\tau +B_1{D}_1{sin}^2\tau +\left(A_1{D}_1+B_1{C}_1\right)sin\tau cos\tau \hfill \\ +cos\tau \left[A_0{C}_1+A_1{C}_0-C_1\left(1+{\mu }^2\right)+2i\mu D_1\right]\\ \hfill +sin\tau \left[A_0{D}_1+B_1{C}_0-2i\mu C_1-D_1\left(1+{\mu }^2\right)\right].\end{array}$$

  (15)
  One uses [82,83,84,85]

  $${sin}^2\tau ={\displaystyle \frac{1-cos2\tau }{2}},{cos}^2\tau ={\displaystyle \frac{1+cos2\tau }{2}},sin2\tau =2sin\tau cos\tau .$$

  (16)

  and finally, by performing the calculus in eq (14) with $k>1$, one infers the following system of equations

  $$\left\{\begin{array}{ccc}-{\mu }^2{C}_0+A_0{C}_0+{\displaystyle \frac{1}{2}}{A}_1{C}_1+{\displaystyle \frac{1}{2}}{B}_1{D}_1=0\hfill & \hfill & \\ B_1{C}_0+A_0{D}_1-D_1\left(1+{\mu }^2\right)-2i\mu C_1=0\hfill & \hfill & \\ A_1{C}_0+A_0{C}_1+2i\mu D_1-C_1\left(1+{\mu }^2\right)=0\hfill & \hfill \end{array}\right.$$

  (17)
  Using this system of equations, one can determine the coefficients $C_0,C_1$, and $D_1$, respectively. Hence, one determines the solution $x_2$ from eq. (10) and implicitly the solution x in eq. (8). We demonstrate that an ion confined within an electrodynamic trap can be treated as a HO and derive the associated MH equation of motion.
• As demonstrated in [58], one chooses a solution of the form $x_2=e^{i\mu \tau }\Phi$. Hence,

  $${\dot{x}}_2=e^{i\mu \tau }\left(i\mu \Phi +\dot{\Phi }\right)$$

  (18)

  $${\ddot{x}}_2=e^{i\mu \tau }\left(\ddot{\Phi }+2i\mu \dot{\Phi }-{\mu }^2\Phi \right)$$

  (19)

We introduce the solution $x=x_2+x_2^{*}$ in the MH equation (1) and derive

$$\begin{array}{c}{e}^{i\mu \tau }\left(\ddot{\Phi }+2i\mu \dot{\Phi }-{\mu }^2\Phi \right)+e^{-i\mu \tau }\left(\ddot{\Phi }-2i\mu \dot{\Phi }-{\mu }^2\Phi \right)\hfill \\ \hfill +f\left(\tau \right)\left(e^{i\mu \tau }+e^{-i\mu \tau }\right)\Phi =0.\end{array}$$

(20)

By performing a Fourier series expansion [86] one can write

$$\Phi =\sum _{k\in \mathbb{Z}}{c}_k{e}^{ik\tau },f=\sum _{k\in \mathbb{Z}}{a}_k{e}^{ik\tau },a_k^{*}=a_{-k},$$

(21)

where $a_k$ are arbitrary constants and $e^{i\mu \tau }=cos\mu \tau +isin\mu \tau$. Consequently, eq. (20) becomes

$$\begin{array}{c}\sum _{k\in \mathbb{Z}}{c}_k\left[-k^2-2\mu k-{\mu }^2\right]e^{i\left(\mu +k\right)\tau }+\sum _{k\in \mathbb{Z}}\left(-k^2+2\mu k-{\mu }^2\right)e^{i\left(k-\mu \right)\tau }\hfill \\ \hfill +\sum _{p\in \mathbb{Z}}{a}_p\sum _{k\in \mathbb{Z}}{c}_k{e}^{i\left(k+\mu +p\right)}+\sum _{p\in \mathbb{Z}}{a}_p\sum _{k\in \mathbb{Z}}{c}_k{e}^{i\left(k+p-\mu \right)\tau }=0.\end{array}$$

(22)

It can be noticed that the equation above contains terms such as $e^{i\mu \tau }\left({\displaystyle \sum _{k+p=n}}{a}_p{c}_k\right)e^{in\tau }$. One denotes $b_n={\displaystyle \sum _{k+p=n}}{a}_p{c}_k={\displaystyle \sum _{p\in \mathbb{Z}}}{a}_p{c}_{n-p}$. As a result, eq. (22) can be recast as

$$e^{i\mu \tau }E+e^{-i\mu \tau }{E}^{*}=0,$$

(23)

with

$$\left\{\begin{array}{c}E=\sum _{k\in \mathbb{Z}}\left[-c_k{\left(k+\mu \right)}^2+b_k\right]e^{ik\tau },\hfill \\ E^{*}=\sum _{k\in \mathbb{Z}}\left[-c_k{\left(k+\mu \right)}^2+b_k\right]e^{-ik\tau }.\hfill \end{array}\right.$$

(24)

By introducing eqs. (24) in eq. (23), one obtains

$$c_k{\left(k+\mu \right)}^2=b_k=\sum _{p\in \mathbb{Z}}{a}_{k-p}{c}_p,k\in \mathbb{Z},$$

(25)

which is amenable to

$$\sum _{p\in \mathbb{Z}}\left[a_{k-p}-{\left(k+\mu \right)}^2{\delta }_{kp}\right]c_p=0,$$

(26)

where ${\delta }_{kp}$ stands for the Kronecker delta function. One further denotes

$$A=\left[a_{k-p}-{\left(k+\mu \right)}^2{\delta }_{kp}\right].$$

(27)

Further on, we introduce

$$C=\left(\begin{array}{c}⋮\\ c_{-2}\\ c_{-1}\\ c_0\\ c_1\\ c_2\\ ⋮\end{array}\right).$$

(28)

Then $A·C=0$, where C is the column matrix introduced above. One further infers

$$detA=0,$$

(29)

which is exactly the MH equation that determines the Floquet exponent $\mu$ as a function of the $a_k$ coefficients (particularly as a function of the eigenvalue a and q parameter). It is assumed that $a_k=0$ for $\left|k\right|>2$ (Mathieu case) and $a_1=a_{-1}=0$. In such case

$$f\left(\tau \right)=a_0+a_2{e}^{2i\tau }+a_{-2}{e}^{-2i\tau }.$$

(30)

As explained above, one also supposes that $a_2=a_{-2}=q$ and $a_0=a$. Then

$$f\left(\tau \right)=a+2qcos2\tau ,$$

(31)

which is exactly the classical MH equation. Based on eq. (25) one derives

$$c_k{\left(k+\mu \right)}^2=a_0{c}_k+c_{k-2}{a}_2+a_{-2}{c}_{k+2}=a{c}_k+q\left(c_{k+2}+c_{k-2}\right),$$

(32)

with $q\ne 0$ and

$$\left\{\begin{array}{c}{c}_{k+2}=\alpha c_k+\beta c_{k-2},ifk>0,\hfill \\ c_{k-2}=\gamma c_k+\zeta c_{k-2},ifk<0.\hfill \end{array}\right.$$

(33)

where $\alpha ,\beta ,\gamma$, and $\zeta$ are known coefficients. When $|c_k|\sim 0$ and $\left|k\right|>2$, the column matrix C [described by eq. (28)] can be recast as

$$C=\left(\begin{array}{c}⋮\\ c_{-2}\\ 0\\ c_0\\ 0\\ c_2\\ ⋮\end{array}\right),$$

(34)

where it is assumed that $c_1=c_{-1}=0$ is an initial condition (hypothesis). In such case, the function $\Phi$ [in eq. (21)] can be expressed as

$$\Phi =c_0-2c_{2}cos2\tau ,$$

(35)

with $c_2=c_{-2}$. Then, for $k=0$ one finds

$$c_0{\mu }^2=a{c}_0+2c_{2}q⟼{\mu }^2=a+2{\displaystyle \frac{c_2}{c_0}}q.$$

(36)

We consider the initial conditions

$$\left\{\begin{array}{c}{c}_0=x_0,\hfill \\ 2{\displaystyle \frac{c_2}{c_0}}={\displaystyle \frac{q}{2}}.\hfill \end{array}\right.$$

(37)

Under these circumstances, eq. (36) is recast as

$${\mu }^2=a+{\displaystyle \frac{q^2}{2}}.$$

(38)

We return to eq. (1) and express f as

$$f={\mu }^2+f_{res}.$$

(39)

Then, one uses eq. (38) and the residual interaction function $f_{res}$ is

$$f_{res}=f-{\mu }^2,f_{res}=2qcos2\tau -{\displaystyle \frac{q^2}{2}},$$

(40)

which we discuss below.

3. Stability Diagram of the MH Equation

Tha standard (canonical) form of the Mathieu equation for an ion confined within an electrodynamic trap, characterized by the parameters a and q, is [2,55,57,76,93,94,95]

$${\displaystyle \frac{d^{2}w}{d{\tau }^2}}+\left[a-2qcos2\tau \right]w=0.$$

(45)

We have demonstrated in [2] this is a Sturm-Liouville system [96] characterized by periodic boundary conditions, while the system exhibits solutions exclusively for particular values of the eigenvalue a which depend upon q. It is also established that other solutions exist for every pair $(a,q)$, but we consider only periodic solutions [86].

In order to characterize ion dynamics, one of the possible approaches is to separate the ion micromotion from the secular motion [1]. The MH equation can be recast as (see eq. (1))

$$\ddot{x}+fx=0,\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{o}\mathrm{n}\mathrm{e}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{e}\mathrm{s}f={\mu }^2+f_{res},$$

(46)

where $f_{res}$ describes the residual interaction which accounts for the ion micromotion. In addition

$${\displaystyle \frac{1}{T_0}}{\int }_0^{T_0}{f}_{res}d\tau =0,T_0=\pi ,$$

(47)

so the micromotion term averages out during a period of the RF drive, which leads us to the case of an autonomous Hamilton function for the trapped ion. By separating the micromotion eq. (46) changes accordingly

$$\ddot{x}+{\mu }^{2}x+f_{res}x=0,$$

(48)

which stands for the equation of the HO [9,24,25,81,97] of period $2\pi /\sqrt{\left|\mu \right|}$ , with $\mu \in \mathbb{R},2\mu \notin \mathbb{Z}$. It can be observed that ion dynamics is decomposed into a slowly oscillating part (the secular motion) at frequency $\mu$ (see eq. 38) and a fast parametric drive (the micromotion) $f_{res}$. The secular ion motion (or the macromotion) is characterized by the equation

$$\ddot{x}+{\mu }^{2}x=0.$$

(49)

where $\mu$ is the Floquet exponent [55,58]. We now revert to eq. (25) and choose $k=0$. Consequently

$${\mu }^2=a_0+{\displaystyle \frac{c_{-2}}{c_0}}{a}_2+{\displaystyle \frac{c_2}{c_0}}{a}_{-2}.$$

(50)

From eq. (46) one derives $f_{res}=f-{\mu }^2$, which we multiply by x. We consider the solution $x=x_2+x_2^{*}$ and we employ eq. (21). Hence

$$f_{res}x=\left(\sum _k{a}_k{e}^{ik\tau }-{\mu }^2\right)\left(e^{i\mu \tau }\sum _p{c}_p{e}^{ip\tau }+e^{-i\mu \tau }\sum _p{c}_p{e}^{-ip\tau }\right).$$

(51)

One denotes

$$\varphi \left(\tau \right)=\sum _p{c}_p{e}^{ip\tau },{\varphi }^{*}\left(\tau \right)=\sum _p{c}_p{e}^{-ip\tau },$$

(52)

where ${\varphi }^{*}\left(\tau \right)$ stands for the complex conjugate of $\varphi \left(\tau \right)$. Because

$${\int }_0^{2\pi }{e}^{ik\tau }d\tau ={\displaystyle \frac{1}{ik}}{e}^{ik\tau }{|}_0^{2\pi }=0,\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}k\ne 0,$$

(53)

$${\int }_0^{\pi }{e}^{ik\tau }d\tau ={\displaystyle \frac{1}{ik}}\left(e^{ik\pi }-1\right)=\left\{\begin{array}{cc}0,\hfill & k\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}\hfill \\ 2i/k,\hfill & k\mathrm{o}\mathrm{d}\mathrm{d}\hfill \end{array}\right.$$

(54)

We revert to eq. (47) and denote

$$\begin{array}{c}I=\underset{0}{\overset{\pi }{\int }}{f}_{res}xd\tau =\sum _{p,k}{a}_k{c}_p\left[{\displaystyle \frac{e^{i\left(k+p+\mu \right)\pi }-1}{i\left(\mu +k+p\right)}}-{\displaystyle \frac{e^{-i\left(k+p+\mu \right)\pi }-1}{i\left(\mu +k+p\right)}}\right]\hfill \\ -{\mu }^2\sum _p{c}_p{\displaystyle \frac{1}{i\left(\mu +p\right)}}\left[e^{i\left(\mu +p\right)\pi }-e^{-i\left(\mu +p\right)\pi }\right]\\ \hfill =\sum _{p,k}{a}_k{c}_p{\displaystyle \frac{2sin\left(\mu +k+p\right)\pi }{\mu +k+p}}-{\mu }^2\sum _p{\displaystyle \frac{2sin\left(\mu +p\right)\pi }{\mu +p}}.\end{array}$$

(55)

As demonstrated in [2], by choosing $a_0=a,a_2=a_{-2}=2q,\left|a_k\right|=0$ for $\left|k\right|>2$, one derives [1]

$$\begin{array}{c}I=\left\{\left(\sum _p{c}_p{\displaystyle \frac{2sin\left(\mu +p\right)\pi }{\mu +p}}\right)\left(a_0-{\mu }^2\right)+2q\sum _p{c}_p\left[{\displaystyle \frac{2sin\left(\mu +p\right)\pi }{\mu +2+p}}+{\displaystyle \frac{2sin\left(\mu +p\right)\pi }{\mu -2+p}}\right]\right\}\hfill \\ \hfill =2\sum _p{c}_{p}sin\left(\mu +p\right)\pi \left[{\displaystyle \frac{a_0}{\mu +p}}+{\displaystyle \frac{2q}{\mu +2+p}}+{\displaystyle \frac{2q}{\mu -2+p}}\right]=0.\end{array}$$

(56)

By assuming $c_1=c_{-1}=0;\left|c_k\right|=0$ for $\left|k\right|>2$ and $c_2=c_{-2}$, for values $p=-2,0,2$, and the micromotion term averages out (see [1])

$$\begin{array}{c}{c}_0\left[{\displaystyle \frac{a_0-{\mu }^2}{\mu }}+q\left({\displaystyle \frac{1}{\mu +2}}+{\displaystyle \frac{1}{\mu -2}}\right)\right]\hfill \\ +c_2\left[{\displaystyle \frac{a_0-{\mu }^2}{\mu }}+q\left({\displaystyle \frac{1}{\mu +4}}+{\displaystyle \frac{1}{\mu }}\right)\right]\\ \hfill +c_{-2}\left[{\displaystyle \frac{a_0-{\mu }^2}{\mu }}+q\left({\displaystyle \frac{1}{\mu -4}}+{\displaystyle \frac{1}{\mu }}\right)\right]=0.\end{array}$$

(57)

which is further amenable to

$${\displaystyle \frac{c_2}{c_0}}=-{\displaystyle \frac{2\left(a_0^2-\mu \right)+2q\left(2+\frac{\mu }{\mu +2}+\frac{\mu }{\mu -2}\right)}{a_0-{\mu }^2+2q\left(\frac{\mu }{\mu +2}+\frac{\mu }{\mu -2}\right)}}.$$

(58)

In case when q is very small [1,2,76]

$$a_0={\mu }^2+{\displaystyle \frac{q^2}{2\left({\mu }^2-1\right)}}+{\displaystyle \frac{\left(5{\mu }^2+7\right)q^4}{32{\left({\mu }^2-1\right)}^3\left({\mu }^2-4\right)}}+\dots ,\mu \ne 1,2,3.$$

(59)

and the even solution can be expressed as

$${\displaystyle \frac{c_2}{c_0}}={\displaystyle \frac{q}{4\left(\mu +1\right)}}+{\displaystyle \frac{\left({\mu }^2+4\mu +7\right)q^3}{128{\left(\mu +1\right)}^3\left(\mu +2\right)\left(\mu -1\right)}}+\dots ,\mu \ne 1,2,$$

(60)

which completely determines the MH equation for an ion confined in a Paul trap, while it stands as the main result of the paper. In Ref. [1] the values used are

$$a={\mu }^2-{\displaystyle \frac{1}{2}}{q}^2,{\displaystyle \frac{c_2}{c_0}}={\displaystyle \frac{q}{4}},$$

(61)

as the value of the Floquet exponent $\mu$ was chosen small, and implicitly the characteristic value (eigenvalue) a and parameter q in the MH equation, both of them real-valued, as previously emphasized in Section 1. Nevertheless, this hypothesis does not coincide with the values used in [1].

It was demonstrated in [2] that for ion trap operating points that lie within the first stability region of the MH equation, and for small values of the eigenvalue a and q parameter, the coefficients $c_k$ rapidly converge to zero. Under such conditions, higher harmonics become insignificant and the fundamental frequency prevails. It was also demonstrated that the solution of the MH equation is the sum between a shift that corresponds to the fundamental frequency and a correction due to the occurrence of higher harmonics associated with ion dynamics. The problem of evaluating the maximum amplitude of stable oscillations under given initial conditions was approached. The maximum is characterized by upper bounds determined by the relative position of the electrodes. For optimum trap operating parameters, it is essential to consider all axes for which the equation of motion is of Mathieu type. A nonzero solution can be expressed as the combination of two linear dependent functions, which form a fundamental solution of the MH eq. Ion trajectories are investigated by building a family of ellipses, while the transfer matrix is supplied along with the maximum amplitude of stable oscillations [2].

The frontiers of the stability diagrams are defined by the equations [56,57,76,93]

$$a_0\left(q\right)=-{\displaystyle \frac{q^2}{2}}+{\displaystyle \frac{7q^4}{128}}-{\displaystyle \frac{29q^6}{2304}}+{\displaystyle \frac{68687q^8}{18874368}}+\dots ,$$

(62)

$$\begin{array}{c}{a}_1\left(-q\right)=b_1\left(q\right)=1-q-{\displaystyle \frac{q^2}{8}}+{\displaystyle \frac{q^3}{64}}-{\displaystyle \frac{q^4}{1536}}-{\displaystyle \frac{11q^5}{36864}}+{\displaystyle \frac{49q^6}{589824}}\hfill \\ \hfill -{\displaystyle \frac{55q^7}{9437184}}-{\displaystyle \frac{83q^8}{35389440}}+\dots ,\end{array}$$

(63)

$$a_2\left(q\right)=4+{\displaystyle \frac{5q^2}{12}}-{\displaystyle \frac{763q^4}{13824}}+{\displaystyle \frac{1002401q^6}{79626240}}-{\displaystyle \frac{1669068401q^8}{458647142400}}+\dots ,$$

(64)

$$b_2\left(q\right)=4-{\displaystyle \frac{q^2}{12}}+{\displaystyle \frac{5q^4}{13824}}-{\displaystyle \frac{289q^6}{79626240}}+{\displaystyle \frac{21391q^8}{458647142400}}+\dots .$$

(65)

When $a,\left|q\right|\ll 1$ and $r\ge 7$ (case when $a_r$ is approximately equal to $b_r$), the characteristic values of the frontiers of the stability region are described by the following power series approximation [9,57]

$$\begin{array}{c}\hfill a_r,b_r=r^2+{\displaystyle \frac{q^2}{2\left(r^2-1\right)}}+{\displaystyle \frac{\left(5r^2+7\right)q^4}{32{\left(r^2-1\right)}^3\left(r^2-4\right)}}+{\displaystyle \frac{\left(9r^4+58r^2+29\right)q^6}{64{\left(r^2-1\right)}^5\left(r^2-4\right)\left(r^2-9\right)}}+\dots .\end{array}$$

(66)

The analytical modeling performed in the paper is valid for all electrodynamic trap designs. It is is essential to establish whether ion trajectories are stable or unstable, depending on the experimental conditions, as well as the boundaries stable and unstable regions associated to the ion dynamics. It is established that these limits correspond to combinations of cosine and sine elliptic series, as shown in [2,9,55,58,75,76]. For such values solutions of the MH equations are periodic but not bound, while they determine the point at which the ion trajectory becomes unbound.

Appendix A.1. Harmonic Oscillator with Damping

It is generally known that a trapped ion can be assimilated with a kicked-damped HO model for an ion trap [24,109] or a parametric oscillator [50,51].

We consider the equation of motion for a damped HO (unidimensional case) [4,110]

$$m\ddot{x}=-kx-\alpha \dot{x},$$

(A1)

where k denotes the elastic constant and $\alpha >0$ is the damping coefficient. Finally, one derives

$$\ddot{x}+2\lambda \dot{x}+{\omega }_0^2=0,with{\displaystyle \frac{k}{m}}={\omega }_0^2,2\lambda ={\displaystyle \frac{\alpha }{m}}.$$

(A2)

A solution of this equation is $x=e^{rt}$, where r is either real $r\in \mathbb{R}$ or complex $r\in \mathbb{C}$. Then, eq. (A2) is cast as

$$r^2+2\lambda r+{\omega }_0^2=0.$$

(A3)

Further on, we will shortly discuss the solutions of eq. (A3).

• Case 1
  We denote

  $$\Delta =4{\lambda }^2-4{\omega }_0^2>0,$$

  (A4)

  which implies that $\lambda$ > $\omega$. The solutions of eq. (A3) are

  $$r_{12}=-\lambda \pm \sqrt{{\lambda }^2-{\omega }_0^2}.$$

  (A5)
  Then, a solution of eq. (A2) is

  $$x=C_1{e}^{r_{1}t}+C_2{e}^{r_{2}t}=C_1{e}^{-\lambda +\sqrt{{\lambda }^2-{\omega }_0^2}}+C_2{e}^{-\left(\lambda +\sqrt{{\lambda }^2-{\omega }_0^2}\right)},$$

  (A6)

  which describes damped oscillations that are non-periodical, with $C_1,C_2$ constants of integration.
• Case 2

  $$\Delta =0\Rightarrow \lambda ={\omega }_0.$$

  (A7)
  In such case the solution is expressed as

  $$r_1=r_2=-\lambda \Rightarrow x=(B+Ct)e^{-\lambda t},$$

  (A8)

  with $B,C$ constants.
• Case 3

  $$\Delta <0\Rightarrow {\omega }_0>\lambda .$$

  (A9)
  Then

  $$r_{1,2}^{\prime }=-\lambda \pm i\sqrt{{\omega }_0^2-{\lambda }^2},$$

  (A10)

  and

  $$x=C_3{e}^{{r_2}^{\prime }t}+C_4{e}^{{r_2}^{\prime }t}=e^{-\lambda t}\left(C_3{e}^{i\sqrt{{\omega }_0^2-{\lambda }^2}t}+C_4{e}^{-i\sqrt{{\omega }_0^2-{\lambda }^2}}t\right),$$

  (A11)

  with $C_3=C_4^{*}$ and $C_4^{*}$ is the complex conjugate of $C_4$. Considering that

  $$e^{ix}=cosx+isinx\Rightarrow cosx={\displaystyle \frac{e^{ix}+e^{-ix}}{2}},sinx={\displaystyle \frac{e^{ix}-e^{-ix}}{2i}},$$

  (A12)

  eq. (A11) can be cast as

  $$x=e^{-\lambda t}\left(D_{3}cos\sqrt{{\omega }_0^2-x^2}t+D_{4}sin\sqrt{{\omega }_0^2-x^2}t\right),$$

  (A13)

  with $D_3,D_4$ constants.

Appendix B.1. Hill’s Method

By using Hill’s method one derives a homogeneous system of linear equations [2]

$$c_n+{\gamma }_n\left(c_{n-1}+c_{n+1}\right)=0,n=0,\pm 1,\pm 2,\dots .$$

(A20)

with

$${\gamma }_n={\displaystyle \frac{q}{\left(2n-i\mu \right)}}-a.$$

(A21)

The infinite determinant of the system is [2,58,94,104]

$$△\left(\mu \right)=\left|\begin{array}{ccccccc}\dots & \dots & \dots & \dots & \dots & \dots & \dots \\ \dots & 1& {\gamma }_{-2}\left(\mu \right)& 0& 0& 0& \dots \\ \dots & {\gamma }_{-1}\left(\mu \right)& 1& {\gamma }_{-1}\left(\mu \right)& 0& 0& \dots \\ \dots & 0& {\gamma }_0\left(\mu \right)& 1& {\gamma }_0\left(\mu \right)& 0& \dots \\ \dots & 0& 0& {\gamma }_1\left(\mu \right)& 1& {\gamma }_1\left(\mu \right)& \dots \\ \dots & 0& 0& 0& {\gamma }_2\left(\mu \right)& 1& \dots \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \end{array}\right|,$$

(A22)

and the Floquet exponent is determined by the equation

$$△\left(\mu \right)=0.$$

(A23)

The infinite determinant described by eq. (A22) is absolutely convergent, and it represents a meromorphic function [111] of the Floquet exponent $\mu$[55,58,75,103], with simple poles for $\mu =\pm i\left(\sqrt{a}+2s\right),;$$s=0,\pm 1,\pm 2,\dots$ [2,58]. Therefore

$$△\left(\mu \right)-{\displaystyle \frac{C}{cosh\left(\pi \mu \right)-cos\left(\pi \sqrt{a}\right)}},$$

(A24)

is an even and periodic meromorphic function. If the constant C is determined so that the function described by eq. (A24) has no pole at $\mu =i\sqrt{a}$, then it will exhibit no other pole, which makes it a constant. Because $△\left(\mu \right)\to 1$ as $\mu \to \infty$, such constant is equal to 1. To derive the expression of C, one chooses $\mu =0$ and infers [2,56,94]

$$\begin{array}{c}△\left(\mu \right)=1-{\displaystyle \frac{\left[1-△\left(0\right)\right]\left[1-cos\left(\pi \sqrt{a}\right)\right]}{cosh\left(\pi \mu \right)-cos\left(\pi \sqrt{a}\right)}}\hfill \\ \hfill ={\displaystyle \frac{cosh\left(\pi \mu \right)-1+△\left(0\right)\left[1-cos\left(\pi \sqrt{a}\right)\right]}{cosh\left(\mu \pi \right)-cos\left(\pi \sqrt{a}\right)}}.\end{array}$$

(A25)

Because the Floquet exponent is determined by the eq. (A23), one obtains

$$cosh\left(\pi \mu \right)=1+2△\left(0\right){sin}^2\left({\displaystyle \frac{1}{2}}\pi \sqrt{a}\right).$$

(A26)

Eq. (A26) defines the stability frontiers of the MH equation when $cosh\left(\pi \mu \right)=\pm 1$.

Ref. [2] analyses the dynamical stability of the MH equation for a trapped ion, when the trap operating points $a,q\in \mathbb{R}$ are located within the first stability region, focused on the associated dynamics is bounded or unbounded.

The issue of stability and instability intervals for the MH equation is very well explained in [56,76]. Further on, we review the analysis performed in [94]. If $a,q\in \mathbb{R}$ are both real, then eq. (A26) shows that $cosh\left(\pi \mu \right)\in \mathbb{R}$. When $-1<cosh\left(\pi \mu \right)<1$, then $\mu \in \mathbb{C}$ is a complex number, $i\mu$ is not an integer and any solution of the MH equation is bounded along the real axis x. The domains of stability located within the parameter plane $\left(a,q\right)$ are characterized by $-1<cosh\left(\pi \mu \right)<1$. When $cosh\left(\pi \mu \right)>1$ one can consider the Floquet exponent $\mu$ as real and non-zero. On the other hand, if $cosh\left(\pi \mu \right)<-1$ one may consider $\mu -i$ as real and non-zero. In both cases the solution of the MH equation is not bounded along the real x axis. Unstable regions of the MH equations are characterized by $cosh\left(\pi \mu \right)>1$ or $cosh\left(\pi \mu \right)<-1$. The stable and unstable regions are separated by curves characterized by $cosh\left(\pi \mu \right)=\pm 1$, where one solution is bounded and periodic, while the general solution stays unbounded [94].

One can distinguish amongst several cases [2,55,75,76]

• $\mu \in \mathbb{C}$ (pure imaginary) and $i\mu \notin \mathbb{Z}$
• the frontiers of the stability domains are defined by $i\mu \notin Z$ (not integer)
• the associated dynamics is unstable when $\mu \in \mathbb{R}$ or $\mu -i\in \mathbb{R}$, with $\mu =i\theta$

Stability of ion few crystals is firstly discussed in [1], based on a well established analytical model that is frequently used in the literature to describe regular and nonlinear ion dynamics in 3D QIT. Numerical modelling of quantum manifestation of order and chaos for ions confined in a Paul trap is explored in [112], where it is established that quasienergy states statistics serves to discriminate between integrable and non-integrable quantum dynamics. Double well dynamics for systems made of two trapped ions (in either a Paul or a Penning trap) is investigated in [72], where the RF trapping voltage effect towards enhancing or altering quantum transport in the chaotic separatrix layer is also examined. Dynamical stability of two-ion crystals in a Paul trap is explored in [113], based on the pseudopotential approximation [60,62,114]. Nonlinear dynamics of ions confined in a RF trap is reported in [115], where a new deterministic melting region is reported, along with crystallization in a secondary Mathieu stability region. Irregular dynamics of single ion dynamics in an EIT with axial symmetry is approached in [116], by employing both analytical and numerical modelling.

The dynamics of a trapped ion in a mass spectrometer, which undergoes the action of both quadrupole RF and dipolar DC excitation is investigated in Ref. [117], for a classical quadrupole 3D trap that exhibits hyperbolic geometry. It is demonstrated that dipolar excitation qualitatively changes the associated dynamics. The equation of motion is shown to be a classical MH equation that is perturbed with a constant inhomogeneous term, along with a small quadratic nonlinearity. An experimental approach implements motional parametric excitation and reports coherent spin-motion coupling of ions obtained by employing a spin-dependent force [118], for 2D crystals consisting of 100 ⁹Be ions confined in a Penning trap.

Appendix B.2. Electrodynamic Ion Traps (EIT) Operating under SATP Conditions. Damping Case

We will shortly discuss the case of an electrodynamic trap that operates under Standard Atmospheric Temperature and Pressure (SATP) conditions [92,119,120,121,122]. The equation of motion is

$${\displaystyle \frac{d^{2}x}{d{t}^2}}+2\lambda {\displaystyle \frac{dx}{dt}}+{\omega }^2\left(t\right)x=0,$$

(A27)

where the second term characterizes friction (damping) in air, whilst $\lambda$ stands for the damping coefficient. A solution to this equation is

$$x=e^{\rho t}f\dot{x}=e^{\rho t}\left(\rho f+\dot{f}\right),\ddot{x}=e^{\rho t}\left({\rho }^{2}f+2\rho f+dotf+\ddot{f}\right),$$

(A28)

where f is a periodic, time dependent function. By introducing eqs. (A28) in eq. (A27) and choosing $\rho =-\lambda$, one derives

$$\ddot{f}+\left({\omega }^2-{\lambda }^2\right)=0,$$

(A29)

which is exactly a MH type equation.

One of the first papers that reports trapping of macroscopic dust particles under SATP conditions and the occurrence of ordered structures is Ref. [123], which emphasizes that friction in air results in an efficient cooling of the trapped particle. Ordered structures are also reported in [119]. The dynamics of damped single charged particles levitated in a Paul trap is investigated in [124] based on an analytical model. The modified stability diagrams in the $\left(a,q\right)$ parameter space are derived, demonstrating that stable regions are not only enlarged but also shifted [2,105].

A multipole linear RF trap [125] is described in [126], which delivers an effective potential that characterizes three additional stable quasi-equilibrium points. The trap is used to levitate a group of charged silicate microspheres under SATP condition. The setup exhibits strong dependence on the RF field modulation and effective potential.