Dynamics of Radiation Belt Particles

This paper reviews basic concepts of particle dynamics underlying theoretical aspect of radiation belt modeling and data analysis. We outline the theory of adiabatic invariants of quasiperiodic Hamiltonian systems and derive the invariants of particle motion trapped in the radiation belts. We discuss how the nonlinearity of resonant interaction of particles with small-amplitude plasma waves, ubiquitous across the inner magnetosphere, can make particle motion stochastic. Long-term evolution of a stochastic system can be described by the Fokker-Plank (diffusion) equation. We derive the kinetic equation of particle diffusion in the invariant space and discuss its limitations and associated challenges which need to be addressed in forthcoming radiation belt models and data analysis.


Introduction
The stability of charged particles trapped in Earth's magnetic field was well established by 1960 (e.g., Northrop andTeller 1960) providing a theoretical basis for the existence of radiation belts discovered by pioneering space missions (Van Allen 1959;Vernov et al. 1959). It was shown that in the approximately dipole field of the inner magnetosphere charged particles undergo three quasiperiodic motions each associated with an adiabatic invariant. A set of three invariants defines a stable drift shell encircling Earth. Subsequent experiments revealed that particle intensities across the belts can vary significantly with time (see review by Roederer 1968), which requires violation of one or more of the adiabatic invariants. Theoretical interpretation of the variability of radiation belt intensities was largely inspired by experiments in particle acceleration by random-phased electrostatic waves in synchrocyclotron devices (e.g., Burshtein et al. 1955;Keller and Schmitter 1958) and the development of quasi-linear theory of weak plasma turbulence (e.g., Drummond and Pines 1961; A.Y. Ukhorskiy ( ) · M.I. Sitnov Applied Physics Laboratory, Johns Hopkins University, 11100 Johns Hopkins Rd, Laurel, MD 20723, USA e-mail: ukhorskiy@jhuapl.edu Romanov and Filippov 1961;Vedenov et al. 1961). It was concluded that in the absence of large-amplitude perturbations in the electric and magnetic fields the adiabatic invariants of trapped particles can be violated by waves, which can resonantly interact with the quasiperiodic particle motions. Since both the density and energy density of radiation belt particles is negligible compared to other plasma populations, their motion does not affect the fields that govern it. Thus, in accordance with the quasi-linear theory it was suggested that the evolution of radiation belt intensities can be described as a diffusion in the adiabatic invariants under the action of prescribed wave fields, with the diffusion coefficients determined by resonant wave-particle interactions (see reviews Dungey 1965;Trakhtengerts 1966;Tverskoy 1969). While the diffusion framework of radiation belt particle acceleration and loss was well developed within the first decade after the discovery of the belts (e.g., Falthammar 1965;Kennel and Engelmann 1966), the micro-physical origins of particle diffusion and the limitations of the diffusion framework were not fully realized until the development of nonlinear dynamics in 1980-90s.
The goal of this paper is to review basic physical concepts of particle dynamics underlying theoretical apparatus of radiation belt modeling and data analysis. The review is intended primarily for graduate students and non-experts in radiation bet physics who wish to have a brief yet systematic introduction into the field. The material for this review is based on classical monographs on radiation belt particle dynamics such as Roederer (1970) and Schulz and Lanzerotti (1974), several monographs on nonlinear dynamics including (Lichtenberg and Lieberman 1983;Sagdeev et al. 1988;Zaslavsky 2005), as well as a number of original research papers referenced in the text.
We start with outlining the theory of adiabatic invariants of quasiperiodic Hamiltonian systems, then we discuss the motion of charged particles trapped in a quasi-dipole magnetic field of the inner magnetosphere and derive the adiabatic invariants for each of the three quasiperiodic motions of trapped particles. In Sect. 3 we discuss resonant interaction of particles with small-amplitude regular wave fields. We show that particles at resonance with a given harmonic of the spectrum can be trapped in the wave potential where they undergo nonlinear oscillations and phase mixing. The overlap of particle populations at resonance with adjacent harmonics of the spectrum results in stochasticity of particle motion. In the space of adiabatic invariants particle dynamics then resembles random motion of Brownian particles due to collisions with gas molecules. In Sect. 4 we derive the equation of quasilinear diffusion in the invariant space, often used in radiation belt models, and discuss its relation to the Fokker-Plank kinetic equation of long-term evolution of stochastic systems governed by Markov processes. In Sect. 5 we focus on some limitations underlying the diffusion approximation and associated challenges which need to be addressed in forthcoming radiation belt models and data analysis. In summary we provide a reference table of the most commonly used formulas discussed in this review.

Quasiperiodic Motion and Adiabatic Invariants
Adiabatic invariants are approximate constants of motion of a slowly changing system. The change of an adiabatic invariant approaches zero asymptotically as some physical parameter approaches zero. Adiabatic invariants are of great importance for the analysis of stability of the quasiperiodic particle motion in radiation bets in the presence of small perturbation forces, such as various plasma waves or slow variation of the ambient magnetic field due to changing solar wind and geomagnetic conditions.

General Considerations
Rigorous theory of adiabatic invariants was developed for Hamiltonian systems (e.g., Landau and Lifshitz 1976;Goldstein 1980;Arnold et al. 2010), which in a one-dimensional case are described by equations: where H is the Hamiltonian function, p and q are the canonically conjugate momentum and coordinate variables.
If the Hamiltonian of a system does not depend on time explicitly, then the energy E = H (p, q) is an invariant of motion, i.e. E = const. For an integrable 1 system with periodic motion the action-angle variables are defined by: where the integration is carried out over one period of motion. And by: (4) From the above equations it follows that the action is an integral of motion, which determines its nonlinear frequency ω(I ): Consider now a slow varying one-dimensional system with a Hamiltonian: where the control parameter λ exhibits slow time dependence: where ω is a characteristic frequency of the periodic motion at λ = const. If in the case when λ = const the system is integrable, then the action (2) is an adiabatic invariant of motion. The integration in this case is carried out over one cycle of motion along unperturbed trajectories specified by λ = const.
To examine adiabatic invariance of action we make a canonical transformation from the original variables p and q to the action-angle space (I, θ ) using the generating function from Eq. (3): where momentum p = p(q, I, λ) is specified at λ = const by: Transformation of the Hamiltonian in this case is given by: Since H = H (I, λ) and S = S(q(I, θ), I, λ), the canonical equations of motion in new variables have the following form: where ω(I, λ) = ∂H (I, λ)/∂I is the frequency of periodic motion at λ = const. Since θ is a cyclic variable, the generating function can be expanded in a Fourier series: To estimate the invariant change over long-term evolution of the system we insert the above expansion into the right hand side of the first equation in system (11) and integrate it over time: whereλ∂S k /∂λ = F (εt) is a slow varying function and the phase θ in the exponent is given by: The frequency, which is also a slow varying function ω = ω(εt), usually has zeros only in the complex plane (e.g., Birmingham 1984). To estimate the integrals in expression (13) in this case they are analytically continued over the complex plane. According to the stationary phase method (e.g., Olver 1974) the only non-zero contributions to the integrals occur in the regions of the stationary phase,θ(t 0 ) = 0, given by ω(εt 0 ) = 0 with εt 0 ≡ τ 0 = O(1). The integrals on the right hand side of expression (13) can then be approximated as: 2π According to the above expression the invariant change is determined by the imaginary part of the phase at the stationary point: whereω is the average frequency and τ 0i is the imaginary part of τ 0 . Consequently the invariant change can be approximated as: This means that the average (over many periods of fast oscillatory motion) change of the adiabatic invariant is exponentially small, i.e. approaches zero faster than any power of ε. As was proved by Kruskal (1962) the action integral is an adiabatic invariant (i.e. is conserved to all orders in ε) of any Hamiltonian system with quasiperiodic solutions. It has to be noted that the change in the invariant is no longer exponentially small if the oscillation frequency goes to zero,ω ∼ ε. This corresponds to the case when the system crosses a phase space separatrix, which can result in non-negligible violation of the invariant (e.g., Cary et al. 1986;Neishtadt 1986). Implications of invariant violation at separatrix crossings to dynamics of the outer belt particles is discussed in some detail in Sect. 5.3, until then we assume that the estimate (17) holds.

Adiabatic Invariants of Radiation Belt Particles
The motion of a charged particle in time-varying electromagnetic field of the inner magnetosphere is described by the Hamiltonian function (e.g., Landau and Lifshitz 1959): where e is the electric charge of the particle, m is its mass, c is the speed of light, A is the vector potential of the magnetic field, ϕ is the electrostatic potential, and P and r are the canonically conjugate variables: where γ is the relativistic factor, and v is the particle's velocity. We assume charge neutrality, so the electromagnetic field is given by: Canonical equations of motion can then be written in the familiar Lorentz form: In a static quasi-dipole field of the inner magnetosphere the quasiperiodic motion of trapped particles (Fig. 1a) is a superposition of three independent motions (e.g., Northrop Fig. 1 Quasiperiodic motion of charged particles in the inner magnetosphere. Panel (a) shows a trajectory of a 1 MeV particle with electron charge and the mass of 20 m e (necessary to resolve the gyromotion) over one drift period around Earth. Three components of the quasiperiodic motions are illustrated by subsequent panels. Panel (b): particle gyration in a homogeneous magnetic field. Panel (c): bounce gyrocenter motion along the field lines. Panel (d): drift-bounce gyrocenter motion around Earth, computed for a 1 MeV electron and Teller 1960): (1) particle gyromotion about its guiding center (Fig. 1b), also referred to as the Larmor motion; (2) the bounce motion of particle guiding center along the field lines between the conjugate reflection points in the northern and southern hemispheres (Fig. 1c); and (3) the longitudinal gradient-curvature drift motion of particle guiding center around Earth (Fig. 1d). Each motion is associated with its own adiabatic invariant.
The first adiabatic invariant is associated with the particle gyromotion. If the characteristic time of field variations (τ ) is slow compared to the particle gyration period (T = 2π/Ω): τ T , and the spatial scales of field variations (L) greatly exceed the Larmor radius (ρ): L ρ, particle gyration ρ (R, t, ψ) can be separated from the motion of its guiding center R(t), which can be considered independent of this gyration: where ψ is the gyration phase: dψ/dt = Ω/γ . The first adiabatic invariant I 1 , in this case, can be estimated from expansion of field quantities about the particle guiding center. After expressing an element of Larmor orbit as: where p ⊥ is the momentum associated with the gyromotion, and expanding the magnetic field vector potential in (19) in Taylor series up to the first order about the guiding center, we obtain the following estimate for the action integral: where · · · = 1 2π dψ · · · . To estimate the second term on the right hand side of (24), we use the identity p ⊥ = mΩρ ×b, whereb = B/B, which yields: After substituting the above expression into (24), and using ρ = p ⊥ /mΩ for the Larmor radius, we obtain the following estimate for the first adiabatic invariant: A charged particle gyrating in strong magnetic field is equivalent to a closed loop of electric current j = eΩ/2πγ with the area πρ 2 . It therefore has a magnetic moment equal to: Thus, the first adiabatic invariant of a guiding-center particle is related to its magnetic moment; they become equal in the limit of particle velocities much smaller than the speed of light (γ = 1). It has to be noted that μ is a constant of motion only in the guiding-center approximation. If the guiding-center approximation does not hold, i.e. when the magnetic field changes over one gyroperiod become non-negligent, particle orbits can be considerably more complex than in a stably trapped example shown in Fig. 1. In stretched magnetic field configurations, such as in the magnetotail, MeV electrons and keV ions exhibit complex trajectories illustrated in Fig. 2, which shows a Speiser orbit (e.g., Speiser 1965) of an energetic proton in the magnetotail. While the first adiabatic invariant still exists for such orbits, as long as their motion remain quasiperiodic, it is no longer related to μ or the magnetic moment (e.g., Büchner and Zelenyi 1989).
Even in the guiding-center approximation the first invariant is a non-local quantity, which can significantly complicate its derivation from in situ particle measurements. While the perpendicular momentum in expression (26) is defined at the location of particle measurements, the magnetic field intensity has to be estimated at the gyrocenter, not sampled by the measurement.
The guiding-center approximation holds well for particles inside the electron and the inner proton belts. The description of particle motion, in this case, can be significantly simplified. Equations for the guiding-center motion were originally derived by expanding the Lorentz equation of motion about the guiding center and then removing fast oscillations by the gyrophase averaging (e.g., Landau and Lifshitz 1959;Sivukhin 1965;Northrop 1963). This procedure, however, has two fundamental shortcomings. First, the equations obtained by the gyrophase averaging do not have the Hamiltonian structure of the original Lorentz equation (e.g., Balescu 1988). As a result they do not conserve the phase space density and therefore are in violation of the Liouville's theorem. Consequently these equations cannot be used for the description of collective phenomena in plasmas. Second,  Tsyganenko (1996) magnetic field model at P dyn = 4 nPa and zero tilt angle from the initial location r = (−15, −5, 0) R E and the equatorial pitch angle of 87°t he obtained equations do not conserve energy in time-independent fields. Nonconservation appears in second order terms in the Larmor radius expansion (e.g., Cary and Brizard 2009) and can present difficulties in modeling long-term effects in particle dynamics.
Both problems are successfully solved in a Hamiltonian theory of the guiding-center motion (see review, Cary and Brizard 2009). The six dimensional guiding-center phase space is given by (R, ψ, p , μ). In the absence of potential electric field and when E B (which is generally true for the inner magnetosphere) and u E v ⊥ (where u E is the E × B drift velocity: u E = cE ×b/B), relativistic guiding-center Hamiltonian function can be written as: H (R, p , μ, t) = mc 2 γ (R, p , μ; t) The following noncanonical guiding-center equations of motion are then derived from the variational principle with the use of the above Hamiltonian function (Cary and Brizard 2009;Ukhorskiy et al. 2011): where the effective electromagnetic field: is defined in terms of the effective electromagnetic potentials: In the absence of large electric currents ∇ × B/B 0 and B * = B + cp eBb × ∇B. For a static magnetic field equations (29) are then reduced to: The first equation in (32), which describes the motion of particle guiding center along the magnetic field lines, can be written as: where s measures the distance along filed lines from the magnetic equator (minimum of B(s) in a dipole-like magnetic field). From conservation of kinetic energy and the first adiabatic invariant it follows that particle motion along magnetic field lines also satisfies: where α is the particle pitch angle: sin α = p ⊥ /p. From Eq. (34) it follows that the particle pitch angle, increases while it moves along a field line from the equator to higher latitudes, where the magnetic field intensity is higher. If at the equator a particle had the pitch-angle α eq , then the parallel component of its velocity will become zero (α = π/2) at the point s m : B(s m ) = B(0)/ sin 2 α eq , and the particle will get reflected back towards the equator. To demonstrate that Eq. (33) describes particle oscillations between the conjugate reflection points, consider particle motion in the vicinity of the magnetic equator, where the field can be approximated by the first two non-zero terms of the Taylor expansion: After substituting expression (35) into Eq. (33) we obtain a harmonic oscillator equation: with ω b angular frequency of particle oscillations across the equator. The second adiabatic invariant is computed by integrating the parallel component of particle canonical momentum along a guiding-center bounce orbit: where following common convention we multiplied definition (8) by 2π . The integral in the second term on the right hand side is the total magnetic flux enclosed by the unperturbed orbit, which, in this case, is zero since bounce motion is along magnetic field lines. Since the second invariant is calculated for a fixed magnetic field configuration, the kinetic energy is conserved and the invariant can be expressed as: where the integration is carried out between the conjugate bounce points s m and s m along a fixed magnetic field line. It has to be noted that some textbooks such as Roederer (1970), for example, and papers use different notation: J for I 2 and I for J , which we chose not to follow to avoid confusion with other notations in the paper. The third equation of (32) describes the longitudinal guiding-center drift across the magnetic field lines around Earth, which is referred to as the gradient-curvature drift. This motion corresponds to the third adiabatic invariant: where Φ is the magnetic flux across the drift orbit. The above expression is dominated by the second term, since the gradient-curvature drift velocity is small p ⊥ /mγ U D . The magnetic flux can be computed by shifting the integration contour to any contour C on the drift-bounce surface closed around Earth: dl is the contour element, and dS is the element of a surface encircled by the contour. In the axisymmetrical dipole magnetic field particles gradient-curvature drift along L = const surfaces, where L is a dipole coordinate which is constant along a given magnetic field line, and is equal to the distance from the dipole center to the field line at the equator measured in Earth radii. The magnetic flux through a drift orbit in this case is equal to: where B 0 31000 nT is the magnetic field intensity on Earth's surface at the equator, and R E 6380 km is the Earth's radius. In a realistic nondipolar field of the inner magnetosphere the third invariant of trapped particles is often quantified with the use of the generalized L-value or L-star based on relation (41) for a dipole field: Physically L * is the radial distance (in Earth radii) to the equatorial points of the drift-bounce shell on which the particle would be found if all nondipolar contributions to the magnetic field would be adiabatically turned off. Energy and spatial dependencies of the gyration, bounce, and the gradient-curvature drift frequencies in a dipole magnetic field is shown in Fig. 3. Characteristic frequencies of relativistic (∼1 MeV) electrons in the center of the outer belt (L ∼ 4-5) are separated by approximately three orders of magnitude: gyration frequency ∼kHz, bounce frequency ∼Hz, Fig. 3 Contours of constant adiabatic gyration, bounce, and gradient-curvature drift frequencies of equatorially mirroring electrons in a dipole field (adapted from Schulz and Lanzerotti 1974) while the gradient-curvature drift frequency ∼mHz. This means that in a dipole or an approximately dipole field the quasiperiodic motions corresponding to the three adiabatic invariants are well decoupled. The gradient-curvature drift does not affect the bounce motion, which does not alter the gyration. It also means that it is possible to violate higher invariants without changing lower invariants. For instance, in the process of resonant wave-particle interaction with the gradient-curvature drift motion an ultra-low frequency (ULF) wave can violate the third invariant without altering either the first or the second invariants.

Particle Dynamics in Wave Fields
In this section we discuss resonant interaction of particles with small amplitude waves. We show that even if wave fields are regular and no external randomness is introduced into the system, the nonlinearity of resonant wave-particle interaction combined with the overlap of particle populations in resonance with adjacent harmonics of the wave spectrum result in a stochastic particle motion. In the space of adiabatic invariants particles exhibit random walk motion similar to the Brownian motion of heavy particles due to collisions with light molecules in gasses. In our consideration we use a specific example of resonant interaction between the drift motion of the outer belt electrons and the ULF waves, resulting in a stochastic radial motion of electrons across the drift shells. However, the discussed properties of the stochastic motion are general and are equally applicable to resonant interaction of waves with the bounce and the gyromotion of trapped particles.
Consider an electron bouncing at the magnetic equator (I 2 = 0) and drifting around Earth due to the gradient of its dipole magnetic field, which in the equatorial plane is given by: B(L) =ẑB 0 /L 3 . According to Eqs. (32), U D = 3φμc/γ eR E L, and the unperturbed motion of the electron guiding center is the rotation around Earth at constant L: where ϕ is the azimuthal angle andφ is the unitary vector in the azimuthal direction. For relativistic particles with γ 2 1, the first invariant can be approximated by: μ = mc 2 γ 2 L 3 /2B 0 , and the drift frequency can then be written as function of L: ω D (L) = aL −5/2 , where a includes all constant terms.
Consider now a perturbation of the periodic drift motion due to a small-amplitude azimuthal ULF wave field of the following form: where E 0 is the wave amplitude, ω is the frequency spacing between adjacent spectral harmonics, M is the number of harmonics, and ψ m are their phase shits. In the wave field, the particle guiding center will also experience radial motion due to the After introducing a new variable I = L −2 , assuming that the change I = I − I 0 is small, and approximating the frequency as: ω D (I ) = aI 5/4 ω 0 + ω D I , where ω 0 = ω D (I 0 ) and ω D = 5ω 0 /4I 0 , the above equations can be written as: Finally, after the following substitutions: dynamical system (46) can be written as: where M 1 + M 2 + 1 = M. The above system corresponds to the following Hamiltonian function: While the above derivation was carried out for electron interaction with azimuthal ULF waves (see also Elkington et al. 1999Elkington et al. , 2003Ukhorskiy et al. 2005;, the final form of Hamiltonian function (49) is quite general and describes a wide class of wave-particle interactions including interactions with particle gyration and the bounce motion (e.g., Southwood et al. 1969;Smith and Kaufman 1978;Jaekel and Schlickeiser 1992;Shklyar and Matsumoto 2009). For a particle with a given value of I , the sum on the right hand side of the first equation in (48) is dominated by the term m 0 closest to the resonance: |I − m 0 | → min. Neglecting for now contributions of non-resonant terms and introducing a new angle variable, Ψ = θ − m 0 t + ψ m0 + π , we obtain the following equation for particle oscillations in resonance vicinity:Ψ where Ω 2 NL = K/4π 2 is the frequency of nonlinear oscillations. The above equation describes dynamics of a nonlinear pendulum with the Hamiltonian: where δI = I − m 0 . Its phase portrait is shown in Fig. 4. Singular points of the motion are defined by the conditions:Ψ = 0 and d dt δI = −∂V /∂Ψ = 0, i.e.: This yields:Ψ s = 0, Ψ s = πn, n = 0, ±1, . . . . At singular points the velocityΨ s is zero, and the potential V (Ψ s ) has minima (even n) or maxima (odd n). Thus the singularities are of the elliptic type at n = 2k and saddles at n = 2k + 1, where k = 0, ±1, . . . . The system has two different types of solutions. When H < Ω 2 NL the solutions correspond to particles trapped in the potential well of the wave and oscillating about the elliptical singular points (blue trajectory in Fig. 4). At H > Ω 2 NL the solutions correspond to untrapped particles with unbounded trajectories (yellow trajectory in Fig. 4). The trajectory separating the phase space regions corresponding to the different types of solutions is called the separatrix. It passes through the pointΨ = 0, Ψ = π and therefore corresponds to H s = Ω 2 NL .  51) it follows that the maximum separatrix width in action space, which is also referred to as the resonance width, is equal to: It specifies the maximum deviation of I from the resonance value m 0 at which particles can still be trapped at resonance with a wave of the amplitude K/4π 2 . Expression (53) also gives the upper limit for the change in action due to nonlinear resonance with one wave harmonic. The period of trapped particle oscillations in the potential well of the wave depends on the proximity to the resonance. Close to the resonance, i.e. at small δI (Ψ = 0) values, particles exhibit harmonic oscillations with the period T 0 = 2π/Ω NL , as directly follows from Eq. (50) for small Ψ values. With increase in δI (Ψ = 0), oscillations become nonlinear and their period grows. As trajectories approach the separatrix the period goes to infinity logarithmically (see panel (a) in Fig. 5). To elucidate it, let us consider a trapped particle oscillating in close proximity of the separatrix such that: /H s 1, = H s − H . From Eq. (51) it follows that the action integral of the trapped particle motion is given by: The oscillation period can then be computed as: As H approaches H s = Ω 2 NL , the turning points (Ψ = 0) approach Ψ = ±π and the expression in the denominator on the right hand side of expression (55) goes to zero. Thus, the largest contributions to the integral in (55) is given by vicinities of the turning points. To evaluate the integral, we can therefore use Taylor expansion of the potential function in the denominator about the reflection points. After introducing a new variable: x = Ψ − π , and expanding up to the first non-vanishing, non-constant term we obtain: where the integration goes from the reflection point: x 2 = √ 2 /H s , to the center of the resonance island: x 1 = π . After integrating and taking the limit → 0 we obtain the following expression for the oscillation period: which goes to infinity as the trajectory approaches the separatrix. As a result of nonlinear dependence of the oscillation period on the resonance proximity, trapped particles undergo phase mixing. Particles who originally had the same phase but slightly different values of δI oscillate about the center of the resonant island at different frequencies. Consequently their phases gradually separate and the motion becomes eventually uncorrelated. Phase mixing is illustrated in Fig. 5. Panels (b)-(g) of the figure show evolution of the phase distribution function computed for an ensemble of 10 5 particles trapped at resonance with a single wave. Initially all particles have the same phase (Ψ = 0) but were evenly distributed over δI from −2Ω NL to 2Ω NL . After first several oscillation periods, T 0 , particles spread over almost entire phase range from −π to π , but still exhibit strong phase bunching indicated by multiple pronounced peaks in the distribution function. Eventually the peaks break down and spread over the phase interval. After a large number of drift periods the process produces a smooth distribution function accept for a singular peak around Ψ = 0 corresponding to the particles at exact resonance Ψ =Ψ = 0.

Resonance Overlap
Wave perturbation in the Hamiltonian (49) consists of multiple wave harmonics. From our previous considerations it follows that there is a layer of trapped particles of the width I centered at each of the spectral harmonics. The larger is the wave amplitude K the wider are these layers (see Eq. (53)). If the amplitude increases to the point when the resonant width I exceeds the spacing ω between spectral harmonics m 0 and m 0 ± 1, the resonant populations trapped by adjacent harmonics overlap, and the particle motion is no longer bounded to a single resonance. From this qualitative consideration it follows that resonances overlap, if K π/2 (Chirikov 1960). Particle population initially at resonance with one wave harmonic can then spread over the entire system (maximum spread restricted only by the width of the spectrum: I = M). Phase mixing in this case results in exponential divergence of particle trajectories with similar initial condition, which is an attribute of chaotic dynamics. Generally speaking, chaotic systems are the systems described by regular dynamical equations (the Lorentz equations (21) in this case) with no stochastic coefficients, but at the same time with solutions that are similar to some stochastic processes.
Transition to stochasticity due to resonance overlap is best illustrated with a special case of a regular broad-band wave field, i.e. when all phase shits in (48) are zero ψ m = 0 and M 1,2 → ∞. The following identity: with ν = 1 can then be applied to the right hand side of the first equation (48), which yields: It means that at time moments t m = 2πm particle experience sharp kicks in action I , while between t m and t m+1 the action is conserved and particles move at constant angular frequencẏ θ = I = const. After defining I m = I (t m − 0) and θ m = θ(t m − 0) Eqs. (59) can be reduced to an algebraic map: where the change in action was computed by integrating over a small vicinity of t m . To illustrate the dynamics described by map (60), which is commonly referred to as the standard map, its phase portrait was computed for different values of the nonlinearity parameter K (see Fig. 6). At K = 0.5 most of phase space trajectories are stable. Primary resonant islands are centered at I = 0, 1 and θ = π . The space between the primary islands has a complex structure. It is populated by chains of smaller islands of various sizes and periodicities associated with higher-order resonances, such as the resonances between trapped particle oscillations and the wave field: kω(I ) − m ω = 0, where ω(I ) is the frequency of particle motion about a primary resonant island. Each island chain is bounded by a separatrix. The area near a separatrix is most susceptible for the onset of chaos, since particle velocity near saddle points approaches zero and their dynamics becomes sensitive to small perturbations. An arbitrary small periodic perturbation destroys the separatrix and creates a stochastic layer of extremely complicated phase space topology where chaotic trajectories are embedded with infinite number of islands (Figs. 6c, d). Even at K 1 there is a thin stochastic layer in the vicinity of each separatrix. With an increase in K the width of stochastic layers grows until stochastic layers connect across all values of I , which results in transition to global stochasticity. A detailed analysis which takes into account higher order resonances (Chirikov 1979) shows that transition to global stochasticity in the standard map corresponds to K = 1, which is known as the modified Chirikov criterion (Fig. 6b). With further increase K the area of stochastic phase space region keeps growing, while the area of stable islands keeps shrinking (Fig. 6c).

Transition to Kinetic Description
The motion of individual charged particles is described by the Lorentz equations of motion (21). Particle distribution function or the phase space density evolves in accordance with the Liouville's theorem. If particle motion becomes stochastic, which in a collision-free plasma can be caused by interactions with waves, then correlations among dynamics of individual particles decay. Consequently, the description of long-term evolution of particle distribution function can be reduced to a Fokker-Planck equation, similar to the description of diffusion in gas, which we discuss in this section.

Phase Space Density
For a system of large number of particles one can introduce the phase space density: where N is the total number of particles, and f (p, q, t)dpdq is the number of particles in the volume dpdq centered at z = (p, q) at time t . If particles are not lost or introduced into the system, their evolution in phase-space satisfies the continuity equation: If p and q are the canonically conjugate momentum and coordinate variables corresponding to a Hamiltonian function (Eq. (1)), then ∇ z ·ż = 0, and the last equation can be written as: which is known as the Liouville's theorem stating that the phase space density is conserved along particle trajectories. The phase space density is directly related to observable quantities such as particle flux and intensity. The intensity j α (E, r) of particles of a given class and kinetic energy is defined as the number of particles coming from a given direction which impinge per unit time, unit solid angle and unit energy, on a surface of unit area oriented perpendicular to their direction of incidence. If δN is the number of particles with kinetic energies between E and E + δE impinging on the area δA with normaln during time interval δt, and whose direction of incidence lie in the solid angle δΩ oriented along p (see Fig. 7), then: δN = j α (E, r)δA cos αδΩδEδt.
After comparing expressions (64) and (65) we obtain: Kinetic energy and momentum of a relativistic particle are related as: from which we obtain that mγ dE = pdp. Therefore the intensity and the phase space density are related as:

Diffusion in Action Space
In the action-angle coordinates the Liouville's equation can be written as: To illustrate derivation of the diffusion equation we as previously use a one-dimensional example. Consider evolution of plasma with small-amplitude waves described by the following Hamiltonian: where H 0 corresponds to the unperturbed system without waves, V is the perturbation due to waves, and the dimensionless parameter ε indicates that wave amplitudes are small. The corresponding equations of motion are: where ω 0 is the frequency of the unperturbed motion.
The above system evolves at two characteristic time scales. It exhibits rapid oscillations in the angle variable θ and slow change in the action I due to resonant wave-particle interactions. An ensemble of particles with initially same values of I but distributed in θ will initially rotate coherently, sinceθ ω(I ). However, if the system is stochastic, then I of different particles will undergo different small variations due to their interactions with the wave field (see previous section). After some time the ensemble will spread in I and will rotate at different frequencies. Consequently the ensemble will exhibit phase mixing, i.e. correlations between particle θ(t) and its initial values will decay, and eventually particle distribution in I will become independent of the initial distribution in θ . On timescales longer than the phase correlation decay time (τ c ), it then become possible to derive a reduced description of long-term evolution of the system by averaging the Liouville's equation over the fast angular variable θ .
We start with expanding the distribution function up to first order in ε: where in accordance to quasi-linear theory the ε 2 factor insures that slow variations in the distribution function F appear only as a second order term: dF /dt = O(ε 2 ). We then expand the Liouville's equation up to second order and find solutions order by order. Assuming that the wave frequency is related to the wave number by dispersion relationship of the plasma (ω m = ω(m)) and f m ∼ e −i(ωmt−mθ) we find in first order: The angular dependence appears in second order, which we remove by averaging: With the use of integration by parts we obtain: After inserting expression (73) we obtain the diffusion equation for slow evolution of phase averaged distribution function: with the quasi-linear diffusion coefficient: where we used the identity: Im(ω m − mω 0 ) −1 = iπδ(x − x 0 ), which indicates that quasilinear diffusion is produced by resonant wave-particle interactions. Alternatively, for systems described by maps, such as standard map (60) discussed in the previous section, the kinetic equation is based on the Fokker-Plank equation of Markov processes. It is assumed that on time scales longer than the correlation decay time τ c the transition probability W t (I − I, t, I, T ), which is the probability that an ensemble of phase points having an action I at time t suffers an increment in action I after a time T , does not depend on the angle variable and that: Long-term evolution of the distribution function F is then described by the Fokker-Plank equation (see for example Lichtenberg and Lieberman 1983): with the coefficients defined as: where T is a characteristic time scale, such that T > τ c . Derivation of (80) also implies that moments I i I j I k and higher are 0. Coefficients of the Fokker-Planck equation are related as: 1 2 Following Landau (1937) let us show it for a one dimensional system by expanding the change in action I up to the second order in time: where:İ With the use of Hamiltonian equations and after regrouping terms we obtain: After averaging over the angle variable θ we obtain: all other terms vanish since H is a periodic function of θ . Similarly, for the second moment we obtain: By comparing (86) and (87) we find the relation (82) in a one dimensional case.
The Fokker-Plank equation (80) can therefore be recast in the form of quasi-linear diffusion equation (77): where the diffusion coefficient D ij = B ij /2. Trapped particle motion in the inner and the outer radiation belts often conserve both the first and the second adiabatic invariants. Equation (88), in this case, reduces to a one dimensional equation in the third invariant, which is often written in terms of L * . Recalling that I 3 ∝ Φ ∝ 1/L * (see Eq. (42)), we obtain: Since L * has the physical meaning of the dimensionless distance to the equatorial points of the drift-bounce shell particles would have in a dipole field, written in this form the diffusion equation in the third invariant is known as the radial diffusion equation. In radiation belt models, Eq. (88) is usually written in terms of L * , energy, and pitch-angle. When all three invariants are violated, the diffusion equation cast in these variables can have a complicated structure, even if additional assumptions are made, such as that the pitch-angle and energy diffusion are uncoupled from radial transport and that the cross-diffusion terms in energy and pitch angles can be neglected (for a detailed discussion see for example Schulz and Lanzerotti 1974). Let us go back to our example Hamiltonian (49), which we derived for radial transport of radiation belt electrons due to interaction with ULF waves. In the case of a standard map (ψ m = 0) the correlation decay time for K 1 can be estimated as (e.g., Zaslavsky 2002): where the time step T of the map defines the characteristic time scales of changes in the action variable due to the resonant wave-particle interaction. It has physical meaning similar to the time interval between random weak collisions experienced by a heavy Brownian particle in a gas. The fact that individual steps of the map are small, i.e. I /I 1, makes this analogy even closer.
From estimate (90) it follows that for large K the diffusion coefficient of map (60) can be estimated from a single step, I = K 2π sin θ , T = 2π . This yields: which is often referred to as the quasi-linear estimate because it does not include higherorder corrections due to subsequent steps (see for example, Lichtenberg and Lieberman 1983). It can be expected that the deviations of the diffusion coefficient D QL (K) from this one-step estimate is the highest for moderate values of K, when according to Eq. (90) it takes more than one effective collision to randomize particle phases. Additionally, in the standard map case, there are islands of regular particle motion embedded into stochastic regions of phase space at any finite value of K, where particle trajectories are stable. Particle trajectories can be trapped in vicinity of the boundary between the stochastic and regular phase space regions, where the action variable changes almost linearly causing deviations of Diffusion coefficient for a system with Hamiltonian (55) with random shifts among frequency harmonics (ψ m ∈ [0, 2π)). The nonlinearity parameter ε corresponds to K in our notations (Cary et al. 1990) transport from pure diffusion (e.g., Zaslavsky 2002). Rechester and White (1980) analyzed transport properties in the standard map at different values of K by calculating the diffusion coefficient both numerically and analytically including higher-order correlations due to finite correlation decay time. Their results are shown in Fig. 8a. The diffusion coefficient oscillates about its quasi-linear value with the maximum value exceeding D QL by more than factor of 2, and minimum around D QL /2. Since the correlation time becomes shorter and the phase space area occupied by stable islands decreases, deviations become smaller with increase in K.
The existence of random phase shits ψ m ∈ [0, 2π) between different harmonics of the wave spectrum in Eq. (49), considerably changes the dynamic properties of the system. If resonances overlap (K > 1), the islands of regular motion are completely destroyed by random shifts, and the system becomes stochastic everywhere across the phase space. The analysis of particle transport at different values of K (Cary et al. 1990;Helander and Kjellberg 1994) show that the diffusion coefficient in this case can still exhibit large deviations from the quasi-linear value, D QL : while it never gets smaller than D QL , at K 18 it reaches the maximum of 2.3D QL (see Fig. 8b).
In reality, additional stochasticity may be introduced into the system due to random nature of the wave fields. Phase shifts ψ m at different harmonics of the spectrum in Eq. (49) may no longer be stationary in this case. Their values can change in some characteristic time intervals T corresponding to the spatial or temporal coherence of the problem. For instance, variations of the solar wind dynamic pressure is one of the dominant drivers of the ULF waves in the inner magnetosphere (e.g., Takahashi and Ukhorskiy 2007). ULF waves can violate the third adiabatic invariant of trapped electrons in the process of resonant interaction with their drift-bounce motion discussed in Sect. 3. Oscillations in dynamic pressure are attributed to the Alfvén turbulence in the solar wind. The phase shifts between different harmonics of the ULF wave spectrum therefore change on the time scales of the autocorrelation time of the solar wind turbulence, ∼3 hr (e.g., Jokipii and Coleman 1968).
Electromagnetic ion cyclotron (EMIC) waves are considered to be one of the dominant local mechanism of electron losses from the outer radiation belt (e.g., Thorne and Kennel 1971;Horne and Thorne 1998;Summers et al. 1998;. Resonant interaction of EMIC waves with electron gyromotion breaks the first adiabatic invariant and can cause electron scattering into the atmospheric loss cone and their subsequent loss via precipitation. Free energy for the EMIC wave growth is supplied by the positive temperature anisotropy (T ⊥ > T ) of energetic (∼10-100 keV) ions (e.g., Cornwall 1965;Kennel and Petscheck 1966). EMIC waves grow to observable amplitudes at frequencies of maximum growth rate out of small-amplitude electromagnetic noise propagating along the field lines through the regions of positive anisotropy (e.g., Gomberoff and Neira 1983;Horne and Thorne 1994). EMIC wave activity can extend over >10°about the magnetic equator (e.g., Erlandson and Ukhorskiy 2001) and last for tens of minutes producing pitch-angle scattering of radiation belt electrons over many bounce periods. Detailed spectral analysis  revealed that wave events consist of many short (∼30 sec) wave packets. Consequently phase shits among the harmonics of EMIC spectra vary at time scales comparable to the duration of individual wave packets.
Numerical simulations  showed that if additional extrinsic stochasticity is introduced into the system by varying phase shifts ψ m among spectral harmonics of the wave perturbation (49) at time intervals comparable to the time T between effective collisions (90), then particle motion becomes stochastic even if resonances do not overlap. The diffusion coefficient in this case agrees well with its quasi-linear estimate (91). At the time scales longer than the correlation decay time τ c the system can then be described by diffusion equation (88) with quasi-linear diffusion coefficients. The correlation decay time τ c in this case depends on both the collision time and wave amplitude similar to expression (90).

Limitations and Challenges
During over five decades since the discovery of radiation belts the concept of diffusion in the invariant space has been successfully applied for the analysis of transport, acceleration, and loss of radiation belt particles. Radial diffusion due to drift-resonant interaction with solar-wind driven ULF fields was the first mechanism proposed to explain acceleration of electrons and protons in radiation belts (Kellog 1959;Tverskoy 1964;Dungey 1965;Falthammar 1965). Subsequent analysis showed that radial diffusion causes not only acceleration but loss of particles from the outer belt (e.g., Bortnik et al. 2006;) and can be driven by variety of plasma waves including waves excited internally by instabilities in ring current plasma such as stormtime Pc5 waves (Lanzerotti et al. 1969;Ukhorskiy et al. 2009). Local resonant interaction of electron gyromotion with whistler waves was initially considered to be primarily responsible for electron losses from the belts (Dungey 1963;Cornwall 1964;Kennel and Petscheck 1966). Local wave-particle interactions are now recognized as both loss and acceleration mechanisms. As was mentioned in the previous section EMIC waves are considered to be one of the primary mechanism of local losses outside of the plasmasphere. Whistler chorus (e.g., Horne and Thorne 1998;Summers et al. 1998) and magnetosonic (e.g., ) waves were identified as mechanisms of local acceleration of radiation belt electron, more efficient than energization due to radial diffusion (e.g., Horne 2007). A number of recent review papers Shprits et al. 2008aShprits et al. , 2008bThorne 2010) provide detailed discussions and reference lists on diffusion theory of the radiation belt processes. In this section we discuss to what extent particle motion in the belt can be described in terms of three adiabatic invariants, some limitations of the diffusion approximation and associated challenges which need to be addressed in forthcoming radiation belt models and data analysis.
Diffusion approximation applies to the situations when in zero order radiation belt particles are stably trapped in quasiperiodic motion associated with three adiabatic invariants. This implies that the magnetic field has a slow-varying quasi-dipole configuration, such that the time scales of the three periodic motions are well separated, and the electric field is small, such that the E × B drift is negligible compared to the gradient-curvature drift. In this case particle invariants can be violated only in the process of resonant wave-particle interaction. Reducing the description from the full Vlasov equation to a Fokker-Plank equation in the invariant space also requires that waves have small enough amplitudes, such that nonlinear phase-dependent effects can be neglected, and the characteristic time scales of the described processes are longer than the phase correlation decay time. Variability of radiation belt intensities do not always satisfy these conditions.

Large Perturbations
The beginning of large geomagnetic storms driven by coronal mass ejections (CMEs) is typically marked by a sudden storm commencement (SSC), a few tens of nT intensification in the low-latitude ground-based magnetic field intensity, lasting typically for some tens of minutes. SSCs are produced by interplanetary shocks on the front end of CMEs, which can compress the magnetopause inside geosynchronous orbit. As an interplanetary shock impacts the magnetosphere it launches a large-amplitude fast magnetosonic wave. The leading portion of the bipolar electric field pulse associated with the wave can exceed 200-300 mV/m (Wygant et al. 1994) and is predominantly westward. According to spacecraft observations large SSC events produce injections of tens of MeV electrons and protons all the way into the inner radiation belt (Blake et al. 1992;Wygant et al. 1994;Looper et al. 2005). A number of test-particle simulations of the effects of shock-induced waves on the radiation belts were conducted with the use of empirical field models (e.g., Li et al. 1993;Gannon et al. 2005) as well self-consistent electromagnetic fields from global MHD models (e.g., Hudson et al. 1997;Kress et al. 2007. Simulations showed that trapped particles can E × B drift inward with the wave front through multiple L-shells, undergoing significant energization in a fraction of a drift period due to conservation of the first adiabatic invariant (see Fig. 9a). Particles are energized as long as they stay in phase with the azimuthally propagating wave front. This process, therefore, depends on the azimuthal phase of particle gradient-curvature drift motion and cannot be described with a phase-averaged Fokker-Planck equation (80). Full Liouville's equation (69) must be solved to model rapid particle energization by shock-induced electric field pulses.
Early observations of depletions of the outer belt intensities during storm main phase (Dessler and Karplus 1960;McIlwain 1966) were attributed to an adiabatic response of relativistic electrons to a slow (compared to electron drift period) increase in ring current intensity, which is referred to as the D st effect. An increase in ring current intensity decreases the magnetic flux Φ enclosed by electron drift-bounce orbits. To conserve Φ electrons move outward to regions of lower magnetic field intensity. Since μ = p 2 ⊥ /2mB is also conserved, the outward motion decreases electron energies. Thus, measurements of electrons within a fixed energy at a fixed radial location after increase in ring current register electrons previously located at lower radial distances where their energy was higher and their phase space density lower so that a lower intensity is measured. In recent years with the development of more quantitative empirical models of storm-time magnetic field (e.g., Tsyganenko andSitnov 2005, 2007;, it became apparent that the ring current has much more profound effect on the outer radiation belt. Test particle simulations (Ukhorskiy et al. 2006b) show that storm-time intensification of highly asymmetric partial ring current produces fast outward expansion of electron gradient-curvature drift orbits leading to their loss through the magnetopause. Depending on the storm magnitude, particles from a broad Lrange of outer belt can be permanently lost. These theoretical predictions were recently  (Kress et al. 2007). Panel (b): An example of large-amplitude (>200 mV/m) whistler waves in the inner magnetosphere  confirmed by the observational analysis of multi-spacecraft data (e.g., Millan et al. 2010;Turner et al. 2012). Due to its rapid nature and dependence on the magnetic local time (azimuthal angle) this effect can be described only with full Liouville's equation (69).
Typically pitch-angle and energy diffusion coefficients in radiation belt models are computed based on statistical properties of waves derived from time-averaged spectral intensity data. For whistler chorus waves characteristic time-averaged wave amplitude is ∼0.5 mV/m (Meredith et al. 2001). Recent analysis of instantaneous wave data (Cattell et al. 2008;Kellog et al. 2011; showed that whistler chorus waves can have very large amplitudes >200 mV/m (Fig. 9b). Such large-amplitude whistler waves can accelerate electrons by more than an MeV in less than a second (Cattell et al. 2008), trap electrons (Kellog et al. 2010), and cause their prompt scattering into the loss cone and consequent precipitation into the atmosphere . While it was suggested that some aspects of particle response to large-amplitude coherent waves can be described with a Fokker-Planck equation (Albert 2010), bounce and gyrophase dependent aspects of wave particle interactions require fully kinetic treatment.

Non-diffusive Transport
In the previous section we showed that the diffusion approximation (80) is valid only on time scales τ D much longer than the correlation decay time τ c . Regardless of whether the stochasticity in the system has an extrinsic (noise) or intrinsic (nonlinearity) nature, for moderate wave amplitudes K < 10, typical for the inner magnetosphere, there is the following hierarchy of time scales: Red line: average of ( L * (t)) 2 over all realization of electron motion. Dashed yellow line: radial diffusion ( L * (t)) 2 = 2D QL t (Adapted from  where T is the time between effective collisions in the system. While τ D always exists for an unbounded system, for a bounded system there is an additional requirement: where v T is the characteristic velocity of stochastic transport and L is the system size. This requirement means that phase correlations must decay before particle distribution spreads over the entire system. Theoretical estimates and detailed numerical simulations (Ukhorskiy et al. 2006a; suggest that condition (93) may never be satisfied for radial transport in the outer electron belt. As a result radial transport always exhibits large deviations from the radial diffusion approximation. The consequences of non-diffusive electron transport are illustrated by Fig. 10 showing results of test-particle simulations of radial transport under typical magnetospheric conditions. Panel (A) in Fig. 10 shows snapshots of an ensemble of particles at various stages of radial transport driven by quiet-time oscillations in the solar wind dynamic pressure. Initially all particles had the same value of the third invariant, but were evenly distributed over the drift phase. Last snapshots (panels (A)4a and (A)4b in Fig. 10) correspond to the time moment when the ensemble expanded up to the magnetopause. Particle distribution function (bottom panels) at this point still exhibits large number of pronounced peaks indicative of persistent phase correlations (compare with Fig. 5). Panel (B) in Fig. 10 shows time evolution of ensemble-averaged ( L * ) 2 (black lines) computed for 15 statistically identical time intervals. Had the diffusion approximation been valid, all transport curves would have been very close to each other and monotonically grow in time. However individual transport curves exhibit large deviations from each other and from the straight line corresponding to diffusion with locally estimated quasi-linear diffusion coefficient, indicating that the diffusion approximating has not been attained over the time by which particle ensembles spread over the entire system.  Tsyganenko and Sitnov (2007) magnetic field at P dyn = 3 nPa. Test particles were launched at r = (8, 0, 0) with equatorial pitch angles of 80°(red), 59°(blue), and 20°(cyan) (adapted from Ukhorskiy et al. 2011). Panel (B): Schematic illustration of bifurcating particle dynamics in a simplified case of a symmetric (north-south and east-west) dayside compressed magnetic field. Field-aligned profiles of magnetic field intensity (panels (a)-(e)), and phase portraits of bounce motion (panels (f)-(j)) at different points of electron drift orbit around Earth (adapted from Ukhorskiy et al. 2011). Panel (C): Earthward boundary (L M is the L value at midnight) of the bifurcating orbit region at different values of the solar wind dynamic pressure P dyn computed with guiding-center simulations in the Tsyganenko (1996) magnetic field model for different values of the equatorial pitch angle and midnight

Drift Orbit Bifurcations
Many observational techniques rely on computing electron phase space density as a function of three adiabatic invariants. In particular, radial (L * ) profiles of electron phase space density computed at constant values of the first and second invariants are used as a diagnostic of relative roles of local and global acceleration mechanisms across the outer electron belt (e.g., Green and Kivelson 2004;Chen et al. 2007). If the phase space density has a local peak at some L * value, much exceeding the phase space density value at the outer boundary of the belt, it is considered to be an indication of additional electron acceleration process operating locally at this L * value. Recent studies (e.g., Öztürk and Wolf 2007;Wan et al. 2010;Ukhorskiy et al. 2011) suggest that this argument should be used with great caution.
In a dayside compressed magnetosphere electrons can exhibit three types of trajectories ( Fig. 11(A)). As was discussed in Sect. 2, stably trapped particles (shown in cyan color) participate in three distinct quasiperiodic motions, timescales of the motion are separated by multiple orders of magnitude, and all three adiabatic invariants are conserved. Particles from the magnetopause loss cone intersect the magnetopause and escape the belt before completing a full circle around Earth (red color). Since the drift trajectories are not closed in this case, only the first and the second invariants exist and are conserved (before particles are lost). The third type of trajectories that undergo bifurcations is shown in blue. The existence of bifurcating orbits has been known for a long time (e.g., Northrop and Teller 1960;Northrop 1963;Roederer 1970;Shabansky 1971). More recently though it was realized that at drift orbit bifurcations particle trajectory crosses a separatrix in the phase space plane associated with particle bounce motion. According to the general theory (Cary et al. 1986;Neishtadt 1986) its second invariant is violated at each of the separatrix crossings. This follows from the fact that at the approach to the separatrix particle frequency goes to zero and according to Eq. (17) the change in the adiabatic invariant is no longer exponentially small.
If a guiding-center particle is trapped in a stationary magnetic field, both its energy and the first invariant are conserved: which means that the magnetic field intensity at the bounce points (B(s m ) = B m ) is also constant. In a compressed geomagnetic field the distribution of B(s) along the field lines can exhibit two qualitatively different profiles (Figs. 11(B)a and 11(B)c). On the night side, B(s) has a single minimum at the equator similar to a dipole field (U profile) (Figs. 11(B)a, 11(B)e). On the dayside, adjacent to the magnetopause, however, B(s) has a local maximum at the equator and two minima below and above the equator (W profile) (Figs. 11(B)b-d).
Consider a particle initially bouncing across the equator from points with some value B m of the magnetic field intensity and gradient curvature drifting from the nightside to the dayside into the region where the magnetic field is compressed and has a local maximum at the equator. At some point of the drift trajectory the B(s) profile changes from the U to the W shape (Figs. 11(B)a and 11(B)b). As the particle drifts further into the dayside, the height of the equatorial maximum in the W profile grows. If the magnetic field intensity at the maximum increases up to B m , the particle can no longer cross the equator and its drift orbit exhibits a bifurcation. To conserve the magnetic field intensity at the bounce points, the particle branches off the equator into one of the local B(s) minima pockets (Figs. 11(B)b and 11(B)c). The trajectory traverses the dayside region either below or above the equator never crossing it until the point where the field at the equator decreases back to the B m value at the bounce points. The trajectory then bifurcates again and the particle resumes bouncing across the equator (Figs. 11(B)c and 11(B)d).
At drift orbit bifurcations the particle phase space trajectory crosses a separatrix (Figs. 11(B)g and 11(B)h), which divides the (p , s) phase plane into three distinct regions. The region outside the separatrix corresponds to the bounce motion across the equator, while two lobes connected at a saddle point correspond to trajectories trapped below and above the equator. As the particle approaches the separatrix, its instantaneous bounce period increases logarithmically (as discussed in Sect. 2) and in some small vicinity of the separatrix becomes comparable to the drift period. In this vicinity the quasiperiodic character of the bounce motion is broken, since the effective potential of the motion there is changing at the time scales of the instantaneous bounce period and can no longer be considered slowly varying (ε in Eq. (8) is no longer small). Close to the separatrix the second invariant is therefore not conserved. At two consecutive separatrix crossings corresponding to bifurcations off the equator and back, the invariant exhibits jumps. As a result by the time the particle resumes its motion across the equator it accumulates a nonzero change in the second invariant. Each bifurcation also leads to radial and pitch angle jumps. Consequently when the particle drifts back to its initial location on the nightside, the drift orbit does not close on itself as in the case of stably trapped particles (Fig. 11(A)).
The range of the second invariant (or equatorial pitch-angle) values affected by bifurcations at given radial locations depends on the degree of the day-night asymmetry in the geomagnetic field, which is mostly controlled by the solar wind dynamic pressure (P dyn ).
To quantify the extent of the phase space region affected by bifurcations, we calculated the Earthward boundary of the bifurcating orbits at three different values of equatorial pitch angle as function of P dyn using guiding-center simulations in the Tsyganenko 96 magnetic field model (Tsyganenko 1996) at moderate values of the dynamic pressure (P dyn < 10 nPa). The radial location of the boundary was quantified by L at midnight, L M . The results are shown in Fig. 11(C). As can be seen from the figure, a broad range of the outer belt trajectories is affected by bifurcations. At geosynchronous orbit, for instance, at P dyn > 6 nT all orbits with the equatorial pitch angles α eq > 50°(which constitutes most of the pitch-angle distribution) exhibit bifurcations.
In the bifurcating region particle drift motion around Earth is no longer quasiperiodic (i.e. does not have three independent integrals of motion): there is no slow varying control parameter λ in the Hamiltonian function (see Eqs. (6) and (7)), which can be adjusted to turn the bifurcations off. For the drift motion, bifurcations are a property of the unperturbed Hamiltonian. The third adiabatic invariant therefore is undefined for bifurcating orbits and particle phase space density cannot be transformed into the invariant space. An alternative methodology is required for the analysis of relative roles of various acceleration mechanisms extending into the bifurcation region of the outer belt phase space.

Summary
In summary we provide a reference table (Table 1) of relativistic formula from this chapter, which are most commonly used in modeling, theory, and the analysis of radiation belt particle data.
Magnetic moment M = μ/γ Magnetic moment is an adiabatic invariant of motion only in non-relativistic limit (γ 1) This approximation assumes that the magnetic field is curl-free, a more general case is treated by (29)-(32); p is the momentum component andb is the unit vector parallel to magnetic field, U D is the guiding-center velocity perpendicular to magnetic field (32)

Second invariant
The integration is carried out along a fixed magnetic line between conjugate bounce points Third invariant Φ = C A · dl C is any contour on the drift-bounce surface closed around Earth (40) Dipole L r = LR E sin 2 ϑ Is constant along a dipole field line, it measures the distance from the dipole center to the field line at the equator in Earth radii (R E 6380 km); r is the radial distance and ϑ is the co-latitude to the equatorial points of the drift-bounce shell on which the particle would be, if all nondipolar contributions to the magnetic field would be adiabatically turned off; B 0 31000 nT is the magnetic field intensity on Earth's surface at the equator (41) Liouville's equation ∂f ∂t +q · ∂f ∂q +ṗ · ∂f ∂p = 0 f = f (p, q) is the density of particles in the phase space (PSD), (p, q), conserved along particle trajectories Intensity j α (E, r) = p 2 f (p, r) Is the number of particles with kinetic energy E coming from a given direction α which impinge per unit time, unit solid angle and unit energy, on a surface of unit area oriented perpendicular to their direction of incidence T > τ c is the characteristic time of diffusion, I i is the adiabatic invariant change over T Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.