Fitting the trajectories of particles in the equatorial plane of a magnetic dipole with epicycloids

Motion of charged particles in geomagnetic field has important applications in plasma and space physics, for example, to dynamics of radiation belts and ring current Lyons84 ; Gombosi . Most commonly, this motion is described using the Guiding Center (GC) approximation northrop , which is one of the most successful approximations in plasma physics. Under the GC approximation, the motion of a charged particle is split into two motions considered separately: fast gyration of the particle around the guiding center and a much slower drift of the guiding center. The accuracy of this approximation decreases as the energy of the particle increases and while it is excellent for ring current electrons, it is hardly applicable to cosmic rays.

Motion of electrons and ions in the magnetosphere is, of course, three dimensional. However, many important insights into physics of various populations of the charged particles and their dynamics can be gained from considering particles in the equatorial plane of the terrestrial dipole li98 ; sarris2005 ; kabin2017 . This is often done to simplify the analysis. Although seldom used in practice, there is an analytical solution for Lorentz equations of motion of charged particles in the equatorial plane a magnetic dipole kusaka ; Juarez1949 ; Avrett1962 . Naturally, for small particle energies this solution approaches the GC expectations, but for higher energies can deviate from it very significantly Avrett1962 .

In a typical application of GC, the particle is considered to be “collapsed" into the guiding center and only the motion of the guiding center itself is discussed. The gyration speed of the particle in the approach only appears through the first adiabatic invariant (the magnetic moment), which is assumed to be conserved exactly Kabin2019 . Thus, all information about the phase of the particle as it gyrates around the gyrocenter is lost. For some applications, however, it may be desirable to keep the phase information even in the GC approximation. For example, results of a particle interaction with high-frequency waves may depend on the relative phases of the particle and the wave. Some numerical codes may switch between GC approximation and integration of full Lorentz equations depending on the magnetic field parameters. In such cases, the phase of the particle transitioning from the region of GC approximation to full Lorentz integration is often assigned randomly.

Full motion of a charged particle in the GC approximation is represented by an epicycloid johnston . In this paper we study how well the particle trajectory can be approximated by an epicycloid curve. In particular, we use the results of the analytical solution for the particle motion in the equatorial plane of a dipole field Avrett1962 to determine the parameters of the epicycloid. We also study the accuracy of the epicycloid approximation based on GC parameters. In this paper we consider two versions of the GC approximation: the simplest possible implementation based on the initial position of the particle, and a more accurate one in which parameters at the guiding center are used. We show that with well-chosen parameters, the epicycloid approximation of the particle trajectory remains accurate for surprisingly large energies, even for those approaching the threshold of confinement in the dipole field Kabin2018 . Such high energies are usually considered to be well outside the area of applicability of the GC approximation.

Here we consider magnetic field aligned with the $z$ axis of a right-handed Cartesian coordinate system $(x,y,z)$:

where ${\cal M}$ is the strength of magnetic dipole moment expressed in units of T$\cdot$ m³. For $z=0$ this magnetic field coincides with that in the equatorial plane of a magnetic dipole. Motion of charged particles in this field is described by Lorentz equations:

where $m$ is the mass of the particle, $e$ is its charge, ${\bf v}$ its velocity and ${\bf r}$ its position. For magnetic field (1) motion along the $z$ axis is decoupled from the motion perpendicular to it, so for the rest of the paper we only consider the trajectories in the $(x,y)$ plane. It is convenient to introduce

to simplify the notation.

For magnetic field (1) there is a class of very simple solutions of (2) which are circles with the center at the origin. If $r$ is the radius of the circle, the corresponding circular speed is

This circular speed provides a convenient scale for the velocity in this problem.

Motion of charged particles in the equatorial plane of a magnetic dipole has many practical applications and has been studied extensively in the past kusaka ; Juarez1949 ; Avrett1962 . This problem has an analytical solution which can be expressed in term of elliptic integrals. Typical particle trajectories are confined to an annulus defined by $r_{\min}$ and $r_{\max}$ which can be calculated as:

Here $v$ is the speed of the particle (it remains constant during the motion), $r_{*}$ is the radius where the particle velocity is pointing in the radial direction, and $v_{*}$ is the speed of the circular motion given by (4) at distance $r_{*}$ Avrett1962 . If $v$ exceeds $v_{*}/4$ then $r_{\max}$ becomes undefined; in this case the particle trajectory is unbounded and eventually the particle moves to infinity. Here we only consider bounded trajectories. Transition from bounded to unbounded trajectories for various axisymmetric fields is discussed in detail by Ref. Kabin2018, .

Period of a single cycle (gyroperiod) can be calculated as Avrett1962 :

where integral $I_{2}$ is a function of $v/v_{*}$ only; the subscript 2 is used for consistency with Ref. Avrett1962, . The integration variable $\xi$ in (6) is $r_{*}/r$. Note, that equation (6) gives the “synodic" period, i.e. the time it takes the particle initially, for example, at $r_{\min}$ reach the $r_{\min}$ position again.

The average drift around the origin of the coordinate system can be calculated as

with

Note, that $\Omega_{D}$ is negative for $A>0$, i.e. positively charged particles in upward magnetic field drift clockwise. And integral $I_{1}$ is also a function of $v/v_{*}$ only. Note that both integrals $I_{1}$ and $I_{2}$ have integrable singularities at both ends of the interval, thus certain care needs to be taken for accurate numerical evaluation.

A typical trajectory of a particle in the equatorial magnetic field is shown with red lines in figures 3 and 1.

Cycloidal motion is, of course, expected in the GC approximation. In this section we seek to approximate the trajectory of a particle making the best use of the analytical parameters of its trajectory discussed in section II. Specifically, we place the center of the cycloid at $(r_{\min}+r_{\max})/2$ and ensure that the amplitude of circular motion is $(r_{\min}-r_{\max})/2$. We also use the exact expressions for the drift and gyration frequencies. Then the equation for the epicycloidal approximation to the particle trajectory can be written as

where $t$ is time and $t_{0}$ is the initial time which defines the initial phase of the particle in its trajectory. If $t_{0}=0$ then the particle is initially $x=r_{\min}$ and $y=0$, which is a convenient initial condition. With this choice of initial conditions, the epicycloid touches the exact trajectory at $r_{\min}$ and $r_{\max}$ points. Note, that in equation (9) $\omega$ is the “sidereal" angular frequency of the particle which is computed as

The function ${\rm sign}(A)$ in this expression accounts for the direction of gyration of the particle which depends on the sign of the charge and the direction of the magnetic field.

Figure 1 shows an example of a comparison between the particle trajectory calculated numerically by integrating equations (2) with that approximated by an epicycloid. In this figure, $A=10$ and $v/v_{*}=0.22$. Red line shows the numerically calculated particle trajectory and blue dashed line is the epicycloid approximation.

A traditional geometric interpretation of equations (9) is a path traced by a point attached to a circle of radius $b$ rolling without slipping around the outside of a fixed circle of radius $a$, where the point is at a distance $R$ from the center of the moving circle johnston . Parameters $a$, $b$, and $R$ can be calculated as follows:

These geometric parameters are also shown in figure 1. Note, that although $a$ is close to $r_{*}$, they are different, see the upper panel of figure LABEL:abR. Figure LABEL:abR shows the variation of the epicycloid parameters, $a$, $b$, and $R$, normalized with $r_{*}$, with dimensionless speed of the particle $v/v_{*}$. It follows from equations (5), (6), and (7) that the curves in figure LABEL:abR for different initial positions of the particle coincide; parameters $a/r_{*}$, $b/r_{*}$, and $R/r_{*}$ depend only on $v/v_{*}$. Parameters $R$ and $a+b$ have obvious physical interpretations as the gyroradius of a particle and the distance of the guiding center from the origin, respectively. The physical meaning of the parameters $a$ and $b$ individually is unclear.

Note, that figure 1 corresponds to a very large speed of the particle, which is nearly 90% of the confinement threshold speed for a dipole. For a proton in the Earth magnetic field with $r_{*}=3R_{E}$ (typical value for the inner radiation belt) this speed would correspond to about 800 MeV, which is much larger than typical energy of a radiation belt proton. Nevertheless, even at such extremely high energies, equation (9) provides a fair approximation for the trajectory of a particle. Furthermore since the exact frequencies are used, the epicycloid approximation does not contain secular error. Of course, in cases of practical interest, analytical expressions for such parameters are unlikely to be available. Nevertheless, it is remarkable how accurate cycloidal approximation (9) with well-chosen parameters can be, even for very high particle energies which would normally be considered to be well outside the GC validity.

A more typical application of the cycloid approximation would be based on GC quantities, which we discuss next. Accuracy of the cycloid approximation is further discussed in section V.

The gyrofrequency and the gyroradius for a particle moving in magnetic field (1) can be calculated as

The gradient drift speed is e.g. Ref. Lyons84, .

giving the drift frequency of

Note, that in the above expressions the gyrofrequency and the drift frequency can be either positive or negative, depending on the sign of $A$. These definitions simplify the expression for the epicycloid approximation for the trajectory below.

An important question arising in applying equations (12-14) to describe the motion of a specific particle is what value for $r$ should be used in these equations burby2013 . The simplest interpretation of $r$ in equations (12) and (14) is to consider it to be the instantaneous location of the particle. Under this interpretation, however, the gyroperiod, gyroradius, and the drift frequency all change as the particle moves. If, for example, we assume that the particle is initially at $r_{\min}$ then, because the magnetic field decreases with $r$, as the particle moves, $\omega_{g}$ decreases while $r_{g}$ and $\Omega_{GD}$ both increase.

It is well known that the accuracy of the GC approximation improves if some sort of an average value is used for the magnetic field when applying equations (12) and (14) burby2013 . One way to improve the accuracy of the GC approximation is by calculating the gyroperiod and gyroraidus of the particle based on the previous estimation of the gyrocenter instead of the current position of the particle. For example, assuming that the particle is initially at $r_{\min}$ we can set up the following iteration for the gyrofrequency, gyroradius, and drift frequency:

with $r_{g}^{(0)}=vr_{\min}^{3}/|A|$. Here, superscript ⁽ⁱ⁾ refers to the iteration number (not the power). Our numerical calculations show that iteration (15) for magnetic field (1) converges for $v<0.204v_{*}$, independent of the initial value of $r_{\min}$. Note, that the velocity corresponding to this convergence threshold is smaller than that corresponding to the confinement threshold.

Using expressions (15), epicycloid equation for the trajectory can be written as:

Here we use the same initial conditions as before: for $t=0$ the particle is at $x=r_{\min}$ and $y=0$.

Figure 3 shows exact trajectory of the particle with red line, the simplest GC approximation ($i=0$ in equation (15)) with pink dotted line, and GC approximation with converged parameters with blue dashed line for $v/v_{*}=0.15$. Although figure 3 corresponds to a smaller speed than figure 1, the approximation to the trajectory is not nearly as good as that in figure 1. This is not surprising, considering that approximation (16) does not include any exact results. It is, however, obvious that iterating parameters for the GC approximation dramatically improved the results. It should also be noted that figure 3 still corresponds to a very high speed which is 60% of the confinement speed. For smaller values of $v$, GC is considerably more accurate, as discussed further in section V. Because even the iterated frequencies still differ from the exact ones (6) and (7), GC approximation suffers from secular error, and over many periods diverges from the exact trajectory. This type of error is also easily visible in figure 3.

Next we analyze the accuracy of the epicycloid approximations introduced in the previous sections. First, we define the time-dependent distance $\delta$ between the exact trajectory and approximations described in sections III and IV:

Here $x(t)$ and $y(t)$ are given by equations (9) or (16) and the exact trajectory is obtained either analytically, as described in Ref. kusaka, ; Avrett1962, or by a high accuracy numerical integration of the equations of motion.

Figure 4 shows the approximation error $\delta$ of equation (9) as a function of time for several values of $v/v_{*}$. The time here is normalized with $T$ given by equation (6), which is different for every value of the velocity. Because of symmetry, only half of the gyroperiod is shown. As expected, the approximation error goes to zero at $r_{\min}$ and $r_{\max}$ and peaks roughly midway between these turning points. The black dashed line in Figure 4 shows that the largest error, $\max(\delta)$, occurs at slightly different points during the cycle; for small $v$ the largest error occurs at quarter of the period. Comparing figures 1 and 4 it is clear that the error shown in figure 4 is larger than the distance between the exact and approximated trajectories in figure 1. This occurs because in equation (6) both the exact and approximated position of the particle are used at the same moment in time; thus $\delta$ is not exactly the geometrical distance between the two curves. If one is interested strictly in the distance between the two trajectories, without regard for the time, then the appropriate estimation of accuracy would be just the geometrical distance between the two curves, calculated as

Equation (18) uses the numerically calculated trajectory as the reference curve and at its every point finds the smallest distance to the approximated trajectory. It is obvious how this definition can be modified to use the approximated trajectory as the reference curve producing a similarly defined function $\delta_{1}^{\prime}$.

Figure 5 shows the comparison between $\delta$ (blue line) defined by (6) and the geometrical distance $\delta_{1}$ (red line) defined by (18). For reference, $\delta_{1}^{\prime}$ is also shown as a red dashed line; the maximum value of $\delta$ is, of course, the same as the maximum of $\delta_{1}^{\prime}$, which is emphasized by a horizontal black dashed line. All plots are made for $v=0.22v_{*}$. Figure 5 shows that for this speed $\delta_{1}$ is roughly a factor of 5 smaller than $\delta$ and peaks somewhat closer to $r_{\min}$. The value of $\delta_{1}$ is consistent with the distance between the real and approximated trajectories seen in figure 1. Thus, the accuracy of approximating the trajectory of the particle with epicycloid (9) is several times higher than approximating the position of the particle along the trajectory. The errors for the trajectory approximation (16) are qualitatively the same but larger.

It is natural to expect for approximations (9) and (16) to work best for small values of the particle velocity, when the GC approximation is most accurate, and to break down for large velocities approaching the confinement limit. This is illustrated by figure 6 which shows the $\max(\delta)$ with a solid blue line and $\max(\delta_{1})$ with a dashed blue line as a function of $v/v_{*}$ for approximation (9); note the logarithmic scale of this figure. Also shown are $\max(\delta_{2})$ where $\delta_{2}$ is computed with equation (6) for the simple GC approximation, based on the instantaneous position of the particle (red dashed line), and $\max(\delta_{3})$ for the iterated GC approximation (red dotted line). For small $v$ the error $\max(\delta)\sim v^{2}$ for both approximations (16) and (9), while the geometric distance between the curves decreases at an even faster rate of $\max(\delta_{1})\sim v^{3}$. Although the error for simplest GC approximation based on the particle position also decreases as $\sim v^{2}$, it is always roughly a factor of 10 larger than that for either the iterated GC or the epicycloid approximation. It is impressive, that for most of the considered range of $v$, the iterated GC approximation works nearly as well as the epicycloid approximation with the exact frequencies. The iterated GC approximation only loses it accuracy for high values of the velocities, and in fact this approximation cannot be computed at all for $v>0.204v_{*}$ when iterations (15) diverge. Figure 6 only shows the accuracy of all the approximations for the first gyro-cycle; because of the secular errors inherent in the GC approximation, over time it becomes progressively worse. The cycloid approximation with the exact frequencies is free of the secular errors. Furthermore, somewhat surprisingly for such a simple approximation, it remains reasonably accurate (on the order of 10%) even at speeds close to the confinement threshold.

Another important way in which the accuracy of approximations (9) and (16) can be quantified is by considering the velocity of the particle. Here we focus on equation (9), but the results for the GC approximation are quite similar. During motion in magnetic field, in the absence of electric field, the velocity of the particle remains constant. Differentiating (9) we arrive at the following equation for the velocity for the epicycloid approximation:

where the average value of $v_{approx}^{2}$ over the period is given by

Unfortunately, $v_{approx}$ always depends on time, so approximation (9) is never exact. Expanding the expressions for $r_{\min}$, $r_{\max}$, etc. into Taylor series for small $v$ (see e.g. Ref. Avrett1962, ) it is easy to show that the second term in the above equation for $\left<v_{approx}^{2}\right>$ dominates and this term approaches $v^{2}$ for small $v$. This is not surprising, since the first term in equation (19) represents the drift speed and the second the gyration speed. One of the requirements for GC approximation applicability is that the drift speed is much smaller than the particle speed northrop .

Figure 7 shows $v_{approx}^{2}/v^{2}$ as a function of time ($t/T$) for one period for several values of $v/v_{*}$. As velocity of the particle deceases, the ratio $v_{approx}^{2}/v^{2}$ approaches horizontal line at value 1; however, for larger values of $v/v_{*}$ the amplitude of the change in $v_{approx}$ is quite large. The dependence of this amplitude $(r_{\max}^{2}-r_{\min}^{2})\Omega_{D}\omega/2v^{2}$ on $v/v_{*}$ is shown in figure 8. For small $v$ this normalized amplitude decreases linearly with $v$.

Finally, figure 9 illustrates the deviation of $\left<v_{approx}^{2}\right>$ from $v^{2}$ as a function of $v/v_{*}$. The quantity $\left|\left<v_{approx}^{2}\right>-v^{2}\right|/v^{2}$ shown in figure 9 is the normalized deviation from the expected value. For small values of $v$ this error decreases as $v^{2}$, the same rate as $\max(\delta)$ in figure 6.

Comparison of figures 8 and 9 suggests that the largest cause of the error of approximation (9) is likely associated with the amplitude of the velocity oscillations during the gyration. This conclusion is also supported be the fact that the trajectory itself is approximated by (9) better than the particle position along the trajectory, as indicated by figure 6. It would be interesting to develop a modification to the epicycloid approximation which represents the velocity of the particle more accurately.

Approximating particle trajectories with cycloidal curves can be applied to other axisymmetric fields as well. A simple generalization of (1) is

Magnetic fields with such spatial dependencies were considered, for example, by Ref. Kabin2018, ; $n=3$, obviously, corresponds to the equatorial plane of a magnetic dipole. As long as $n>0$ the shape of the trajectories remains very similar to those discussed in this paper so far. The epicycloid approximations for these trajectories can be calculated in exactly the same way, with the only difference being that there are no known analytical solutions for the values of $n$ other than $3$ and $1$ kabin2015 . Therefore, the epicycloid parameters have to be calculated either from direct numerical integration of the particle motion for one period or from GC approximation. The situation changes somewhat for $n<0$, when magnetic field intensity increases, rather than decreases, with distance. This behavior of the magnetic field changes the sign of the magnetic field gradient in (13), and, therefore, the direction of the drift velocity. Therefore, for $n<0$, the trajectory is represented by a hypocycloid Gordon1984 . In fact, equations (11) still apply in these cases, however, parameter $b$ becomes negative in which case equation (9) describes a hypocycloid. Figure 10 shows representative examples of such trajectories for $n=-3$. The curves in both upper and lower panes are hypocycloids, but the orbits in the lower one enclose the origin, while the ones on the upper do not. The accuracy of the trajectory approximations remains qualitatively similar in all cases.

In this paper we described several approximations for trajectories of particles moving in the equatorial plane of a magnetic dipole. These approximations, inspired by the Guiding Center approach, are based on cycloids. The best results are obtained, not surprisingly, when exact periods and gyration amplitudes are used. Approximation in which the gyroradius, gyrofrequency, and the drift period of a particle are based on its initial position is by far the least accurate. Intermediate accuracy is achieved by using GC approximation in which the values of the gyroradius, gyrofrequency, and the drift period are iterated. For small values of the velocity, this approximation performs virtually as well as that based on the exact parameters, but for larger velocities it loses accuracy faster. Furthermore, since the radius of convergence of the gyrofrequency iterations is found to be smaller than the threshold velocity for particle confinement by the dipole field, this iteration procedure cannot be carried out for all values of initial velocity of the particle. A detailed analysis of accuracy of these approximations shows that the trajectory of the particle is reproduced more accurately than its velocity. Although only motion in the equatorial field of a dipole is analyzed in detail, the results for other similar axisymmetric magnetic fields are virtually identical, except that epicycloids are replaced with hypocycloids in some cases. In the future, it might be worthwhile to try to improve the representation of the velocity in cycloid-based approximations which would allow to design a more accurate yet still simple approximation.