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 on the stability of systems subjected to 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 and resonances arising from parametric excitation [12]. What is more, 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 achieve mass-selective removal of ions from a Paul trap [16,17,18].

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

Hamiltonian dynamics of a single quantum harmonic oscillator (QHO) [13,46,47] in presence of dissipation and parametric driving is investigated in Ref. [4], where it is demonstrated that time-dependent parametric frequency drives the system perpetually out of equilibrium. What is more, fine-tune of the system parameters enables one to control the competition between dissipation and parametric driving. To achieve control of QHO, a key issue lies in characterizing and suppressing the inherent noise [48,49]. The intrinsic noise spectrum of a trapped ion can be characterized and experimental demonstrations of QHO control are reported [25].

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 further investigates the issue of the solutions of the MH equation that describes 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 solutions in order to solve the MH equation, based on separating the micromotion from the secular motion of the ion. It is well known that the solution of the MH equation can be expressed as a Hill series [50,51,52,53]. We now demonstrate that the solution of the MH can be cast as a Fourier series, while 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 from 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 infer a recursive relationship between the coefficients that we introduce. By using our method we reproduce the MH equation which determines the Floquet exponent as a function of these coefficients that also depend on the $\left(a,q\right)$ parameters in the MH equation. Section 3 uses the results obtained in Section 2 and discusses ion dynamics by separating the ion micromotion from the secular motion. By considering the residual interaction we supply the solutions for the MH equation and the corresponding coefficients along with the value of $a_0$, as functions of the Floquet exponent. The MH equation can be regarded as a HO equation.

In Section A.1 we demonstrate the 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 has. In Section B we investigate solutions of the MH equation for an ion in a Paul (RF) trap, treated as a parametric oscillator (PO). Based on previous results from [2], we review the Hill method employed to determine the Floquet coefficient. We demonstrate the solutions can be bounded or unbounded.

The results in the paper are applicable in ultrahigh-resolution spectroscopy, mass spectrometry (MS), and in the domain of quantum technologies (QT) based on ion traps [54,55,56,57], with an emphasize on quantum sensors and quantum metrology.

2. The Solution 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 on the fact that the corresponding transition can not be classified as an order → chaos transition which occurs at a critical value of a control parameter. What is more, in the vicinity of this limit the system exhibits sensitivity to perturbations, which enables the experimenter to investigate crystal melting well ahead the instability limit. The influence of the ion micromotion on ion crystal stability is also discussed in [1], while numerical modelling is used to characterize dynamical stability and discriminate between four possible regimes. Images of large Coulomb crystals consisting of trapped ions are reported in [42], where a Molecular Dynamics (MD) numerical modelling is used to account for the cold elastic collisions between Coulomb crystals and virtual very light atoms. A review of strongly coupled Coulomb systems consisting of ultracold trapped ions is performed in [58], while 2D crystals in zig-zag configurations and structural transitions are discussed in [59,60,61]. Recent approaches achieve direct observation of micromotion for multiple ions in a laser-cooled trapped ion crystal, while novel techniques to measure the micromotion amplitude are implemented [62]. A numerical approach to determine phase transitions in small-ion Coulomb crystals confined in an electrodynamic ion trap (EIT) is reported in [63], where the phase transitions points are introduced as extremes of the proposed interpolated functions. The novel concept of fractal quasi-Coulomb crystals is introduced in [64], for surface electrodynamic traps (SET) with the Cantor Dust electrode configuration. New experimental evidence of orientational melting occurring in a 2D crystal of up to 15 ions is reported in [65].

A classical study of the dynamical stability for trapped ion systems is performed in [66], based on the model introduced in [1,67]. The associated dynamics is shown to be quasiperiodic or periodic. The time evolution of two coupled oscillators in an RF trap is described using a model that depends on the chosen control parameters, and it is demonstrated that ion dynamics is integrable only for discrete values of the ratio between the axial and the radial frequencies of the secular motion. Then, the Morse theory is employed to qualitatively discuss system stability. The results are further extended to many body strongly coupled trapped ion systems, locally studied in the vicinity of equilibrium configurations that identify ordered structures. These equilibrium configurations exhibit a large interest for ion crystals or quantum logic.

For example, the Kibble-Zurek (KZ) mechanism plays a fundamental role in defect formation with universal scaling laws in nonequilibrium phase transitions. A generalized KZ mechanism that describes defect formation in trapped ion systems is discussed in [68,69], where the analytical models proposed enable investigation of KZ physics in inhomogeneous systems, with applications of interest in other nonequilibrium platforms.

We begin our treatment by considering the Mathieu-Hill (MH) equation for a trapped particle [1]

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

(1)

Equation (1) is a linear differential equation with a periodic coefficient, as one assumes that f is a periodic function $f\left(\tau +T_0\right)=f\left(\tau \right)$, with the particular case $T_0=\pi$. As we discuss a Paul (RF) trap, $\tau =\Omega t/2$ is the dimensionless time [1,54], where $\Omega$ is the radiofrequency of the oscillating trapping voltage supplied between the trap electrodes. What is more, the function $f:\mathbb{R}\to \mathbb{R}$ is real and continuous, of period $T_0>0$. Furthermore, the f function can be expressed as a Fourier series [2,50,70,71]

$$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 [53] explains that the general solution (which is a complex function) can be expressed as [1]

$$x_1=Q\left(\tau \right)\mathsf{\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 $\mathsf{\Phi }\left(\tau \right)=\mathsf{\Phi }\left(\tau +T_0\right)$ represents a complex, periodical and twice differentiable function, that characterizes the micromotion [2,72]. In case of a MH equation, ion dynamics is stable in case when the Floquet exponent is complex $\mu \in \mathbb{C}$ [73,74]. Hence, the solution can be cast as

$$x=e^{i\mu \tau }\mathsf{\Phi }\left(\tau \right)+e^{-i\mu \tau }{\mathsf{\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$. One can also expand $\mathsf{\Phi }\left(\tau \right)$ as a Fourier series

$$\mathsf{\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 Equation (1) can be expressed as the equation of a (quantum) harmonic oscillator (HO) [26,27,28,75,76]

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

(8)

where $A_0$ is a known constant and

$$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 as

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

  (10)
  By differentiating twice with respect to $\tau$ one finds:

  $${\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 Equation (1), where we introduce the expression of $f\left(\tau \right)$ supplied by Equation (2)

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

  (13)
  Equation (13) is also verified for the complex conjugate $x_2^{*}$. Then, one uses Equations (11) and (12) which are introduced in Equation (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)
  Equation (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 [77,78,79,80]

  $${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 Equation (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.
• As demonstrated in [53], one chooses a solution of the form $x_2=e^{i\mu \tau }\mathsf{\Phi }$. Hence,

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

  (18)

  $${\ddot{x}}_2=e^{i\mu \tau }\left(\ddot{\mathsf{\Phi }}+2i\mu \dot{\mathsf{\Phi }}-{\mu }^2\mathsf{\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{\mathsf{\Phi }}+2i\mu \dot{\mathsf{\Phi }}-{\mu }^2\mathsf{\Phi }\right)+e^{-i\mu \tau }\left(\ddot{\mathsf{\Phi }}-2i\mu \dot{\mathsf{\Phi }}-{\mu }^2\mathsf{\Phi }\right)\hfill \\ \hfill +f\left(\tau \right)\left(e^{i\mu \tau }+e^{-i\mu \tau }\right)\mathsf{\Phi }=0.\end{array}$$

  (20)

  By performing a series expansion one can write

  $$\mathsf{\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, Equation (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, Equation (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 Equation (24) in Equation (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 a and q adimensional coefficients). 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)

  and we have the classical MH equation. Based on Equation (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 Equation (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. In such case, the function $\mathsf{\Phi }$ (in Equation (21)) can be expressed as

  $$\mathsf{\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, Equation (36) is recast as

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

  (38)

  We return to Equation (1) and express f as

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

  (39)

  Then, one uses Equation (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.

• Discussion

When $\left|k\right|>1$, then $a_k=0$. In case of the HO [2,8,9,27]$a_0>0$ and from Equation (32) one infers

$$c_k{\left(k+\mu \right)}^2=a_0{c}_k,$$

(41)

and we distinguish two distinct cases

i)

$k=0\Rightarrow a_0={\mu }^2$ with $c_0\ne 0$.

ii)

$c_k\ne 0\Rightarrow {\left(k+\mu \right)}^2={\mu }^2=a_0\Rightarrow k+\mu =\pm \mu$ and the two solutions of the equations are

$$\left\{\begin{array}{c}{k}_1=0,\hfill \\ k_2=-2\mu \notin \mathbb{Z}.\hfill \end{array}\right.$$

(42)

where the latter case is discarded. Then, $\mathsf{\Phi }=c_0$ and

$$x=e^{i\mu \tau }{c}_0+e^{-i\mu \tau }{c}_0^{*},$$

(43)

which stands for the equation of the HO for $\mu \in \mathbb{R}$, whilst $c_0^{*}$ stands for the complex conjugate of $c_0$. Hence, the MH equation turns into the equation of a classical HO

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

(44)

By using the time dependent variational principle (TDVP) [81,82] and by expressing the Hamilton equations of motion for a degree of freedom, it was demonstrated that the quantum equation of motion in the Husimi $\left(Q\right)$ representation for a boson confined in a nonlinear ion trap results [23,83]. This equation of motion is shown to be fully consistent with the equation describing a perturbed classical oscillator. Hence, it is demonstrated that a phase-space representation of quantum mechanics, such as the Husimi or Wigner representation, reveals the structures of the corresponding phase-space [22] while also establishing a correspondence between the classical and quantum dynamics for such mesoscopic systems [84,85].

3. Stability Diagram of the MH Equation

Tha standard form of the Mathieu equation with parameters a and q is [2,50,52,71,86,87,88]

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

(45)

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 Equation (1))

$$\ddot{x}+fx=0,whereoneintroducesf={\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)

which means the micromotion term averages out during a period of the RF drive. By separating the micromotion Equation (46) changes accordingly

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

(48)

which stands for the equation of the HO [9,22,23,76,89] 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 }_x$ 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 [50,53]. We now revert to Equation (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 Equation (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 Equation (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{when}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{even}\hfill \\ 2i/k,\hfill & k \mathrm{odd}\hfill \end{array}\right.$$

(54)

We revert to Equation (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], if one considers $a_0=a,a_2=a_{-2}=2q,\left|a_k\right|=0$ for $\left|k\right|>2$, then

$$\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)

If $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$, one derives

$$\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,71]

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

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 values of the parameters a and q in the MH equation. Nevertheless, this hypothesis does not coincide with the values used in [1]. The frontiers of the stability diagrams are defined by the equations [51,52,71,86]

$$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,52]

$$\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)

4. Discussion

The paper investigates particular solutions of the classical Mathieu-Hill equation for ions confined in electrodynamic ion traps (EITs), by employing the Floquet theory [50,53]. It is well established that single trapped ions are an excellent tool to test a variety of important physical models realized as time-dependent (quantum) harmonic oscillators (QHO) [13,21,22,23,25,76,90,91]. In addition, investigations of dynamical stability for systems subjected to periodic parametric driving, such as trapped ions, is an issue of large scientific interest [12,92,93,94]. The time evolution for such systems is better explained within an integrated approach based on the pseudopotential approximation and the resonances induced by parametric excitation.

The paper suggests two methods to solve the MH equation. It is demonstrated that the solution of the MH equation can be expressed as a Fourier series. We apply the Floquet theory, separate the micromotion from the secular ion motion and then supply the general solution. Thus, the MH equation is shown to be similar to the equation of a HO. The paper introduces two different methods to solve the MH equation with complex solutions. In the first case we show that a linear system of equations results, which allows one to determine the coefficients of the MH equation solutions.

The second method suggested delivers the MH equation which yields the Floquet coefficient as a function of the $a_k$ coefficients in the series expansions we have used. The associated ion dynamics in a Paul trap is shown to be consistent with the equation of motion for a classical HO. Finally, we obtain the expression for $a_0$ and the coefficients of the even solution as a funcion of the Floquet coefficient, for very low values of the q parameter.

It has been established that a phase-space representation of quantum mechanics, such as the Husimi or Wigner representation, reveals the structures of the corresponding phase-space [22]. In particular, for the case of regular classical dynamics, the Husimi function of an eigenstate (or of a Floquet state in the case of a driven system) is localized in phase-space along the corresponding quantizing torus. By employing the TDVP and the Hamilton equations of motion, it was demonstrated that the law of dynamics for a boson confined in a nonlinear ion trap is equivalent to the equation of motion for the perturbed classical oscillator [23].

The paper suggests new methods to derive the solutions of the MH equation of motion for systems of trapped ions, demonstrating they are equivalent with the dynamics of the classical harmonic oscillator. The HO with damping is discussed in Section A, where the number of solutions depends on the sign of the discriminator associated to the 2nd order equation of motion. For a positive sign of the discriminator, the solution of the equation of motion describes non-periodic, damped oscillations. The parametric HO is discussed in Section B, based on the Floquet theory. It is demonstrated the solution are linearly independent and they establish a fundamental system of solutions. Then, we shortly discuss the Hill method to determine the Floquet coefficient [2,50,53,95,96] and supply the equation that defines the stability frontiers for the MH equation, followed by a discussion on the stability and instability intervals. A very short comment on EITs that operate under SATP conditions ends the paper.

A future paper will discuss the quantum mechanical dynamics of a charged particle in an oscillating quadrupole field [47,97], where the time-dependent wave functions are in a simple one-to-one correspondence with the wave functions of an ordinary static-field HO.

Applications of the methods used in the paper span MS, high-resolution spectroscopy and trapping of ion crystals, where the latter is of extreme importance for optical atomic clocks, quantum metrology and quantum information processing (QIP) with ultracold ions [28,32,98,99,100].

Finally, the paper is both an add-on and a follow-up to the paper recently published in Photonics [2].

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 [22,101] or a parametric oscillator [47,102]. The technique of single-frequency ion parking is a largely used soft ionization method employed in electrospray mass spectrometry (ESI-MS), especially intended for biopolymer analysis. It is based on gas-phase charge-reduction ion/ion reactions in an electrodynamic ion trap (EIT), in tandem with supplying a supplementary oscillatory voltage which selectively inhibits the reaction rate of an ion species of interest [103]. In addition, a forced-damped HO model is adapted for the digital 3D ion trap (DIT) geometry. What is more, Quadrupole Ion Trap - Mass Spectrometry (QIT-MS) [104] is also employed to perform analysis of micro-sized dry inorganic and bioparticles including red blood cells (RBCs) and different sizes of MCF-7 breast cancer cells [105]. Other applications include real-time tandem MS of secondary organic aerosol (SOA) components [106] or enhanced resolution MS [107]. Parametric excitation is also employed to achieve mass-selective removal of ions from an EIT (Paul trap) [17,18].

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

$$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, Equation (A2) is cast as

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

(A3)

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

1

Case 1

We denote

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

(A4)

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

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

(A5)

Then, a solution of Equation (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.

2

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.

3

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)

Equation (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,53,87,96]

$$△\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 Equation (A22) is absolutely convergent, and it represents a meromorphic function [109] of the Floquet exponent $\mu$[50,53,70,95], with simple poles for $\mu =\pm i\left(\sqrt{a}+2s\right),;s=0,\pm 1,\pm 2,\dots$ [2,53]. 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 Equation (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,51,87]

$$\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 Equation (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)

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

A discussion on the stability of the MH equation for a trapped ion, in case of trap operating points $a,q\in \mathbb{R}$ that lie within the first stability region, is performed in [2], pointing out when the associated dynamics is bounded or unbounded.

The issue of stability and instability intervals for the MH equation is very well explained in [51,71]. Further on, we present the analysis performed in [87]. If $a,q\in \mathbb{R}$ are both real, then Equation (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 stable regions in 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 along which $cosh\left(\pi \mu \right)=\pm 1$, where one of the solutions is bounded and periodic, while the general solution remains unbounded [87].

One can distinguish amongst several cases [2,50,70,71]

• $\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], where an analytical model is introduced that is frequently used in the literature to characterize regular and nonlinear dynamics for systems of trapped ions in 3D QIT. Numerical modelling of quantum manifestation of order and chaos for ions confined in a Paul trap is explored in [110], where it is demonstrated one can use the quasienergy states statistics to discriminate between integrable and chaotic regimes of quantum dynamics. Double well dynamics for two trapped ions (in a Paul or Penning trap) is investigated in [67], where the RF drive influence in enhancing or modifying quantum transport in the chaotic separatrix layer is also discussed. The stability of two-ion crystals in a Paul trap is explored in [111], based on the pseudopotential approximation [55,57,112]. The nonlinear dynamics of ions confined in a Paul trap is reported in [113], where the existence of a new deterministic melting region and crystallization in a secondary Mathieu stability region are observed. Irregular dynamics of single ion dynamics in an EIT with axial symmetry is approached in [114], by employing both analytical and numerical modelling. It is also established that period-doubling bifurcations represent the preferred route to chaos.

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. [115], 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 of motional parametric excitation, which demonstrates coherent spin-motion coupling of ions obtained with a spin-dependent force is reported in [116], for the case of 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 [85,117,118,119,120]. 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 Equations (A28) in Equation (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. [121], which emphasizes that friction in air results in an efficient cooling of the trapped particle. Ordered structures are also reported in [117]. The dynamics of single charged particles in a Paul trap in presence of the damping force is investigated in [122] based on an analytical model. The modified stability diagrams in the $\left(a,q\right)$ parameter space are derived, demonstrating that the stable regions are not only enlarged but also shifted [2,123].

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