Equations of Motion Near Cyclotron Resonance

This work compares several versions of the equations of motion for a test particle encountering cyclotron resonance with a single, field-aligned whistler mode wave. The gyro-averaged Lorentz equation produces both widespread phase trapping (PT) and “positive phase bunching” of low pitch angle electrons by large amplitude waves. Approximations allow a Hamiltonian description to be reduced to a single pair of conjugate variables, which can account for PT as well as phase bunching at moderate pitch angle, and has recently been used to investigate this unexpected bahavior at low pitch angle. Here, numerical simulations using the Lorentz equation and several versions of Hamiltonian-based equations of motion are compared. Similar behavior at low pitch angle is found in each case.


INTRODUCTION
Cyclotron-resonant wave-particle interactions are a crucial ingredient in magnetospheric dynamics, especially in the radiation belts, and there is a vast tradition of simulating the process as quasi-linear diffusion of phase space density by a broad-band spectrum of small, incoherent waves (Thorne, 2010;Thorne et al., 2013), following the pioneering work of Lyons et al. (1971) and Lyons et al. (1972). A complementary approach is that of test particle simulation, most often in the presence of a single, coherent wave whose amplitude need not be small. Inan et al. (1978) noted both quasi-linear and nonlinear behavior, including the "loss cone reflection effect" whereby low pitch angles increase rather than decrease below zero. In the quasi-linear regime, connections between the two perspectives have been provided by Lemons et al. (2009), Lemons (2012), Allanson et al. (2022), and a unified picture of quasi-linear and nonlinear behavior was obtained by Albert (2001), Albert (2010). These studies all used specified and idealized models of the waves, while Liu et al. (2010), Liu et al. (2012) examined test particles driven by waves from self-consistent particle-in-cell simulations.
This work compares several versions of the equations of motion for a test particle encountering cyclotron resonance with a single, field-aligned whistler mode plane wave. The Lorentz force law, resolved into components parallel and perpendicular to the background magnetic field and gyroaveraged, is commonly used for such simulations. Hamiltonian descriptions are in principle equivalent, and with several approximations they allow the reduction to a one-dimensional (1D) system (one action-angle pair, plus the independent variable playing the role of time). If the time dependence is slow enough, particle motion is nearly along instantaneously drawn contours, with invariant breaking at separatrix crossings. There is a rich literature of work based on these concepts, which has been exploited in this context to some degree. Among others, Shklyar (1986) Albert (1993), Albert (2000), Artemyev et al. (2018) further approximated the Hamiltonian as equivalent to that of a time-dependent pendulum and obtained quantitative estimates of energy and pitch angle changes, which have proved useful and reliable.
Recently, using the gyro-averaged Lorentz equation, Kitahara and Katoh (2019), Gan et al. (2020) found both widespread (or "anomalous") phase trapping (APT) and "positive phase bunching (PPB)" of low pitch angle electrons by large amplitude waves. Both phenomena lead to pitch angle increase, in contrast to the phase bunching behavior that is the usual alternative to phase trapping, and are associated with low pitch angle, which violates a certain approximation made in obtaining the pendulum Hamiltonian. Albert et al. (2021), Artemyev et al. (2021) presented generalizations of the pendulum Hamiltonian which avoid that specific approximation, but still relied on several others. In particular, differences in the first-order (in wave amplitude) term of the phase evolution equation are present among several versions of the equations of motion. This work shows numerically that, despite these differences, the generalized 1D Hamiltonian reproduces the behavior at low pitch angle, and is therefore an appropriate framework for the future development of refined analytical estimates.

GYRO-AVERAGED EQUATIONS OF MOTION
Starting with the Lorentz equation for a charged particle in a background magnetic field and a single whistler-mode wave, where p = mvγ is the mechanical momentum, γ is the relativistic factor, B 0 is the local geomagnetic field strength with equatorial value B eq , and E w and B w are the electric and magnetic fields of the wave. Gyro-averaged equations of motion valid near a single resonance have been obtained by many authors, including (Chang and Inan, 1983;Bell, 1984;Albert et al., 2012;Li et al., 2015;Kitahara and Katoh, 2019). For primary resonance (ℓ = −1) between an electron (charge q = −e) and a parallel-propagating whistler wave, equation 3 of Albert et al. (2012) The angle ξ is a combination of wave phase and gyrophase, Ω is the local nonrelativistic electron gyrofrequency eB 0 /mc, and η is the refractive index kc/ω. The standard resonance condition is just dξ/dt = 0, neglecting the term proportional to B w .
Equations 3-9 of Kitahara and Katoh (2019) are very similar after shifting ξ by π/2, using ηE w = B w (in Gaussian units), and accounting for the opposite sign convention in wave phase: These two versions are brought into agreement by invoking the lowest-order resonance condition, which consists of setting the bracketed expression in the equation for dξ/dt to zero. Ginet and Heinemann (1990), Ginet and Albert (1991) used a Hamiltonian version of the equations of motion near resonance with a constant-frequency wave propagating obliquely to a constant background magnetic field B 0 , The Hamiltonan formulation uses canonical momentum P = p + qA/c, where c is the speed of light, and A is the vector potential that describes both B 0 and the wave electromagnetic field. A canonical transformation was made from (x, P x , y, P y , z, P z ) to variables (I, ϕ, X, P X , z, P z ), with z the distance along B 0 in slab geometry. I and ϕ correspond to standard first adiabatic invariant and gyrophase but have modifications proportional to the wave amplitude. After gyro-averaging, and specializing to the case of a parallel-propagating wave, the variables (ϕ, z, t) appeared only in the combination kdz − ωt − ϕ (equation 19 of Ginet and Heinemann (1990) with k x = 0 and sℓ = 1). Albert (1993) generalized the treatment to include slow dependence of Ω and η on z, obtaining the Hamiltonian

TIME-DEPENDENT HAMILTONIAN EQUATIONS
and using normalized variables (ωz/c, ωt, ωI/mc 2 , P z /mc) as in Albert (1993). Appropriate partial derivatives of H give equations of motion for (I, ϕ, P z , z), e.g., dI/dt = −zH/zϕ and dϕ/dt = zH/zI, from which It is also found that dH/dt = zH/zt equals dI/dt, so that I − H is a constant, denoted c 2 : Frontiers in Astronomy and Space Sciences | www.frontiersin.org June 2022 | Volume 9 | Article 910224 2 I − H c 2 .
( 8 ) Following Shklyar (1986), Albert (1993) solved this for P z after approximating H by ϒ, obtaining These can be used to eliminate P z in the equations of motion, giving as a closed set of equations in (I, ξ, z, t). Since P 0 is defined as always positive, explicit minus signs account for the motion of the particle toward the equator. The bracketed expression in the equation for dξ/dt gives the lowest order resonance condition. Retaining the wave term in H to first order gives again allowing P z to be eliminated. The correction to P z /ϒ significantly affects Eq. 7, giving which is also a closed set of equations in (I, ξ, z, t). Ginet and Heinemann (1990) and Ginet and Albert (1991) proceeded to transform to variables (ξ, P ξ , μ, P μ ,φ,Ĩ), with P ξ canonically conjugate to ξ. However, doing so in an inhomogeneous setting reintroduces explicit time dependence in place of z dependence (see equation 68 of Ginet and Albert, 1991). Instead, following Shklyar (1986), Albert (1993) divided the equations for dI/dt and dξ/dt by the equation for dz/dt and attempted to write the results in Hamiltonian form using z as the independent variable. With a Hamiltonian K of the form

POSITION-DEPENDENT HAMILTONIAN EQUATIONS
FIGURE 1 | Evolution of 24 electrons starting at z/R e = 1 and interacting with a whistler mode wave, with particle and wave parameters as given in the text. Red curves show results for the equations of motion given in Eq. 2, and blue curves used Eq. 3.
FIGURE 2 | Evolution of 24 electrons starting at z/R e = 1 and interacting with a whistler mode wave, according to equations of motion based on K(I, ξ, z, t). Blue curves show results for the equations of motion given in Eq. 10, and red curves used Eq. 12.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org June 2022 | Volume 9 | Article 910224 the choice gives dI dz −K 1 cos ξ, which agrees with (dI/dt)/(dz/dt) from Eq. 10. Using then gives dξ dz or, once more using the lowest-order resonance condition, dξ dz It is clear that the first-order term in dξ/dz obtained by this procedure, which enforces the form of Eq. 13, is not the same as that of (dξ/dt)/(dz/dt) from either Eq. 10 or Eq. 12. The analogous disagreement is evident between equations 3.8 and 3.10 of Shklyar (1986), who treated the simpler case of an electrostatic wave and nonrelativistic protons. (Both equations give versions of dξ/ds, the typesetting error in equation 3.8 notwithstanding.)

SIMULATIONS AND DISCUSSION
The consequences of the disagreement in the first-order terms of the various ξ evolution equations is studied here numerically. We choose wave and particle parameters following Kitahara and Katoh (2019); Gan et al. (2020). A Taylor expansion of the geodipole magnetic field about the equator gives the variation along a field line as B/B eq 1 + 4.5z 2 /(LR e ) 2 , with L = 4, where R e is the radius of the Earth and LR e is the field line equatorial crossing distance. The cold electron density is constant, and chosen to give the ratio of plasma frequency to gyrofrequency as f pe /f ce = 4 at the equator. The field-aligned whistler mode wave has frequency such that ω/Ω e = 0.3 at the equator. We consider ensembles of 24 electrons, with energy 20 keV, uniformly distributed in initial gyrophase. We take equatorial pitch angle α 0 = 5°a nd B w /B eq = 3 × 10 -4 (with B w fixed), since this case seems particularly complex, exhibiting a mixture of conventional phase trapping and "anomalous" phase trapping (as opposed to the oppositely directed change associated with phase bunching for larger pitch angles). The particles are launched towards the equator (z = 0) from a distance of 1 R e , and the equations of motion are advanced with a standard Runge-Kutta integrator with variable step size. Figure 1 shows results using Eq. 2 (in red) and Eq. 3 (in blue). The sets of trajectories are not expected to be identical because of accumulated phase differences far from resonance. Nevertheless the overall behavior is very similar, showing no significant change until reaching resonance around z/R e = 0.35, after which the  equatorial pitch angle increases either over a sustained period (conventional phase trapping, PT) or transiently. The long-time behavior of the phase angle ξ is oscillatory for PT but monotonic otherwise. This corresponds to the NL1 regime of Gan et al. (2020), also referred to as positive phase bunching (PPB). Numerically, PT was identified by a change of sign in dξ/dt from one time step to the next after crossing below z/R e = 0.1. Of 24 simulated particles, 10 became PT using either Eq. 2 or Eq. 3. Figure 2 shows results using Eq. 10 (blue) or Eq. 12 (red). The equatorial pitch angle α 0 obtained from the normalized variables (I, z) via B B eq sin 2 α 0 sin 2 α p 2 The behavior turns out to be very similar to the previous run, with 9 instances of PT, leading to α 0 ≈ 25°at z = 0, with the rest of the particles ending up with α 0 spread between about 4°and 14°. Finally, Figure 3 shows results using Eqs. 15, 18. Again the results are very similar in the final α 0 values reached by PT or PPB particles, and in the number of each. The number of PT particles in this run is 8, which does not deviate much from the previous values given the small number (24) of particles in each simulation.
We conclude that the reduced Hamiltonian K(I, ξ, z) of Eq. 13 captures the nature of the particle dynamics, including APT and PPB, with fidelity comparable to the other models. This is propitious because it allows access to a rich body of work on invariant breaking at separatrix crossings (e.g., Cary et al., 1986), enabling both qualitative understanding and quantitative analytical estimates.
Some steps have already been taken in that direction. Figure 4 shows the results of Figure 3 in the (I, ξ) plane, with PT trajectories (identified as above) over the interval 0.4 > z > 0 shown in red, and become limited in ξ while reaching large values of I. The remaining paths, shown in blue (over the interval 0.4 > z > 0.22, for clarity), do not reach such large values of I but are less FIGURE 5 | Contours of K(I, ξ, z), at several values of z shortly before and after resonance crossing, according to motion based on K(I, ξ, z). O-points are shown as diamonds, and X-points (if present) are shown with an X symbol, with the contour through them is in cyan. The red, dashed curve shows the initial value of I.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org June 2022 | Volume 9 | Article 910224 5 restricted in ξ. Figure 5 shows contours of K(I, ξ, z) at several fixed values of z chosen during the trapping process, based on Figure 3. They indicate that at early times (large values of z) there is only a single, O-type fixed point, while an X-point and separatrix, as well as another O-point, form around the time of the trapping process. Contours circling the O-point at ξ = π/2 correspond to the (red) PT trajectories of Figure 4, and PPB trajectories (in blue) are connected to the development of the O-point at low I and ξ = −π/2. Similar contours, developed from Eq. 13 with further approximation, were obtained and studied by Albert et al. (2021), Artemyev et al. (2021). Quantitative analysis of separatrix formation and crossing, invariant breaking, and energy and pitch angle change will be the subject of future work.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
JA conceived this work, and wrote and ran the simulations. AA helped analyze approaches to reducing the dimensionality of the system of equations. WL, QM, and LG consulted on the work and made several useful suggestions on the simulations and presentation.

FUNDING
JA was supported by NASA grant 80NSSC19K0845 and the Space Vehicles Directorate of the Air Force Research Laboratory. WL, LG, and QM also acknowledge the NASA grants 80NSSC20K1506 and 80NSSC20K0698, NSF grant AGS-1847818, and the Alfred P. Sloan Research Fellowship FG-2018-10936.

ACKNOWLEDGMENTS
The views expressed are those of the author and do not reflect the official guidance or position of the United States Government, the Department of Defense or of the United States Air Force. The appearance of external hyperlinks does not constitute endorsement by the United States Department of Defense (DoD) of the linked websites, or the information, products, or services contained therein. The DoD does not exercise any editorial, security, or other control over the information you may find at these locations.