On the properties of a class of higher-order Mathieu equations originating from a parametric quantum oscillator

Evolution in time of an arbitrary initial state for a parametrically driven quantum oscillator is an interesting problem since there exist regions in parameter space (defined by the amplitude and frequency of the driving) where the moments of the probability distribution can diverge in time. While the first moment satisfies a Mathieu equation, the higher-order moments follow Mathieu like equations of order greater than two. It is not very often that a physical problem gives rise to higher-order Mathieu equations. Hence, we give a detailed study of the different stability zones associated with the parametric quantum oscillator, using perturbative techniques traditionally associated with the Mathieu equation. We verify our results by numerical analysis, thus demonstrating that for the higher-order Mathieu equations, the traditional perturbation theory methods give a consistent account of the stability zones.


Introduction
Mathieu equations describe parametrically driven oscillators. They have been investigated for decades using a variety of techniques like harmonic balance, Lindstedt-Poincare perturbation theory, multiple timescale techniques and Krylov-Bogoliubov amplitude and phase equations. An important feature of the Mathieu equation is that systems described by it show instability in a set of tongue-shaped regions in the parameter space spanned by the driving frequency and the amplitude of the drive. The tongues emanate from certain points on the frequency axis where the natural frequency is n 2 times the driving frequency with n = 1, 2, 3 . . . being the set of integers. The motion is strictly periodic only on the lines separating the stable from the unstable zones. In the stable zone, the motion is quasiperiodic in general. The techniques referred to above are used to study the boundaries where the motion is periodic. These techniques have been covered comprehensively in various texts and lecture notes (e.g. [1][2][3]). The study of coupled Mathieu equations became popular about three decades ago and can be found in the treatment of elastic pendula (two coupled degrees of freedom) [4], in the study of colliding particle beams in a two-dimensional accelerator [5] or in the description of two-mode dynamics of rotating rings with variable spin-speed [6]. The fact that the usual techniques for obtaining the boundaries of the single-variable Mathieu equation can be used to calculate the unstable zones of the coupled system was demonstrated clearly by Mahmoud [7].
More recently, a different perspective on this class of problems was discussed by Landa et al. [8]. While the usual emphasis had been on finding the instability zones and the periodic solutions on the boundary of the instability zones, these authors developed a technique for finding the classical as well as the quantum solutions in the region where the motion is bounded. One of the important problems addressed in this work is the one-dimensional quantum parametric oscillator which describes the motion of a single ion confined in a radio-frequency trap [9]. A comprehensive account of the existence of periodic orbits for a variety of generalised forms of the Mathieu equation was provided by Younesian et al. [10]. This included nonlinear Mathieu equations treated by harmonic balance and averaging [11], bifurcations in cubic Mathieu equations [12] and resonances in quasi-periodic Mathieu equations [13]. Another direction in which the Mathieu equation was extended was the introduction of time delay [14]. That the usual perturbative techniques associated with the Mathieu equation could be used for a nonlinear Mathieu equation with delay was established by Morrison and Rand [15]. Yet another extension of these driven systems has been the introduction of systems with fractional derivatives. The appropriateness of the usual perturbative techniques for exploring the regions of bounded motion has recently been established [16][17][18].
While the second-order (recently even fractional order) non-autonomous differential equations, which is the class to which the Mathieu equation belongs, have been studied from a formal mathematical standpoint as well as from the applicational viewpoint of perturbation theory, the higher-order equations have been, to the best of our knowledge, studied only from a formal perspective [19][20][21][22]. It turns out that the one-dimensional quantum parametric oscillator discussed above gives rise to an infinite set of higher-order Mathieu like equations. Our aim in this paper will be to explore the effectiveness of the usual perturbative techniques to study the unstable zones of these higher-order Mathieu-like equations. We will use known perturbative techniques to find the boundaries between the regions where the solution is bounded and where they are unbounded. It will turn out that the boundaries thus obtained will agree with the boundaries obtained numerically.
A combination of Heisenberg's equation of motion and Ehrenfest's theorem in quantum mechanics describes the time development of the expectation value of any operator O(x) under a given Hamiltonian H. If Ψ (x, t) is the wave function at any time t (we consider only one-dimensional space), then the expecta- which is a complex valued function generally written as Ψ (x, t) = √ P(x, t)e iφ(x,t) , where P(x, t) = |Ψ (x, t)| 2 is the probability distribution associated with the state at any time t. Ehrenfest's theorem has rarely been used except to study classical trajectories. However, it can provide useful information about various moments of the probability distribution and as such can yield information regarding the form of the distribution. A particularly interesting situation is the motion of a particle in a time-dependent harmonic trapping function represented by the potential Ω . The Hamiltonian operator for the quantum parametric oscillator is [23] In the above, p = −ih d dx is the momentum operator. In this case, for an initial wave function which is Gaussian, the time development of Ψ (x, t) can be found if the dynamics is in the part of ω − space where the classical motion is bounded [24][25][26][27]. For arbitrary initial wave functions, it is virtually impossible to arrive at the time development exactly. In this situation, Ehrenfest's theorem is very useful since it provides exact equations for all the moments. Each moment gives rise to a nonautonomous dynamical [28] system, and this very physical problem gives rise to the whole set of interesting mathematical systems which are worth studying in their own right. In this work, we will focus on the mathematical aspects of the dynamical systems arising from Ehrenfest's theorem for the moments.
The process of arriving at the equations of motion is straightforward. One needs to use Eq. (1) and carry out evaluation of [O, H] with the H of Eq. (2). As an example, we explicitly work out the first moment x which represents the mean position of the quantum particle. From Eq. (1) and the basic commutation relation [x, p] = ih, we see that and One move derivative of Eq. (3) and use of Eq. (4) leads to This is the dynamics of the mean position of the particle, and it follows the classical trajectory obtained from the Hamiltonian of Eq. (2). If f (t) = cos(Ωt), then we have the familiar Mathieu equation.
The dynamics of all moments x n have to be obtained by a process of repeated differentiation and use of Eq. (1). The moments x n are associated with different features of the probability distribution. In particular for n = 2, is the variance. This is the most important moment as its nonzero value signals a complete departure from classical physics. For n = 3, we have which is the skewness, and for n = 4 which is the kurtosis of the distribution. We write down the equations of motion (the derivation of the dynamics of V is shown in "Appendix A", the others follow in a similar manner) For the discussions in this paper, we will focus entirely on f (t) = cos(Ωt).
The properties of Mathieu equation (Eq. 6) are well known [1][2][3]. We will recall them in Sect. 2 to carefully explain the perturbation theory technique that we will use and the numerical procedure that we will follow in the rest of the paper.
If we consider the dynamics of the variance V and redefine 4ω 2 asω 2 , then Eq. (10) becomes for f (t) = cos(Ωt) We further simplify by redefining the combination ω 2 as¯ , so that The above equation corresponds to the dynamical system with the structure The trace of the matrix A i j (t) is clearly zero, and hence, the three Floquet multipliers μ 1 , μ 2 and μ 3 satisfy μ 1 μ 2 μ 3 = 1. Consequently, the periodic solutions can have the periods T, 2T and 3T . For a subharmonic response of period 3T (frequency Ω 3 ), we would have periodic solutions originating fromω = Ω 3 , 4Ω 3 . . . , etc. For a periodic response of period 2T , the solutions start atω = Ω 2 , 3Ω 2 . . . , etc and for a response of period T , the solutions start atω = Ω, 2Ω . . . , etc.
In a similar fashion, the equations of motion of S and K are in principle capable of showing periodicities T, 2T, 3T, 4T and T, 2T, 3T, 4T, 5T , respectively.
Using different perturbation techniques [2], we will study the properties of the dynamical systems shown in Eqs. (10)- (12). The techniques (harmonic balance and Krylov-Bogoliubov) will be geared to finding the periodic trajectories. In the case of Eq. (13) specifically, we obtain the following: For Eq. (11), the conclusions are qualitatively similar with the unstable zone emanating fromω = Ω alone. For Eq. (12), however, we have instability zones emanating fromω = Ω 4 , Ω 2 and Ω. In this case, the instability zone emanating fromω = Ω is far wider than those for Eqs. (6), (10) and (11).
We briefly recapitulate the primary results of the Mathieu equation in Sect. 2. In Sect. 3, we deal with the dynamics of the variance in detail, and in Sect. 4, we study the skewness and the kurtosis. We conclude with a short summary in Sect. 5.

Properties of Mathieu equation
In this section, we review the properties of the Mathieu equation which is a second-order differential equation It has the dynamical system formẊ = Y andẎ = −(ω 2 + cos(Ωt)X ). Since the dynamical system is traceless, the two Floquet multipliers μ 1,2 satisfy μ 1 μ 2 = 1. This means that the system can have responses of period T and 2T . Since T = Ω 2π , harmonic balance indicates that period 2T solutions should emerge from ω = Ω 2 and period T solution should emerge from ω = Ω. We first explore the straightforward perturbative method for obtaining periodic solutions of period 2T , around ω = Ω 2 . Accordingly, we write and expand Substituting into Eq. (16) and equating the coefficient of each power of to zero, we havë and inserting this solution in Eq. (19b), For X 1 to be finite, we cannot have cos( Ωt 2 ) or sin( Ωt 2 ) terms on the right-hand side. These terms are known as secular terms and have to be removed order by order in perturbation theory. This requires B = 0 and δ 1 = − 1 2Ω if the solutions are of the cos( Ωt 2 ) variety and A = 0 and δ 1 = + 1 2Ω if the solutions are of the form sin( Ωt 2 ). Thus, the curves on which the solution is periodic are given by That the solution is divergent between the two lines of periodicity above and quasi-periodic outside is shown by the technique of Krylov-Bogoliubov. It proceeds by writing the solution X (t) as near the point ω = Ω 2 . We will assume A(t) and B(t) to be slowly varying (so thatȦ <<Ä and similarly for B) and find the dynamical systems for A and B. The fixed point will give the periodic solutions. Taking derivatives and keeping the relevant termṡ With ω 2 = Ω 2 4 + Ωδ to the lowest order in δ, we insert Eqs. (22) and (24) in Eq. (16) Setting the coefficients of cos( Ωt 2 ) and sin( Ωt 2 ) to zero, we geṫ leading to the solution The answer reveals that the slow variation of A(t) comes from the fact the coefficient of t in the exponential is small. It also shows that if δ 2 > 2 4Ω 2 , the solutions are periodic with frequency δ 2 − 2 4Ω 2 and hence the solution for X (t) [as seen from Eq. (22)] is quasi-periodic. For δ = ± 2Ω , the solution X (t) is periodic as A(t) and B(t) become constants. This agree exactly with the curves we found from the perturbation theory for periodic solutions. For δ < | 2Ω | (the region bounded by the curves where the solution is periodic), X (t) diverges. Hence, we have a situation where quasi-periodic solutions are separated from diverging solutions by curves on which the solution is periodic. Searching for periodic solutions is a good strategy for finding zones where the solutions may diverge. It should be noted that the Krylov-Bogoliubov technique is extremely versatile as it gives a complete picture of what is going on, in the vicinity of the periodic solution.
It should be noted that the multiple timescale technique (see Ref. [29]) yields identical answers. The numerically obtained picture of solutions of the Mathieu equation are shown in Fig. 1.
To solve the higher-order differential equations, we used the Runge-Kutta (RK) recursion formula. Higherorder differential equations can be represented by a set of first-order differential equations. Runge-Kutta algorithm has the form y i+1 = y i +hφ(x i , y i , h), to numerically solve the first-order differential equation, where φ is an increment function of the approximation of f (x, y) in the interval of x i ≤ x ≤ x i+1 . We have used fourth-order RK formulas, where φ is the weighted average of four derivative evaluations k 1 , k 2 , k 3 and k 4 on the interval x i ≤ x ≤ x i+1 , that is, where, We have used a step size of h = 0.01. Smaller values of h do not change the computed boundaries.
With this introduction to the techniques to be used, we now turn to a fairly elaborate discussion of the thirdorder equation Eq. (14) satisfied by V (t). The derivation of the equation is given in "Appendix A".

The properties of the third-order non-autonomous equation for the variance
The general discussions of Eqs. (10-12) from the mathematical point of view can be found in Refs. [19][20][21][22]. We begin by noting the fact that Eq. (13) can have periodic motion with periods 3T, 2T and T and by asking for the region in (ω − ) space where the solution of period 3T should be obtained. This has to be the primary resonance nearω = Ω 3 and we writē Perturbation theory is what we want to begin with and write Substituting the expansion in Eq. (13), we have The analysis of the above equation up to O( 2 ) is presented in "Appendix B". We find that there is no zone of instability to O( 2 ) and the periodic orbits appear on the curveω = − 27 2 128Ω 3 . We now focus on the orbits of period 2T . A harmonic balance argument suggests that these will exist in the vicinity ofω = Ω 2 . We writeω = Ω 2 + δ, where δ is small and try a Krylov-Bogoliubov technique with the trial solution We provide the expressions forV ,V and ... V in "Appendix B".
From the coefficients of sin( Ωt 2 ), cos( Ωt 2 ), sin( 3Ωt 2 ) and cos( 3Ωt 2 ) in Eq. (B.10) of "Appendix B", The eigenvalues λ of this set of equations are obtained from the condition which gives the following quadratic equation in λ 2 : The two roots for λ 2 are obtained as We note that neither of the roots can be positive, and hence, there is no instability zone aroundω = Ω 2 . The periodic orbit is obtained for λ = 0 which corresponds to δ = − 2 6Ω 3 . This is interesting since the anticipated instability zone is not there. It can be traced to an accidental cancellation of terms in cos( Ωt 2 ) and sin( Ωt 2 ) coming from the third and fourth terms of Eq. (14).
We now try to find out the solutions of period T which are expected to occur nearω = Ω. Settingω = Ω + δ, we carry out the Krylov-Bogoliubov method with the ansatz, where A 0 , A 1 , B 1 , A 2 , B 2 are all slowly varying function of time. Our calculation technique is exactly the same as detailed in the derivation of Eqs. (29)-(32d) and finally yieldṡ The eigenvalue λ satisfies a fifth-order equation which always has a zero eigenvalue. Of the four remaining eigenvalues, two are purely imaginary and the other two have the form ± 2 Ω 2 16 − 4δ 2 , and hence, there is a zone of divergent solution bounded by δ = ± 8 exactly as for the Mathieu equation near ω = Ω 2 . We have solved Eq. (14) numerically using Runge-Kutta methods to obtained the instability zones as shown in Fig. 2. The instability zones agree exactly with those found in Fig. 1 for small values of . As seen from our calculation, the first instability zone appears atω = 2ω = Ω and its width is exactly the same as the instability zone in Fig. 1 around ω = Ω 2 . This explains why the recent numerical work of Hashemloo et al. [27] finds that the width increases indefinitely only in the parameters zones where the mean position increases. We end this section by noting that the absence of an instability zone aroundω = Ω 3 (which corresponds to ω = Ω 6 ) is presumably due to the fact that in the absence of the periodic forcing, the system exhibits a zero mode. The absence of a zone around ω = Ω 4 is accidental because of the coefficients of f andḟ in Eq. (10) differing exactly by a factor of two.

The properties of the fourth-and fifth-order equations for the skewness and kurtosis
We begin our discussion with the skewness S, rewriting Eq. (11) to O( ). This leads to with f = cos(Ωt). Writing this as a fourth-order dynamical systemẊ i = A i j X j with periodic coefficients, we find that M i j is traceless and hence the solutions X (t) can have periodicities 4T, 3T, 2T and T . There are solutions with frequency Ω 4 , Ω 3 , Ω 2 and Ω. The difference between Eq. (37) and the system studied in Sects. 2 and 3 is that the unperturbed solution S 0 of Eq. (37) [i.e. solution for = 0] is a two-frequency solution (this is a significant difference that can occur with the higher-order Mathieu equation) S 0 (t) = A 1 cos(ωt) + B 1 sin(ωt) + A 2 cos(3ωt) +B 2 sin(3ωt). The expected resonant response can happen when either ω or 3ω matches Ω 4 , Ω 3 , Ω 2 and Ω. Harmonic balance indicates which matching would lead to dominant effect and inspection shows that the conditions ω = Ω 6 or Ω 2 are the ones capable of producing response at O( ). The details of the calculations are presented in "Appendix C". Here we simply report that there is no divergent region near ω = Ω 6 , while in the vicinity of ω = Ω 2 , we have the boundaries given by the curves shown in Eq. (C.8a).
We show the instability zones obtained from our computational work and the boundaries obtained from perturbation theory in Fig. 3.
We now turn to the dynamics of the kurtosis K Eq. (12) and writing it as a dynamical system find that once again the matrix A is traceless and the possible time periods of the solution can be 5T, 4T, .., T with frequencies Ω 5 , Ω 4 , .., Ω. Noting that the = 0 system has a solution with frequencies 2ω and 4ω, we find from a visual inspection of harmonic balance that at O( ), the possibilities are the matches at ω = Ω 8 and ω = Ω 4 and ω = Ω 2 . We rewrite Eq. (12) in a more streamlined form by writing 2ω asω and ω 2 as¯ . With there new redefinition, we have (keeping terms to O(¯ ) only as our calculations will be done to the order) .....
The relevant resonance-inducing terms on the r.h.s. of Eq. (41) are of the two varieties: one set proportional to A 1 cos( Ω 4 ) and B 1 sin( Ω 4 ) which have δ 1 as a prefactor and the other set proportional to A 2 cos( Ω 2 ) and B 2 sin( Ω 2 ) where coefficients contain terms with δ 1 and also without. Setting δ 1 = 0 to remove resonance is not an option. Hence, the resonance is removed by setting A 1 = B 1 = 0 and This produces a zone starting at ω = Ω 8 where the solutions diverge in time. On the boundaries of this zone, the solution is periodic with a period 2T .
An identical calculation holds near ω = Ω 4 , but now the width of the zone is exactly double that the zone around ω = Ω 8 . We finally consider the case ω = Ω 2 +δ. In this case, as for V (t), the perturbative technique used above does not work. With ω = Ω 2 + δ, we write Eq. (12) to the lowest order in δ and as We use the Krylov-Bogoliubov method with Straightforward algebra now leads tȯ As expected there is a zero mode, as there was for the variance. The other eigenvalues are λ = ±i 60 137 + O( ) and the remaining four obtained from the solution of λ 4 + λ 2 20δ 2 − 9 2 32 + 64δ 4 + 153 1024 4 = 0 | which is very close to 8 . However, this is not the complete story since this is true only when λ 2 has real roots. If λ has complex roots with λ 2 having the structure λ 2 = −α ± iβ, where α is positive being in the stable zone δ > | 3 8 √ 10 | then λ = 4 α 2 + β 2 e ±iθ/2 where cos(θ ) = α √ α 2 +β 2 and sin(θ ) = β √ α 2 +β 2 . A negative cos(θ ) puts θ in the range π 2 to 3π 2 and cos( θ 2 ) is in the range π 4 to 3π 4 . This implies cos( θ 2 ) is positive for π 2 < θ < 3π 4 and this is what magnifies the instability zone around ω = Ω 2 , way beyond what happens for x , V and S. As can be seen from Fig. 4, our perturbative calculations are in accordance with the numerical data. It should be seen that at this order of Mathieu equation, the main feature of instabilities arising at lower frequencies than the second-order equation is clearly borne out as the accidental cancellation no longer happen. However, in this case as increases, the region in parameter space where the kurtosis diverges becomes very large as seen from our numerical work and this makes the long time quantum dynamics in such systems of particular interest. The long time shape of the probability distribution needs to have a power law fall off in the regions where the kurtosis diverges.

Conclusions
Linear differential equations with periodic coefficients and of order greater than two are a rarity physical sciences. The quantum parametric oscillator provides a real physical system [24][25][26][27] where such equations make a natural appearance when the dynamics of higher moments is considered. Since description of the quantum state is in terms of probabilities one studies the dynamics of the mean value x instead of the position of x itself. But the mean is not sufficient to picture the distribution and hence moments like the variance are essential. The fact that the moments follow Mathieu equations of higher order imply that while the mean may remain bounded, the fluctuation around the mean or the spread parameter kurtosis may diverge with time leading a drastic change in shape of the probability distribution. Hence, the quantum parametric oscillator provides the physical justification for investigating higher-order Mathieu equations.
The discussion on higher-order Mathieu-like equations have till now been restricted to mathematical analysis. We wanted to explore the practical issue of whether the perturbation techniques that work so well for finding the instability zones of the Mathieu equations would be equally useful for higher-order Mathieu equation. As seen from Figs. 2, 3 and 4, that is indeed the case. It will be noted that we never found any instability zone around lowest admissible frequency ω = Ω 6 in our study of the variance in Sect. 3. This we attribute to the existence of a zero mode in the system. The fact that there was no instability zone around ω = Ω 4 was due to an accidental cancellation. This did not occur for the kurtosis and instability zones can be seen emanating from ω = Ω 4 and ω = Ω 8 in Fig. 4. From a purely technical point of view, the fact that traditional technique like harmonic balance and Krylov-Bogoliubov provides results which could be verified in their zones of validity by numerical analysis indicates the robustness of these methods and should lead to interesting future investigations. ...
... For the solution V 1 to exist, we cannot have resonating terms on the right-hand side and hence δ 1 = 0. At this order ...
A 0 cos(Ωt) + 3 16Ω 2 A 1 cos 4Ωt 3 We need to pick out the resonance inducing secular terms from the right-hand side of Eq. (B.6) and set them equal to zero. This gives 2Ω 2 9 and hence we obtain This indicates that the period 3T orbit can exist only along the lineω = − 27 We now turn to Eq. (31) of the text which is the input for the Krylov-Bogoliubov technique for study of Eq. (14) near the frequencyω = Ω 2 . We consequently need the expression forV (t),V (t) and ... V (t) keeping in mind that A 1 (t), B 1 (t), A 2 (t) and B 2 (t) of Eq. (31) are slowly varying quantities, hence the above mentioned derivatives of V (t) will not contain any time derivatives of A 1 , B 1 , A 2 and B 2 higher than the second. With this, we geṫ The periodic solution S(t) = A 1 cos( Ω 2 t) + A 2 cos( 3Ω 2 t) occurs on the line obtained from consistency of 3 4 A 2 = ( 2δ 1 Ω + 1 2 )A 1 and δ 1 Ω A 2 + A 1 24 = 0 which leads to δ 1 = − Ω 8 and the solution S(t) = B 1 sin( Ω 2 t) + B 2 sin( 3Ω 2 t) exists on the line obtained from consistency of 3 4 B 2 = (− 2δ 1 Ω + 1 2 )B 1 and δ 1 Ω B 2 − B 1 24 = 0, which leads to δ 1 = Ω 8 . Around ω = Ω/2, we have the periodic boundaries at Fig. 3 ω These boundaries are shown blue shaded region in Fig. 3.