The Earth’s magnetosphere can support various electromagnetic wave modes that play a vital role in near-Earth space dynamics. These electromagnetic waves of which the whistler mode is of particular importance, concurrently interact with higher energy radiation belt particles which are trapped in a magnetic mirror configuration of the geomagnetic field (Bell and Buneman, 1964; Helliwell, 1965; Kennel and Petschek, 1966; Lyons et al., 1972; Gendrin, 1981; Omura et al., 1991; Bortnik et al., 2008). Two primary aspects of whistler mode wave-particle interactions are the amplification of the waves by an unstable hot particle distribution, and the precipitation and/or acceleration of these particles by the waves. In general, modeling wave amplification along with particle scattering and acceleration is a difficult problem that requires a self-consistent solution to the Vlasov-Maxwell system of equations. Several researchers have approached the self-consistent problem using particle-in-cell (PIC) or hybrid simulations (Katoh and Omura, 2007; Gibby et al., 2008; Hikishima et al., 2010; Omura and Nunn, 2011; Wang et al., 2024); however, such computationally intensive self-consistent simulations can be unnecessary in scenarios where particle induced modifications to the wave are small and particle scattering by the waves is of greater interest. In other words, to investigate particle dynamics while minimizing computational costs, the test particle method can be employed, where the feedback of particles on the waves is neglected.

A common approach is to specify the wave-fields and neglect feedback from the particles which results in a simpler but only approximate system of equations. For small amplitude and incoherent signals such as plasmaspheric hiss, quasi-linear diffusion theory is a common method to track the time-evolution of a particle distribution (Kennel and Petschek, 1966; Inan et al., 1978; Abel and Thorne, 1998; Albert, 1999; Summers, 2005). For arbitrary wave-fields (which can go beyond the scope of linear theory), however, the dynamics of the particle distribution may need to be investigated using test particle simulations (Maldonado et al., 2016; Fu et al., 2019) or other methods (Hu and Krommes, 1994). Although test particle simulations have been successful in previous work, simulations with large-amplitude and coherent waves can still require many millions of particles to accurately describe the nonlinear phase-space dynamics of the particle distribution function (Dysthe, 1971; Matsumoto and Omura, 1981; Bell, 1984; Gibby et al., 2008; Nunn, 1974; Albert et al., 2012).

In order to alleviate the computational cost of traditional test-particle methods, a more efficient backward test particle solution to the Vlasov equation is employed here to evaluate the nonlinear effects of large-amplitude coherent waves on the energetic particle distribution function. The backward test particle approach exploits the conservation of phase space as articulated in Louisville’s theorem and permits determination of the prior state of a particle distribution if the wave-fields are known. Practically, this technique allows regions of interest in phase-space to be defined a priori which thus greatly reduces the number of particles that need to be tracked in the simulations and avoids complications from potential under sampling. The method has been previously utilized to efficiently model particle precipitation by coherent whistler mode waves since the loss-cone is well defined (Harid et al., 2014; Nunn and Omura, 2015). In addition, the technique is well-suited for analyzing acceleration to high energies and second-scale variations to the particle distribution function. In this work, the detailed temporal dynamics of the nonlinear wave-induced trap in phase-space are shown at a higher resolution than has been shown in previous works by full PIC (Figure 7 in Hikishima and Omura, 2012) or hybrid (Figure 4.9 in Harid, 2015) simulations. This high resolution allows for the accurate determination of “scattered” fields to investigate salient features of amplified and triggered magnetospheric waves.

The mathematical basis of modeling wave-particle interactions is via the Vlasov-Maxwell system of equations. The Vlasov equation governs the evolution of electron (with mass of m, and charge of q) phase space density f r , v in a collision-free plasma, as: ∂ f ∂ t + v ∂ f ∂ r − q m E w + v × B ∂ f ∂ v = 0 . (1) where the quantities r and v correspond to the position and velocity coordinates of phase-space. E w corresponds to the wave electric field while B is the total magnetic field. The total magnetic field can be decomposed into B = B w + B 0 where B w is the wave magnetic field and B 0 represents the background geomagnetic field. The Earth’s magnetic field retains a dipole shape in the inner magnetosphere, and forces radiation belt electrons to experience helical motion around the background field. The electrons can constantly interact with the circularly polarized whistler mode waves that are propagating along the field line. Here, the waves are assumed to propagate parallel to the magnetic field lines and other (non-whistler mode) wave modes are ignored.

For a circularly polarized whistler mode wave with frequency ω and wavenumber k that is propagating parallel to the background magnetic field lines ( − z direction), the expression for the wave magnetic field , B w , in Cartesian coordinates is: B w = R x ^ − j y ^ B w e j ϕ w + ω t + ∫ k d z (2)

Counter streaming electrons that travel at the appropriate positive velocity will experience an approximately static wave fields and significant energy exchange. This is referred to as Doppler shifted cyclotron resonance or gyro-resonance where resonance velocity ( v r ) is given by Equation 3: v r = ω c γ − ω k (3)

The terms ω c is the electron cyclotron frequency ( ω c = q B 0 m , where B 0 = B 0 ) and γ is known as the Lorentz factor ( γ = 1 + p 2 m 2 c 2 , where c = 1 μ 0 ε 0 ) which is included to take relativistic corrections into account.

Maxwell’s equations govern the evolution of the wave electric and magnetic fields, but the wave equations in a magneto-plasma can be simplified under the assumption of a narrowband modulating wavepacket, given the coherence of the signals observed in the data. The term ω t + ∫ k d z in Equation 2 corresponds to the phase variation of a monochromatic plane wave and can be thought of as a feature of an injected (incident) carrier wave. It is worth mentioning, although the focus of the current study is on a monochromatic wave with a singular frequency, waves with a finite but narrow bandwidth are also called narrowband waves and can lead to nonlinear wave-particle interactions. The quantity B w e j ϕ w corresponds to the complex wavepacket that modulates the carrier whistler wave. Both the amplitude ( B w ) and phase ( ϕ w ) of the modulating wavepacket will be slowly varying functions of position and time, even if the quantity B w e j ϕ w may be varying rapidly. Under the slowly-varying or narrowband assumption the evolution equations for the amplitude and phase of the modulating wavepacket can be simplified as (Nunn, 1974; Matsumoto and Omura, 1981; Omura and Nunn, 2011; Gibby et al., 2008): ∂ ∂ t − v g ∂ ∂ z B w = − μ 0 v g 2 J E (4) ∂ ∂ t − v g ∂ ∂ z ϕ w = − μ 0 v g 2 J B B w (5)

These narrowband wave equations have the advantage of separately quantifying the effect of the wave growth and frequency change and can provide useful physical intuition behind the evolution of a wavepacket that is propagating at the group velocity ( v g ) of the whistler wave. Specifically, Equation 4 shows that the wave amplitude is driven by the component of the resonant current that is in the direction of the wave electric field ( J E ). For a coherent wavepacket, it is assumed that the wave electric field, E w , is given by, E w = v p B w , where v p is the phase-velocity of the wave.

The geometry of the resonant currents ( J R ) with respect to the wave fields due to the energetic electrons are delineated in Figure 1. The variables v ∥ ( p ∥ ) and v ⊥ ( p ⊥ ) are parallel and perpendicular components of electron velocity (momentum) relative to the background geomagnetic field ( B 0 ). The angle between v ⊥ and B w ( − B w ) is referred to as the gyrophase ζ ( φ ). The quantities J E and J B are orthogonal components of resonant currents ( J R = J B + j J B ) and requires computing the appropriate integral directly over the phase space distribution function ( f ) in cylindrical ( v ⊥ , v ∥ , φ ) coordinate coordinate system, which is given by Equations 6, 7, respectively: J E = q ∫ f v ⊥ 2 ⁡ sin ⁡ φ d v ∥ d v ⊥ d φ (6) J B = q ∫ f v ⊥ 2 ⁡ cos ⁡ φ d v ∥ d v ⊥ d φ (7)

The dynamics of resonant electrons in a monochromatic whistler mode wave field ( B w , E w ) immersed in a background inhomogeneous magnetic field ( ∂ ω c ∂ z ), is in general governed by the Lorentz force. d z d t = v ∥ (8) d p ∥ d t = q m γ B w p ⊥ ⁡ sin ⁡ ζ − p ⊥ 2 2 m γ ω c ∂ ω c ∂ z (9) d p ⊥ d t = − q ⁡ sin ⁡ ζ E w + v ∥ B w + p ∥ p ⊥ 2 m γ ω c ∂ ω c ∂ z (10) d ζ d t = k v r − v ∥ − qcos ζ p ⊥ E w + v ∥ B w (11) where p ∥ ( p ⊥ ) are parallel (perpendicular) components of electron momentum relative to the background geomagnetic field. If only electrons that are close to resonance are examined and the small contribution of centripetal acceleration due to the wave is neglected, the equations of motion can be written as Gołkowski et al. (2019): d ζ d t = θ (12) d θ d t = ω t r 2 sin ⁡ ζ + S (13)

Here the variable θ = k v ∥ − v r represents a normalized change of the electron’s parallel velocity from resonance. The quantity ω t r = q k v ⊥ B w m is known as the trapping frequency. The quantity S is called the “S-parameter” and is a collective inhomogeneity factor based on Equation 14 which quantifies the effect of background inhomogeneity ( ∂ ω c ∂ z ) as well as the frequency sweep rate observed by the particle ( d ω d t ): S = 1 ω t r 2 k v ⊥ 2 2 ω c + 3 2 v r ∂ ω c ∂ z + 2 ω + ω c ω d ω d t (14)

Differentiating Equation 12 with respect to time and plugging into Equation 13, results in a nonlinear ordinary differential equation which represents a forced pendulum equation where the forcing term is proportional to S : d 2 ζ d t 2 = ω t r 2 sin ⁡ ζ + S (15)

For S = 0 , Equation 15 takes the form of a conventional pendulum equation and the particle will oscillate around ζ = π at the trapping frequency ω t r in a manner similar to which a pendulum oscillates in a constant gravitational field. For values in the range − 1 < S < 1 , the central phase angle around which the particle oscillates is moved to   ζ 0 = − sin − 1 ⁡ S .

Of particular importance is the formation of a wave-induced trap in phase-space which changes structure depending on the location along the geomagnetic field line. Figure 2 shows the shape of the phase-space trap at several positions (for a monochromatic whistler mode wave) (Gołkowski et al., 2019). The variable ξ on the vertical axis is defined by θ 2 ω t r and is essentially a normalized deviation of v ∥ from resonance.

The trapped trajectories correspond to closed curves in phase-space while the untrapped particles follow open curves. The trapped and untrapped electron populations are separated by a boundary known as a separatrix and are shown by the red contours. Formally, the separatrix can only exist when − 1 < S < 1 . If the background magnetic field were homogenous ( S = 0 ), then the initially trapped electrons stay trapped and untrapped particles stay untrapped indefinitely. However, for a spatially dependent S , this is not necessarily the case. Electrons can either swing around the trap or be trapped for many trapping periods before being detrapped. The exact dynamics are strongly dependent on the initial phase angle, energy, and pitch angle of the electrons. Therefore, by setting S = 1 , the minimum amplitude required for phase trapping, B t r = m q k v ⊥ k v ⊥ 2 2 ω c + 3 2 v ∥ ∂ ω c ∂ z , can be obtained.

The backward test particle numerical model essentially solves the Vlasov Equation 1 for a given wave at a particular location along the geomagnetic field line. Since the Vlasov equation is an advective-type partial differential equation (PDE), information propagates around phase space in a complicated manner. An accurate method of computing the distribution is by using the method of characteristics, which in the context of Vlasov equation is equivalent to Liouville’s theorem. By neglecting the transverse spatial motion of electrons and only considering spatial variation along the field line, the Vlasov Equation 1, can be written as Equation 16: ∂ f ∂ t + d z d t ∂ f ∂ z + d p ∥ d t ∂ f ∂ p ∥ + d p ⊥ d t ∂ f ∂ p ⊥ + d φ d t ∂ f ∂ φ = 0 (16) where the terms in parenthesis corresponds to equations of motion in Equations 8–11. This is done by considering characteristic curves, which are curves along which the distribution function is advected. This turns the PDE into a set of ODEs. More specifically, consider a general advection equation as: ∂ f ∂ t + c → r ∂ f ∂ r = 0 (17)

This type of equation describes advection of the quantity f t , r at “speed” c → at “position” r . To find the characteristics, we find the trajectories, r t for which the total derivative vanishes as shown in Equation 18: d f t , r d t = ∂ f ∂ t + d r d t ∂ f ∂ r = 0 (18)

The original advection Equation 17 can only be satisfied if d r d t = c → is satisfied. Another interpretation is that in the frame of reference moving at speed c → , the quantity f does not change. Thus, by solving for the trajectories, we find the curves along which f is advected (Harid et al., 2014). In the case of the Vlasov equation, the characteristic curves are found by solving Equations 8-11. This means the value of the distribution function at any particular point can be determined by tracing the characteristic curves back until time zero. This method requires a grid generated over ( φ , v ∥ , α ) in phase space at any position of interest, z , within the interaction region. The characteristics are traced backward (formally d t → − d t for the equations of motion and narrowband wave equations) until time zero.

The simulations use an interaction region that is approximately ± 2000  km around the geomagnetic equator. The phase space grid has N v ∥ × N α × N φ = 200 × 50 × 32 = 320000 grid points. All simulations are performed using a centered dipole geomagnetic field model. We use the cold density model from Carpenter and Anderson (1992) to determine the cold plasma parameters under quiet geomagnetic conditions. At L = 4.9 , the equatorial gyrofrequency is f c = 6.8 kHz; the simulations use an input wave frequency of f 0 = f c 2 = 3.4 kHz and a cold plasma density of N c = 400   e l c m 3 . A subtracted bi-Maxwellian distribution of particles is considered for the energetic electrons, which is given by Equation 19: f p ∥ , p ⊥ = C b e − p ∥ 2 2 p t h ∥ 2 e − p ⊥ 2 2 p t h ⊥ 2 − e − p ⊥ 2 2 β p t h ⊥ 2 (19) where the parameter β determines depletion level of the loss cone distribution, such that β = 0 referrs to a pure bi-Maxwellian and larger values of β give a more depleted loss cone. The quantities p ∥ and p ⊥ are the particle momenta parallel and perpendicular to the geomagnetic field, respectively. Additionally, C b = N h 1 − β p t h ∥ p t h ⊥ 2 2 π 3 , where N h is the hot plasma density, and p t h ∥ ( p t h ⊥ ) represents the average hot plasma momenta parallel (perpendicular) to the geomagnetic field. All simulations use p t h ∥ = 0.17 c and p t h ⊥ = 0.4 c where c is the speed of light in free space. Such a distribution represents a Maxwellian distribution without the loss cone particles that would be absent in the magnetosphere. Figure 3A shows the hot electron density as a function of position along the field line. Figures 3B, C shows a color-map of the phase-space subtracted bi-Maxwellian distribution function at a distance −10,000 km and 101 km from the magnetic equator, respectively. As shown, the distribution function is higher at high pitch angles which is characteristic of the radiation belt population. The rapid decay of particle density (away from z = 0) due to the magnetic bottle configuration is clearly demonstrated in Figure 3A and is consistent with the interaction region being dominated by the near-equatorial region.

The wave fields are taken into account by illuminating the entrance of the interaction region with an input signal and using this as a boundary condition for the wave equations. A fourth-order Runge-Kutta (RK4) time-stepping scheme is used to evolve Equations 8–11. The wave equation is time stepped using a semi-Lagrangian scheme with cubic spline interpolation.

For a distribution function of this form, the pitch angle anisotropy is given by A = β p t h ⊥ 2 p t h ∥ 2 + p t h ⊥ 2 p t h ∥ 2 − 1 . For the case of a classic bi-Maxwellian ( β = 0 ), the expression simplifies to A = p t h ⊥ 2 p t h ∥ 2 − 1 = T ⊥ T ∥ − 1 , which readily conveys the concept of anisotropy or directional dependence of the distribution. For a bi-Maxwellian, without the loss cone, if the perpendicular and parallel temperatures (thermal velocities) are equal, the anisotropy is identically zero. Thus, waves are unstable to the plasma if the electron temperature is higher in the direction perpendicular to the magnetic field than parallel to it. For this reason, the term “temperature anisotropy” is used since it provides simple physical intuition behind why the waves are unstable. In the general case of a subtracted bi-Maxwellian, the pitch angle anisotropy will always be higher than that of a classic bi-Maxwellian and one can have anisotropy just from the presence of the loss cone even if thermal temperatures are in equilibrium.

We first use BTP model to explore a small number of particle trajectories as has been done in other works using forward time stepping methods (Inan, 1977; Albert, 2002; Tao et al., 2012) to illustrate some basic physical phenomena. The test particle trajectories will be compared to the linearized equations of motion to emphasize the need for a scattering model which includes both linear and nonlinear effects. Afterwards, the backward test particle model is used to investigate particle scattering and dynamics of the distribution function.

As mentioned before, phase trapping is a nonlinear process and therefore is not in the scope of linear scattering theory. The distribution function is calculated when particles are linearly scattered by the wave. The linearly scattered distribution can be reconstructed by running test particle trajectories backwards in time over linearized equations of motion and employing Liouville’s theorem. The goal was to be able to compare the distribution dynamic in the linear regime with the results when full equations of motion are considered in constructing the distribution function.

Figure 9 shows the Linear scattered particle distribution (top a-e panels) for five different locations (z = 1750, 750, −250, −1,250, and −2,250 km) along the field line at t = 0.30   s along with the nonlinear scattering case (bottom g-k panels) for a monochromatic half a second pulse with 1 pT wave amplitude. As is shown in Figures 9A–E, the linearized equations of motion give the same results as the nonlinear model (Figures 9G–K). This is expected when the wave amplitude shown in Figure 9F is below the phase trapping threshold. In general, since the wave amplitude can grow to values larger than the phase trapping threshold (Bell and Inan, 1981), the scattering model should solve the full equations of motion. There is no formation of phase-space hole and both models show features that can be characterized as linear striations or ripples in phase-space.

Following the previous Figure setup, Figure 10 shows the distribution function for a 20 pT wave amplitude while other simulation parameters remain unchanged. No electron hole or hill shows up in the linear case (panels a–e) since phase trapping cannot be considered within the scope of linear scattering theory (linearized equations of motion).

We developed a backward test particle numerical model that calculates scattering of the energetic electrons distribution by coherent whistler mode waves. This model requires specifying the wave fields a priori and is quite useful at evaluating the effect of waves on the particles in terms of scattering. This study provides important insights into the nonlinear dynamics of wave-particle interactions in Earth’s radiation belts, particularly through the lens of coherent whistler wave instability and whistler mode wave amplification. By utilizing a backward test particle model, we have demonstrated the formation of distinct phase space structures, such as electron phase space holes, in the nonlinear regime, which are not present under linear conditions. These findings underscore the limitations of quasi-linear diffusion theory in capturing the full scope of wave-induced perturbations. The differences observed between the linear and nonlinear perturbations confirm the significant role of nonlinear phase trapping in driving the amplification of whistler mode waves. Although, the presented model is not able to reproduce the effect of the particles on the wave’s amplitude and phase, the formation of nonlinear structures in phase-space can be obtained and analyzed and the expected radiation currents from the perturbed distribution can be evaluated. These results contribute to a more comprehensive understanding of radiation belt dynamics and the complex processes involved in particle acceleration and scattering, with implications for future studies in space weather modeling and the prediction of radiation belt behavior.