Mathieu-Hill Equation Stability Analysis for Trapped Ions. Anharmonic Corrections for Nonlinear Electrodynamic Traps

1. Introduction

Linear differential equations (LDEs) with variable (periodic) coefficients are ubiquitous in both physics and engineering, but their solutions are generally identified by means of numerical simulations. In an effort to identify a solution for such system, it is essential to infer the so called Floquet or characteristic exponents that define a fundamental matrix associated to the system [1]. One of the most prevalent approaches used in case of LDEs with periodic coefficients is based on truncated Fourier series [2], whose coefficients are derived by means of the Harmonic Balance (HB) approximate method [3,4] that is employed to investigate non-linear oscillating systems described by ordinary non-linear differential equations (NLDEs) [5].

Mathieu functions of period $\pi$ or $2\pi$, also known as elliptic cylinder functions, were introduced in 1868 by Mathieu [6] together with the so-called modified Mathieu functions [7], in order to characterize the vibrations of an elastic membrane placed in a fixed elliptical hoop [8,9]. As the Mathieu equation does not possess closed-form analytic solution, its applications are affected with analytical [10] and numerical approximation schemes [11], along with nonlinear analysis of the associated stability charts [12,13]. Therefore, analytical solutions of Mathieu (Hill) or Duffing equations are generally investigated by means of various perturbation techniques [14,15,16], while other approaches also use non-perturbative techniques to characterize the parametric damping nonlinear Mathieu–Duffing oscillator [17].

During the last 150 years, Mathieu-type linear or nonlinear parametric differential equations have been the subject of large scientific interest, because of their frequent use in applied mathematics [18,19], quantum physics [20,21], quantum optics [22], engineering [12], analytical chemistry, mechanics [13,14], as well as in general relativity [23] and astrophysics [24]. For instance, the Mathieu-Hill equation is employed to characterize periodic orbits in a classical model of the magnetic hydrogen atom in Ref. [25]. An infinite Hill determinant may be computed, which enables one to infer the discriminant of the relevant Whittaker-Hill equation [26] that is used to characterize the stability of the H atom periodic orbits. Such method is less elaborate in comparison with the numerical integrations of the associated orbits.

Analytic periodic approximations for differential equations of Hill type are investigated in [3], where two different methods are employed in the straightforward case of a Mathieu equation. The first method is triggered by the HB method [4,27], with an aim to identify analytic approximations for the critical values and focus on periodic solutions of the Mathieu equation. It is demonstrated that kind of solutions are valid for any values of the parameter q. The second method exploits truncations of Fourier series [2].

1.1. Applications of Mathieu (Hill) Equation in Case of Electrodynamic Ion Traps. Nonlinear (Anharmonic) Traps. Kicked Mathieu-Duffing (Parametric) Oscillator

Scientists have always focused on developing new experimental techniques intended for confining single particles (neutral or electrically charged atoms and molecules, photons, electrons, antimatter, etc.) in a sharply defined region in space, under conditions of minimal perturbations (by creating an almost interaction free environment) [28,29,30,31,32]. This pursuit led to the development of radically new methods aimed at confining atomic and subatomic particles, which resulted in the maturation of optical tweezers, ion trap or magneto-optical (MOT) trap based techniques (including optical lattices) [20,33,34,35]. Prof. W. Paul had the idea of using a rotating or vibrating saddle–like electrostatic potential to confine ions or electrically charged particles [28,36] within a well defined area, under conditions of dynamical equilibrium [37,38]. Ion dynamics in a Paul (or RF) trap is characterized by a Mathieu-Hill type equation [20,39]. The Paul trap apparatus [in case of both 2D linear ion traps (LIT) and 3D versions] has been developed and refined for high finesse quantum engineering experiments, high precision spectroscopy [40,41], along with classical mass spectrometry (MS) [42,43,44,45,46,47] or chemical analysis [48], including the detection of aerosols and chemical warfare [49,50,51,52,53,54,55,56,57]. Besides, ion traps also enable exceptional control in preparing and manipulating atomic quantum states [58,59,60,61,62,63], which is why their wide area of applications also includes quantum logic [64,65,66,67,68], quantum sensing [69,70,71,72], quantum metrology [73,74] and even time fractals [75] or time crystals [76]. To these one adds high accuracy optical frequency standards [77,78,79,80], which are amongst the most sensitive quantum sensors [81,82] used to perform searches for physics beyond the Standard Model (BSM) [83,84,85] or to disseminate atomic time scales and redefine the SI unit of time, the second [86,87,88,89,90].

An electrodynamic trap is highly susceptible to geometric imperfections and electrode misalignment [44,91,92], that lead to the occurrence of local minima in the trapping potential and nonlinear resonances in the ion dynamics [47,93,94,95,96]. In case of 2D traps, motional coupling between axial and radial directions is reported [97,98]. Deviations of the trap potential from an ideal quadrupole lead to parasitic effects [99] and impose limitations on the number of ions that can be confined, and implicitly on the signal to noise ratio (SNR) [70,100,101]. The effects produced by the presence of higher order anharmonic terms [65] of the trapping potential (for linear strings of trapped ions) is investigated in [102], where two distinct effects are emphasized: a) an alteration of the oscillation frequencies and amplitudes of the ions’ normal modes of vibration in case of many-ion crystals, because each ion experiences a different curvature in the trap potential, along with b) significant amplitude-dependent shifts of the normal-mode frequencies, triggered by increased anharmonicity or higher excitation amplitude. As the ratio between the anharmonic and harmonic terms typically increases when the ion-electrode distance is reduced, it is demonstrated that anharmonic effects become more critical by lessening the ion trap scale. This requires special care, particularly in quantum computing with trapped ions [70,103] or in case of single ion [104,105,106] and multi-ion optical clocks [79,107,108,109]. On the other hand, it was demonstrated that an increased dodecapole potential generated in an asymmetric linear ion trap (LIT) [97] opens new pathways for exciting applications in mass spectrometry or Coulomb crystals. Such a setup also enables dark ion (inappropriate confined species) elimination from a trapped ion-chain consisting of physical qubits [47], or performing searches for time and parity violations that could constrain sources of new physics BSM [83,110,111].

Ref. [112] introduces a flow equation to approximate the solution of a weakly nonlinear Mathieu equation (NME), employed to characterize ion dynamics in the neighbourhood of the stability boundary of ideal traps. The HB method [4] is employed in [113] to explore the coupling effects of hexapole and octopole fields on ion dynamics in a quadrupole ion trap (QIT). Ion motion characteristics, e.g. motion center displacement, secular frequency shift, nonlinear resonance curve and buffer gas damping effects have been investigated. It turns out that hexapole fields have larger impact on ion motion center displacement, while octopole fields are mainly responsible for the ion secular frequency shift. In addition, the nonlinear features induced by hexapole and octopole fields may enhance or cancel each other. Ion dynamics in nonlinear Paul traps is investigated in [114] using the theoretical HB method, which allows one to derive the analytical ion trajectory and ion motion frequency in the superimposed octopole field by solving the NME [99]. The HB method is then validated by means of the numerical fourth-order Runge-Kutta (4th RK) method. The incremental HB method alongside with the path-following technique is applied in Ref. [115], where the Mathieu-Duffing nonlinear oscillator model (a Mathieu equation with cubic nonlinear term [4,116]) is employed to explore the steady-state response of a nonlinear piezoelectric energy harvester which undergoes external and parametric excitation. Ref. [115] also establishes that for some specific combination of the system parameters, vibration amplitudes and harvested power can be amplified up to three or five times in comparison to the classical broadband nonlinear energy harvester based on the forced (kicked) Duffing oscillator (DO) [15,116,117]. A method to generate stability plots for the Mathieu equation in case of a toroidal ion trap mass analyser is presented in [118].

The response of a DO to a harmonic excitation in the presence of (viscous) damping is known to exhibit, among other features, hysteretic and chaotic behaviors [1,13,98,116,119]. It is demonstrated that in particular cases the damping force may induce instability in the ion dynamics [120,121]. For instance, Ref. [122] explores the role of field inhomogeneities in altering the stability boundaries in nonlinear Paul traps, taking into account higher order terms in the equation of motion. The Poincaré-Lighthill-Kuo (PLK) method is employed in [123] to derive an analytical expression on the stability boundary and characterize the ion trajectory within a nonlinear Paul trap. The paper shows that a multipole superposition model (which essentially involves the octopole component) explains quite well how the field inhomogeneities shift the stable trapping region.

Particle dynamics in a nonlinear quadrupole Paul trap with octopole anharmonicity (which is a well-known dissipative system [38,124]), is described by a NME [4,125] which includes all perturbing contributions such as: damping, multipole terms of the potential and a harmonic excitation force (periodic kicking). To complicate the picture even more, the ion also undergoes interaction with a laser field [126], in presence of contributions from the hexapole and octopole fields that superpose over the harmonic trap potential. Hence, it was demonstrated that a damped parametric oscillator (PO) levitated in a Paul trap exhibits fractal properties and complex chaotic orbits, along with the emergence of strange attractors (which are often called fractal sets) [126,127]. Ion motion on the strange attractor exhibits sensible dependence on the initial conditions. Several recent papers also demonstrate that ion dynamics within a Paul trap can be assimilated with the equation that describes a damped driven DO [98,99]. The equation of motion for such system is similar to a perturbed Duffing-type equation [15], which is a generalization of the linear differential equation that describes damped and forced harmonic motion [128,129].

Quantum dynamics of general time-dependent three coupled oscillators based on using the unitary transformation method is explored in Ref. [130], which also applies to particle trapping in ion traps and electromagnetic fields. Ref. [131] explores classical and quantum integrability for a trapped ion Hamiltonian in 3D space, under conditions of a non-quadratic potential with superpositions of the hexapole and octopole terms in the series expansion of the electric potential [132]. The Painlevé analysis is employed to determine new integrable cases. It is also demonstrated that the 3D perturbed Hamiltonian is completely integrable in the sense of Liouville [131]. Late results also show the effective potential of an ion in 3D multipole fields in the mixed-mode is characterized by multiple isolated local minima. In addition, the specific number and spatial positions of the stable quasi-equilibrium points depend on the $a/q$ parameters ratio and the polarity of the electrodes DC voltage component [133]. An iterative method is employed to derive the stability parameters from measured secular frequencies, which is presented in Ref. [104].

1.2. Mass Spectrometry with Ion Traps. Late Developments

As the Mathieu equation describes the associated ion dynamics, the dependence of an electrodynamic (Paul) trap on both the ion mass and charge makes confinement of two ion species quite intricate. This is an outcome of the fact that heavier ion species exhibit a lower specific charge $Q/M$ ratio, which results in weaker confinement by the RF trapping field [134]. In contrast, confining electrons in a Paul trap requires drive frequencies around 1 GHz, while confinement of different ion species (Al, Ca, Ba, Sr, Yb) requires oscillating electric quadrupole fields in the range of tens of MHz. Commercial Off-The-Shelf (COTS) Quadrupole Mass Spectrometers (QMS) available on the market use oscillating quadrupole fields that operate at frequencies starting from the MHz range (in case of light ions) down to tens of kHz, when investigating heavy molecules characterized by low $Q/M$ ratio. Probability plots illustrate that heavy ions are more probably distributed near the maximum of their oscillation amplitude, while lighter ions position themselves at the center of their oscillation [135].

What is more, it is estimated that COTS MS and modified COTS hardware represent versatile tools that are capable of lessening the costs associated with exploration missions [136]. Hence, evaluation and optimization of COTS systems for space mission requirements is presumed to provide affordable analytical tools that enable a wider science community to bring a contribution in space science missions. An example would be the Mass Spectrometer Observing Lunar Operations (MSOLO) [137] that is intended to help NASA’s Volatiles Investigating Polar Exploration Rover (VIPER) mission science team to investigate the chemical composition of the lunar soil and search for water on the surface of the Moon. Another example is the Development and Advancement of Lunar Instrumentation Program (DALI) implemented by NASA. The DALI program has funded the development of the Environmental Analysis of the Bounded Lunar Exosphere (ENABLE) project. One of the tasks of ENABLE- performed at the Southwest Research Institute (sWRI) - lies in adapting a COTS QMS to operate as a space-qualified prototype instrument that could serve aboard multiple mission platforms to monitor environmental conditions under multiple lunar mission scenarios [138,139].

1.3. Ultraprecise Optical Atomic Clocks Based on Ultracold Ions. Current Directions of Action

Different ion trap geometries and electrode space arrangements are generally used (especially hyperbolic and cylindrical traps), in order to achieve a harmonic (electric) trapping potential around the trap centre. Whilst hyperbolic geometries surround the trap centre almost entirely, cylindrical ones usually exhibit open-endcap designs that allow better axial access to the centre and enhanced interaction between levitated particles and laser beams used for cooling or manipulation [140], which makes them suited for quantum logic, quantum optics and quantum metrology applications. For this reason, classical (hyperbolic) endcap electrode geometries are frequently employed when building single ion optical clocks, as they provide certain benefits such as a saddle shaped trapping potential which is essentially quadrupole, along with an open electrode structure that enables effective fluorescence detection [104,106]. Furthermore, a hyperbolic trap potential exhibits true azimuthal angle independence due to the electrode symmetry [141]. An open endcap trap design also exhibits a smaller quadrupole field component in the series expansion of the trap potential, while higher order terms must be considered when calculating the potential apart from the origin. Consequently, a careful trap design must consider and minimize issues such as potential frequency shifts and systematic errors in the atomic clock operation.

The sensitivity of optical atomic clocks [77,142] as quantum sensors is limited by the Standard Quantum Limit (SQL), imposed by the inherent projection noise of a quantum measurement [143]. Recent developments in atomic physics have enabled the experimental generation of many-body entangled states to boost the performance of quantum sensors beyond the SQL [60,71,144,145]. Furthermore, it is anticipated that the use of compound atomic clocks can enhance the stability of single ion clocks with long clock transition lifetimes to levels comparable to that of optical lattice clocks [80,106]. For ion species with shorter lifetimes, the stability can be improved directly by increasing the number of ions, but this approach requires special care in the selection of the atomic transition and offers the potential for a stability beyond the SQL , which could be a viable method to further improve the stability of optical clocks [109] and provide quantum-limited optical time transfer [87,146] to achieve intercontinental clock comparisons through a common-view node in geostationary orbit (GEO). Several missions such as NASA’s Deep Space Atomic Clock (DSAC) or several European Space Agency (ESA) missions are based on deploying optical clocks in space.

1.4. Structure of the Paper

All the aspects presented above motivate the interest towards further investigations on the Hill (Mathieu) class of equations, in order to identify new methods and techniques to mitigate anharmonic effects in case of custom made Paul traps used in laboratories, including anharmonic terms in the potential energy [147]. As the technology matures and the trap dimensions continuously scale down, studies are focused on exploring the anharmonic terms of the RF trapping potential for different implementations (geometries) of a quadrupole Paul trap [101,102]. Even if ion traps are intensely investigated both theoretically and experimentally, the effect of small perturbations from ideal conditions due to higher order (multipole) terms of the electric potential, trap asymmetry and geometric imperfections [92,98], along with patch potentials or other perturbing effects, is an issue of high interest. The paper brings new contributions towards this direction.

Section 1 is also conceived as a review of the applications of ion traps in quantum technologies, mass spectrometry and searches for physics BSM, where among the References the reader will discover some of the most recent advances in the domain. As traps used in experiments across world laboratories are far from being ideal, parasitic effects alter the dynamics of trapped ions and intense efforts are performed to mitigate and minimize their influence. Section 2 investigates the properties of the Mathieu-Hill equations. Based on the Floquet theorem use Floquet’s theorem and Hill’s Method solution, a recurrence relation is inferred. The Brillouin method is also explained. Section 3 investigates stability properties of the solutions of the Mathieu-Hill equation for a trapped ion. It is shown that the solution can be expressed as a Hill series, which introduced into the Mathieu equation yields a recurrence relationship that helps finding the characteristic (Floquet) exponent by employing Sträng’s technique, which is explicitly detailed in Appendix A.1. The stable solution of motion corresponds to a superposition of harmonic motions. The stability diagram for an electrodynamic (Paul) trap is illustrated as a function of the Mathieu function eigenvalues, along with the stability diagram for a combined (Paul and Penning) trap, where both axial and radial dynamics is considered. Ion trajectories in the phase space are discussed and we infer the transfer matrix along with the maximum amplitude of stable oscillations for trapped ions.

Section 4 approaches an issue of large interest by exploring the competition between the ion micromotion and the multipole anharmonicities of the trap electric potential, with an aim to characterize ion dynamics and eventually discriminate between ordered and irregular (chaotic) motion [126]. We supply the modified frontiers of the stability diagram in case of an electrodynamic trap with nonlinear octopole term of the electric potential by using the perturbation theory, as a function of the perturbation parameter we choose.

2. Mathieu-Hill Equations

The Hill differential equation is similar to Mathieu’s equation, but exhibits a more general nature. It arises in Hill’s method [148,149] to establish the motion of the Lunar Perigee and it represents a generalisation of the Mathieu equation [18,26,150]. The Mathieu-Hill equation is a homogeneous, second-order ordinary differential equation (ODE) with variable (periodic) coefficients [4,13]. It characterizes dynamical systems that exhibit intrinsic periodicity and parametric behaviour, such as modulation of radio carrier waves, transverse vibrations of a tense elastic membrane, stability of periodic motion in a non-linear system, the focus and defocus of beams in particle accelerators [151,152], or it is employed to explore planet dynamics [24]. The Hill equation takes the canonical or standard form of a reduced linear equation of the second order [8,10,19,153]:

$$\frac{d^{2}w}{d{\tau }^2}+J\left(\tau \right)w=0,$$

(1)

where J stands for a continuous and even function, periodic in $\tau$. The period is usually taken equal to $\pi$ for historical reasons [153]. The differential equation (1) was discussed by Mathieu in 1868 in connexion with the problem of vibrations of an elliptic membrane [6]. It is also assumed that the function $J\left(\tau \right)$ can be expressed as a Fourier series [26]

$$J\left(\tau \right)={\theta }_0+\sum _{s=1}^{\infty }2{\theta }_{s}cos2\tau ,$$

(2)

which converges within an infinite band in the $\tau$ plan that includes the real axis. If ${\theta }_s=0$ for $s\ge 2$, then eq. (1) turns into the Mathieu equation. What is more, eq (1) can also be arranged as:

$$\frac{d^{2}w}{d{\tau }^2}+\left[a-2q\psi \left(\tau \right)\right]w=0,$$

(3)

where $\psi \left(\tau \right)$ represents a continuous and even function of period $\pi$, while a and q denote adimensional parameters [154]. In the particular case of Mathieu’s general equation or Hill’s equation, a fundamental system of solutions consists of $e^{\mu \tau }\varphi \left(\tau \right)$ and $e^{-\mu \tau }\varphi \left(-\tau \right)$, as the equation is invariable to the change $-\tau \to \tau$. Consequently, Floquet demonstrated that the complete solution of Mathieu’s general equation (1) can be expressed as [26,150,153,155,156]:

$$w=A{e}^{\mu \tau }\varphi \left(\tau \right)+B{e}^{-\mu \tau }\varphi \left(-\tau \right),$$

(4)

where $\mu \in \mathbb{C}$ is called Floquet or characteristic exponent and it is a definite function of the a and q parameters, $\varphi$ is a periodic class $C^2$ function (twice continuously differentiable) [157], while $A,B$ stand for arbitrary constants. Then, based on the Floquet theorem and according to Hill’s method one may assume a series solution

$$w=e^{\mu \tau }\sum _{s=-\infty }^{\infty }{c}_{2s}{e}^{2is\tau }=\sum _{s=-\infty }^{\infty }{c}_{2s}{e}^{\left(\mu +2is\right)\tau }.$$

(5)

By introducing eq. (5) into eq. (1) one infers a recurrence relation:

$${\left(\mu +2is\right)}^2{c}_{2s}+\sum _{m=-\infty }^{\infty }{\theta }_{2m}{c}_{2s+m}=0,s\in \mathbb{Z},$$

(6)

where $\mathbb{Z}$ is the set of integer numbers, while ${\theta }_{-2m}={\theta }_{2m}$. If one tries to eliminate the term $c_{2s}$ in eq. (6), then a nonconvergent infinite determinant [26] would result. To avoid such an outcome, every expression in eq. (6) is divided to its central term ${\left(\mu +2is\right)}^2+{\theta }_0$. To determine the characteristic exponent, one multiplies the system matrix described by eq. (6) with a diagonal matrix. Then, a matrix results whose diagonal entries are each one equal to unity, assuming that none of the terms ${\left(\mu +2is\right)}^2+{\theta }_0$ vanishes. By denoting the determinant of this matrix as $▵\left(i\mu \right)$, the equation that determines the Floquet exponent $\mu$ is [26]

$$▵\left(i\mu \right)=0.$$

(7)

What is more, eq. (7) can be further arranged as:

$$cosh\pi \mu =1-2▵\left(0\right){sin}^2\frac{1}{2}\pi {\sqrt{\theta }}_0.$$

(8)

When $\mu$ has been determined, the $c_{2s}$ coefficients can be inferred in terms of $c_0$ and co-factors of $▵\left(i\mu \right)=0$. Thus, the solution of the Hill differential equation is complete. In case of a fairly rapid convergence of the determinant $▵\left(i\mu \right)$, eq. (7) can be used in its algebraic as well as recursive and explicit form, which represents Hill’s original method [9,26]. Evaluation of the determinant $▵\left(0\right)$ and use of eq. (8) represents an alternative method when this determinant converges quite well. A practical method to solve the Hill equation was suggested by Brillouin [153], based on the expression

$${sinh}^2\left(\frac{1}{2}\pi \mu \right)=w_1^{\prime }\left(\frac{\pi }{2}\right)w_2\left(\frac{\pi }{2}\right),$$

(9)

where $w_1$ and $w_2$ stand for two fundamental solutions. If a periodic fundamental solution exists for the Mathieu equation, then the existence of the other fundamental solution is forbidden and the characteristic exponent of the periodic solution corresponds to the a and q parameters of the characteristic curves that separate the stability domains. In case of the Hill equation, such property is not generally satisfied. Eq. (9) is useful to compare with Whittaker’s theory on the Hill equation [26]. The advantage of applying such method lies in the fact that it works with periodic $J\left(\tau \right)$ functions that exhibit a finite number of discontinuities.

For example, Ref. [158] investigates the modes which diagonalize the dynamical problem for linearly coupled Mathieu equations, which leads to the Floquet–Lyapunov transformation where the motion is associated with decoupled linear oscillators. The method is then used to solve the Heisenberg equations of the corresponding quantum-mechanical problem, and to determine the quantum wavefunctions for stable oscillations in the configuration (Hilbert) space. Such transformation and solution can be applied to more generic linear systems with periodic coefficients, such as coupled Hill equations and periodically driven parametric oscillators.

3. Stability of the Solutions of the Mathieu-Hill Equation for a Trapped Ion

The equation of motion for an ion confined within a Paul (electrodynamic) trap exhibits the standard form of the Mathieu equation [18,159,160,161]

$$\frac{d^{2}w}{d{\tau }^2}+\left[a-2qcos\left(2\tau \right)\right]w=0,$$

(10)

with $\tau =\mathsf{\Omega }t/2$ dimensionless, where $\mathsf{\Omega }$ is the RF of the trapping voltage $V_0$ applied between the electrodes. To investigate the stability of ion trajectories one uses the stability properties of the Mathieu equation solutions. As shown in eqs. (4) and (5), the solution of eq. (10) can be expressed as a Hill series

$$w\left(\tau \right)=A{e}^{\mu \tau }\sum _{s=-\infty }^{\infty }{c}_{2s}{e}^{2is\tau }+B{e}^{-\mu \tau }\sum _{s=-\infty }^{\infty }{c}_{2s}{e}^{-2is\tau },$$

(11)

where the A and B constants are determined from the initial conditions. If $i\mu \notin \mathbb{Z}$, then using eq. (11) one infers a fundamental system of solutions associated to eq. (10) for $A=1,B=0$ and $A=0,B=1$ (the Floquet theorem [18,26]). The characteristic exponent $\mu$ and the $c_{2s}$ coefficients are functions of the parameters a and q [159].

By introducing the solution supplied by eq. (5) in the Mathieu equation (10), the latter is cast into

$$\sum _{s=-\infty }^{\infty }{c}_{2s}\left[{\left(\mu +2is\right)}^2+a-2q\frac{e^{2i\tau }+e^{-2i\tau }}{2}\right]e^{\left(\mu +2is\right)\tau }=0,$$

(12)

which after matching the terms according to the powers of s results into a recurrence relation

$$-q{c}_{2s-2}+\left[{\left(\mu +2is\right)}^2+a\right]-q{c}_{2s+2}=0.$$

(13)

One further multiplies with $i^2=-1$ and then divides with the middle term in eq. (13), which leads to an infinite system of linear equations which are homogeneous in the $c_{2s}$ coefficients [155]:

$$\frac{q}{{\left(2s-\mu i\right)}^2-a}{c}_{2s-2}+c_{2s}+\frac{q}{{\left(2s-\mu i\right)}^2-a}{c}_{2s+2}=0,$$

(14)

which can be written as

$${\gamma }_{2s}{c}_{2s-2}+c_{2s}+{\gamma }_{2s}{c}_{2s+2}=0,s=0,\pm 1,\pm 2,\dots ,$$

(15)

with

$${\gamma }_{2s}=q{\left[{\left(2s-i\mu \right)}^2-a\right]}^{-1}.$$

(16)

The system of equations (15) admits a nontrivial solution if and only if the corresponding infinite determinant $▵\left(i\mu \right)$ vanishes for s noninfinite [26,162]

$$△\left(i\mu \right)=\left|\begin{array}{ccccc}\dots & & & & \\ & {\gamma }_{-2}& 1& {\gamma }_{-2}& & \\ & & {\gamma }_0& 1& {\gamma }_0& \\ & & & {\gamma }_2& 1& {\gamma }_2\\ & & & & & \dots \end{array}\right|.$$

(17)

Due to the symmetry of $△\left(0\right)$, ${\gamma }_{-j}={\gamma }_j$. Solving this determinant is not a simple issue, which is why one uses the approach of Whittaker [26]. The method is described in Appendix A, and it was introduced by Sträng [162,163]. The equation

$$▵\left(\mu \right)=0,$$

(18)

supplies the characteristic exponent $\mu$, while eqs. (15) enable one to recursively determine the coefficients $c_{2s}$. The infinite order determinant $▵\left(\mu \right)$ is absolutely convergent and it represents a meromorphic function [164,165] of $\mu$, with simple poles for $\mu =\pm i\left(\sqrt{a}+2s\right),;s=0,\pm 1,\dots$. Hence, eq. (17) is equivalent with [18]

$$cosh\left(\pi \mu \right)=1+2▵\left(0\right){\left[sin\left(\pi \sqrt{a}/2\right)\right]}^2.$$

(19)

Hereinafter we discuss the stability of the solutions of eq. (10) for $\tau ,a,q\in \mathbb{R}$. The solution of eq. (11) is stable and respectively unstable, if $\left|w\left(\tau \right)\right|$ is bounded, respectively unbounded along the $\tau >0$ semi-axis. From eq. (19) one infers $cosh\left(\pi \mu \right)\in \mathbb{R}$, which renders the solution of eq. (11) as stable when $\left|cosh\left(\pi \mu \right)\right|<1$, so for $i\mu$ real non-integer. If $\left|cosh\left(\pi \mu \right)\right|>1$, then $\alpha =\Re e\mu \ne 0$, while $\left|w\left(\tau \right)\right|$ is bounded along the $\tau$ axis. In such case $\left|w\left(\tau \right)\right|$ is bounded along the $\tau >0$ semiaxis only when $A=0,\alpha >0$ or $B=0,\alpha <0$. As an outcome, the $\left(a,q\right)$ plan is divided into stability regions with $\left|cosh\left(\pi \mu \right)\right|<1$ and instability regions characterized by $\left|cosh\left(\pi \mu \right)\right|>1$, separated by the curves $\left|cosh\left(\pi \mu \right)\right|=1$ for which a stable and periodic solution of eq. (10) exists, although the general solution is unbounded. The curves characterized by integer values of $i\mu$ are called characteristic curves. They are naturally introduced by means of the Sturm-Liouville (eigenvalue) problem [26,154,166] for the Mathieu functions $c{e}_m\left(z,q\right)$ and $s{e}_m\left(z,q\right)$ [19,161], which are treated as characteristic functions of eq. (10) with the limit conditions

$$\frac{dw\left(0\right)}{d\tau }=\frac{dw\left(\pi \right)}{d\tau }=0,w\left(0\right)=w\left(\pi \right)=0.$$

(20)

The $c{e}_m\left(z,q\right)$ and $s{e}_m\left(z,q\right)$ functions [7,167] are even and odd respectively, known up to a constant factor, with period $\pi$ when m is even, and with period $2\pi$ for m odd [18,26]. The characteristic curves $a_m\left(q\right)$ and $b_m\left(q\right)$, correspond to the Mathieu functions $c{e}_m\left(z,q\right)$ with $m=0,1,\dots$, respectively $s{e}_m\left(z,q\right)$, with $m=1,\dots$. In addition, the characteristic curves are analytical in q, while being characterized by the subsequent properties [26]

$$a_{2m}\left(q\right)=a_{2m}\left(-q\right),b_{2m+2}\left(q\right)=b_{2m+2}\left(-q\right),a_{2m+1}\left(q\right)=a_{2m+1}\left(-q\right),$$

(21)

$$a_{2m+1}\left(-q\right)=b_{2m+1}\left(q\right),a_m\left(q\right)<b_{m+1}\left(q\right)<a_{m+1}\left(q\right),$$

(22)

with $m=0,1,\dots$, and $q>0$.

The family of curves characterized by $\left|cosh\left(\pi \mu \right)\right|=1$ coincides with the family made up by the characteristic curves $a_m\left(q\right)$ with $m=0,1,\dots$, and $b_m\left(q\right)$ with $m=1,\dots$, which divide the $\left(a,q\right)$ plan into stability regions [ranging from $a_m\left(q\right)$ to $b_{m+1}\left(q\right)$ for $q\ge 0$, and from $a_m\left(q\right)$ to $a_{m+1}\left(q\right)$ or from $b_m\left(q\right)$ to $b_{m+1}\left(q\right)$ for $q\le 0$], and instability regions located below $a_0\left(q\right)$ or ranging from $b_m\left(q\right)$ to $a_m\left(q\right)$.

Ref. [168] shows that in case of microparticles levitated within a 2D linear Paul trap (LPT) operating under Standard Temperature and Pressure (STP) Conditions (in air), the Mathieu equations describing the trapping process are homogeneous along the x and y axes, with the well known solutions and stability domains. On the other hand, the z-axis motion is described by an inhomogeneous Mathieu equation. In this case, stability is obtained for (a) $i\mu \in \mathbb{R}$ ; (b) $\mu \in \mathbb{R}$ and $\mu <\mathsf{\Lambda }$ (c) $\mu -i\in \mathbb{R}$ and $|\mu -i|<\mathsf{\Lambda }$, where $\mu$ is the Floquet exponent, $\mathsf{\Lambda }=K/m\mathsf{\Omega }$, K stands for the coefficient which describes the drag aerodynamic force and m represents the particle mass. Hence, stability regions for solutions of the inhomogeneous equations of motion in presence of drag forces include both the stability regions of the homogeneous equations along with a part of the instability regions, which extends them considerably.

The Mathieu equation (10) exhibits $\pi -$ and $2\pi -$ periodic solutions on continuous stability curves $a=a\left(q\right)$, starting from points $a=n^2,n=1,2,3,\dots$. What is more, the periodic solutions and the boundaries between stability and instability regions in the $(a,q)$ parameter plane can be found by means of the Lindstedt-Poincaré method [16]. A cosine elliptic $c{e}_n$ and a sine elliptic $s{e}_n$ functions are associated to any value of n, where each one of these functions has its own characteristic number $a_{c{e}_n}$ and $a_{s{e}_n}$. The Mathieu functions and their eigenvalues (characteristic numbers) in power series of q, for $a,\left|q\right|\ll 1\left(q\simeq 0\right)$, are expressed as:

$$a_0\left(q\right)=a_{c{e}_0}\left(q\right)=-\frac{q^2}{2}+\frac{7q^4}{128}-\frac{29q^6}{2304}+\frac{68687q^8}{18874368}+\frac{123707q^{10}}{104857600}+\dots .$$

(23)

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

(24)

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

(25)

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

(26)

$$b_2\left(q\right)=a_{c{e}_2}\left(q\right)=4-\frac{q^2}{12}+\frac{5q^4}{13824}-\frac{289q^6}{79626240}+\frac{21391q^8}{458647142400}+\dots ,$$

(27)

The eigenvalue associated with the even solutions of the Mathieu functions, $c{e}_k(z,q)$ is labelled by $a_k\left(q\right),k=0,1,2,\dots$, while the one associated with the odd Mathieu function, $s{e}_k(z,q)$, is denoted as $b_k\left(q\right),k=1,2,\dots$, as illutrated in Figure 1. The higher rank terms that describe the frontiers of the stability regions are given in Appendix B. To graphically illustrate the stability diagram (associated to the Mathieu equation) that characterizes the dynamics of an ion confined within a combined quadrupole trap (a combination between a Penning and a Paul trap), we have used eqs. (23 - 27) and eqs. (A20 - A28). The result is shown in Figure 1.

The stable solution of the equation of motion (10) (for $\mu =i\beta ,\beta \in \mathbb{R}$) corresponds to a superposition of harmonic motions [20,31,101]

$$w\left(t\right)={\alpha }_1\sum _{s=-\infty }^{\infty }{c}_{2s}cos\left[\left(s+\frac{\beta }{2}\right)\mathsf{\Omega }t\right]+{\alpha }_2\sum _{s=-\infty }^{\infty }{c}_{2s}sin\left[\left(s+\frac{\beta }{2}\right)\mathsf{\Omega }t\right],$$

(28)

with the frequency spectrum

$${\nu }_s=\left(2s\pm \beta \right)\frac{\mathsf{\Omega }}{4\pi },s=0,1,\dots .$$

(29)

In Figure 3 we have illustrated the stability diagram for a combined (Paul and Penning trap). One notices that $a_0,b_1,a_1,b_2,a_2$ represent the frontiers of the stability diagram associated with the canonical Mathieu equation that describes axial motion. We also introduce

$$\begin{array}{c}{c}_0\left(q\right)=c-a_0\left(-q/2\right)/2=c+\frac{1}{16}{q}^2-\frac{7}{4096}{q}^4+\frac{29}{2304}\frac{q^6}{64}\hfill \\ \hfill -\frac{68687}{18874368}\frac{q^8}{512}-\frac{123707}{104857600}\frac{q^{10}}{2048},\end{array}$$

(30)

$$\begin{array}{c}{c}_1\left(q\right)=c-a_1\left(-q/2\right)/2=c-\frac{1}{2}+\frac{q}{4}+\frac{q^2}{64}-\frac{q^3}{512}+\frac{1}{1536}\frac{q^4}{32}+\frac{11}{36864}\frac{q^5}{64}\hfill \\ \hfill -\frac{49}{589824}\frac{q^6}{128}+\frac{55}{9437184}\frac{q^7}{256}+\frac{83}{35389440}\frac{q^8}{512}-\frac{12121}{15099494400}\frac{q^9}{1024}\end{array}$$

(31)

$$\begin{array}{c}{d}_1\left(q\right)=c-b_1\left(-q/2\right)/2=c-\frac{1}{2}-\frac{q}{4}+\frac{q^2}{64}+\frac{q^3}{1024}+\frac{1}{1536}\frac{q^4}{32}-\frac{11}{36864}\frac{q^5}{64}\hfill \\ \hfill -\frac{49}{589824}\frac{q^6}{128}-\frac{55}{9437184}\frac{q^7}{256}+\frac{83}{35389440}\frac{q^8}{512}+\frac{12121}{15099494400}\frac{q^9}{1024}\end{array}$$

(32)

Further on, we will consider trap operating points of $\left(a,q\right)$ that lie within the first stability region, with $0<\beta <1$ and $q\ge 0$. By minimizing the parameters a and $q\left(a/q\ll 1,q<0.4\right)$, the coefficients $c_{2s}$ (with $s\ne 0$) rapidly converge to zero, in such a way that higher harmonics are practically insignificant and the fundamental frequency ${\nu }_0=\beta \mathsf{\Omega }/\left(4\pi \right)$ prevails [169]. In such event $w\left(t\right)=w_0\left(t\right)+w_1\left(t\right)$, where $w_0\left(t\right)$ stands for the average of the $w\left(t\right)$ shift during the interval $\left[t,t+2\pi /\mathsf{\Omega }\right]$ and $\left|w_1\left(t\right)\right|\ll \left|w_0\left(t\right)\right|,\left|d{w}_1/dt\right|\gg \left|d{w}_0/dt\right|$. Under such circumstances, the solution of the Mathieu equation can be approxi-

mated as the solution of the equation

$$\frac{d^{2}w}{d{t}^2}=\frac{{\mathsf{\Omega }}^2}{4}\left[-a+q\left(\frac{a}{2}+2\right)cos\mathsf{\Omega }t-q^2{cos}^2\mathsf{\Omega }t\right]w.$$

(33)

By averaging eq. (33) over a period of the RF voltage (which is exactly the case of the pseudopotential approximation, when the electric forces that are time-averaged over an RF period generate a harmonic potential [20,31,170,171]), one derives

$$\frac{d^2{w}_0}{d{t}^2}=-{\beta }^2{w}_0,\beta =\sqrt{a+\frac{q^2}{2}}.$$

(34)

The shift $w_0$ in eq. (34) corresponds to the harmonic motion of fundamental frequency ${\nu }_0$, while the correction $w_1$ is determined by the higher harmonics $\left(2-\beta \right)\mathsf{\Omega }/4\pi$, $\left(2+\beta \right)\mathsf{\Omega }/4\pi ,\dots$. The axial frequencies are twice the radial frequencies: $a_z=-2a_r$ and $q_z=-2q_r$, while the trapped ions describe approximate Lissajoux trajectories when the ratio between the frequencies along the axes is $1:1:2$. The initial conditions have no effect whatsoever on the stability of ion trajectories, but they establish their position with respect to the area located between the trap electrodes for every operating point $\left(a,q\right)$ fixed within a stability domain. In this regard, the problem of evaluating the maximum amplitude of stable oscillations under given initial conditions is approached. The maximum admits upper bounds determined by the relative position of the electrodes. For optimum operating parameters it is mandatory to take into account all axes for whom the equation of motion is of Mathieu type.

One considers the stable solution of eq. (28) expressed as

$$w\left(\tau \right)={\alpha }_1{w}_1\left(\tau \right)+{\alpha }_2{w}_2\left(\tau \right),$$

(35)

with

$$w_1\left(\tau \right)=\sum _{s=-\infty }^{\infty }{c}_{2s}cos\left(2s+\beta \right)\tau ,w_2\left(\tau \right)=\sum _{s=-\infty }^{\infty }{c}_{2s}sin\left(2s+\beta \right)\tau ,$$

(36)

where both $w_1\left(\tau \right)$ and $w_2\left(\tau \right)$ are differentiable functions. If $w\left(\tau \right)\ne 0$ (which also implies that the coefficients ${\alpha }_1,{\alpha }_2\ne 0$), then $w_1$ and $w_2$ are linearly dependent functions. Hence, $w_1$ and $w_2$ form a fundamental system of solutions of the Mathieu eq. (10) and the Wronskian determinant [7] is

$$W=w_1\left(\tau \right){\dot{w}}_2\left(\tau \right)-{\dot{w}}_1\left(\tau \right)w_2\left(\tau \right)\ne 0.$$

(37)

Furthermore, the derivative of the Wronskian determinant is

$$\dot{W}=w_1\left(\tau \right){\ddot{w}}_2\left(\tau \right)-{\ddot{w}}_1\left(\tau \right)w_2\left(\tau \right)\ne 0.$$

(38)

The modulo of the solution of eq. (35) admits the following upper limit

$$\mathsf{\Lambda }=\sum _{-\infty }^{\infty }\left|c_{2s}\right|\sqrt{{\alpha }_1^2+{\alpha }_2^2}.$$

(39)

One denotes

$$d\left(\tau \right)=-W^{-1}\left[w_1\left(\tau \right){\dot{w}}_1\left(\tau \right)+w_2\left(\tau \right){\dot{w}}_2\left(\tau \right)\right],$$

(40)

$$b\left(\tau \right)=W^{-1}\left[w_1^2\left(\tau \right)+w_2^2\left(\tau \right)\right],c\left(\tau \right)=W^{-1}\left[{\dot{w}}_1^2\left(\tau \right)+{\dot{w}}_2^2\left(\tau \right)\right],$$

(41)

$$\epsilon =W{\mathsf{\Lambda }}^2{\left[\sum _{s=-\infty }^{\infty }\left|c_{2s}\right|\right]}^{-2}.$$

(42)

By eliminating the parameters ${\alpha }_1$ and ${\alpha }_2$ in eqs. (35) and (39), it follows that

$$c\left(\tau \right)w^2\left(\tau \right)+2d\left(\tau \right)w\left(\tau \right)\dot{w}\left(\tau \right)+b\left(\tau \right){\dot{w}}^2\left(\tau \right)=\epsilon .$$

(43)

Out of eq. (40) - (42) one infers

$$b\left(\tau \right)c\left(\tau \right)-d^2\left(\tau \right)=1.$$

(44)

Hence, for an initial given phase $\tau ={\tau }_0$, eq. (43) represents the equation of an ellipse of area $\pi \epsilon$ in the plan $\left(w,\dot{w}\right)$. For initial conditions given such that the point $\left(w\left({\tau }_0\right),\dot{w}\left({\tau }_0\right)\right)$ is located on the ellipse described by eq. (43) or within it, one finds $\left|w\left(\tau \right)\right|<M$ at any moment in time. The ions are not captured at the electrodes if all points w (with $\left|w\right|<M$) belong to the projection along the w axis of the region located within the electrodes.

To investigate ion trajectories in the phase space it will suffice to build the family of ellipses described by eq. (43), for $\tau$ ranging between $0\le \tau <\pi$. Indeed

$$\left(\begin{array}{c}w\left(\tau +\pi \right)\\ \dot{w}\left(\tau +\pi \right)\end{array}\right)=M\left(\tau \right)\left(\begin{array}{c}w\left(\tau \right)\\ \dot{w}\left(\tau \right)\end{array}\right),$$

(45)

with the transfer matrix expressed as

$$M\left(\tau \right)=\left(\begin{array}{cc}cos\pi \theta +d\left(\tau \right)sin\pi \theta & b\left(\tau \right)sin\pi \theta \\ -c\left(\tau \right)sin\pi \theta & cos\pi \theta -d\left(\tau \right)sin\pi \theta \end{array}\right).$$

(46)

The maximum amplitude of stable oscillations for trapped ions results from eq. (43)

$$w_{max}=\underset{0\le \tau \le \pi }{max}\sqrt{\epsilon b\left(\tau \right)}.$$

(47)

3.1. the Kicked Damped Parametric Oscillator

One considers the differential equation (DE) which describes a damped, kicked parametric oscillator (PO)

$$\ddot{u}+f\left(\tau \right)\dot{u}+g\left(\tau \right)u=h\left(\tau \right),$$

(48)

where $f,g$ and h are continuous functions of $\tau$, while $h\left(\tau \right)$ stands for the kicking term, which is usually a periodic. The homogeneous equation

$$\ddot{u}+f\left(\tau \right)\dot{u}+g\left(\tau \right)u=0,$$

(49)

exhibits the fundamental solutions which we denote as ${\phi }_1$ and ${\phi }_2$. Then, the general solution of eq. (48) is expressed as

$$\begin{array}{c}u\left(\tau \right)=c_1{\phi }_1\left(\tau \right)+c_2{\phi }_2\left(\tau \right)+{\phi }_2\left(\tau \right)\underset{{\tau }_0}{\overset{\tau }{\int }}\frac{{\phi }_1\left({\tau }^{\prime }\right)h\left({\tau }^{\prime }\right)}{W\left({\tau }^{\prime }\right)}d{\tau }^{\prime }\hfill \\ \hfill -{\phi }_1\left(\tau \right)\underset{{\tau }_0}{\overset{\tau }{\int }}\frac{{\phi }_2\left({\tau }^{\prime }\right)h\left({\tau }^{\prime }\right)}{W\left({\tau }^{\prime }\right)}d{\tau }^{\prime },\end{array}$$

(50)

with $W\left(\tau \right)={\phi }_1\left(\tau \right){\dot{\phi }}_2\left(\tau \right)-{\dot{\phi }}_1\left(\tau \right){\phi }_2\left(\tau \right)\ne 0$ is the Wronskian determinant. Eq. (49) changes into the normal form

$$\ddot{w}+J\left(\tau \right)w=0,$$

(51)

with

$$w\left(\tau \right)=u\left(\tau \right)exp\frac{1}{2}\underset{{\tau }_0}{\overset{\tau }{\int }}f\left({\tau }^{\prime }\right)d{\tau }^{\prime }.$$

(52)

In particular, if the J function is periodic and continuous, the normal form described by eq. (51) is a Hill equation. In case when f is a constant function $f\left(\tau \right)=\lambda >0$, one derives

$$u\left(\tau \right)=w\left(\tau \right)exp\left[-\frac{\lambda }{2}\left(\tau -{\tau }_0\right)\right],J\left(\tau \right)=g\left(\tau \right)-\frac{{\lambda }^2}{4}.$$

(53)

It can be noticed that the solution of eq. (48) is stable if $2\left|\Re e\mu \right|<\lambda$, with the Floquet exponent $\mu$ given by eq. (4). If w is stable, then definitely u is also stable. Nevertheless, there exist unstable solutions w for which the corresponding solutions u are stable. The frontiers of the stability domains are supplied by the equation $2\left|\Re e\mu \right|=\lambda$.

4. Anharmonic Corrections for Electrodynamic (Paul) Traps. Perturbation Method Analysis

An issue of large interest lies in exploring the competition between the ion micromotion and the multipole anharmonicities of the trap electric potential, in order to discriminate between stable and unstable (chaotic) dynamics [172,173,174,175]. Ref. [126] reports on the dynamics of an ion confined in a nonlinear Paul trap [61,122,176], assimilated with a time-periodic differential dynamical system. As such systems are characterized by low dissipation, exotic phenomena can be observed such as strange attractors, limit cycles, period doubling bifurcation and fractal basin boundaries [75,127,175]. We investigate the following equation of motion, which describes the axial dynamics for a particle of mass M and electric charge Q, confined within a quadrupole Paul trap, with anharmonicity derived from an octopole (quartic) electric potential $\lambda z^4/4$ [126,177]:

$$\frac{d^{2}z}{d{\tau }^2}+\left[a-2qcos\left(2\tau \right)\right]z+\lambda z^3=0,$$

(54)

where the adimensional parameters a and q are expressed as

$$\left\{\begin{array}{cc}a=-\frac{16Q{U}_0}{M{\mathsf{\Omega }}^2\left(r_0^2+2z_0^2\right)}\hfill & \\ q=\frac{8Q{V}_0}{M{\mathsf{\Omega }}^2\left(r_0^2+2z_0^2\right)}\hfill & ,\hfill \end{array}\right.$$

(55)

while $r_0$ and $z_0$ stand for the Paul trap radial and axial dimensions. Generally, in case of a typical 3D Paul trap, $a\le 0.05$ and $q\le 0.3$. The frequency of the applied a.c. voltage (micromotion) is denoted by $\mathsf{\Omega }$, $U_0$ and $V_0$ stand for the static and time-varying trapping voltages, respectively. Analytical modelling of the dimensionless equation (54) is based on employing techniques characteristic to the global bifurcation theory [14,128]. Numerical modelling is used to explore the associated dynamics and discuss chaos in such a nonlinear system [172,173,174,175].

Further on we introduce the perturbation parameter $\epsilon =2q$ and choose $\lambda =\epsilon \beta$, with the anharmonicity parameter $\beta$ fixed. Then, one can perform a series expansion of the a parameter and z coordinate as a function of the $\epsilon$ parameter:

$$a=\sum _{k=0}^{\infty }{\epsilon }^k{a}_k,z=\sum _{k=0}^{\infty }{\epsilon }^k{z}_k.$$

(56)

4.1. Solutions of the Mathieu Equation

The frontiers of the stability domains in the plan of the control parameters $q-a$ are determined by the periodic solutions of eq. (54) [104,126]. There is a requirement that z is a periodic solution, of period $2\pi$ and known parity, Within the limit $\epsilon \to 0$ the z solution is expressed as $cosn\tau$ or $sinn\tau$, with $n\in \mathbb{Z}$. One considers the case when $\epsilon =0(\epsilon =2q,\lambda =\epsilon \beta )$, with $\beta$ an anharmonicity parameter. Then, eq. (54) writes as

$$\frac{d^{2}z}{d{\tau }^2}+az=0,$$

(57)

whose solution can be expressed as $z=cosn\tau$. Then, $\ddot{z}=-n^{2}cosn\tau$ and eq. (57) changes into

$$-n^2+a=0\Rightarrow a_0=n^2.$$

(58)

From eq. (56) one infers

$$\ddot{z}=\sum _{k=0}^{\infty }{\epsilon }^k{\ddot{z}}_k,$$

(59)

$$z^3=\sum _{k=0}^{\infty }\sum _{k^{\prime }=0}^{\infty }\sum _{k"=0}^{\infty }{\epsilon }^{k+k^{\prime }+k"}{z}_k{z}_{k^{\prime }}{z}_{k"}=\sum _{p=0}^{\infty }{\epsilon }^p\sum _{\genfrac{}{}{0pt}{}{k,k^{\prime }\ge 0}{k+k^{\prime }\le p}}{z}_k{z}_{k^{\prime }}{z}_{p-k-k^{\prime }},p=k+k^{\prime }+k^{\prime \prime }.$$

(60)

We use the expression for $z^3$ derived above and revert to eq. (54)

$$\sum _{k=0}^{\infty }{\epsilon }^k{\ddot{z}}_k+a\sum _{k=0}^{\infty }{\epsilon }^k{z}_k-cos\left(2\tau \right)\sum _{k=0}^{\infty }{\epsilon }^{k+1}{z}_k+\beta \sum _{p=0}^{\infty }{\epsilon }^{p+1}\sum _{\genfrac{}{}{0pt}{}{k,k^{\prime }\ge 0}{k+k^{\prime }\le p}}{z}_k{z}_{k^{\prime }}{z}_{p-k-k^{\prime }}=0.$$

(61)

Then, eq. (61) can be cast into

$$\sum _{p=0}^{\infty }{\epsilon }^p{\ddot{z}}_p+a\sum _{p=0}^{\infty }{\epsilon }^p{z}_p-cos\left(2\tau \right)\sum _{p=0}^{\infty }{\epsilon }^p{z}_{p-1}+\beta \sum _{p=0}^{\infty }{\epsilon }^p\sum _{\genfrac{}{}{0pt}{}{k,k^{\prime }\ge 0}{k+k^{\prime }\le p-1}}{z}_k{z}_{k^{\prime }}{z}_{p-k-k^{\prime }-1}=0.$$

(62)

Therefore, by employing eqs. (56) and (54), after identifying the rank of the $\epsilon$ parameter, one infers a system of differential equations that recursively determines the $a_k$ parameters and $z_k$ coordinates:Eq. (62) can be expressed as

$$\frac{d^2{z}_p}{d{\tau }^2}+a_0{z}_p=z_{p-1}cos\left(2\tau \right)-\beta \sum _{\genfrac{}{}{0pt}{}{k,k^{\prime }\ge 0}{k+k^{\prime }\le p-1}}{z}_k{z}_{k^{\prime }}{z}_{p-k-k^{\prime }-1}-\sum _{k=0}^{p-1}{a}_{p-k}{z}_k,p\ge 1,$$

(63)

where $z_0=cos\left(n\tau \right)$ or $z_0=sin\left(n\tau \right)$, with $a_0=n^2,n=0,1,\dots$. Therefore, the frontiers of the Mathieu stability diagram are either even functions $a_{+}^{\left(k\right)}$ or odd functions $a_{-}^{\left(j\right)}$ that depend on the parameter q, which within the limit $q\to 0$ approaches $k^2$, respectively $j^2$. The k and j indexes are integer, with cu $k\ge 0$ and $j>0$. By performing a series expansion of these functions depending on q parameter $\left(\left|q\right|\ll 1\right)$, one derives the first frontiers of the stability diagram [178]

$$a_0\left(q\right)=\epsilon +2q\beta -\frac{q^2}{2}\left(1+4\beta \right)-\frac{q^3}{8}\beta \left(21+16\beta \right)+\dots ,$$

(64)

$$\begin{array}{c}{b}_1\left(q\right)=a_{-}^{\left(1\right)}=1+\epsilon -\frac{q}{2}\left(2+3\beta \right)-\frac{q^2}{32}\left(4+20\beta +9{\beta }^2\right)\hfill \\ \hfill -\frac{q^3}{256}\left(12+46\beta +60{\beta }^2+21{\beta }^3\right)+\dots ,\end{array}$$

(65)

$$\begin{array}{c}{a}_1\left(q\right)=a_{-}^{\left(1\right)}=1+\epsilon +\frac{q}{2}\left(2-3\beta \right)-\frac{q^2}{32}\left(4+4\beta -3{\beta }^2\right)\hfill \\ \hfill -\frac{q^3}{256}\left(4+20\beta -11{\beta }^2\right)+\dots ,\end{array}$$

(66)

$$\begin{array}{c}{b}_2\left(q\right)=4+\epsilon -\frac{3}{2}q\beta -q^2\left(\frac{1}{12}+\frac{1}{4}\beta +\frac{9}{128}{\beta }^2\right)\hfill \\ \hfill -\frac{q^3}{64}\beta \left(\frac{19}{6}+\frac{61}{20}\beta +\frac{15}{32}{\beta }^2\right)+\dots ,\end{array}$$

(67)

where $\beta =\left(\chi -\epsilon /6\right)/\left(2q\right)$ and $\chi$ is a constant.

Hence, the frontiers of the stability diagrams for the Mathieu equation with anharmonic perturbation are characterized by eqs. (64) - (67). It is evident how the frontiers of the stability diagram are shifted towards negative values along the a axis, in the control parameters plan $a-q$. Figure 4 and Figure 5 illustrate the stability diagram for the dynamics of an ion in the anharmonic trap we have investigated.

It can be noticed that the first stability region is bounded by the curves $a_0\left(q\right)$ a̧nd $b_1\left(q\right)$, while the second stability region is bounded by the curves $a_1\left(q\right)$ and $b_2\left(q\right)$, respectively. As opposed to the pseudopotential approximation when the electric RF potential is described by a polynomial of rank 2 in $z^2$ [20,30,31,114,156], consideration of the micromotion leads to stability regions that are qualitatively similar to those of the Mathieu equation for low enough values of c, but generally with frontiers sensibly altered as a function oi the anharmonicity parameter $\beta$. Hence, explicit analitic calculus which is numerically illustrated, enables one to establish the differences between the stability domains of the nonlinear dynamical system that is periodic in time (described by eq. (Figure 5)) and those of the autonomous dynamical system that is associated by means of the pseudopotential approximation [171].

4.2. the Frontiers of the Stability Diagram for the Mathieu Equation with Nonlinear Term

This section of the paper presents the technique we have developed to derive $a_0$. The frontiers of the stability diagrams of the Mathieu equation with nonlinear term (anharmonic perturbation) are given by eqs. (64 - 67), which correspond to the following limit solutions for $q=0$:

$$a_0=0;a_1=1and{z}_0=cos\tau ;a_1=-1and{z}_0=sin\tau .$$

(68)

The motion of a single charged particle in a nonlinear Paul trap, in presence of the damping force, is investigated theoretically in [120] and the modified stability diagrams in the parameter space are calculated. The results show that the stable regions in the $a--q$ parameter plane are not only enlarged but also shifted. Our results also show that the frontiers of the stability diagrams for the Mathieu equation are shifted towards negative values of the a parameter (within the plan defined by the control parameters $a-q$). In disagreement with the pseudopotential approximation for a Paul (RF) trap [114,171], when the electric potential is described by a polynomial of rank 2 in $z^2$, by taking into account the micromotion (when not in the pseudopotential approximation) stability regions result that are qualitatively similar to those of the Mathieu equation, but with frontiers that are shifted depending on the anharmonicity parameter $\lambda \left(\beta \right)$. Further on we express eq. (63) in case $p=1$ as

$$\frac{d^2{z}_1}{d{\tau }^2}+a_0{z}_1=z_{0}cos\left(2\tau \right)-a_1{z}_0-\beta z_0^3,$$

(69)

and we compute the right term for $z_0=cos\tau$. Then, the right hand of eq. (69) is cast into

$$cos\left(2\tau \right)cos\tau -a_{1}cos\tau -\beta {cos}^3\tau ,$$

(70)

which can be further expressed as (See Appendix C)

$$\left(2{cos}^2\tau -1\right)cos\tau -a_{1}cos\tau -\beta {cos}^3\tau .$$

(71)

One writes

$$cos3\tau =cos2\tau cos\tau -sin2\tau sin\tau =2{cos}^3\tau -cos\tau -2\left(1-{cos}^2\tau \right)cos\tau ,$$

(72)

as

$$cosx=\frac{e^{ix}+e^{-ix}}{2}.$$

(73)

Then, one derives

$$\begin{array}{c}{\left(cosx\right)}^n=\frac{1}{2^n}\sum _{k=0}^n{C}_n^k{e}^{ix\left(n-k\right)}{e}^{-ixk}=\frac{1}{2^n}\sum _{k=0}^n{C}_n^k{e}^{ix\left(n-2k\right)}\hfill \\ \hfill =\frac{1}{2^n}\sum _{k=0}^n{C}_n^k\left[cos\left(n-2k\right)x+isin\left(n-2k\right)x\right].\end{array}$$

(74)

Hence, eq. (74) becomes

$${\left(cosx\right)}^n=\frac{1}{2^n}\sum _{k=0}^n{C}_n^{k}cos\left(n-2k\right)x,$$

(75)

Using eq. (75) we infer

$$cos3\tau =4{cos}^3\tau -3cos\tau \mapsto co{s}^3\tau =\frac{1}{4}\left(cos3\tau +3cos\tau \right)$$

(76)

Then, eq. (69) becomes

$$\frac{d^2{z}_1}{d{\tau }^2}+a_0{z}_1=\left(2-\beta \right){cos}^3\tau -\left(a_1+1\right)cos\tau$$

(77)

By using eq. (76), eq. (77) changes accordingly

$$\frac{d^2{z}_1}{d{\tau }^2}+a_0{z}_1=\left(2-\beta \right)\frac{1}{4}cos3\tau +\left(\frac{1}{2}-\frac{3}{4}\beta -a_1\right)cos\tau .$$

(78)

The general solution of the harmonic oscillator (HO) writes as

$$z_{1g}=\alpha cos\tau +\beta sin\tau .$$

(79)

We choose an even solution, such as $z_{1g}=\alpha cos\tau$. Then, one tries a particular solution such as

$$z_{1part}=Acos3\tau +Bcos\tau ,$$

(80)

with

$${\ddot{z}}_{1part}=-9Acos3\tau -Bcos\tau .$$

(81)

Then, eq. (78) changes into

$$-8Acos3\tau =\frac{1}{4}\left(2-\beta \right)cos3\tau +\left(\frac{1}{2}-\frac{3}{4}\beta -a_1\right)cos\tau ,$$

(82)

which enables us to infer

$$a_1=\frac{1}{4}\left(2-3\beta \right),A=\frac{1}{32}\left(\beta -2\right)\Rightarrow z_1=cos\tau +\frac{1}{32}\left(\beta -2\right)cos3\tau$$

(83)

Hence, the solution of the Mathieu equation can be expressed as

$$z=z_0+z_1+\dots =1+cos\tau +\frac{1}{32}\left(\beta -2\right)cos3\tau +\dots .$$

(84)

More explicit details with respect to solving the Mathieu equation by means of the perturbation theory and determining approximations to all other Mathieu functions and eigenvalues are supplied in Appendix C.

5. Discussion

We use Floquet’s theorem, where any fundamental solution of Mathieu’s equation is compelled to satisfy the boundary value equations. In addition, the Wronskian determinant of the system does not vanish. Mathieu’s equations will admit an even an odd pair solution, which will define any other solution. We then assume Hill’s Method solution which we introduce into Mathieu’s equation, and infer a recursive relation as a function of the a and q parameters, as well as the coefficients of the Hill’s solution $c_s$. We show these coefficients $c_s$ exhibit non-trivial solutions if the infinite determinant satisfies $△\left(i\mu \right)=0$ (vanishes) for noninfinite s. Then, a holomorphic function with a single pole (similar to the determinant) is introduced, and we use the Liouville theorem for complex calculus and derive the expression of the Floquet characteristic exponent. It is assumed the Floquet characteristic exponent $\mu$ is chosen to satisfy the prerequisite such that the determinant described by eq. (17) vanishes. The solution will be unbounded unless the Floquet exponent is imaginary, and we supply this solution as a function of $△\left(0\right)$. To calculate this last determinant, one uses the Sträng recursion formula (method) [163], which is described in Appendix A.1.

We then discuss the stability of the Mathieu-Hill equation solution for a trapped ion and characterize ion dynamics as stable, respectively unstable, depending if the solution is bounded, respectively unbounded. We show the $\left(a,q\right)$ plan is divided into stability regions and instability regions, separated by the curves $\left|cosh\left(\pi \mu \right)\right|=1$ for which a stable and periodic solution of eq. (10) exists, although the general solution is unbounded. The curves characterized by integer values of $i\mu$ are called characteristic curves. They are naturally introduced by means of the Sturm-Liouville (eigenvalue) problem [26,166] for the Mathieu functions $c{e}_m\left(z,q\right)$ and $s{e}_m\left(z,q\right)$ [19], which are treated as characteristic functions of eq. (10) with limit conditions. We also show the periodic solutions and the boundaries between stability and instability regions in the $\left(a,q\right)$ parameter plane can be found by means of the Lindstedt-Poincaré method. We illustrate the stability diagram for the Mathieu equation in case of an electrodynamic trap, as a function of the associated eigenvalues and show the stable solution to be a superposition of harmonic motions.

We then represent the stability diagram for the Mathieu equation in case of a combined (Paul and Penning) trap, illustrating the frontiers for axial and radial motion, taking into account an extra parameter that is proportional to the cyclotronic frequency in a Penning trap.

Then, we consider trap operating points of that lie within the first stability region. We show that by minimizing the parameters a and q, the coefficients of the Hill series solution rapidly converge to zero, in such a way that higher harmonics are practically insignificant and only the fundamental frequency prevails. We infer the solution of the Hill equation in such case and find the solution by means of the pseudopotential approximation. We show the initial conditions have no effect whatsoever on the stability of ion trajectories, but they establish their position with respect to the area located between the trap electrodes for every operating point located within a stability domain. In this regard, the problem of evaluating the maximum amplitude of stable oscillations under given initial conditions is approached. The maximum admits upper bounds determined by the relative position of the electrodes. To investigate ion trajectories in the phase space we build a family of ellipses, while we supply the transfer matrix and find the maximum amplitude of stable oscillations. We exemplify our analytical model for the case of a damped, kicked PO and characterize the solutions and discuss the frontiers of the stability domains.

Axial stability for ion dynamics in a nonlinear Paul trap is investigated. In case of the octopole trap (with an anharmonicity of order 4) considered, the frontiers of the stability diagram have been explicitly determined as a series expansion of the control parameters chosen. We illustrate graphically that the frontiers of the stability diagrams for the Mathieu equation with anharmonic perturbation are shifted towards negative values along the a axis, in the control parameters plan $a-q$.

We consider the results to be of interest for 2D and 3D ion traps used for different applications in high-resolution spectroscopy, MS, and especially in the domain quantum technologies (QT) based on ion traps, among which optical clocks (as extremely high precision quantum sensors), quantum logic and quantum metrology.

Appendix A.1. Sträng’s Recursion Formula for △0

According to Sträng’s method [163], one defines

$$A_j=∥\begin{array}{ccccccccc}1& {\gamma }_{2j}& 0& & & & & & \\ {\gamma }_{2(j-1)}& 1& {\gamma }_{2(j-1)}& & & & & \\ 0& {\gamma }_{2(j-2)}& 0& & & & \\ & & & \dots & \dots & & & & \\ & & & & \dots & \dots & & & \\ & & & & & 1& {\gamma }_{2(j-2)}& 0\\ & & & & & {\gamma }_{2(j-1)}& 1& {\gamma }_{2(j-1)}\\ & & & & & 0& {\gamma }_{2j}& 1\end{array}∥,$$

(A9)

with ${△}_j=det\left(A_j\right)$ and $△\left(0\right)={lim}_{j⟶\infty }{△}_j$. In addition, $△\left(0\right)\in {△}_j$. Then $A_j$ can be decomposed in terms of $A_{j-1}$

$$A_j=∥\begin{array}{ccccc}1& {\gamma }_{2j}& & & \\ {\gamma }_{2(j-1)}& .& .& .& \\ & .& A_{j-1}& .& \\ & .& .& .& {\gamma }_{2(j-1)}\\ & & & & {\gamma }_{2j}& 1\end{array}∥.$$

(A10)

Using the Laplace decomposition method [180], the determinant of the above matrix is expressed as

$$det\left(A_j\right)=∥\begin{array}{cccc}.& .& .& \\ .& A_{j-1}& .& \\ .& .& .& {\gamma }_{2(j-1)}\\ & & {\gamma }_{2j}& 1\end{array}∥-{\gamma }_{2j}∥\begin{array}{ccccc}{\gamma }_{2(j-1)}& .& .& .& \\ & .& r{A}_{j-1}& .& \\ & .& .& .& {\gamma }_{2(j-1)}\\ & & & {\gamma }_{2j}& 1\end{array}∥.$$

(A11)

In the equation above $r{A}_{j-1}$ represents $A_{j-1}$ with its most right column discarded. According to Sträng [163], $lA$ stands for the matrix A with its leftmost column removed, $uA$ denotes the matrix A with its lowest row removed, while $dA$ is the matrix A with its upper most row removed. Finally, $uldr\left(A_{j-1}\right)=A_{j-2}$, and owing to the symmetry involved, $det\left(rd\left(A_{j-1}\right)\right)=det\left(ul\left(A_{j-1}\right)\right)$. Following this technique one finds[162,163]

$${△}_j={△}_{j-1}-2{\gamma }_{2j}{\gamma }_{2(j-1)}det\left(rd\left(A_{2(j-1)}\right)\right)+{\left[{\gamma }_{2j}{\gamma }_{2(j-1)}\right]}^2{△}_{j-2}.$$

(A12)

By using the Laplace decomposition method once more we infer

$${\mathsf{\Omega }}_j=det\left(ul\left(A_j\right)\right)=det\left(rd\left(A_j\right)\right),$$

(A13)

which allows us to obtain

$${\mathsf{\Omega }}_j=det\left(A_{j-1}\right)-{\gamma }_{2j}{\gamma }_{2(j-1)}{\mathsf{\Omega }}_{j-2},$$

(A14)

so as to

$$\frac{{△}_{j-1}-{\mathsf{\Omega }}_j}{{\gamma }_{2j}{\gamma }_{2(j-1)}}={\mathsf{\Omega }}_{j-1}=det\left(rd\left(A_{j-1}\right)\right),$$

(A15)

and

$${△}_j={△}_{j-1}+2\left({\mathsf{\Omega }}_j-{△}_{j-1}\right)+{\left[{\gamma }_{2j}{\gamma }_{2(j-1)}\right]}^2{△}_{j-2}.$$

(A16)

Out of eq. (A16) one infers

$${\mathsf{\Omega }}_j=\frac{{△}_j+{△}_{j-1}-{\left[{\gamma }_{2j}{\gamma }_{2(j-1)}\right]}^2{△}_{j-2}}{2}.$$

(A17)

By introducing eq. (A17) into eq. (A12) one derives

$$\begin{array}{c}{△}_j=\left[1-{\gamma }_{2j}{\gamma }_{2(j-1)}\right]{△}_{j-1}+\left[{\left({\gamma }_{2j}{\gamma }_{2(j-1)}\right)}^2-{\gamma }_{2j}{\gamma }_{2(j-1)}\right]{△}_{j-2}\hfill \\ \hfill +{\gamma }_{2j}{\gamma }_{2(j-1)}{\left[{\gamma }_{2(j-1)}{\gamma }_{2(j-2)}\right]}^2{△}_{j-3}.\end{array}$$

(A18)

Further on, one defines ${\alpha }_{2j}={\gamma }_{2j}{\gamma }_{2(j-1)}$ and $1-{\alpha }_{2j}={\delta }_{2}j$, which finally leads to

$${△}_j={\delta }_{2j}{△}_{j-1}-{\alpha }_{2j}{\delta }_{2j}{△}_{j-2}+{\alpha }_{2j}{\alpha }_{2(j-1)}^2{△}_{j-3}.$$

(A19)

It can be observed that one can recursively solve the equation above for $△\left(0\right)={lim}_{j⟶\infty }{△}_j$ to an accuracy as fine as desired. In addition, computer algebra programs such as Axiom, Maxima, Maple, Mathematica and SageMath, can be used to identify all stable values of $\mu$, which satisfy eq. (A8) as real values (i.e. all iso-$\mu$ values for which $\mu$ is exclusively imaginary).

Appendix C.1. Perturbation Theory

We revert to eq. (A32). If $n=0$, the only non-trivial, periodic solution of the equation, apart from an arbitrary multiplication constant, is $y_0\left(z+\right)=1$ [181]. Then, eq. (A32) turns into

$$\frac{d^2{y}_1}{d{z}^2}=2cos2z-a_1.$$

(A38)

By integrating eq. (A38) one obtains

$$y_1\left(z\right)=-\frac{1}{2}{a}_1{z}^2-\frac{1}{2}cos2z+c_{1}z+c_2$$

(A39)

which exhibits a periodic solution only if $a_1=0$ (implicitly the constants $c_1,c_2=0$). Therefore, the solution is

$$y_1\left(z\right)=-\frac{1}{2}cos2z.$$

(A40)

Hence. eq. (A33) can be expressed as

$$\frac{d^2{y}_2}{d{z}^2}+a_0{y}_2=2y_1\left(z\right)cos2z-a_2{y}_0\left(z\right)-a_1{y}_1\left(z\right).$$

(A41)

where

$$2y_{1}cos2z=-cos\left(2z\right)cos\left(2z\right)=-\frac{1}{2}\left[cos\left(4z\right)+1\right],$$

(A42)

and we have used

$$cosxcosy=\frac{1}{2}\left[cos\left(x+y\right)+cos\left(x-y\right)\right].$$

(A43)

As a result, eq. (A41) can be finally written as

$$\frac{d^2{y}_2}{d{z}^2}=-\frac{1}{2}\left[cos\left(4z\right)+1\right]-a_2.$$

(A44)

After integration one infers

$$y_2\left(z\right)=-\frac{1}{4}\left(1-2a_2\right)z^2+\frac{1}{32}cos\left(4z\right)+c_{3}z+c_4,$$

(A45)

which is periodic only if $a_2=-\frac{1}{2}$ (along with the requirement that $c_3,c_4=0$), and one derives

$$y_2\left(z\right)=\frac{1}{32}cos\left(4z\right).$$

(A46)

We revert to eq. (A35) and focus on the first term in the right hand part of the equation

$$2y_2\left(z\right)cos\left(2z\right)=\frac{1}{16}cos\left(4z\right)cos\left(2z\right),$$

(A47)

and use eq. (A43) to express it as

$$cos\left(4z\right)cos\left(2z\right)=\frac{1}{2}\left[cos\left(6z\right)+cos\left(2z\right)\right].$$

(A48)

Hence, eq. (A34) is cast into

$$\frac{d^2{y}_3}{d{z}^2}=\frac{1}{32}cos\left(6z\right)-\frac{7}{32}cos\left(2z\right)-a_3,$$

(A49)

which we integrate into

$$y_3\left(z\right)=-\frac{1}{2}{a}_3{z}^2+\frac{7}{128}cos\left(2z\right)-\frac{1}{1152}cos\left(6z\right)+c_{5}z+c_6.$$

(A50)

The equation above is periodic if and only if $a_3=0$ (and of course the coefficients $c_5,c_6=0$), which simplifies to

$$y_3\left(z\right)=\frac{7}{128}cos\left(2z\right)-\frac{1}{1152}cos\left(6z\right).$$

(A51)

The method can be further applied to any order based on the algorithm presented above, which leads to [181]

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

(A52)

and

$$c{e}_0\left(z,q\right)=1-\frac{q}{2}cos\left(2z\right)+\frac{q^2}{32}cos\left(4z\right)-\frac{q^3}{128}\left(7cos2z-\frac{1}{9}cos6z\right)+\dots .$$

(A53)

The same technique can be employed to identify approximations to all other Mathieu funcitons and eigenvalues, but taking into account that things are lightly dissimilar in case when $n/neq0$. The example below illustrates the expansions of $s{e}_3\left(z,q\right)$ and $b_3\left(q\right)$, for which $a_0=3$ and $y_0\left(z\right)=sin3z$. The equation for $y_1\left(z\right)$ appears as

$$\frac{d^2{y}_1}{d{z}^2}+9y_1=2sin3zcos2z-a_{1}sin3z$$

(A54)

One uses

$$sinacosb=\frac{1}{2}\left[sin\left(a+b\right)+sin\left(a-b\right)\right].$$

(A55)

Then eq. (A54) modifies accordingly

$$\frac{d^2{y}_1}{d{z}^2}+9y_1=sin5z+sinz-a_{1}sin3z,$$

(A56)

which exhibits a periodic solution provided that $a_1=0$. On the other hand, when $a_1\ne 0$ the solution would contain the non-periodic term proportional to $a_{1}zcos3z$. Therefore, one must choose $a_1=0$ to derive

$$y_1\left(z\right)=\frac{1}{8}sinz-\frac{1}{16}sin5z+c_{7}cos3z+c_{8}sin3z,$$

(A57)

where $c_7$ and $c_8$ represent two constants. As $s{e}_3\left(z,q\right)$ is an odd function and $cos3z$ is an even function, one infers $c_7=0$. What is more, it is redundant to include the $sin3z$ term, as this harmonic is already included in $y_0\left(z\right)$. Incorporation of this term is analogous with multiplying the solution with a constant, which means that we also choose $c_8=0$.

We now turn our attention to the equation that characterizes $y_2\left(z\right)$

$$\frac{d^2{y}_2}{d{z}^2}+9y_2=2y_1\left(z\right)cos2z-a_{2}sin3z$$

(A58)

which can be expressed as

$$\frac{d^2{y}_2}{d{z}^2}+9y_2=-\frac{1}{8}sinz+\left(\frac{1}{16}-a_2\right)sin3z-\frac{1}{16}sin7z,$$

(A59)

which exhibits a periodic solution when the parenthesis vanishes, namely $a_2=1/16$, which allows one to infer the solution as

$$y_2\left(z\right)=-\frac{1}{64}sinz+\frac{1}{640}sin7z.$$

(A60)

By further iterating one infers

$$b_3\left(q\right)=9+\frac{q^2}{16}-\frac{q^3}{64}+\frac{13}{20480}{q}^4+\frac{5}{16384}{q}^5-\frac{1961}{23592960}{q}^6+\frac{609}{104857600}{q}^7+\dots ,$$

(A61)

and

$$s{e}_3\left(z,q\right)\simeq =sin3z+q\left(\frac{1}{8}sinz-\frac{1}{16}sin5z\right)-q^2\left(\frac{1}{64}sinz-\frac{1}{640}sin7z\right)+\dots .$$

(A62)

The normalised function is

$$\begin{array}{c}s{e}_3\left(z,q\right)\simeq =sin3z+q\left(\frac{1}{8}sinz-\frac{1}{16}sin5z\right)\hfill \\ \hfill -q^2\left(\frac{1}{64}sinz+\frac{5}{312}sin3z-\frac{1}{640}sin7z\right)+\dots .\end{array}$$

(A63)

This approach can be employed for any specific value of n and to any order. However, for common n values it yields a series in $q^2$ where the coefficient of $q^{2r}$ contains the factor $(n-r)$ in the denominator. Accordingly, when using this series for a specific value of n it must be truncated at the $q^2(n-1)$ term, while the last few terms of this series are erroneous.