Introduction

Magnetized plasmas are ionized gases in which the ambient magnetic field significantly alters particle trajectories. They are known to play crucial roles in various fields, including fusion research, solar-terrestrial and astrophysical environments, and plasma-based industries. Physically, the fundamental challenge in understanding magnetized plasmas is the vast disparity between the very fast gyromotion time scale of a charged particle and the much slower characteristic time scales for collective instabilities. Gyrokinetics provides a unified framework for describing the long-term spatiotemporal evolution of magnetized plasmas (see refs. 1,2,3 and references therein). As one of the major achievements in modern plasma physics, this theory was established by an asymptotic construction of the magnetic moment adiabatic invariant \(\bar{\mu }\), initially through the gyro-averaging method4,5,6,7,8 and later through the Lie-transform perturbation theory9,10,11,12,13. The details of the charged particle’s gyromotion are not of dynamical importance in the resulting gyrokinetic equation, which thereby reduces the kinetic problem from six dimensions to five. After half a century of intense pursuit, gyrokinetics now is the basis of numerous simulation codes and theoretical models used to study plasma instabilities, turbulence, and transport processes1,2,14,15,16. A notable example in this regard is the recent launch of ambitious simulation projects that strive to deliver a high-fidelity whole-device model of magnetic fusion devices within the gyrokinetic framework (see, e.g., ref. 17).

Despite its practical success, the validity regime of gyrokinetic theory is typically stated by the nonlinear gyrokinetic ordering in terms of the formal expansion parameters1,8,18 and often assumed. The gyrokinetic theory is a reduced kinetic theory derived from adiabaticity. To ensure its generality, it is essential to establish the robustness of adiabatic invariant across a wide range of plasma conditions. Previous studies, however, mainly focus on the construction of adiabatic invariant4,19,20,21, a clear physics picture of the destruction mechanism is still lacking. The aim of this study is to illuminate what kind of physics sets the validity limits of gyrokinetic theory. Below, we will perform a thorough analysis of the nonlinear particle dynamics to investigate (i) under which conditions the adiabaticity would be broken and, accordingly the gyrokinetic theory would become invalid; (ii) why the high-frequency gyrokinetic theory, which was developed in the early 1980s18,22,23,24,25,26, has not gained much popularity in implementation; (iii) whether the superadiabaticity27 could also be utilized to construct a reduced kinetic theory.

In this work, we show the conventional ordering assumption is not precise enough for the strict validity of the gyrokinetic theory. The particle dynamics must constitute a boundary layer problem to guarantee the persistent existence of adiabatic invariant over a reasonably wide range of conditions. Consequently, in contrast to previous ordering arguments which are semi-quantitative in nature, we find a quantitative, frequency-independent threshold in the normalized amplitude below which the gyrokinetic theory is generally valid for low-frequency perturbations. Given the ordering of spatiotemporal scales and fluctuation strength in the standard low-frequency gyrokinetic theory1, one could conclude that the corresponding normalized fluctuation level is considerably lower than the threshold. Therefore, the existence of such a threshold exhibits the robustness of the low-frequency gyrokinetic theory. The adiabaticity in high-frequency regimes, however, is sensitive to wave parameters, which raises concerns about the fundamental hypotheses of high-frequency gyrokinetic theory. Further analyses suggest that it is not possible to construct a reduced kinetic equation from superadiabaticity.

Results

The model

To elucidate the key physics involved in these highly complex problems, we restrict our attention to Taylor’s model4. This paradigm concerns the charged particle motion in a uniform magnetic field B = Bez with transverse perturbations, which, after a straightforward derivation, is described by the Hamiltonian

$$H=\frac{1}{2}({p}^{2}+{q}^{2})+A\cos \omega t\cos q.$$
(1)

Here ω is the wave frequency and time is normalized to the gyroperiod Ω−1 = mc/eB, with e being the charge and m the mass. Though simple, Taylor’s model is of fundamental importance in magnetized plasmas and was one of the seminal bases for the modern development of gyrokinetic theory1. It represents the simplest paradigm for wave-particle dynamics in both the electrostatic drift wave (with \(A=-{k}_{\perp }^{2}e\delta \phi /(m{\Omega }^{2})\))4 and shear Alfvén wave (with A = (k⊥vA/Ω)(δB⊥/B))28 (see Methods, section ‘Derivation of the model Hamiltonian’), where k⊥ is the perpendicular wavenumber, vA denotes the Alfvén velocity, and δϕ and δB⊥ are, respectively, amplitudes of the fluctuating electrostatic potential and magnetic field. Note that here, the nonlinearity parameter A is determined by the spatial scale and amplitude of the perturbation, but is independent of the wave frequency. The particle dynamics in Eq. (1) can be conceptually separated into two parts: the gyromotion about the background magnetic field (the q2/2 term), and the wave-particle trapping due to the perturbation (the \(A\cos \omega t\cos q\) term). Furthermore, one can easily verify that the phase space flow arising from Eq. (1) exhibits time-reversal symmetry. Whenever (q(t), p(t)) is an orbit of the system with the initial condition (q0, p0), both ( − q( − t), p( − t)) and (q( − t), − p( − t)) will be physically allowable orbits, with the initial conditions ( −q0, p0) and (q0, −p0), respectively. For periodic orbits, the associated dynamical invariant can be formulated as the action integral

$$I=\oint _{\gamma }pdq,$$
(2)

where γ is a closed orbit at constant time.

Further progress is possible if one introduces a variable h conjugate to t, yielding the extended phase space Hamiltonian

$${{\mathcal{H}}}=\frac{1}{2}({p}^{2}+{q}^{2})+A\cos \omega t\cos q-h.$$
(3)

In this way, the time-dependent one-dimensional Hamiltonian system is replaced by a two-dimensional Hamiltonian system where time plays a role analogous to that of an angle variable. Technically, to illustrate the long time scale (ωt ≫ 1) dynamical complexity of the system (which determines the existence of invariants), we construct the Poincaré map in extended phase space: \({{{\boldsymbol{x}}}}_{n+1}={{{\mathcal{M}}}}_{T}{{{\boldsymbol{x}}}}_{n}\), where x = (q, p), and the mapping \({{{\mathcal{M}}}}_{T}\) is defined such that the point xn at tn is advanced by Eq. (3) to the next crossing point xn+1 at tn+1 = tn + T, with T = 2π/ω.

Stability analysis

The global property of a Hamiltonian system can be understood by examining the orbits close to fixed points, which are either elliptic (stable) or hyperbolic (unstable). As depicted in Fig. 1, \({{{\mathcal{M}}}}_{T}\) exhibits two different types of fixed points in accordance with the significance of nonlinearity. The primary fixed point xp  = (0, 0) is also the fixed point of the original system (1). It becomes unstable when the linear parametric resonance29 occurs. Secondary fixed points xs, as assured by Poincaré–Birkhoff theorem30, are nonlinearly generated by a finite amplitude perturbation, half of which are elliptic and the other half hyperbolic in an alternating sequence. In general, the stability of secondary fixed points is not amenable to analytical analysis, numerical computations are mandatory.

Fig. 1: Examples of (primary and secondary) fixed points and (adiabatic, superadiabatic, and chaotic) orbits.
figure 1

Colors indicate distinct particles in the phase space, and x = (q, p) are the canonical coordinates. In this case, first, we evolve 180 particle orbits in time according to the Hamiltonian (1) for ω = 0.5 and A = 1.3, with initial conditions (q0, p0) selected randomly in [ − 5, 5] × [ − 4, 4]. Then the phase space portrait of the 2D Hamiltonian map \({{{\boldsymbol{x}}}}_{n+1}={{{\mathcal{M}}}}_{T}{{{\boldsymbol{x}}}}_{n}\) is constructed such that the point xn at tn is advanced to the next crossing point xn+1 at tn+1 = tn + T, with the wave period T = 2π/ω.

Full size image

To clarify the nature of the primary fixed point, linear stability analysis has been performed analytically via a perturbative treatment [Methods, section “Stability of primary fixed point”]. It follows that different unstable domains can be labeled by the marginal stability condition (ω = 2/m, A = 0), with m a positive integer. Notably, a numerical evaluation of the linear growth rate \({{\rm{Im}}}\lambda\) (see Fig. 2) identifies two distinct regimes. In the high-frequency regime with \(\omega \gtrsim {{\mathcal{O}}}(1{0}^{-1})\), the threshold amplitude Ac varies rapidly with ω. However, the low-frequency regime (\(\omega \lesssim {{\mathcal{O}}}(1{0}^{-1})\)) is characterized by a dense spectrum corresponding to high- order resonances with m ≫ 1, where \({{\rm{Im}}}\lambda\) is insensitive to m and Ac piles up at a constant value Ac = 1.

Fig. 2: Stability diagram for the Taylor’s model.
figure 2

The color contours show \({{\rm{Im}}}\lambda\) of Eq. (17), which measures the growth rate of separation of initially close orbits near the primary fixed point x = (0, 0). ω and A are the wave frequency and amplitude, respectively. The color bar represents the magnitude of \({{\rm{Im}}}\lambda\). The white (magenta) line is the threshold numerically evaluated for the onset of superadiabatic (chaotic) orbits.

Full size image

The dynamical features of the system are intrinsically linked to the type of fixed points and, as illustrated in Fig. 1, can be put into three categories. The adiabatic orbit rotates with respect to the primary fixed point in the Poincaré map. The associated Kolmogorov–Arnold–Moser (KAM) torus30 is deformed but remains intact topologically. It possesses the adiabatic invariant \(\bar{\mu }\), which can be derived from Eq. (2) via a perturbative treatment, as in ref. 4. In regions of phase space where resonances take place, the original KAM torus will be destroyed, and the perturbation series for \(\bar{\mu }\) fails to converge20. Nonetheless, the emergence of a new KAM torus gives rise to a superadiabatic invariant27, leading to the orbit rotating around a secondary fixed point. We refer to this type of orbit as the superadiabatic orbit. The distinguishing feature of chaotic orbits is the destruction of separating KAM tori. Quantitatively, the chaotic orbit in the present study is determined by the finite size Lyapunov exponent λL (for details, see “Methods: Lyapunov analysis of orbit topology”). Examples in Fig. 3 show the topological differences between adiabatic, superadiabatic, and chaotic orbits in extended phase space. The orbit in Fig. 3c displays a chaotic behavior with λL ≃0.014.

Fig. 3: Comparison of the adiabatic, superadiabatic, and chaotic orbits in the extended phase space.
figure 3

The orbits are plotted for A = 1.3 and ω = 0.5 with different initial conditions (q0, p0): (q0 = 3.2, p0 = 0) for the adiabatic orbit in (a), (q0 = 2.5, p0 = 0) for the superadiabatic orbit in (b), and (q0 = 1.0, p0 = 0) for the chaotic orbit in (c).

Full size image

We now examine the transition from adiabatic to superadiabatic/chaotic orbits by scanning A and ω. A result, shown by the white line in Fig. 2, is the identification of a quantitative, frequency-independent threshold Ac ≃ 1 for \(\omega \lesssim {{\mathcal{O}}}(1{0}^{-1})\), which ensures the general applicability of adiabaticity in the low-frequency regime. To understand the physics of this threshold, we notice that Eq. (1) can be rendered into a boundary layer problem in the low-frequency limit, yielding

$${\omega }^{2}\frac{{d}^{2}q}{d{\theta }^{2}}+q-A\cos \theta \sin q=0,$$
(4)

with θ = ωt. Equation (4) displays manifestly the disparity between the short gyromotion time scale (∣∂θ∣ ~ ω−1), the wave-particle trapping time scale (\(| {\partial }_{\theta }| \sim \sqrt{A}/\omega\)), and the long time scale of the wave period (∣∂θ∣ ~ 1). On the gyromotion time scale, \(\cos \theta\) can be treated as a constant, and one can easily show that the primary fixed point becomes unstable when A > 1, in agreement with Fig. 2. Here, the adiabaticity threshold Ac = 1 could be interpreted as the condition that the wave-particle trapping frequency (\(\simeq \sqrt{A}\)) of resonant particles is comparable to the gyrofrequency (=1). On the long time scale ∂θ ~1, we have \(q-A\cos \theta \sin q\simeq 0\), then the secondary fixed points, which appear on time scales long compared to the wave period, are only possible if \(A\simeq | q/\sin q| \ge 1\). Additionally, numerical analyses suggest that Ac  = 1 also gives a quantitative estimate for the chaotic threshold in low-frequency regime, as seen in Fig. 2. There are no superadiabatic or chaotic orbits for A < 1 in the low-frequency limit (see for example, Fig. 4a, b). Figure 4c shows that a chaotic layer arises as the homoclinic tangle30 formed in the proximity of the hyperbolic primary fixed point. But the onset of chaos via the heteroclinic tangle30, originating from unstable hyperbolic secondary fixed points, is observed as well (see Fig. 4d).

Fig. 4: Phase space portrait of the 2D Hamiltonian map in low-frequency regime.
figure 4

Colors indicate distinct particles. In this case, we evolve 180 particle orbits in time according to the Hamiltonian (1) for a ω = 0.06, A = 0.9, b ω = 0.059, A = 0.9, c ω = 0.06, A = 1.01, and d ω = 0.059, A = 1.01, where initial conditions (q0, p0) are selected randomly in [ − 3, 3] × [ − 1.5, 1.5]. Then the 2D Hamiltonian map \({{{\boldsymbol{x}}}}_{n+1}={{{\mathcal{M}}}}_{T}{{{\boldsymbol{x}}}}_{n}\) is constructed such that the point xn at tn is advanced to the next crossing point xn+1 at tn+1 = tn + T, with the wave period T = 2π/ω.

Full size image

Validity of gyrokinetic theory

From the discussion above, it follows that the gyrokinetic theory is valid in the low-frequency regime. Even though this is consistent with the gyrokinetic ordering assumptions1, the underlying physical interpretation is substantially different. The formal validity of gyrokinetic theory relies on the particle dynamics constituting a boundary layer problem in the present paradigm. This requirement establishes a quantitative, frequency-independent adiabaticity threshold, ensuring the persistent existence of the adiabatic invariant across a wide range of conditions. Thus, here, the validity of the gyrokinetic theory is a quantitative issue. Once below the threshold, gyrokinetics exhibits equal accuracy for all low-frequency waves. In the conventional qualitative paradigm, however, the low-frequency gyrokinetic theory is expected to be more accurate for the wave with lower frequency. Therefore, the inherent nature of particle dynamics as a boundary layer problem is essential for the validity of gyrokinetic theory. Specifically, recalling the spatiotemporal scales and relative fluctuation levels assumed in the low-frequency nonlinear gyrokinetic theory1: \({k}_{\perp }v/\Omega \sim {{\mathcal{O}}}(1)\), ∣ω/Ω∣ ≪ 1 and ∣eδϕ/mv2∣, ∣δB⊥/B∣ ≪ 1, one can readily show that the associated normalized fluctuation amplitude is typically lower than the threshold Ac = 1. The existence of such a quantitative and frequency-independent adiabaticity threshold thus reveals the robustness of the standard low-frequency gyrokinetic theory.

The separation of time scales does not apply for high-frequency perturbations. Accordingly, Fig. 2 shows that both the superadiabatic and chaotic thresholds are highly sensitive to the wave frequency. Thus the adiabaticity is not universally preserved, as opposed to the low-frequency case. This difference calls into question the validity of high-frequency gyrokinetic theory, which is derived from the assumption of adiabatic orbits18,24,25,26. In addition, the breakdown of adiabaticity also vitiates the fundamental hypothesis of high-frequency gyrokinetics that the existence of gyrocenter coordinates is independent of ω. More specifically, noting that the perpendicular particle velocity in the present uniform plasma model satisfies vx ∝ p and vy ∝ q4, the phase space portraits of \({{{\mathcal{M}}}}_{T}\) (shown in Fig. 5) indicate that the near identity transformation between the particle and gyrocenter coordinates, which is expressed as an asymptotic expansion in powers of the wave amplitude in Lie perturbation theory1,2,18,24, will be invalidated by the presence of superadiabatic/chaotic orbits. The high-frequency gyrokinetic theory is, therefore, strictly valid only in the linear limit with A = 0+22,23. Nevertheless, when considering small but finite amplitude perturbations, the linear high-frequency gyrokinetic theory becomes inadequate for describing the particle dynamics near resonant points (ω = 2/m), due to the neglect of parametric resonance effect. Consequently, it cannot be utilized to develop quasilinear or weak turbulence theories.

Fig. 5: Phase space portrait of the 2D Hamiltonian map in a high-frequency regime.
figure 5

Colors indicate distinct particles. In this case, we evolve 180 particle orbits in time according to the Hamiltonian (1) for a ω = 0.50, A = 0.9, b ω = 0.52, A = 0.9, c ω = 1.00, A = 0.9, d ω = 0.98, A = 0.9, e ω = 0.50, A = 1.3, f ω = 0.52, A = 1.3, g ω = 1.00, A = 1.3, and h ω = 0.98, A = 1.3, where initial conditions (q0, p0) are selected randomly in [ −5, 5] × [ − 4, 4]. Then the 2D Hamiltonian map \({{{\boldsymbol{x}}}}_{n+1}={{{\mathcal{M}}}}_{T}{{{\boldsymbol{x}}}}_{n}\) is constructed such that the point xn at tn is advanced to the next crossing point xn+1 at tn+1 = tn + T, with the wave period T = 2π/ω.

Full size image

Unlike in the low-frequency regime, Fig. 2 shows that while high-frequency perturbations break the adiabaticity more easily, the onset of chaos deviates significantly from the superadiabatic threshold. In this situation, one may attempt to construct a new reduced kinetic theory from the superadiabatic invariant. Unfortunately, this possibility is excluded. First, the intrinsically nonperturbative nature of superadiabatic orbits and the associated complexity of phase space structures (cf. Fig. 5) make it difficult, if not impossible, to derive a comprehensive analytic description for the superadiabatic invariant. Secondly, a general reduced kinetic theory requires that all particles possess the same invariant over a reasonably broad range of conditions. However, Fig. 5 demonstrates that different particles could have distinct superadiabatic invariants at a fixed ω and, for each specific particle, the existence of corresponding superadiabatic invariant also depends sensitively on physical parameters.

Whilst it is well established that the Lie series20 of adiabatic/superadiabatic invariant will not converge for chaotic orbits due to the presence of nonlinear resonance28,31, intuition suggests that it might be possible to construct a gyrokinetic theory valid for a sufficiently long time via the gyro-averaging approach, if the orbit is only marginally chaotic. However, we find that once the chaotic threshold is exceeded, homoclinic and heteroclinic tangles can scatter the particle motion into different types of orbits without regularity, as shown in Fig. 6. In this scenario, the gyro-averaging (which refers to the integration along unperturbed adiabatic/superadiabatic orbits here) is ill-defined. It is thus impossible to properly account for the chaotic motion in both gyro-averaging and Lie perturbation approaches.

Fig. 6: Typical evolution of the 2D Hamiltonian map for marginally chaotic orbits.
figure 6

The particle is randomly scattered into different types of orbits by homoclinic (a) or heteroclinic (b) tangles.

Full size image

Discussion

We investigate the validity of gyrokinetic theory by employing the simplest paradigm. It is demonstrated that the conventional ordering assumption may not be precise enough for the strict validity of the gyrokinetic theory. The particle dynamics must constitute a boundary layer problem to guarantee the persistent existence of adiabatic invariance over a wide range of conditions. Accordingly, contrary to previous semi-quantitative ordering arguments, we find a quantitative, frequency-independent adiabaticity threshold Ac below which the gyrokinetic theory is universally valid in the low-frequency regime. In particular, using the spatiotemporal scales and fluctuation levels in the low-frequency nonlinear gyrokinetic theory, it can be estimated that the corresponding normalized fluctuation amplitude is considerably lower than the threshold Ac. In this sense, the existence of such a threshold reveals the robustness of the standard low-frequency gyrokinetic theory. In the high-frequency regime, however, the adiabaticity displays sensitive dependence on the wave parameters, which thereby raises concerns about the fundamental hypotheses of high-frequency gyrokinetic theory. Furthermore, unlike the situation with adiabaticity, it is not feasible to develop a new general reduced kinetic theory based on superadiabaticity.

From these considerations, we conclude that once the wave amplitude satisfies A < Ac, the gyrokinetic theory is reasonably applicable to the low-frequency drift-Alfvénic turbulence, whose wave frequency is typically \({{\mathcal{O}}}(1{0}^{-2})\) smaller than the gyrofrequency16. Meanwhile, recalling the definition of A, it is worth mentioning the condition A < Ac implies that the drift kinetic theory, which assumes k⊥ → 0+ and ω ≪ Ω, applies for arbitrary fluctuation amplitudes, consistent with previous results32,33. But still, it is of great concern whether one can employ the gyrokinetic theory, for example, in the vicinity of the separatrix of modern divertor tokamaks, where the radial wavenumber tends to diverge due to the magnetic topology change34. Since there exist global nonlinear gyrokinetic codes including the divertor X-point region, such as XGC35, it would be important to compare the predictions of turbulent fluxes and plasma profile evolution across the separatrix resulting from nonlinear gyrokinetic simulations with those of first-principle kinetic analyses at fixed boundary conditions (by this, we mean decoupling the well known complications of extending kinetic simulations into the plasma scrape-off layer.).

Moreover, we remark that while this study only deals with a single-frequency fluctuation in order to simplify the presentation, the present analysis can be readily generalized to turbulence with a broad spectrum, and the conclusions remain valid. In addition to the time scale considered here, spatial variations may also break the adiabatic invariant, as observed numerically in spherical tokamaks36,37. Though a similar approach could also apply here, a systematic and analytically treatable formalism incorporating both the spatial and time scales on an equal footing is beyond the intended scope of this study.

Methods

Derivation of the model Hamiltonian

Following4, we consider the motion of a charged particle immersed in a uniform magnetic field B = Bez with the electrostatic perturbation

$${{\boldsymbol{E}}}=-\nabla \Phi \equiv \nabla [\delta \phi \cos (\omega t)\cos ({k}_{\perp }x)].$$
(5)

The associated particle Hamiltonian can be written

$$H=\frac{1}{2m}[{p}_{x}^{2}+({{p}_{y}-\frac{eBx}{c}})^{2}+{p}_{z}^{2}]+e\Phi ,$$
(6)

where the momenta py = mvy + eBx/c and pz are constant. Then, the equations of motion are

$${\dot{p}}_{x} = \frac{eB}{mc}({p}_{y}-\frac{eBx}{c})-e{k}_{\perp }\delta \phi \cos (\omega t)\sin ({k}_{\perp }x),\\ \dot{x} = \frac{{p}_{x}}{m}.$$
(7)

By taking py = 0, Eq. (7) becomes an oscillator under the influence of external force

$$\ddot{x}+{\Omega }^{2}x=-\frac{e{k}_{\perp }\delta \phi }{m}\cos (\omega t)\sin ({k}_{\perp }x).$$
(8)

Normalizing the space to \({k}_{\perp }^{-1}\) and time to Ω−1, Eq. (8) is readily cast as

$$\ddot{q}+q=A\cos \omega t\sin q,$$
(9)

which is derivable from the Hamiltonian (1) with \(A=-{k}_{\perp }^{2}e\delta \phi /(m{\Omega }^{2})\), q = k⊥x, and \(p=\dot{q}\). Therefore, it can be recognized that parameter A is the normalized amplitude of the perturbed force and quantifies the strength of nonlinearity.

Equation (1) offers a simple paradigm for the wave-particle dynamics in shear Alfvén wave (SAW)28 as well. To see this, let us consider an obliquely propagating (k = kzez + kxex) shear Alfvén wave, whose electromagnetic fluctuation, by employing the Walén relation38, can be expressed as δB = δByey and δE = (vA/c)δByex. Then the Lorentz equations of motion are

$$m{\dot{v}}_{x} =\frac{e}{c}[{v}_{y}B+({v}_{A}-{v}_{z})\delta {B}_{y}],\\ {v}_{y} =-\Omega x,\\ m{\dot{v}}_{z} =\frac{e}{c}{v}_{x}\delta {B}_{y},$$
(10)

and the energy \({v}^{2}={v}_{x}^{2}+{v}_{y}^{2}+{({v}_{z}-{v}_{A})}^{2}\) is a constant of motion.

Although Eq. (10) states that SAW can influence the parallel particle dynamics and bring relevant extra-dimensions to the system, it was shown39 that the particle motion is non-chaotic for parallel propagation SAWs (k = kzez). The perpendicular dynamics with a finite kx is essential for the chaotic behavior of the system40. To provide analytical insights and, thus, clarity in illustration, we focus on the particles initially cold in the laboratory frame (\({v}_{z}^{2}(0)\ll {v}_{A}^{2}\)). To the first order in δBy/B, one arrives at the formula

$$\ddot{x}+{\Omega }^{2}x=({v}_{A}-{v}_{z}(0))\Omega \frac{\delta {B}_{y}}{B}.$$
(11)

Again, if we normalize x to \({k}_{\perp }^{-1}\) and time to Ω−1, and take the SAW frequency ωA ≃ k∥vA/Ω, Eq. (11) can be rendered into

$$\ddot{q}+q=\frac{{k}_{\perp }({v}_{A}-{v}_{z}(0))}{\Omega }\frac{\delta {B}_{y}}{B}.$$
(12)

Writing \(\delta {B}_{y}=\delta {B}_{\perp }\cos (({\omega }_{A}-{k}_{\parallel }{v}_{z})t)\sin ({k}_{\perp }x)\) and ω = ωA − k∥vz(0) ≃ ωA, Eq. (12) becomes identical to Eq. (9) with A = (k⊥vA/Ω)(δB⊥/B) and q = k⊥x.

Furthermore, it should be stressed that the arguments presented in the current study can also be applied to ascertain the existence of adiabaticity when retaining parallel dynamics. In particular, using \({({v}_{z}-{v}_{A})}^{2}={v}^{2}-{v}_{x}^{2}-{v}_{y}^{2}\), Eq. (12) can be cast as

$$\ddot{q}+q=\frac{\sigma \sqrt{{k}_{\perp }^{2}{v}^{2}-{\Omega }^{2}({\dot{q}}^{2}+{q}^{2})}}{\Omega }\frac{\delta {B}_{\perp }}{B}\cos ({k}_{\parallel }{\int}^{t}\sigma \sqrt{{v}^{2}-\frac{{\Omega }^{2}}{{k}_{\perp }^{2}}({\dot{q}}^{2}+{q}^{2})}\frac{dt^{\prime} }{\Omega })\sin q,$$
(13)

in which σ denotes the sign of vA − vz. For the primary fixed point, linearizing Eq. (13) around (q = 0, p = 0) yields the same expression as Eq. (15) with \({v}^{2}={({v}_{z}(0)-{v}_{A})}^{2}\). Hence, Fig. 2 also characterizes the primary fixed point stability of the system (10). On the other hand, noting that Eq. (13) is a boundary layer problem in the low-frequency limit, we find the secondary fixed points appearing on long time scales are only possible if

$$\frac{{k}_{\perp }v}{\Omega }\frac{\delta {B}_{\perp }}{B}\simeq | \frac{q}{\sqrt{1-\frac{{\Omega }^{2}}{{k}_{\perp }^{2}{v}^{2}}({\dot{q}}^{2}+{q}^{2})}\cos ({k}_{\parallel }{\int}^{t}\sigma \sqrt{{v}^{2}-\frac{{\Omega }^{2}}{{k}_{\perp }^{2}}({\dot{q}}^{2}+{q}^{2})}\frac{dt^{\prime} }{\Omega })\sin q}| > 1.$$
(14)

As a result, the system (10) still possesses the quantitative and frequency-independent adiabaticity threshold Ac = (k⊥v/Ω)(δB⊥/B) = 1 in the low-frequency regime.

Stability of primary fixed point

The derivation of the stability condition of the primary fixed point starts from the Hamiltonian in Eq. (1). By approximating \(\cos q\simeq 1-{q}^{2}/2\), the corresponding equation of motion reduces to a Mathieu equation:

$$\ddot{q}+(1-A\cos \omega t)q=0.$$
(15)

According to the Floquet theory, the formal solution of Eq. (15) can be expressed as

$$q={e}^{i\lambda t}\sum\limits_{n}{Q}_{n}{e}^{in\omega t},$$
(16)

where the Floquet characteristic exponent λ is a measure of the rate of separation of orbits close to the primary fixed point. xp is elliptic (hyperbolic) if λ is real (complex). Substituting Eq. (16) into Eq. (15) and following the procedure in ref. 41, one obtains

$$\lambda =\left\{\begin{array}{ll}\frac{\omega }{\pi }\arcsin \sqrt{D(0){\sin }^{2}(\frac{\pi }{\omega })},&\omega \ne \frac{1}{l},\\ \frac{\omega }{2\pi }\arccos [2D(\frac{\omega }{2})-1],&\omega =\frac{1}{l},\end{array}\right.$$
(17)

where D is the determinant of an infinite tridiagonal matrix D = (djk) with

$${d}_{jk}=\left\{\begin{array}{ll}1\hfill &{{\rm{if}}}\,k=j,\\ \frac{A}{2[{(\lambda +j\omega )}^{2}-1]}\quad &{{\rm{if}}}\,k=j\pm 1,\\ 0\hfill &{{\rm{otherwise.}}}\end{array}\right.$$
(18)

l, j, k are integers. For sufficiently small A ≪ ω2/2, an analytic expression for D(0) can be concisely given as D(0) ≃ 1 − RA2 with \(R\equiv (\pi /\omega )\cot (\pi /\omega )/2(4-{\omega }^{2})\)41. By means of Eq. (17), one can derive the threshold for hyperbolic xp

$${A}_{c}^{2}=\left\{\begin{array}{ll}{R}^{-1}\hfill &{{\rm{if}}}\,R > 0,\hfill\\ {R}^{-1}(1-{\sin }^{-2}\frac{\pi }{\omega })\quad &{{\rm{otherwise.}}}\end{array}\right.$$
(19)

Equation (19) allows one to label different unstable domains using the marginal stability condition (ω = 2/m, A = 0), with m a positive integer.

Lyapunov analysis of orbit topology

Chaotic orbits in this work are determined in terms of a modified version of the finite size Lyapunov exponent (FSLE) λL.

Consider at t = 0 a reference trajectory x(0) and a perturbed trajectory \({{\boldsymbol{x}}}^{\prime} (0)={{\boldsymbol{x}}}(0)+\delta {{\boldsymbol{x}}}(0)\), the basic idea behind FSLE is to quantify the average rate of error growth at different scales when subject to non-infinitesimal perturbations42. Specifically, let us first take an infinitesimal initial perturbation δmin ≪ 1 and choose a set of thresholds at different scales δn = δ0ρn. Here, δmin ≪ δ0 ≪ 1, n = 0, ⋯, N, ρ > 1, and δN accounts for the maximum allowed separation. The parameters we have adopted for the current study are δmin = 10−6, d0 = ∣δx(0)∣ + 10−5 and \(\rho =\sqrt{2}\). The trajectories are advanced using a fourth-order Runge–Kutta algorithm. After the perturbation has grown from δmin to δ0, we measure the time ti(δn) from δn to δn+1. A single error-doubling experiment completes when the error reaches δN. Then this procedure is repeated M times to obtain the set of doubling times {ti(δn)} for i = 1, ⋯, M error-doubling experiments. The growth rate γi is estimated through the average of τi across different scales,

$${\gamma }_{i}=\frac{\ln \rho }{{\langle {\tau }_{i}({\delta }_{n})\rangle }_{\delta }}.$$
(20)

We define the averaged FSLE as the averaging over doubling experiments \({\lambda }_{L}={\langle {\gamma }_{i}\rangle }_{M}\). Numerically, we find Eq. (20) converges much faster than the FSLE originally given in ref. 42.