New method to determine proton trajectories in the equatorial plane of a dipole magnetic field

A parametric description of proton trajectories in the equatorial plane of Earth’s dipole magnetic field has been derived. The exact expression of the angular coordinate contains an integral to be performed numerically. The radial coordinate results from the initial conditions by basic mathematical operations and by using trigonometric functions. With the approximate angular coordinate formula, applicable for a wide variety of cases of protons trapped in Earth’s radiation belts, no numerical integration is needed. The results of exact and approximate expressions were compared for a specific case and small differences were found.


Background
The differential equations of the motion of charged particles in a magnetic dipole field, as approximation of the Earth's magnetic field, were derived already by Stőrmer (1907). Their solutions were obtained only numerically in an attempt to explain the behaviour of aurora borealis. The interest for charged particles trapped inside the dipole magnetic field increased by the discovery of Van Allen belts. A deep discussion of this topic can be found in Elliot (1963). Alfven (1950) developed an approximation splitting particle motion into guiding centre motion and gyration. The method has been enriched by Northrop (1961). Detailed studies of the motion of charged particles in the meridian plane of magnetic dipole fields are due to Markellos et al. (1978a) and (Markellos et al. 1978b).
The specific case of charged particle motion in the equatorial plane of the dipole magnet allowed De Vogelaere (1950) to obtain general periodic solutions of a Hill type differential equation from two numerically found particular solutions. Dragt (1964) in a comprehensive paper, describes a power series solution originating from a double power series solution of Stőrmer (1955). Graef and Kusaka (1938) obtained solutions for the motion in that plane in terms of elliptic integrals of first kind. These include a step of numerical calculation, no matter if the elliptic integral is obtained by a subroutine of the trajectory determining programme or if its value is picked up from tables.
Here parametric expressions of the proton coordinates in the equatorial plane were derived. Both exact radial and approximated angular coordinate expressions are strictly speaking, analytical formulas including only basic arithmetic operations and trigonometric functions. A prescribed accuracy for the angular coordinate is obtained from an exact formula by numerical integration, performed with the desired degree of precision.

Radial coordinate-exact expression
We use a cylindrical coordinate system, Figure 1, with the polar coordinates r, θ in the equatorial plane of the dipole magnetic field (that of paper). Therefore, the z axis is directed along the dipole, the origin being located in its middle. A proton, moving with the velocity v 0 in the equatorial plane will be pushed towards the origin if the magnetic field lines are directed to enter inside the paper, from outside the Earth's body. In other words, we are looking at the figure from the Earth's North magnetic pole. To determine the coordinate r we use two motion constants: the velocity v 0 and the angular momentum, Eq. 3.35 from Miyamoto (2000), p.26.
"Point" means derivative with respect to the time. c t is a constant. q and m are the proton charge and mass respectively. A θ is the only non vanishing vector potential component. For z = 0 this is: with μ the Earth's magnetic dipole moment. Next we use Richardson's (1947) parameter ψ, the angle between the velocity and the theta velocity component, by the substitution v θ = v 0 cosψ, earlier applied to describe ion paths in "wedge" magnetic fields and ion optical studies, see also Ioanoviciu (1989).
After substitution of A θ and r θ : in (2) we have: next we put: and we obtain: The constant c t results from the initial conditions: when the parameter takes the value ψ = ψ i .
Eliminating c l , between (6) and (8), next substituting r by ρ i = r/ri we have: where: Two solutions for ρ i result: Only the solution with "plus" sign offers physically acceptable values.
From the above expressions we obtain the maxima and minima of the proton radial coordinate by derivation with respect to ψ. As ρ i is a function ρ i (cosψ) the derivative is dρ i /dψ = dρ i /d(cosψ)(−sinψ) It vanishes for ψ = π and ψ = 2π values (if we limit to the first loop) that substituted in the solutions for ρ i give the searched extreme values: for cosψ = −1 , and Figure 1 Proton trajectory in the dipole magnet equatorial plane. Defined are the radial r, and angular θ, coordinates (polar coordinates), the parameter ψ, the initial quantities v 0 , r C and ψ i = π/2.
for cosψ = 1 Further radial coordinate expression simplification can be obtained if we select the observation starting point for ψ i = π/2: Here the index "i" has been removed. Now ρ = r/r C , η = μq/(p r r C 2 ) with r C the radial distance for ψ i = π/2, p r the relativistic particle momentum The extreme values of ρ are then: We can use the simplified expressions by looking for an equivalent equation connecting r i and η i to r C and η. By equating the maximum and the minimum of the equation (12), (13) to those of eq. (15), (16) we obtain the following equivalence: The eq. (11) with "minus" sign, is identical with (14) if:. and The stability condition for the proton trajectories results from the equation of ρ. The condition to keep the charged particle trapped inside the dipole magnetic field: ρ must be real. This happens always if: or by detailing η: Angular exact coordinate expression The angular coordinate, θ is obtained by an integration to be performed numerically with the desired accuracy degree.
For the protons of 1 MeV located between 2.5 and 8R E (second belt) the values of η=215.147 and 21.01 while for those of 0.065 MeV η=844.087 and 82.43.
The enumerated η values suggest when the integral giving θ (as r values result straightforward) can be obtained with accuracy. We take the equation (9) and divide it by η. Accounting for 1/η as small quantity, we substitute 1/η = α and ρ = 1 + ε , both α and ε being assumed to be small quantities. Then we obtain: That after the substitutions gives: Parametric calculated exact coordinates ρ=r/r C and θ, as well as the error, (θ-θ app )/θ by using the approximate angular coordinate formula θ app , were tabulated.
By successive approximation, we obtain the following expression for, ε after the third approximation step: ρ obtained as 1 + ε has been substituted in the expression of dθ : As dε has the form: to obtain θ we have to integrate ψ from π/2 to the current ψ value.
After integration it results: Application to a specific case The case of a trapped 60 MeV proton oscillating around the characteristic radial distance of 1.5 Earth radii was considered. The motion characterizing parameter results then to be η = 75.97. The exact coordinates ρ and θ, as well as the relative difference between the exact and the approximated angular coordinates (θ-θ app )/θ in percents were presented, for a complete loop, in Table 1. The error is of the order of 0.002% to 0.12%. In Figure 2 the reduced radial distance r/r C exact values were represented as function of θ app in abscissa, for three consecutive complete loops, corresponding to the intervals, ψ = π/2 to 5π/2, ψ = 5π/2 to 9π/2 and ψ = 9π/2 to 13π/2 respectively.

Conclusion
The coordinates of the protons moving in the Earth's equatorial plane were derived as functions of a parameter. The use of the exact expressions assumes only basic operations and trigonometric functions to be involved for the radial coordinate calculation, while for the angular coordinate a numerical integration is necessary. Combining the exact radial formula with the approximated angular coordinate expression, an entirely analytic set of formulas has been obtained. The accuracy of this description has been shown on an illustrative example. The amazing simplicity of the coordinate expressions suggests possible developments by accounting for field perturbations. Figure 2 Three loops of a 60 MeV proton inside the Earth's dipole magnetic field equatorial plane. The ratio ρ = r/r C is represented in ordinate as function of the angle θ as abscissa. The three loops correspond to ψ parameter variation intervals: π/2 to 5π.2 the first, 5π/2 to 9π/2 the second and 9π/2 to 13π/2 the third respectively.