2. The exact adiabatic invariant

Consider the non-relativistic motion of a particle of mass $m$ and charge $q = |q| \sigma$ in a uniform, time-dependent magnetic field of the form $\boldsymbol{B}(t) = B_{0} B(t/T) \hat {z}$ . Here, $B(t/T)$ is a dimensionless function, $\sigma = \pm 1$ is the sign of the particle charge, $B_{0}$ is the typical magnetic field strength and $T$ is the characteristic time scale of magnetic field variation. Since we are only interested in non-vanishing magnetic fields, we assume without loss of generality that $B$ is strictly positive. We define the dimensionless ordering parameter $\epsilon ^{-1} = \omega _{c}T$ , where $\omega _{c} = {|q|B_{0}}/{m}$ is the typical cyclotron frequency; $\epsilon$ is assumed to be a small parameter representing the ratio between the typical gyro-period and the time scale of magnetic field variation. We adopt units of time such that $T = 1$ . The single-particle equations of motion on the fast time scale $\tau = \epsilon ^{-1}t$ are

(2.1) \begin{equation} \frac {\text{d}^2 \boldsymbol{x}}{\text{d}\tau ^2} = \sigma \left [B(\epsilon \tau ) \frac {\text{d}\boldsymbol{x}}{\text{d}\tau } \times \hat {z} + \epsilon \boldsymbol{E}(\boldsymbol{x},\epsilon \tau ) \right ]. \end{equation}

We assume $\boldsymbol{B}$ is generated by a rotationally symmetric configuration, so the induced electric field is given in cylindrical coordinates by $\boldsymbol{E}(\boldsymbol{x},t) = -{1}/{2} r B'(t) \hat {\theta }$ . Since $z'' = 0$ , we ignore the trivial motion in the $z$ direction and focus only on the transverse particle motion.

To study invariants of this time-dependent system, we must place the equations of motion in autonomous Hamiltonian form. Let $(\tau , p_{\tau })$ be a canonically conjugate position–momentum pair obeying $\dot {\tau } = 1$ . Using polar coordinates, we also define the canonical momenta $p_{r} = {\text{d} r}/{\text{d}\tau }$ and $p_{\theta } = r^2 ({\text{d}\theta }/{\text{d}\tau }) + ( {r^2 \sigma B(\epsilon \tau )}/{2})$ . The Hamiltonian in the extended phase space can be taken to be

(2.2) \begin{equation} \tilde {H}_{\epsilon } = \frac {\epsilon }{2}\left (p_{r}^2 + \left ( \frac {p_{\theta }}{r} - r \sigma \varOmega (t) \right )^2 \right ) + \epsilon p_{\tau }. \end{equation}

The symplectic form is $-\text{d}\vartheta _\epsilon$ with the Liouville one-form $\vartheta _{\epsilon }$ given by

(2.3) \begin{equation} \vartheta _{\epsilon } = \epsilon p_{r}\text{d}r + \epsilon p_{\theta }\text{d}\theta +p_{\tau } \text{d}t . \end{equation}

The equations of motion are

(2.4) \begin{eqnarray} \dot {p_{\theta }} &\ =\ & 0, \nonumber \\ \dot {\theta } &\ =\ & \left (\frac {p_{\theta }}{r^2} - \sigma \varOmega (t) \right)\!, \nonumber \\ \dot {p_{\tau }} &\ =\ & \epsilon \sigma \varOmega '(t) \left (p_{\theta } - r^2 \sigma \varOmega (t) \right)\!, \nonumber \\ \dot {t} &\ =\ & \epsilon, \nonumber \\ \dot {p_{r}} &\ =\ & \frac {p_{\theta }^2}{r^3} - \varOmega ^2(t) r, \nonumber \\ \dot {r} &\ =\ & p_{r}. \end{eqnarray}

In the variables $(r, \theta , t, p_{r}, p_{\theta }, p_{\tau })$ , (2.4) takes the desired form of a nearly periodic Hamiltonian system. As shown in the seminal work of Kruskal (Reference Kruskal1962), nearly periodic Hamiltonian systems encompass a broad class of systems for which a perturbative adiabatic invariant $\mu _{\epsilon } = \mu _{0} + \epsilon \mu _{1} + \cdots$ can be constructed, and therefore provide a general context for defining the notion of a non-perturbative adiabatic invariant. For the sake of completeness, we recall the generalised definition of a nearly periodic system as defined by Burby and others in Burby & Squire (Reference Hyouguchi, Adachi and Ueda2020), Burby & Hirvijoki (Reference Burby, Hirvijoki and Leok2021) and Burby et al. (Reference Burby, Hirvijoki and Leok2023).

Critical to the construction of an adiabatic invariant $\mu _{\epsilon }$ in nearly periodic Hamiltonian systems is the roto-rate vector $R_{\epsilon } = R_{0} + \epsilon R_{1} + \cdots$ which is the unique formal vector field satisfying, in the sense of formal power series in $\epsilon$ .

1. (i) $\lim\nolimits _{\epsilon \to 0} R_{\epsilon } = R_{0}$ .

2. (ii) $\exp (2 \pi \mathcal{L}_{R_{\epsilon }}) = \text{Id}$ .

3. (iii) $[X_{\epsilon }, R_{\epsilon }] =0$ .

It is not hard to show that the conditions imposed on the roto-rate vector ensure $R_{\epsilon }$ is a formal Noether symmetry. Essentially, $R_{\epsilon }$ generates a formally ignorable angle-like coordinate, often the gyrophase-angle in the context of guiding-centre theory. The asymptotic invariant $\mu _{\epsilon } = \mu _{0} + \epsilon \mu _{1} + \cdots$ is then the associated formally conserved quantity obeying $i_{R_{\epsilon }} \text{d}\vartheta _{\epsilon } = -\text{d}\mu _{\epsilon }$ . As argued in Gjaja & Bhattacharjee (Reference Hyouguchi, Adachi and Ueda1992), even under ideal conditions, the series for $\mu _{\epsilon }$ may not converge for any $\epsilon \gt 0$ . In general, there is no coordinate change separating, say, the fast motion of gyrating particles from the slow drift motion. The results of Neishtadt (Reference Neishtadt1984), Ramis & Schäfke (Reference Ramis and Schäfke1996) and Simó (Reference Simó1994) show that such a split may only be possible up to an exponentially small coupling between the fast and slow variables.

While it is too much to ask that the formal series for $R_{\epsilon }$ converges, it may nevertheless be true that a $U(1)$ symmetry persists for $\epsilon \gt 0$ . In such a case, we say the symmetry generator $\mathcal{R}_{\epsilon }$ is a non-perturbative roto-rate vector if $\mathcal{R}_{\epsilon }$ depends continuously on $0\leqslant \epsilon \ll 1$ and admits an asymptotic power series as $\epsilon \to 0^{+}$ . Uniqueness of the roto rate then implies $\mathcal{R}_{\epsilon }$ is asymptotic at all orders to $R_{\epsilon }$ . When a non-perturbative roto-rate vector exists in a nearly periodic Hamiltonian system, the Hamiltonian $\mathcal{J}_{\epsilon }$ generating $\mathcal{R}_{\epsilon }$ is a non-perturbative adiabatic invariant. Indeed, since $\mathcal{R}_{\epsilon }$ is asymptotic to $R_{\epsilon }$ , it follows that $\mathcal{J}_{\epsilon }$ is asymptotic to Kruskal’s adiabatic invariant $\mu _{\epsilon }$ . This leads to the following definition.

In the present case, $\tilde {H}_{\epsilon }$ and $\vartheta _{\epsilon }$ are given by (2.2) and (2.3) respectively, and $X_{\epsilon }$ is the Hamiltonian vector field (2.4). Taking the formal limit $\epsilon \to 0^{+}$ , the limiting dynamics $X_{0}$ consists of periodic motion at an angular frequency of $2\varOmega (t)$ . We therefore have a nearly periodic Hamiltonian system with $\omega _{0}(t) = 2 \varOmega (t)$ and $R_{0} = ( {X_{0}}/{2 \varOmega (t)})$ .

In Qin & Davidson (Reference Qin and Davidson2006), it was shown that

(2.5) \begin{equation} 2\epsilon ^{-1} \mathcal{M}_{\epsilon } = \frac { p_{\theta }^2}{2} \left (\frac {w(t; \epsilon )^2}{r^2} \right ) + \frac {1}{2} \left (\frac {r^2}{w(t;\epsilon )^2} \right ) + \frac {1}{2}(p_{r}w(t;\epsilon ) - \epsilon w'(t;\epsilon )r)^2 - p_{\theta }\sigma , \end{equation}

is an exact invariant of the dynamics when $w(t;\epsilon )$ solves the singularly perturbed Ermakov–Pinney equation

(2.6) \begin{equation} \epsilon ^2 w''(t;\epsilon ) + \varOmega ^2(t) w(t;\epsilon ) = \frac {1}{w^3(t;\epsilon )} . \end{equation}

Taking $\epsilon$ to be zero, one computes that $w(t;0) = \varOmega (t)^{-1/2}$ . Thus, $\epsilon ^{-1}\mathcal{M}_{\epsilon }$ formally limits to magnetic moment $\epsilon ^{-1} \mathcal{M}_{0} = ( {1}/{4 \varOmega })( r^2 \dot {\theta }^2 + p_{r}^2) = \mu _{1}$ , and the Hamiltonian vector field of $\mathcal{M}_{\epsilon }$ limits to $R_{0}$ . To prove $\mathcal{M}_{\epsilon }$ is a non-perturbative adiabatic invariant, it must be shown that $\mathcal{M}_{\epsilon }$ generates a $2 \pi$ -periodic $U(1)$ action, and that there exists a branch of solutions $w(t; \epsilon )$ to the Ermakov–Pinney equation such that both $w(t;\epsilon )$ and $w'(t;\epsilon )$ admit asymptotic power series in $\epsilon$ . Once both of these statements have been shown, one verifies all the conditions of Definition 2.2 are met.

We first assume $w(t;\epsilon )$ is given and show that $\mathcal{M}_{\epsilon }$ generates a $U(1)$ action. In the next section, we show a method by which $w(t;\epsilon )$ may be smoothly constructed. The Hamiltonian vector field of $2\mathcal{M}_{\epsilon }$ generates the flow

(2.7) \begin{eqnarray} \frac {\text{d} t}{\text{d}\xi } & \ =\ & 0, \nonumber \\\frac {\text{d} p_{\tau }}{\text{d}\xi } & \ =\ & \epsilon \left [-ww'\left (\frac {p_{\theta }^2}{r^2} + p_{r}^2 \right ) + r^2 w'\left (\frac {1}{w^3} - \epsilon ^2 w''\right ) + rp_{r}\epsilon \left ({w'}'w + {w'}^2\right )\right ], \nonumber \\ \frac {\text{d} r}{\text{d}\xi } & \ =\ & p_{r}w^2 - \epsilon w w' r, \nonumber \\ \frac {\text{d} p_{r}}{\text{d}\xi } & \ =\ & \frac {p_{\theta }^2 w^2}{r^3} - \frac {r}{w^2} + \epsilon w' \left (p_{r} w - \epsilon w' r \right ), \nonumber \\ \frac {\text{d}\theta }{\text{d}\xi } & \ =\ & \frac {p_{\theta }w^2}{r^2} - \sigma , \nonumber \\ \frac {\text{d} p_{\theta }}{\text{d}\xi } & \ =\ & 0. \end{eqnarray}

The equation for the canonical pair $(r,p_{r})$ closes to give

(2.8) \begin{equation} r'' + r = \frac {p_{\theta }^2 w^4}{r^3}, \end{equation}

which has the general solution $r^2(\xi ; E, \phi ) = p_{\theta }w^2 [E + \sqrt {E^2 - 1}\sin (2\xi + \phi )]$ . Given that $p_{r}(\xi ) ={1}/{w^2}({\text{d} r}/{\text{d}\xi }(\xi ) + \epsilon ww'r(\xi ))$ , $p_{r}(\xi )$ is a $\pi$ -periodic function for any constants $E \geqslant 1, \phi$ . One computes that

(2.9) \begin{align} &\int _{0}^{\pi } \frac {\text{d}\xi }{r^2(\xi )} = \frac {\pi }{p_{\theta } w^2}, \nonumber \\ &\int _{0}^{\pi } r^2(\xi ) \text{d}\xi = (\pi p_{\theta } E) w^2, \nonumber \\ &\int _{0}^{\pi } \left [ \frac {p_{\theta }^2}{r^2(\xi )} + (p_{r}(\xi ))^2\right ]\text{d}\xi = (\pi p_{\theta } E) \left ( \frac {1}{w^2}+\epsilon ^2 {w'}^2\right)\!, \nonumber \\ &\int _{0}^{\pi } r(\xi ) p_{r}(\xi ) \text{d}\xi = \epsilon (\pi p_{\theta }E) ww'. \end{align}

It follows that $\epsilon ^{-1}\mathcal{M}_{\epsilon }$ generates a $2\pi$ periodic $U(1)$ action since

(2.10) \begin{align} \Delta \theta & \equiv \int _{0}^{\pi } \frac {\text{d}\theta }{\text{d}\xi }(\xi ) \text{d}\xi = \pi - \sigma \pi \equiv 0 \, (\text{mod } 2\pi ), \nonumber \\ \Delta p_{\tau } & \equiv \int _{0}^{\pi } \frac {\text{d} p_{\tau }}{\text{d}\xi }(\xi ) \text{d}\xi = 0. \end{align}

We have thus shown that the non-perturbative invariant $\mathcal{M}_{\epsilon }$ generates an angle-like, time-dependent ignorable coordinate. The physical meaning of $\mathcal{M}_{\epsilon }$ is somewhat obscure for $\epsilon \neq 0$ beyond the fact that $\mathcal{M}_{\epsilon }$ is linear combination of $p_{\theta }$ and the approximate constants of motion used in accelerator physics, as shown Qin & Davidson (Reference Qin and Davidson2006).

3. Construction of the auxiliary function

We now construct a smooth function $w(t;\epsilon )$ satisfying the singularly perturbed Ermakov–Pinney equation (2.6) admitting an asymptotic power series in $\epsilon$ limiting to $\varOmega (t)^{-1/2}$ . We show that $w'(t;\epsilon )$ also admits an asymptotic power series obtained by term-by-term differentiation of the series for $w(t;\epsilon )$ . Provided this choice of $w(t;\epsilon )$ is used in the construction of the exact invariant, $\mathcal{M}_{\epsilon }$ admits an asymptotic series and is hence a non-perturbative adiabatic invariant.

It is known by the general results of Sibuya (Reference Sibuya2000, Reference Sibuya2003), and Canalis-Durand et al. (Reference Hyouguchi, Adachi and Ueda2000) that solutions to singularly perturbed analytic ODEs asymptotic to their formal series solution exist for short times. To obtain all-time results for merely smooth $\varOmega (t)$ , use the well-known fact that the Ermakov–Pinney equation is controlled by the singularly perturbed Hill’s equation (Pinney Reference Pinney1950).

(3.1) \begin{equation} \epsilon ^2 \psi ''(t;\epsilon ) + \varOmega ^2(t) \psi (t;\epsilon ) = 0. \end{equation}

Indeed, let $u(t),v(t)$ be any two linearly independent solutions to Hill’s equation. Then there exists smooth, real functions $R(t) \gt 0$ and $\phi (t)\gt 0$ such that

(3.2) \begin{equation} \psi \equiv u(t) + i v(t) = R(t) \exp \left (\mp \frac {i}{\epsilon } \int ^{t}_{0} \phi (t')\text{d}t'\right)\!, \end{equation}

is a never-vanishing solution to Hill’s equation; $R(t)$ can never vanish since that would imply the Wronskian matrix formed from $u(t)$ and $v(t)$ is singular at some point in time, contradicting the linear independence of $u$ and $v$ . The choice of the $\mp$ sign is determined by the condition that $\phi \gt 0$ . $R$ and $\phi$ define a solution to Hill’s equation if and only if

(3.3) \begin{align} \epsilon ^2R'' + \varOmega ^2(t)R &= \phi ^{2}R, \nonumber \\ 2\phi R' + \phi ' R &= 0. \end{align}

The latter equation implies $R = C \phi ^{-1/2}$ , justifying that $\phi$ is also non-vanishing. Provided the initial conditions are such that $C = 1$ , which can always be ensured by rescaling $u$ and $v$ , we have that

(3.4) \begin{equation} \epsilon ^2 R'' + \varOmega ^2(t) R = \frac {1}{R^3}, \end{equation}

which is the desired result. With this basic fact in mind, it is clear that the formal expansion of $w(t;\epsilon )$ into powers of $\epsilon$ is obtained from the usual WKB series of (3.1). Indeed, the first-order WKB approximation yields that $w(t;\epsilon \ll 1) \approx \varOmega ^{-1/2}$ . Results about the accuracy of WKB (Olver Reference Olver1961, Reference Olver1974; Taylor Reference Taylor1982; Costin Reference Hyouguchi, Adachi and Ueda2008) therefore justify that the formal solutions $w^{(n)}(t) = w_{0}(t) + \epsilon ^{2}w_{2}(t) + \cdots + \epsilon ^{n} w_{n}(t)$ to (2.6) approximate some true solution for any fixed $n$ .

To construct an exact solution $\psi (t;\epsilon )$ to Hill’s equation asymptotic to the WKB series at all orders, we use the Bremmer series. Constructed independently in Bremmer (Reference Bremmer1951), Schelkunoff (Reference Schelkunoff1951) and Landauer (Reference Landauer1951) and studied further in Bellman & Kalaba (Reference Bellman and Kalaba1959), Atkinson (Reference Atkinson1960) and Kay (Reference Kay1961), the Bremmer series organises the solution of (3.1) into a convergent sum of waves undergoing a particular number of reflections. Importantly, the Bremmer series agrees with WKB at leading order. In Winitzki (Reference Winitzki2005), it was shown that the Bremmer-series solutions admit an asymptotic $\epsilon$ -series provided there exists a point in time such that all derivatives of $\varOmega (t)$ vanish. By temporarily introducing a cutoff to $\varOmega (t)$ in the far past, we show that this assumption can be dropped. We repeat parts of Winitzski for the sake of a self-contained construction, but refer the reader Winitzki (Reference Winitzki2005) for an interesting discussion on the precision of optimally truncated WKB. We first write $\varOmega (t) = \varOmega _{0} + \phi (t)$ where $\phi (t)$ is smooth and the strict inequality $|\phi (t)| \lt \varOmega _{0}$ holds uniformly. Let $0 \leqslant H(t) \leqslant 1$ be a smooth bump function supported on a domain of the form $[-K,\infty )$ such that $H(t\geqslant 0) \equiv 1$ . We define the deformed magnetic field $\tilde {\varOmega }(t) = \varOmega _{0} + H(t) \phi (t)$ . Clearly, for $t \geqslant 0$ , we have not changed $\varOmega (t)$ so any linearly independent solutions to the deformed equation

(3.5) \begin{equation} \epsilon ^2\tilde {\psi }''(t;\epsilon ) + \tilde {\varOmega }^2(t) \tilde {\psi }(t;\epsilon ) = 0, \end{equation}

will yield a solution to the undeformed Hill’s equation, hence the undeformed Ermakov–Pinney equation, for $t \geqslant 0$ . By deforming the magnetic field, we are able to assume $\tilde {\psi }$ is a positive-frequency WKB solution in the far past, i.e.

(3.6) \begin{equation} \tilde {\psi }(t \to -\infty ;\epsilon ) = \tilde {\varOmega }^{-1/2}(t)\exp \left (-\frac {i}{\epsilon } \int _{0}^{t} \tilde {\varOmega }(t')\text{d}t'\right)\!. \end{equation}

Making this notion more precise, we define the WKB solutions

(3.7) \begin{equation} X_{\pm }(t;\epsilon ) = \tilde {\varOmega }^{-1/2}(t)\exp \left (\mp \frac {i}{\epsilon } \int ^{t}_{0} \tilde {\varOmega }(t')\text{d}t'\right)\!. \end{equation}

Then there exists complex-valued functions $p(t;\epsilon )$ and $q(t;\epsilon )$ such that

(3.8) \begin{equation} \tilde {\psi }(t;\epsilon ) = p(t;\epsilon )X_{+}(t;\epsilon ) + q(t;\epsilon )X_{-}(t;\epsilon ). \end{equation}

The condition that $\tilde {\psi }$ is a WKB solution in the far past implies initial conditions $p(-\infty ;\epsilon ) = 1$ and $q(-\infty ;\epsilon ) = 0$ . Since $p(t;\epsilon )$ and $q(t;\epsilon )$ define four degrees of freedom, we may remove two of these degrees of freedom by demanding that

(3.9) \begin{equation} \frac {\text{d}\tilde {\psi }}{\text{d}t}(t;\epsilon ) = -\frac {i \varOmega (t)}{\epsilon }\big[p(t;\epsilon ) X_{+}(t;\epsilon ) - q(t;\epsilon )X_{-}(t;\epsilon )\big]. \end{equation}

Here, $\tilde {\psi }$ satisfies the deformed Hill’s equation and has the desired form in the far past if and only if

(3.10) \begin{align} p(t;\epsilon ) &= 1 + \int _{-\infty }^{t} \frac {q(t';\epsilon )}{2\tilde {\varOmega }(t')} \frac {\text{d}\tilde {\varOmega }}{\text{d}t}(t')\exp \left (\frac {2i}{\epsilon } \int _{0}^{t'} \tilde {\varOmega }(t'')\text{d}t''\right ) \text{d}t', \nonumber \\ q(t;\epsilon ) &= \int _{-\infty }^{t} \frac {p(t';\epsilon )}{2 \tilde {\varOmega }(t')} \frac {\text{d}\tilde {\varOmega }}{\text{d}t}(t') \exp \left (-\frac {2i}{\epsilon } \int _{0}^{t'} \tilde {\varOmega }(t'')\text{d}t''\right ) \text{d}t' . \end{align}

Equation (3.10) yields a convergent Picard iteration scheme when the initial guess $(p_{0},q_{0}) = (1,0)$ is assumed. Equivalent to the Piccard iteration, we may write $p(t;\epsilon ) = \sum _{n = 0}^{\infty } p_{n}(t;\epsilon )$ and $q(t;\epsilon ) = \sum _{n = 1}^{\infty } q_{n}(t;\epsilon )$ with $p_{0} = 1$ and with

(3.11) \begin{align} q_{n+1}(t;\epsilon ) &= \int _{-\infty }^{t} \frac {p_{n}(t';\epsilon )}{2\tilde {\varOmega }(t')} \frac {\text{d}\tilde {\varOmega }}{\text{d}t}(t')\exp \left (-\frac {2i}{\epsilon } \int _{0}^{t'} \tilde {\varOmega }(t'')\text{d}t''\right ) \text{d}t', \nonumber \\ p_{n+1}(t;\epsilon ) &= \int _{-\infty }^{t} \frac {q_{n+1}(t';\epsilon )}{2\tilde {\varOmega }(t')} \frac {\text{d}\tilde {\varOmega }}{\text{d}t}(t')\exp \left (\frac {2i}{\epsilon } \int _{0}^{t'} \tilde {\varOmega }(t'')\text{d}t''\right ) \text{d}t'. \end{align}

The $\epsilon$ -independent bound

(3.12) \begin{equation} \sum _{n =0 }^{\infty } |p_{n}(t;\epsilon )| + \sum _{n = 0}^{\infty } |q_{n}(t;\epsilon )| \leqslant \exp \left (\int _{-\infty }^{t} \left |\frac {1}{2 \tilde {\varOmega }(t')} \frac {\text{d}\tilde {\varOmega }}{\text{d}t}(t') \right | \text{d}t' \right ) \lt \infty , \end{equation}

proves that the series for $p$ and $q$ converges uniformly on time intervals of the form $(-\infty ,T]$ . To show that that $\tilde {\psi }$ admits an asymptotic series of the form $\tilde {\psi } \sim (1 + \epsilon \tilde {\psi }_{1}(t) + \cdots ) X_{+}$ , we define the variable $\zeta (t) = \int _{0}^{t} \tilde {\varOmega }(t') \text{d}t'$ and the function $\alpha (\zeta ) = ({1}/{2 \tilde {\varOmega }^2}) ({\text{d}\tilde {\varOmega }}/{\text{d}t})(\zeta )$ . Then

(3.13) \begin{align} q_{n+1}(\zeta ;\epsilon ) &= \int _{-\infty }^{\zeta } p_{n}(\zeta ';\epsilon ) \alpha (\zeta ') \exp \left (-\frac {2i}{\epsilon } \zeta ' \right ) \text{d}\zeta ', \nonumber \\ p_{n+1}(\zeta ;\epsilon ) &= \int _{-\infty }^{\zeta } q_{n+1}(\zeta ';\epsilon ) \alpha (\zeta ') \exp \left (\frac {2i}{\epsilon } \zeta ' \right ) \text{d}\zeta '. \end{align}

Equation (3.13) is the advertised Bremmer series, as seen by integrating equations (9a) and (9b) in Bremmer (Reference Bremmer1951). Repeatedly integrating by parts, the standard result that $ \lim\nolimits _{\epsilon \to 0}\int _{D} \, f(\zeta ') \exp (\pm {i \zeta '}/{\epsilon })\text{d}\zeta ' = 0$ for any $f \in L^{1}(D\subset \mathbb{R})$ shows that

(3.14) \begin{equation} q_{1}(\zeta ;\epsilon ) = \int _{-\infty }^{\zeta } \alpha (\zeta ') \exp \left (- \frac {2 i}{\epsilon } \zeta '\right )\text{d}\zeta ' \sim -\frac {\epsilon }{2i} \exp \left (-\frac {2 i}{\epsilon } \zeta \right )\sum _{n = 0}^{\infty } \frac {\epsilon ^{n}}{(2i)^{n}} \frac {\text{d}^{n}\alpha }{\text{d}\zeta ^{n}} (\zeta ). \end{equation}

Since term-by-term integration of an asymptotic series yields an asymptotic series, we have that

(3.15) \begin{equation} p_{1}(\zeta ;\epsilon ) = \int _{-\infty }^{\zeta } q_{1}(\zeta ';\epsilon ) \alpha (\zeta ') \exp \left (\!\frac {2 i}{\epsilon } \zeta '\!\right ) \text{d}\zeta ' \sim -\frac {\epsilon }{2i}\sum _{n = 0}^{\infty } \frac {\epsilon ^{n}}{(2i)^{n}} \!\int _{-\infty }^{\zeta } \!\alpha (\zeta ') \!\left [ \frac {\text{d}^{n} \alpha }{\text{d}\zeta ^{n}} (\zeta ')\right ]\!\text{d}\zeta '. \end{equation}

Continuing in this fashion, its not hard to see that $q_{n}$ admits an asymptotic series of the form $q_{n} \sim \epsilon ^{n} \exp (-({2i}/{\epsilon })\zeta )( q_{n}^{0} + \epsilon q_{n}^{1} + \cdots )$ and $p_{n}$ admits an asymptotic powers series of the form $p_{n} \sim \epsilon ^n (p_{n}^{0} + \epsilon p_{n}^{1} + \cdots )$ . For any $N$ , the finite sums $\sum _{n = 0}^{N} p_{n}$ and $\exp (({2 i}/{\epsilon }) \zeta ) \sum _{n = 1}^{N} q_{n}$ admit asymptotic power series. To see that the limiting functions $p = \sum _{n = 0}^{N} p_{n} + \sum _{n \gt N} p_{n}$ and $\exp (({2 i}/{\epsilon }) \zeta ) q = \exp (({2 i}/{\epsilon }) \zeta ) \sum _{n = 1}^{N} q_{n} + \exp (({2 i}/{\epsilon }) \zeta ) \sum _{n\gt N} q_{n}$ admit asymptotic power series as well, we write $|q_{N+1}(t;\epsilon )|$ $ = \epsilon ^{N+1} C_{N+1}(\zeta , \epsilon )$ . One verifies that $C_{N+1}(t,\epsilon )$ is bounded on any domain of the sort $(\zeta , \epsilon ) \in (-\infty ,\varTheta ] \times (0, \infty ) \equiv \mathcal{D}(\varTheta )$ . This yields the uniform bound

(3.16) \begin{equation} \frac {1}{\epsilon ^{N+1}} \left ( \sum _{n \gt N} |p_{n}(\zeta ;\epsilon )| + \sum _{n \gt N}|q_{n}(\zeta ;\epsilon )| \right ) \leqslant \exp \left (||C||_{L^{\infty }(\mathcal{D}(\zeta ))} \int _{-\infty }^{\zeta } |\alpha (\zeta ')| \text{d}\zeta '\right )\lt \infty . \end{equation}

Hence, $\sum _{n \gt N} p_{n}, \sum _{n \gt N}q_{n} \in O(\epsilon ^{N+1})$ . Summing the asymptotic expression for $p_{n}$ and $q_{n}$ therefore shows that $p$ and $\exp (({2i}/{\epsilon }) \zeta )q$ admit asymptotic power series. Using the form of $\tilde {\psi }$ in the far past, we verify the real and imaginary parts of $\tilde {\psi }$ are linearly independent and normalised so that $\tilde {w}(t;\epsilon ) = |\tilde {\psi }(t;\epsilon )|$ solves the deformed Ermakov–Pinney equation. Defining $w(t;\epsilon ) = \tilde {w}(t;\epsilon )$ for $t \geqslant 0$ and uniquely extending $w(t;\epsilon )$ to $t \lt 0$ , we have a solution to the Ermakov–Pinney equation (2.6) admitting an asymptotic power series of the form $w(t;\epsilon ) \sim \varOmega ^{-1/2}(t) + \epsilon ^2 w_{2}(t) + \cdots$ when $t \geqslant 0$ . To show that $w(t;\epsilon )$ admits the same asymptotic series for $t \lt 0$ , one can again use the Bremmer series with $\tilde {\varOmega }(t)$ replaced by $\varOmega (t)$ and with $p(0;\epsilon )$ and $q(0;\epsilon )$ giving the initial conditions at $t = 0$ . The proof essentially proceeds as above. We remark that the cutoff function $H$ is only used to find suitable initial conditions for $w(t;\epsilon )$ .

To complete the proof that $w(t;\epsilon )$ is the desired function, we now show that $w'(t;\epsilon )$ also admits an asymptotic series. The proof is simple and leverages that $\tilde {w}(t;\epsilon )$ solves an initial value problem. Indeed, we have that

(3.17) \begin{equation} \frac {\text{d}}{\text{d}t}\tilde {w}(t;\epsilon ) = \epsilon ^{-2} \int _{-K}^{t} \left (\frac {1}{\tilde {w}^3(t',\epsilon )} - \tilde {\varOmega }^2 \tilde {w}(t';\epsilon ) \right )\text{d}t'. \end{equation}

Using that $1/z^3$ is analytic away from $z = 0$ , the integrand admits an asymptotic power series. Integrating term by term, we therefore have that

(3.18) \begin{align} \frac {\text{d}}{\text{d}t}\tilde {w}(t;\epsilon ) &\sim \epsilon ^{-2} \!\!\int _{-K}^{t} \!\left [\tilde {\varOmega }^{3/2}(t')(- 3 \tilde {w}_{2}(t') \tilde {\varOmega }^{1/2}(t')\epsilon ^2 + \cdots ) - \tilde {\varOmega }^2(t') (\tilde {w}_{2}(t') \epsilon ^2 + \cdots ) \!\right ] \text{d}t' \end{align}

(3.19) \begin{align} &= a_{0}(t) + \epsilon ^2 a_{2}(t) + \epsilon ^4 a_{4}(t) + \cdots . \end{align}

Since $\tilde {w}(t;\epsilon ) = \varOmega ^{-1/2}_{0} + \int _{-K}^{t} {\text{d}}/{\text{d}t}\tilde {w}(t';\epsilon )\text{d}t'$ , the uniqueness of asymptotic series implies that

(3.20) \begin{equation} \tilde {w}_{0}(t) + \epsilon ^2 \tilde {w}_{2}(t) + \cdots = \varOmega ^{-1/2}_{0} + \int _{-K}^{t} \left ( a_{0}(t') + \epsilon ^2 a_{2}(t') + \cdots \right )\text{d}t'. \end{equation}

Since $\tilde {w}_{n}(t)$ are smooth, differentiating (3.20) gives for $t \geqslant 0$ that

(3.21) \begin{equation} \frac {\text{d}}{\text{d}t}w(t;\epsilon ) \sim \frac {\text{d}}{\text{d}t} w_{0}(t) + \epsilon ^2 \frac {\text{d}}{\text{d}t} w_{2}(t) + \cdots . \end{equation}

Equation (3.21) easily extends to $t \lt 0$ by only a trivial modification of the proof above. Hence, the asymptotic series for $w'(t;\epsilon )$ is obtained by term-by-term differentiation of the asymptotic series for $w(t;\epsilon )$ .

We note that if $\int _{-\infty }^{\infty } |({\text{d}^{n}}/{\text{d}t^{n}})({\varOmega '}/{2 \varOmega })| \lt \infty$ for all $n \geqslant 0$ , then setting $H$ to unity shows $|w(\infty ,\epsilon ) - w(-\infty ,\epsilon )| \sim 1 + 0 \epsilon + 0 \epsilon ^2 + \cdots$ . Hence $\mathcal{M}_{\epsilon }$ is conserved to all orders. This is essentially the result of Littlewood (Reference Littlewood1963). More refined error bounds for the adiabatic invariance of the simple harmonic oscillator are discussed in Wasow (Reference Wasow1973), Stengle (Reference Stengle1977), Neishtadt (Reference Neishtadt1981) and Bonet et al. (Reference Bonet, Sauzin, Seara and València1998).

The Bremmer series is remarkable in that it reorders the formal power series of $\psi (t;\epsilon )$ to give a convergent solution $\psi = \psi _{0} + \psi _{1} + \cdots$ to Hill’s equation with $|\psi _{n}| \in O(\epsilon ^{n})$ . Other methods of constructing convergent WKB-like solutions can be found in Hyouguchi et al. (Reference Hyouguchi, Adachi and Ueda2002), Krivec et al. (Reference Krivec, Mandelzweig and Tabakin2004) and Winitzki (Reference Winitzki2005). When $\varOmega (t)$ is analytic, more refined methods based on Borel summation can be used to give exact solutions $\psi (t;\epsilon )$ with stronger asymptotic and analytic properties. We refer the interested reader to Costin et al. (Reference Hyouguchi, Adachi and Ueda2004), Voros (Reference Voros2004), Costin (Reference Hyouguchi, Adachi and Ueda2008), Takei (Reference Takei2017) and Nikolaev (Reference Nikolaev2023).

4. Non-perturbative dynamics

In this section, we fix $\epsilon$ and study the non-perturbative dynamics of the magnetic system when $\varOmega ^2(t)$ is $T$ -periodic. Using Floquet theory, we identify a form for $w(t)$ which greatly simplifies and elucidates the long-time dynamics of particles. This section illuminates some possible pitfalls of using a perturbative theory when $\epsilon$ is not so small or when long time scales are considered.

Consider Hill’s equation in first-order form with

(4.1) \begin{align} \frac {\text{d}}{\text{d}t} \begin{pmatrix} \psi \\ \varphi \end{pmatrix} = \epsilon ^{-1} \begin{pmatrix} 0 &\quad 1 \\ -\varOmega ^2(t) &\quad 0 \end{pmatrix}\begin{pmatrix} \psi \\ \varphi \end{pmatrix} = \epsilon ^{-1} M(t) \begin{pmatrix} \psi \\ \varphi \end{pmatrix}. \end{align}

Let $\varPhi (t;\epsilon ^{-1}) = \mathcal{T} \exp (\epsilon ^{-1} \int _{0}^{t} M(t')\text{d}t')$ be the fundamental matrix of Hill’s equation, and $B(\epsilon ^{-1}) = \varPhi (T;\epsilon ^{-1})$ the monodromy matrix. Let $\lambda _{1}, \lambda _{2}$ be the eigenvalues of $B(\epsilon ^{-1})$ , and suppose $(\boldsymbol{v}_{1}, \boldsymbol{v}_{2}) \in \mathbb{C}^2 \times \mathbb{C}^2$ is a basis putting $B$ in Jordan normal form. We define $\boldsymbol{z}_{1}(t) = \varPhi (t;\epsilon ^{-1})\boldsymbol{v}_{1}$ and $\boldsymbol{z}_{2}(t) = \varPhi (t;\epsilon ^{-1}) \boldsymbol{v}_{2}$ . Trivially, $\boldsymbol{z}_{1}(t)$ and $\boldsymbol{z}_{2}(t)$ are linearly independent solutions to Hill’s equations obeying $\boldsymbol{z}_{i}(t+T) = B \boldsymbol{z}_{i}(t)$ .

By Abel’s theorem, $\det (B(\epsilon ^{-1})) = 1$ . If $\lambda _{1}= \lambda _{2}^*$ are non-real, then $\lambda _{i}$ are distinct and have unit modulus. It follows that all solutions to Hill’s equation are bounded. Further, there exists a constant $c\gt 0$ such that

(4.2) \begin{equation} w_{\text{per}}(t) = c \big|\boldsymbol{z}^{1}_{1}(t)\big|, \end{equation}

is the unique periodic solution to the Ermakov–Pinney equation. One easily computes that

(4.3) \begin{equation} c = \left |\boldsymbol{z}_{1}^{1}(0)\right |^{-1}\left |\text{Im}\left (\frac {\boldsymbol{z}_{1}^{2}(0)}{\boldsymbol{z}_{1}^{1}(0)}\right )\right |^{-1/2}. \end{equation}

In the case when a periodic solution of Ermakov–Pinney exists, the exact invariant $\mathcal{M}$ Footnote ¹ can be constructed to be time periodic. Partially compactifying the particle phase space by quotienting $t \sim t + T$ , one can verify the invariance of $(\mathcal{M}, \tilde {H}_{\epsilon }, p_{\theta })$ implies the existence of action-angle variables by the Arnold–Liouville theorem. This opens the door to a renewed study of more complex magnetic geometries through, for example, Hamiltonian perturbation theory (Nekhoroshev Reference Nekhoroshev1971, Reference Nekhoroshev1977; Arnol’d Reference Arnol’d1989).

When the $\lambda _{i}$ are strictly real, there is a constant $a\gt 0$ such that

(4.4) \begin{equation} w_{\text{can}}(t) = a |\boldsymbol{z}_{1}(t) + i \boldsymbol{z}_{2}(t)| \end{equation}

solves the Ermakov–Pinney equation. If $\lambda _{1} = \lambda _{2} = \pm 1$ and $B$ is diagonalisable, then $w_{\text{can}}(t)$ is periodic and all solutions to Hill’s equation are bounded. Otherwise, $w(t)$ grows unbounded for any choice of solution to the Ermakov–Pinney equation. It follows all particle trajectories become unbounded at a typically exponential rate. This divergence may be attributed to the particle cyclotron motion resonating with the induced electric field. In such cases, it is not possible to define action-angle variables on the partially compactified particle phase space.

Letting $\Delta (\epsilon ^{-1}) = \text{tr}(B(\epsilon ^{-1}))$ , we have proved the long-term confinement of particles requires $|\Delta (\epsilon ^{-1})| \leqslant 2$ , with stability only ensured when the strict inequality holds. This fact could not have been uncovered at any order in perturbation theory. Perturbative theory predicts that $\mu _{\epsilon }$ is time periodic to all orders. If the series for $\mu _{\epsilon }$ converged, the limiting invariant would be time periodic as well. This would imply all particle trajectories are bounded. However, even small deviations of $\varOmega (t)$ from a constant function may cause particles to become de-confined after some time. Indeed, suppose that $\varOmega ^2(t) = \varOmega _{0}^2 - \epsilon ^2 \phi (t)$ . Then Hill’s equation becomes the one-dimensional Schrödinger equation

(4.5) \begin{equation} -\psi '' + \phi (t) \psi = E(\epsilon ^{-1}) \psi , \end{equation}

with $E(\epsilon ^{-1}) = \varOmega _{0}^2\epsilon ^{-2}$ . For any fixed $\phi (t)$ , the oscillation theory of Haupt (Magnus & Shenitzer Reference Magnus and Shenitzer1957) shows there exists two strictly increasing sequences $\rho _{0} \leqslant \rho _{2} \leqslant \ldots \to \infty$ and $\rho '_{1} \leqslant \rho '_{2} \leqslant \ldots \to \infty$ obeying the inequalities $\rho _{0} \lt \rho _{1}' \leqslant \rho _{2}' \lt \rho _{1} \leqslant \rho _{2} \lt \rho _{3}'\leqslant \rho _{4}' \lt \rho _{2} \leqslant \ldots$ such that the eigenvalues of $B(\epsilon ^{-1})$ are non-real if and only if

(4.6) \begin{equation} E(\epsilon ^{-1}) \in (\rho _{0},\rho _{1}') \cup (\rho _{2}', \rho _{1}) \cup (\rho _{2}, \rho _{3}') \cup \ldots . \end{equation}

Equation (4.6) shows that, in general, infinitely many regions of instability are incurred as $\epsilon ^{-1} \to \infty$ . Regarding $\epsilon ^2 \phi (t)$ as a perturbation, one should expect regions of instability to form around $\epsilon ^{-1}$ values resonating with the unperturbed gyrofrequency, such as in Mathieu’s equation (Arnol’d Reference Arnol’d1983; Kovacic et al. Reference Kovacic, Rand and Sah2018; Butikov Reference Hyouguchi, Adachi and Ueda2018).

The stability properties of Hill’s equation in the form of (4.5) is a classical and well-studied subject (Magnus & Shenitzer Reference Magnus and Shenitzer1957; Magnus & Winkler Reference Magnus and Winkler1961; Levy & Keller Reference Levy and Keller1963; Weinstein & Keller Reference Weinstein and Keller1985; Broer & Levi Reference Broer and Levi1995). Work on the asymptotic stability of Hill’s equation in the more relevant form (4.1) has also been conducted by Weinstein & Keller (Reference Weinstein and Keller1987). Weinstein and Keller have shown that for a certain class of analytic $\varOmega (t)$ , the regions where $|\Delta (\epsilon ^{-1})| \geqslant 2$ grow exponentially small as $\epsilon ^{-1} \to \infty$ . Hence, particle orbits are bounded for most $\epsilon \ll 1$ .

We remark that when $|\Delta (\epsilon ^{-1})| \lt 2$ , the existence of an exact, time-periodic invariant proves the global error in the magnetic moment is uniformly bounded. This condition can be used to correct some of the statements in Qin & Davidson (Reference Qin and Davidson2006) when $\varOmega (t)$ is periodic.

5. Numerical studies

We now use the results of the previous sections to compute and explore the accuracy of the adiabatic invariant $\mu _{\epsilon }$ when $\epsilon$ is not so small. By direct numerical computation, and by comparison with $\mathcal{M}_{\epsilon }$ , we find that $\mu _{\epsilon }$ may not always be useful for studying the particle dynamics. We demonstrate how, with minimal computation, the methods developed in this work present a better alternative to the perturbative theory of confinement in time-varying magnetic fields.

For our example, we choose the $2\pi$ -periodic magnetic field

(5.1) \begin{equation} \varOmega ^{-1/2}(t) = (1 + \exp (\sin (t)). \end{equation}

As we have fully justified, the formal solution to the Ermakov–Pinney equation,

(5.2) \begin{equation} w(t;\epsilon ) = w_{0}(t) + \epsilon ^2 w_{2}(t) + O(\epsilon ^4) = \varOmega ^{-1/2}(t) - \frac {\epsilon ^2 }{4 \varOmega ^2}\frac {\text{d}^2}{\text{d}t^2}\big(\varOmega ^{-1/2}\big)(t) + O(\epsilon ^4), \end{equation}

allows us to compute Kruskal’s adiabatic invariant $\epsilon ^{-1} \mu _{\epsilon } = \mu _{1} + \epsilon \mu _{2} + \epsilon ^2 \mu _{3} + \epsilon ^3 \mu _{4} + O(\epsilon ^4)$ . Namely, we have that

(5.3) \begin{align} \epsilon ^{-1} \mu _{1} &= \frac { p_{\theta }^2}{4} \left (\frac {w^2_{0}(t)}{r^2} \right ) + \frac {1}{4} \left (\frac {r^2}{w_{0}^2(t)} \right ) + \frac {1}{4}(p_{r}w_{0}(t))^2 - \frac {p_{\theta }\sigma }{2}, \nonumber \\ \epsilon ^{-1} \mu _{2} &= -\frac {p_{r}r }{4} \frac {\text{d} w_{0}^2}{\text{d}t}(t),\nonumber \\ \epsilon ^{-1} \mu _{3} &= \frac {p_{\theta }^2}{2} \left ( \frac {w_{0}(t)w_{2}(t)}{r^2}\right ) - \frac {r^2 w_{2}(t)}{2 w_{0}(t)^3} + \frac {r^2}{4}\left ( \frac {\text{d} w_{0}}{\text{d}t}(t)\right )^2 + \frac {p_{r}^2}{2} w_{0}(t) w_{2}(t), \nonumber \\ \epsilon ^{-1} \mu _{4} &= - \frac {p_{r}r}{2} \frac {\text{d}}{\text{d}t}(w_{0} w_{2})(t). \end{align}

To numerically explore the error of the conservation of the $n$ th-order invariant $\epsilon ^{-1} \mu ^{(n)} = \mu _{1} + \cdots + \epsilon ^{n-1}\mu _{n}$ , we set $p_{\theta } = \sigma = 1$ . We define the average $L^{2}$ error in the conservation of $\epsilon ^{-1}\mu ^{(n)}$ over two periods $\Delta \mu _{\epsilon }^{(n)} \equiv {\epsilon ^{-1}}/{4 \pi }||\mu ^{(n)}_{\epsilon }(r(t),p_{r}(t)) - \mu ^{(n)}_{\epsilon }(r(0),p_{r}(0))||_{L^{2}(0,4\pi )}$ . Using the initial conditions $(r(0),p_{r}(0)) = (2,0)$ , corresponding to initially unmoving particles, and the initial conditions $(r,p_{r}) = (2,5)$ , corresponding to initially gyrating particles, we simulate the particle orbits $(r(t),p_{r}(t))$ using (2.1) and the optimal fourth-order symplectic integrator of McLachlan & Atela (Reference McLachlan and Atela1992). We then compute $\Delta \mu ^{(n)}_{\epsilon }$ and plot the result against $\epsilon ^{-1}$ . The results are shown in figures 1 and 2. The numerical step size is chosen proportional to $\epsilon ^{-1}$ so that the integration error is asymptotically smaller than the error in the conservation of $\Delta \mu ^{(n)}$ for $n \in \{1,2,3\}$ . The constant of proportionality was ensured to be sufficiently small by recomputing the figures with significantly smaller step-sizes and observing no appreciable changes to the results.

Figure 1.

Value of $\Delta \mu ^{(n)}$ vs. $\epsilon ^{-1}$ for the initial conditions $(r(0),p_{r}(0)) = (2,0)$ . It is expected that $\Delta \mu ^{(1)}_{\epsilon }$ and $\Delta \mu ^ {(2)}_{\epsilon }$ are both $O(\epsilon ^2)$ since the particle initially has no velocity.

Figure 2.

Value of $\Delta \mu ^{(n)}_{\epsilon }$ vs. $\epsilon ^{-1}$ for the initial conditions $(r(0),p_{r}(0)) = (2,5)$ . The typical asymptotic behaviour as $\epsilon \to 0$ is seen for $\Delta \mu ^{(n)}_{\epsilon }$ .

It is interesting that $\Delta \mu _{\epsilon }^{(n)}$ oscillates as it tends to zero with $\epsilon$ . For any fixed, narrow range of $\epsilon$ values, the relative amplitude of these oscillations increases sharply as the accuracy of $\epsilon ^{-1}\mu _{\epsilon }^{(n)}$ is computed over longer time intervals. This is not unexpected. When $|\Delta (\epsilon ^{-1})| \gt 2$ , the resonance between the electric field and the cyclotron motion causes particle orbits to grow unbounded exponentially fast. This implies $\lim\nolimits _{t \to \infty }|\mu _{\epsilon }(t) - \mu _{\epsilon }(0)| = \infty$ . For these resonant $\epsilon$ values, one therefore expects that $\epsilon ^{-1}\mu _{\epsilon }^{(n)}$ is poorly conserved relative to the smallness of $\epsilon$ . The locations of these ill-behaved $\epsilon$ values can be estimated using the zeroth-order WKB approximation. As a zeroth-order approximation, one has that

(5.4) \begin{equation} \Delta \big(\epsilon ^{-1}\big) = 2 A\big(\epsilon ^{-1}\big) \cos \left (\epsilon ^{-1} \int _{0}^{T} \varOmega (t')\text{d}t' \right ) \equiv A\big(\epsilon ^{-1}\big)\Delta _{\text{as}}\big(\epsilon ^{-1}\big), \end{equation}

for some function $A(\epsilon ^{-1})$ limiting to unity as $\epsilon ^{-1} \to \infty$ . Arbitrarily long-term particle confinement may therefore fail when $\epsilon ^{-1}\int _{0}^{T} \varOmega (t')\text{d}t' \approx \pi N$ for some $N \in \mathbb{Z}$ . For the present example, we compute

(5.5) \begin{equation} \int _0^{2\pi }\varOmega (t^\prime )\text{d}t^\prime = \int _{0}^{2 \pi } ( 1 + \exp (\sin (t'))^{-2}\text{d}t' \approx 1.7452. \end{equation}

We numerically compute $\varDelta (\epsilon ^{-1})$ by finding the fundamental solutions to Hill’s equation and compare the result with $\varDelta _{\text{as}}(\epsilon ^{-1})$ . The results are shown in figure 3. The oscillations in this plot can be directly compared with the oscillations in figure 1.

Figure 3.

Value of $\Delta (\epsilon ^{-1})$ obtained numerically and $\Delta _{\text{as}}(\epsilon ^{-1})$ against $\epsilon ^{-1}$ for relatively large $\epsilon$ . The particle dynamics is stable when $|\Delta (\epsilon ^{-1})|\lt 2$ .

When $\epsilon ^{-1}$ is not much bigger than $20$ , figures 1 and 2 show that $\mu ^{(3)}$ does not improve $\mu ^{(2)}$ by much. Up to $O(\epsilon ^2)$ terms, the approximate conservation of $\mu ^{(2)}$ amounts to the conservation of

(5.6) \begin{align} \epsilon ^{-1} \tilde {\mu }^{(2)} = \frac { p_{\theta }^2}{4} \left (\frac {\varOmega ^{-1/2}(t)}{r} \right )^2 + \frac {1}{4} \left (\frac {r}{\varOmega ^{-1/2}(t)} \right )^2 + \frac {1}{4}\left (p_{r}\varOmega ^{-1/2}(t) - \epsilon r\frac {\text{d}}{\text{d}t} \varOmega ^{-1/2}(t)\right )^2.\nonumber\\[-1pt] \end{align}

For the present example, this should imply the $2\pi$ -flowmap $\exp (2 \pi X_{\epsilon })$ has the approximate invariant

(5.7) \begin{equation} \epsilon ^{-1} \tilde {\mu }_{\text{dis}}^{(2)}(\epsilon ) = \frac {1}{r^2} + \frac {r^2}{16} + \frac {1}{4}\left (2 p_{r}- \epsilon r \right )^2. \end{equation}

In general, non-perturbative theory implies that when $|\Delta (\epsilon ^{-1})| \lt 2$ , there exists an exact invariant of the $2\pi$ -flowmap of the form

(5.8) \begin{equation} \epsilon ^{-1} \mathcal{M}_{\text{dis}}(\epsilon ) = \frac {1}{4} \left (\frac {a^2}{r^2} \right ) + \frac {1}{4} \left (\frac {r^2}{a^2} \right ) + \frac {1}{4}(p_{r}a - \epsilon b r)^2. \end{equation}

Here, $a = w_{\text{per}}(0)$ and $b = w'_{\text{per}}(0)$ are the initial conditions of the periodic solution to the Ermakov–Pinney equation. In the present example, perturbative theory gives the approximation $a \approx 2$ and $ b \approx 1$ . To obtain $a$ and $b$ exactly, one computes an eigenvector of the monodromy matrix $B(\epsilon ^{-1})$ to (4.1). For example, if $\epsilon ^{-1} = 20$ then

(5.9) \begin{equation} B(20) =\begin{bmatrix}-0.8221 & \quad{-}0.9786 \\ 0.0464 & \quad{-}1.1613 \end{bmatrix}. \end{equation}

Here, $B(20)$ has the eigenvector $ \boldsymbol{v} = ( 0.9771, 0.1694 + 0.1286 i )^{t}$ and strictly complex eigenvalues. Using (4.3) to find the correct scaling for $\boldsymbol{v}$ , we have exactly that $a = 2.756$ and $b = 9.554$ . This differs considerably from the perturbative prediction. Evidently, $\mu _{\epsilon }$ fails to be well conserved. This becomes more evident when one looks at $\mu ^{(n)}_{\epsilon }(t)$ in figure 4. The large growth of $\mu _{\epsilon }^{(n)}$ reflects that the eigenvalues of $B(20)$ lie near the real axis, implying $\epsilon ^{-1} = 20$ is nearly resonant. Since $\epsilon ^{-1} = 20$ is not resonant, the Poincaré return theorem ensures that $\mu ^{(n)}_{\epsilon }(t)$ returns close to its initial value.

Figure 4.

Value of $\epsilon ^{-1} \mu _{\epsilon }^{(n)}(t)$ vs. $t$ for the initial condition $(r(0),p_{r}(0)) = (2,0)$ and $\epsilon ^{-1} = 20$ .

In contrast to the perturbative invariant (5.7), the non-perturbative invariant (5.8) fully characterises the orbit structure of the $2\pi$ -flowmap. We demonstrate this by running a full-orbit simulation of $(r(t),p_{r}(t))$ and plotting the resulting Poincaré plot. The result is shown in figure 5 with the level sets of (5.8) overlayed. We note that a full-orbit simulation is unnecessary beyond an independent validation of our results since $(r(2 \pi N), p(2 \pi N))_{N \in \mathbb{Z}}$ can be readily computed from $B(20)$ .

Figure 5.

Poincaré plot of $2\pi$ -flowmap near the fixed point. Level sets of $\mathcal{M}_{\text{dis}}(\epsilon )$ are overlayed in black. Particle orbits are obtained by a full-orbit simulation over 3000 periods.