L-shell and energy dependence of magnetic mirror point of charged particles trapped in Earth’s magnetosphere

Abstract

In the Earth’s inner magnetosphere, there exist regions like plasmasphere, ring current, and radiation belts, where the population of charged particles trapped along the magnetic field lines is more. These particles keep performing gyration, bounce and drift motions until they enter the loss cone and get precipitated to the neutral atmosphere. Theoretically, the mirror point latitude of a particle performing bounce motion is decided only by its equatorial pitch angle. This theoretical manifestation is based on the conservation of the first adiabatic invariant, which assumes that the magnetic field varies slowly relative to the gyro-period and gyro-radius. However, the effects of gyro-motion cannot be neglected when gyro-period and gyro-radius are large. In such a scenario, the theoretically estimated mirror point latitudes of electrons are likely to be in agreement with the actual trajectories due to their small gyro-radius. Nevertheless, for protons and other heavier charged particles like oxygen, the gyro-radius is relatively large, and the actual latitude of the mirror point may not be the same as estimated from the theory. In this context, we have carried out test particle simulations and found that the L-shell, energy, and gyro-phase of the particles do affect their mirror points. Our simulations demonstrate that the existing theoretical expression sometimes overestimates or underestimates the magnetic mirror point latitude depending on the value of L-shell, energy and gyro-phase due to underlying guiding centre approximation. For heavier particles like proton and oxygen, the location of the mirror point obtained from the simulation deviates considerably (∼ 10°–16°) from their theoretical values when energy and L-shell of the particle are higher. Furthermore, the simulations show that the particles with lower equatorial pitch angles have their mirror points inside the high or mid-latitude ionosphere.

Introduction

Some of the solar wind charged particles enter the Earth’s magnetosphere and get trapped along the magnetic field lines. Protons, electrons, helium ions, and ionospheric oxygen ions are commonly seen species in the trapped regions of the Earth’s magnetosphere viz. plasmasphere, ring current, radiation belts and so on. Satellite observations suggest that the energies of these trapped particles range from few $\sim \,\hbox {eV}$ to $\sim 100\,\hbox {MeV}$. In general, $\sim \,\hbox {eV}$ range particles are seen in the plasmasphere, $\sim 1-100\,\hbox {keV}$ in the ring current region, and the most energetic particles, $> 100\,\hbox {keV}$ are dominated in the radiation belt regions (Ebihara and Miyoshi 2011; Millan and Baker 2012). Northrop and Teller (1960) established the stability of charged particles trapped in the Earth’s magnetic field. It was shown that in the Earth’s inner magnetosphere, where the magnetic field is assumed to be nearly dipolar, the charged particles undergo three quasi-periodic motions. These are, gyro-motion around its guiding centre, bounce motion along the magnetic field lines between the conjugate mirror points in the northern and southern hemispheres, and longitudinal gradient-curvature drift motion of particle’s guiding centre around the Earth (Williams 1971). A set of three invariants associated with each of these motions define a nearly stable drift shell encircling the Earth.

Under the guiding centre approximation, the theory provides a set of reduced dynamical equations for the motion of the charged particles in a stable or slowly varying magnetic field, which are averaged over gyro-radius in space or gyro-period in time (Baumjohann and Treumann 2012; Li et al. 2011). However, often it is necessary to carry out these calculations on much longer scales due to the multi-scale nature of magnetized plasmas. It may be noted that the particle is likely to experience the varying magnetic field over larger spatial scales as in the case of larger gyro-radius. Hence, the theoretically estimated mirror point latitude obtained under the guiding centre approximation can deviate from the actual particle trajectory. However, so far, there are no attempts to verify or quantify such effects. The information about the latitude of the magnetic mirror point of the particle during its bounce motion is essential because if the mirror point lies in the neutral atmosphere, the particle can lose its energy through collision and cause heating.

The theoretical expression for the mirror point latitude ($\lambda _\text{m}$) of a bouncing particle indicates that the mirror point latitude is decided only by the equatorial pitch angle ($\alpha _{\text{eq}}$) of the particle (Tsurutani and Lakhina 1997). This theoretical expression is based on the guiding centre approximation, where it is assumed that the magnetic field varies slowly relative to the gyro-period and gyro-radius. For electrons, as the gyro-radius is small, they gyrate very close to the magnetic field line such that the ambient magnetic field over a gyration can be assumed to be nearly constant. However, for protons and other heavier particles like helium and oxygen, gyro-radius is large, and it increases with energy and L-shells. In such a scenario, the ambient magnetic field may not be constant over a gyration, and the guiding centre approximation may not be valid. The presence of such heavier atoms in the Earth’s inner magnetosphere is common (Daglis et al. 1999). In this context, the dependency of magnetic mirror point latitude on the L-shell and energy needs to be verified.

In the present study, we have developed a three-dimensional relativistic test particle simulation model to visualize the trajectories of charged particles trapped in the inner magnetosphere. The magnetic field of the Earth is assumed to be dipolar. The relativistic equation of motion is solved using the sixth-order Runge–Kutta method in order to accomplish the numerically stable computations. We have tested the model by verifying the conservation of energy and all three adiabatic invariants from the fourth- and sixth-order Runge–Kutta methods. We found that the sixth-order Runge–Kutta method has significantly less numerical dissipation and can simulate the physically reliable trajectories of oxygen ions, protons, and electrons of energy range 5 keV to 5 MeV at L = 3–6. We have used this simulation to find the latitudes of the magnetic mirror point of charge particles, having a wide range of energy and L-shells. It may be noted that we have not ignored the gyro-motion of the particles. The particles are allowed to gyrate, bounce, and drift, self-consistently. For charged particles, different equatorial pitch angles, L-shells, energies and gyro-phase are considered to estimate their latitudes of magnetic mirror points through the simulation. Also, their deviation from the theoretical relation between $\lambda _\text{m}$ and $\alpha _{\text{eq}}$ is quantified. Besides, we have obtained the ranges of energy, L-shell, and pitch angle of the particles, for which their mirror points exist in the high- or mid-latitude ionosphere.

The paper is structured as follows. The model equations and numerical schemes used in the simulation are discussed in “Simulation model” section. In “Theoretical background” section, we have introduced the theoretical concepts for the estimates of the latitude of the magnetic mirror point of trapped particles and their dependence on the equatorial pitch angle. The results obtained from the test particle simulations are elaborated in “Results” section. Finally, the present study is discussed and concluded in the sections “Discussion” and “Conclusions”, respectively.

Simulation model

The relativistic equation of motion of a charged particle having charge q and mass m in the presence of magnetic field ${\varvec{B}}$ can be written in the form of Lorentz equation,

$$\begin{aligned} \gamma m_{0}\frac{\text{d}{\varvec{v}}}{\text{d}t} = q{\varvec{v}} \times {\varvec{B}}. \end{aligned}$$

(1)

Once velocity ${\varvec{v}}$ is estimated using the above equation, the position of a charged particle can be obtained using the following equation:

$$\begin{aligned} {\varvec{v}} = \frac{\text{d}{\varvec{r}}}{\text{d}t}. \end{aligned}$$

(2)

Here, $\gamma = (1-v^2/c^2)^{-1/2}$ is the relativistic factor, mass of particle is $m = \gamma m_{0}$, where $m_0$ is the rest mass. The velocity vector is ${\varvec{v}} = [v_{x}, v_{y}, v_{z}]$ and its magnitude is estimated from the kinetic energy, $E_\text{k}$ of particle using the following equation:

$$\begin{aligned} v = c\sqrt{1 - \left( \frac{m_0c^2}{m_0c^2 + E_\text{k}}\right) ^2}. \end{aligned}$$

(3)

In the Earth’s inner magnetosphere, the ambient magnetic field lines are closed and can be assumed to be nearly dipolar (Bittencourt 2011). In this region, the terrestrial dipolar magnetic field $\mathbf{B }_{\text{dip}}(\mathbf{r })$ in the Cartesian coordinate system is expressed as follows (Griffiths 2017; Öztürk 2012):

$$\begin{aligned} {\varvec{B}}_{\text{dip}}({\varvec{r}}) = - \frac{B_{0}{R_{e}}^{3}}{r^5}[3xz\hat{{\varvec{x}}} + 3yz\hat{{\varvec{y}}} + (2z^2 - x^2 - y^2)\hat{{\varvec{z}}}]. \end{aligned}$$

(4)

Here, we have used the magnetic coordinate system, such that the xy-plane represents the magnetic equator. The positive x is radially outward direction from the centre of the Earth, and positive y and z, respectively, represent the magnetic east and the magnetic north of the Earth. At the magnetic equator on the surface of the Earth, i.e., $\mathbf{r }=[x=R_e,y=0, z=0]$, the magnetic field strength is considered as $B_{0} = 3.07 \times 10^{-5}$ T. In this simulation model, we have used the fourth- and sixth-order Runge–Kutta method to solve Eqs. (1) and (2). We have noticed that the numerical dissipation is considerably less in the sixth-order Runge–Kutta method, and it gives fairly good numerical stability to simulation code. The expression for the sixth-order Runge–Kutta method is given below (Luther 1968):

$$\begin{aligned} v_{x}^{t +\bigtriangleup t} = v_{x}^{t} + \frac{\bigtriangleup t}{5}\left[ \frac{16k_{1}}{27} + \frac{6656k_{3}}{2565} + \frac{28561k_{4}}{11286} - \frac{9k_{5}}{10} + \frac{2k_{6}}{11}\right] . \end{aligned}$$

(5)

In discretized form, Eqs. (1) and (2) give a total of six equations, which gives the rate of change of $v_{x}$, $v_{y}$, $v_{z}$, x, y, and z with respect to time. At initial time, $t = 0$, the velocity and position components are [$v\sin(\alpha _{\text{eq}})\cos(\psi )$, $v\sin(\alpha _{\text{eq}})\sin(\psi )$, $v\cos(\alpha _{\text{eq}})$] and [$x_0=L\pm r_\text{L}$, 0, 0], respectively. The magnitude of velocity, v for ions and electrons of same kinetic energy will be different due to their mass. Here, $r_\text{L}$ is gyro-radius and $\psi$ is gyro-phase. The gyro-phase, $\psi$ is the angle made by perpendicular velocity, $v_\perp$ with positive x-direction. It decides the particle entry in the horizontal xy-plane and can vary from 0 to $2\pi$. We have taken $\pm r_\text{L}$ in the initial position coordinate so that particle’s (electron/ion) guiding centre in the equatorial plane is at $x=L$. For electron and ion, we took $x_0=L+r_\text{L}$ and $x_0=L-r_\text{L}$, respectively. We have varied L-shell from $3 - 6R_{e}$. In equation (5), $v_{x}^{t +\bigtriangleup t}$ is the Runge–Kutta approximation of $v_{x}$ at time ${t +\bigtriangleup t}$, which is determined by the present value $v_{x}^{t}$ plus the weighted average of six increments $k_{1}$ to $k_{6}$. The expressions of these increments are given in Luther (1968). The size of time interval ($\bigtriangleup t$) is taken as $\sim 1/50$ times the gyro-period. In a similar way of $v_x$, we computed $v_y$, $v_z$, x, y, z. In this model, first, we computed the velocity components [$v_x$, $v_y$, $v_z$] and then the position components [x, y, z] using the estimates of velocity components.

Fig. 1

Simulated three-dimensional trajectories of proton of energy 10 MeV and $L=4$ (i.e., $x=L-r_\text{L}$) with equatorial pitch angle a $\alpha _{\text{eq}}=15^{\circ }$ and b $50^{\circ }$ in the Earth’s dipolar magnetic field for 60 s. The dipole moment is in the negative z-direction. The black arrow shows the westward motion of the proton due to $\varvec{{\nabla }}B \times {\varvec{B}}$ drift. The values of $\lambda _m$ for proton having $\alpha _{\text{eq}}=15^{\circ }$ and $\alpha _{\text{eq}}=50^{\circ }$ are $40.9^{\circ }$ and $18.6^{\circ }$, respectively. The corresponding energy ($E_\text{k}$) as a function of time is shown in c and d. Conservation of energy with time demonstrates the stability of the numerical scheme

To verify the stability of the simulation model, we have checked the conservation of energy during all simulation runs. As an example, the trajectories of protons having (i) energy $E_\text{k}= 10 \,\text {MeV}$, at $L = 4$, $\alpha _{\text{eq}}=15^{\circ }$, $\psi =90^{\circ }$ and (ii) energy $E_\text{k}= 10 \,\text {MeV}$, at $L = 4$, $\alpha _{\text{eq}}=50^{\circ }$, $\psi =90^{\circ }$ are shown in Fig. 1a, b, respectively. Also, their respective energies, $E_\text{k}$, estimated from the simulation at each time step are plotted as a function of time in Fig. 1c, d, respectively. The three-dimensional trajectories of both protons in the dipolar magnetic field are evident in Fig. 1a, b. As expected, its motion is helical around the magnetic field lines. In addition to this, it performs bounce and azimuthal drift motions due to the gradient and curvature of the magnetic field lines. Also, due to the charge dependency of ${\nabla }{B} \times {\varvec{B}}$ drift, the proton moves westward as shown by the black arrow. It is noticed that the latitudinal extent of the proton during its bounce motion is wider for the proton with $\alpha _{\text{eq}}=15^{\circ }$ (see Fig. 1a) as compared to the proton with $\alpha _{\text{eq}}=50^{\circ }$ (see Fig. 1b). Here, $\lambda _{\text{m}1}$ and $\lambda _{\text{m}2}$ can be considered as magnetic latitudes of two mirror points of the proton in southern and northern hemispheres. The mirror point latitudes ($\lambda _{\text{m}1}=\lambda _{\text{m}2}=\lambda _{\text{m}}$) of the two cases obtained from the simulation are mentioned in the respective subplots. Hence we can conclude that the equatorial pitch angle affects the mirror point latitude. In Fig. 1c, d, we have shown the total kinetic energy computed from the simulation, which is constant. This suggests that the chosen numerical scheme is appropriate to analyse the locations of the magnetic mirror points of the trapped particles.

Fig. 2

Magnetic mirror point latitude ($\lambda _{\text{m}}$) as a function of initial gyro-phase ($\psi$) for (a) oxygen (blue), and proton (magenta) ( c) electron (red) of energy 5 MeV, $L= 4$ and 5, and $\alpha _{\text{eq}} = 30^{\circ }$. Maximum radial distance at equator ($r_{\text{max}}$) as a function of initial gyro-phase ($\psi$) for (b) oxygen, and proton of energy 5 MeV initially placed at $x_0=L-r_\text{L}$ (d) electron of energy 5 MeV initially placed at $x_0=L+r_\text{L}$ with $L=5$, and $\alpha _{\text{eq}}=90^{\circ }$

The initial gyro-phase ($\psi$) of the particle can influence the mirror point of the particle (Liu and Qin 2011; Shalchi 2016). To verify this aspect, we run the simulation model for oxygen, proton, and electron of fixed energy (5 MeV), initial position $x_0=L$ ($L = 4$ and 5), and equatorial pitch angle($\alpha _{\text{eq}} = 30^{\circ }$ and $90^{\circ }$) with $\psi$ varying between $0^{\circ }$ to $360^{\circ }$ with an interval of $10^{\circ }$. It may be noted that we have fixed the initial position of the particles in the equatorial plane with $x_0=L$, as we want to understand the effect of gyro-phase on their mirror point latitudes. The variation of magnetic mirror point latitudes for ions (oxygen and proton) and electron is shown as a function of $\psi$ in Fig. 2a, c, respectively. The dependency of initial velocities [$v_{x0}$, $v_{y0}$, $v_{z0}$] on $\psi$ influences the position of particle during gyration. Therefore, the guiding centre of the particle in the equatorial plane changes with $\psi$ even if the initial position of particle, i.e., $x_0=L$ is same. This results in a sinusoidal type variations of magnetic mirror point latitudes with the initial $\psi$ as seen in Fig. 2a, c. The variation of $\lambda _{\text{m}}$ with $\psi$ for the electron is smaller as compare to the ions because the gyro-radius of the electron is small. Also, the variation of $\lambda _{\text{m}}$ with $\psi$ for the electron has an opposite phase to that of ions as they gyrate in the opposite direction.

For all further runs, we took $\psi$ as $90^{\circ }$ to obtain $\lambda _\text{m}$ from the simulation. The reason for taking $\psi = 90^{\circ }$ can be understood from Fig. 2b, d, which show the maximum radial distance of the particle from the centre of Earth ($r_{\text{max}}$) as a function of initial $\psi$ for the ions ($x_0 = L-r_\text{L}$) and electrons ($x_0 = L+r_\text{L}$), respectively. We have chosen $\alpha _{\text{eq}} = 90^{\circ }$ to see the gyration of the particle in xy-plane as in this case the particle does not perform bounce motion. The initial position ($x_0 = L\pm r_\text{L}$) and corresponding maximum radial distance ($x_0 = L\pm r_\text{L}$) are marked by horizontal blue (oxygen), magenta (proton), and red (electron) lines. It can be noted that, all the particles with initial position ($x_0 = L\mp r_\text{L}$) do not have their guiding centre at (L, 0, 0) as $\psi$ varies from $0^\circ$ to $360^\circ$. The results corresponding to $\psi = 90^{\circ }$ are marked by vertical black lines. It is important to note that when $\psi = 90^{\circ }$, ions and electron with initial position ($x_0 = L\mp r_\text{L}$) will have maximum radial distance as ($x_0 = L\pm r_\text{L}$), respectively. Hence from Fig. 2b, d, it can be concluded that a particles (ion/electron) with initial position $x_0 = L\mp r_\text{L}$ has its guiding centre at L only when $\psi = 90^{\circ }$. Therefore, for the further comparison of simulation and theoretical results, we are taking $\psi =90^\circ$ in the simulation to get the particle guiding centre at (L, 0, 0). The theoretical expression to obtain the magnetic mirror point latitude is explained in next section.

Theoretical background

If we inject a charged particle in the Earth’s dipolar magnetic field with an equatorial pitch angle $\alpha _{\text{eq}}$, it will follow the magnetic field lines and bounce between magnetic mirror points. The latitudes of the magnetic mirror points can be determined by the equatorial pitch angle using the following theoretical expression (Roederer and Zhang 2016),

$$\begin{aligned} \sin^2\alpha _{\text{eq}} = \frac{B_{\text{eq}}}{B_\text{m}} = \frac{\cos^6\lambda _\text{m}}{[1 + 3\sin^2\lambda _\text{m}]^{\frac{1}{2}}}. \end{aligned}$$

(6)

Here, $B_{\text{eq}}$ and $B_\text{m}$ are the strength of the magnetic field at the equator and the mirror point, respectively. The above theoretical expression is derived under the guiding centre approximation (Tsurutani and Lakhina 1997). This equation implies that for a given $\alpha _{\text{eq}}$, the magnetic latitude of the mirror point is independent of the L-shell and energy. In other words, Eq. (6) suggests that all charged particles of a given equatorial pitch angle ($\alpha _{\text{eq}}$) bounce back from the same magnetic latitude in the dipolar magnetic field, irrespective of their radial position of guiding magnetic field line and energy (Baumjohann and Treumann 2012). As this equation is derived under the guiding centre approximation, it ignores the effect of gyration and only traces the trajectory of the guiding centre along the magnetic field line. However, in reality, particle gyrate, and it can influence the mirror points of the particles. For the particles gyrating very close to the magnetic field line, this theoretical relation in Eq. (6) can be still applicable as the magnetic field over a gyration varies slowly. However, if the gyro-radius of a particle is large, the solution of $\lambda _\text{m}$ obtained by tracing guiding centre may deviate from the actual mirror point latitude of the particle. In order to verify this hypothesis, we have carried out the test particle simulations to obtain the magnetic mirror point latitudes of electrons, protons, and oxygen ($\text{O}^+$) ions. In Fig. 3, the Earth’s dipolar magnetic field lines are shown in the yz-plane. In this case, the equation of dipolar magnetic field lines can be written as:

$$\begin{aligned} r/R_e = L\cos^2\lambda , \end{aligned}$$

(7)

where $\lambda$ is the magnetic latitude (which can vary from $-90^{\circ }$ to $90^{\circ }$). As the magnetic field line forms a closed loop, it can pass through the Earth as well. Hence, in the case of the Earth’s dipolar magnetic field, there is an upper limit on magnetic latitude $\lambda _l$ ($\cos^2\lambda _l = \frac{1}{L}$) or lower limit on equatorial pitch angle $\alpha _{l}$ ($\sin^2\alpha _{l} = [4L^6 - 3L^5]^{-1/2}$), which are associated with the footprint of the magnetic field line on the surface of the Earth. For a charged particle to remain trapped in the Earth’s magnetic field and not get lost in the Earth’s atmosphere, it is necessary to have $\lambda _\text{m} < \lambda _\text{l}$ or $\alpha _{\text{eq}} > \alpha _\text{l}$. These theoretical limits are called loss cone boundaries. Particles, which have equatorial pitch angle inside the loss cone boundary would be lost into the dense neutral atmosphere of the Earth before they can reach their mirror points.

Fig. 3

Illustration of the dipolar magnetic field lines of the Earth in the Geocentric Solar Ecliptic coordinate system. The Earth’s centre is at the origin [0 0 0], and the magnetic field lines corresponding to $L = 3$ to 6 are shown. The solid red and magenta lines correspond to the theoretical magnetic mirror points with the equatorial pitch angle of $30^{\circ }$ and $10^{\circ }$, respectively. The dotted lines correspond to the maximum deviation observed in the simulation with theoretical estimates of magnetic mirror point latitude. The dotted black circle represents the upper boundary of the ionosphere at 1000 km above the surface of the Earth

Figure 3 shows the two-dimensional magnetic field lines from $L =3$ to $L=6$ obtained from Eq. (7). The Earth’s centre is placed at the origin, and the Earth is shown with the solid blue circle. The dotted black circle represents the upper ionospheric boundary at 1000 km above the surface of the Earth. The solid red and magenta lines represent the latitudes of magnetic mirror points estimated from theoretical Eq. (6) for equatorial pitch angles of $30^{\circ }$ and $10^{\circ }$, respectively. These pitch angle values are chosen such that they are greater than the loss cone angle, $\alpha _\text{l}$. It is evident that as the equatorial pitch angle decreases, the particle comes close to the Earth and for lower L-shells particle’s mirror point lies inside the ionospheric boundary. Another feature is that since the theoretical expression of $\lambda _\text{m}$ is independent of L-shell, we obtain straight solid lines (i.e., same magnetic mirror point) corresponding to the fixed value of $\alpha _{\text{eq}}$ for all L-shells. Besides, theoretical Eq. (6) suggests that the electron, or proton or oxygen ion should have the same magnetic mirror points if $\alpha _{\text{eq}}$ is the same. However, the guiding centre approximation used in theory is not valid for the protons and other heavier particles, whose gyro-radius is large. To demonstrate this, we have performed the simulations to get the trajectories of electron, proton and oxygen for different L-shells and energies, and obtained their deviation from the theoretical results. As an example, we have shown the maximum deviation of $\lambda _\text{m}$ from its theoretical estimate for the proton having two different equatorial pitch angles in Fig. 3. The dotted red ($\alpha _{\text{eq}} = 30^{\circ }$) and magenta ($\alpha _{\text{eq}} = 10^{\circ }$) lines correspond to the maximum deviation ($\pm \delta _{\max}$) in the magnetic latitude of the mirror points obtained from the simulation for these two fixed equatorial pitch angles but having different L-shell and energy. The actual mirror points of proton lie inside the region bounded by the dotted lines corresponding to a fixed equatorial pitch angle. In the case of the electron, this gap between the solid and dotted lines decreases, and for heavier particles like oxygen, this gap will increase. In Fig. 3, it is evident that the magnetic mirror point latitude of the charged particle depends on its L-shell and energy. In the next section, we have examined this feature in detail and quantified the deviation in the magnetic latitude of the mirror point of charged particles from the theory by varying their L-shell and energy.

Results

In order to investigate the effects of L-shell, energy, and gyro-phase of the charged particles on their magnetic mirror point latitudes, we have performed the simulation for electrons, protons and oxygen ions with energy 5 keV to 5 MeV at $L = 3$ to 6 and $\psi = 0^{\circ }$ to $360^{\circ }$ (with an interval of $10^{\circ }$), and tracked their trajectories. As an input to the simulation code, we shoot a particle of kinetic energy, $E_\text{k}$ at the magnetic equator at position [$x_0$, $y_0$, $z_0$] with a fixed equatorial pitch angle ($\alpha _{\text{eq}}$). The equatorial pitch angle is chosen in the range of $5^{\circ }$ to $85^{\circ }$ with an interval of $5^{\circ }$ such that, for a given L-shell it is always higher than the loss cone angle (i.e., $\alpha _{\text{eq}} > \alpha _\text{l}$).

From simulations, one can estimate the position (x, y, z) and velocity ($v_x$, $v_y$, $v_z$) of the particle at each time step. We have converted the position from the Cartesian coordinate system (x, y, z) to the spherical coordinate system (r, $\lambda$, $\phi$) using coordinate transformation. The simulation output of the magnetic latitude ($\lambda$) of the particle is obtained for one bounce period, and its maximum value is considered as the magnetic mirror point latitude. First, we obtained the magnetic mirror point latitude of particles (electron/ion) for (i) fixed energy and varying L-shells and (ii) varying energy and fixed L-shells by taking $x_0=L\pm r_\text{L}$ and $\psi =90^\circ$. Later, we estimated the penetration depth of particle (electron/ion) in the Earth’s upper atmosphere by varying both L-shell and energy by taking $x_0=L\pm r_\text{L}$ and $\psi =90^\circ$. At last, we quantified the effect of gyro-phase, by varying $\psi$ between $0^\circ$–$360^\circ$ and $\alpha _{\text{eq}}$ between $5^\circ$–$85^\circ$ for fixed L-shell ($L= 5$) and energy (5 MeV) with $x_0=L$. These three cases are analysed separately and presented in the following subsections.

Variation of mirror point latitude with L-shell and energy

Here, we have performed the simulations to trace the trajectories of the electron, proton and oxygen ion having four energies, 5 MeV, 500 keV, 50 keV, and 5 keV, and L-shell $L =3$, 4, 5 and 6 with an equatorial pitch angle of $5^\circ$ to $85^\circ$ with an interval of $5^\circ$. The initial position of electron/ion in the equatorial plane is taken as $x_0\pm r_\text{L}$ with $\psi =90^\circ$. Thus, for a given particle with a fixed value of L-shell and energy, we have a total of 17 simulation runs corresponding to the chosen range of $\alpha _{\text{eq}}$. From each simulation output, we have estimated the magnetic mirror point latitude, $\lambda _\text{m}$ and plotted it as a function of $\alpha _{\text{eq}}$ for different L-shell and energy in Fig. 4. Here, panels (a, d, g, j) are for oxygen ions, (b, e, h, k) are for protons, and (c, f, i, l) are for electrons. Each subplot corresponds to the fixed energy of particles with different L-shells, which are depicted by different lines. The theoretical relation of $\lambda _\text{m}$ and $\alpha _{\text{eq}}$ obtained from Eq. (6) is plotted with the black dashed line in each subplot. It is evident from Fig. 4 that the mirror point latitudes estimated from the simulation deviate from the corresponding theoretical curves. The deviation is larger for the oxygen ions, moderate for the protons, and the lowest for the electrons. Another evident feature is that for the fixed energy, the deviation between the theoretically estimated and the simulated magnetic mirror point latitude shows an increase with L-shell.

Fig. 4

The magnetic mirror point latitude as a function of equatorial pitch angle of the particle with different energy. a, d, g, j are for oxygen ions, b, e, h, k are for protons, and c, f, i, l are for electrons. The dashed black line is drawn from the theoretical Eq. (6) and solid lines correspond to the simulation results of $\lambda _\text{m}$ for fixed $\psi = 90^{\circ }$, fixed energy and different L-shells ($L = 3, 4, 5, 6$)

Fig. 5

The deviation of the theoretical and simulation results of magnetic mirror points as a function of the equatorial pitch angle for different L-shells. a, d, g, j are for oxygen ions, b, e, h, k are for proton, and c, f, i, l are for electrons. The magnitude of the absolute deviation is shown by the colour bar

In order to quantify the deviation in theoretical and simulated magnetic mirror points, we have computed the difference ($\delta =\lambda _\text{m}^{\text{theory}}-\lambda _\text{m}^{\text{simulation}}$) between theoretical and simulated latitudes of the magnetic mirror points. Figure 5 shows the difference between the magnetic mirror point latitudes, $\delta$ for $L = 3$–6 as a function of $\alpha _{\text{eq}}$. The colour bar represents the difference between the theoretical and simulated mirror point latitudes. For oxygen ions (see Fig. 5a, d, g, j) the $\delta$ is found to vary between $0^{\circ }$ to $16^{\circ }$, for protons (see Fig. 5b, e, h, k), it varies between $0^{\circ }$ to $10^{\circ }$ and for electrons (see Fig. 5c, f, i, l) it is between 0° to 0.3°. The maximum deviation for each energy is mentioned in the respective subplots. The deviation for all the particles increases with L-shells. Figure 5 indicates that the deviation between the theoretically estimated and simulated magnetic mirror point latitudes increases with the L-shell, and this deviation becomes more evident for the higher energy particles.

Fig. 6

The magnetic mirror point latitude as a function of the equatorial pitch angle of the particle at different L-shells. a, d, g, j are for oxygen ions, b, e, h, k are for protons, and c, f, i, l are for electrons. The dashed black line is obtained from the theoretical Eq. (6) and the solid lines correspond to the simulation results of $\lambda _\text{m}$ for fixed $\psi = 90^{\circ }$, fixed L-shell and different energy ($E_\text{k} = 5\,\text {MeV}$, 500 keV, 50 keV, 5keV)

Fig. 7

The deviation of theoretical and simulation results of magnetic mirror points as a function of the equatorial pitch angle for different energies. a, d, g, j are for oxygen ions, b, e, h, k are for protons, and c, f, i, l are for electrons. The magnitude of the absolute deviation is shown by the colour bar

Furthermore, we have investigated the dependence of magnetic mirror point latitude on the particle’s energy. For this purpose, we have obtained the simulation results for electrons, protons and, oxygen ions of different energies, L-shell and fixed $\psi =90^{\circ}$. We simulated the magnetic mirror points of the particles for $L = 3$ to 6, corresponding to four energy values. Figure 6 depicts the magnetic mirror points as a function of the equatorial pitch angle for different L-shells. In this figure, each of the subplots is obtained for the fixed L-shell and with four energy values (shown by different colours). Panels (a, d, g, j) are for oxygen ions, (b, e, h, k) are for protons, and (c, f, i, l) are for electrons. The corresponding deviation of magnetic mirror point latitude from the theoretically estimated $\lambda _\text{m}$ for each case is shown in Fig. 7. We noticed that for any fixed L-shell, as the energy of particle increases the simulated magnetic latitude of the mirror point deviates significantly from its theoretical estimate. This deviation is more significant for the oxygen ions and protons, as compared to the electrons. Also, this deviation is more significant for the higher L-shell. Overall, this analysis shows that the magnetic mirror point latitude of a particle is indeed dependent on L-shell, energy and mass, apart from its equatorial pitch angle.

Minimum radial distance of particle

We have investigated the penetration depth of the trapped particles in the Earth’s atmosphere. Our simulation model does not have any atmosphere, ionosphere, and waves in the system to affect the particle’s motion. Hence, the particle’s penetration depth is decided by the location of the latitude of the magnetic mirror point above the surface of the Earth. As mentioned earlier, we used different L-shells, energy, and equatorial pitch angle for electrons, protons, and oxygen ions to estimate the latitudes of the magnetic mirror points. The penetration depth of the charged particle is measured in terms of its minimum radial distance ($r_{\text{min}}$) during its bounce motion, which is the position, r of the particle at mirror point latitude, $\lambda _\text{m}$.

Fig. 8

The minimum vertical height, $H_{\min}$ of oxygen ions (a, d, g, j), protons (b, e, h, k) and electrons (c, f, i, l) above the surface of the Earth as a function of equatorial pitch angle. Here, $H_{\min}=r_{\min}-R_e$. Each subplot is obtained for the fixed L-shell and four different values of energy. The dashed-dotted black and solid lines represent the theoretical and simulation results, respectively. The vertical dashed line corresponds to the loss cone angle, and the horizontal dashed line represents the upper ionospheric boundary at 1000 km. $H_{\min}=0$ indicates the surface of the Earth

We have defined the minimum vertical height above the surface of the Earth corresponding to the magnetic mirror point such that $H_{\min}=r_{\min}-R_e$. The minimum vertical height of the magnetic mirror point above the surface of the Earth is shown in Fig. 8 as a function of the equatorial pitch angle for different energies and L-shells. In Fig. 8, panels (a, d, g, j) are for oxygen ions, panels (b, e, h, k) are for protons and panels (c, f, i, l) are for electrons. Each subplot is for fixed L-shell, and different curves in it correspond to different energies. The dashed-dotted black line in each panel represents the theoretically estimated position of $H_{\min}$ for a given value of $\alpha _{\text{eq}}$. The vertical dashed line indicates a loss cone angle, $\alpha _\text{l}$ corresponding to the L-shell. The particles having an equatorial pitch angle less than the loss cone angle get lost in the atmosphere, and theoretically, their mirror points should lie inside the Earth, i.e., $H_{\min}<0$. Here, $H_{\min}=0$ represents the Earth’s surface. The horizontal dashed line in each panel starting from top, respectively, represents the upper ionospheric boundary at 1000 km above the Earth’s surface. The particles having $r_{\min}< R_e$, i.e., $H_{\min}<0$ does not exist practically because these particles have mirror points inside the Earth’s surface. We are getting such particles in the simulations as the expression of the dipolar magnetic field defines the closed field lines.

Fig. 9

The minimum vertical height of a protons and b electrons above the surface of the Earth as a function of the energy at $L = 3$ with the equatorial pitch angle of $10^{\circ }$ and $15^{\circ }$. The horizontal dashed lines correspond to the ionospheric boundaries at 1000 km and 100 km, respectively, above the surface of the Earth

It is noticed that as the equatorial pitch angle decreases the mirror point of particle moves to lower altitudes in the Earth’s upper atmosphere, and hence, the particle can come closer to the Earth. The charged particles with lower equatorial pitch angle get lost in the loss cone at lower L-shells ($L = 3$ and 4). Also, as the energy of the particle decreases, their $H_{\min}$ decreases, i.e., particles with lower energy have their mirror points close to the Earth. Figure 8 indicates that the particles with equatorial pitch angle between $5^{\circ }< \alpha _{\text{eq}} < 10^{\circ }$ have their mirror points inside the ionosphere, i.e., within 1000 km above the Earth’s surface. In Fig. 9, the minimum vertical height above the surface of the Earth, $H_{\min}$ for proton (panel-a) and electron (panel-b) are shown as a function of their energy for equatorial pitch angle of $10^{\circ }$ and $15^{\circ }$ at $L= 3$. The horizontal dashed lines starting from top, respectively, indicate the ionospheric boundary at 1000 km and 100 km above the Earth’s surface. This figure demonstrates that the protons of energy 5 keV to 300 keV and electrons of energy 5 keV to 5 MeV have their mirror points inside the ionosphere if the equatorial pitch angle is less than or equal to $10^{\circ }$. Such particles entered in the ionosphere and can get lost after their interaction with the ionospheric plasma or neutral species present at that altitudes.

Generally, the particle precipitation occurs at around 100 km. While reaching this altitude, the particle motion gets affected through their collisions with ionized particles and neutrals present below 1000 km. The effect of collisions dominates over the magnetic field effects below 100 km and the particle dynamics in this region mainly governed by the atmospheric chemistry. As we are using a static magnetic field with no atmosphere or ionosphere, there is no possibility to explore the particle loss due to wave–particle interaction, particle–particle interaction, and precipitation in the present model. Also, theoretically, magnetic field lines form a closed loop that can take charged particles up to or inside the Earth’s surface. However, in a real situation before reaching the Earth’s surface, the particle will interact with the atmospheric or ionospheric particles and get lost. Therefore, we have used a limit of $H_{\min}$ between 100 to 1000 km.

Gyro-phase effect on mirror point latitude

In the last two subsections, we have taken an initial $\psi =90^\circ$ to understand the effect of L-shell and energy of particle on its mirror point latitude. However, from Fig. 2 it is clear that magnetic mirror point latitude varies with initial gyro-phase ($\psi$). Here, we have quantified deviation in simulated $\lambda _{\text{m}}$ from its theoretical value for different $\psi$. For this purpose, we perform simulation runs for oxygen, proton, and electron of energy $E = 5\,\text {MeV}$ at $L=5$ by varying $\alpha _{\text{eq}}$ between $5^\circ$ to $85^\circ$ (with an interval of $5^\circ$) and $\psi$ between $0^\circ$ to $360^\circ$ (with an interval of $10^\circ$). The initial position of the particle is the same in all simulation runs, i.e., $x_0=L$. Here, we chose the higher energy, i.e., $E=5 \,\text {MeV}$ and L-shell, $L=5$ because in the previous two subsections we have seen that the deviation $\delta$ is large for higher L-shell and energy. We have a total of 629 simulation runs, and from each run, we estimated the deviation $\delta$ in magnetic mirror point latitude for different combinations of $\alpha _{\text{eq}}$ and $\psi$. In Fig. 10, the deviation in simulated magnetic mirror point latitude from its theoretical value is shown for (a) oxygen, (b) proton and (c) electron. We noticed that $\delta$ is maximum for the oxygen ($\pm 15^\circ$), moderate for the proton ($\pm 10^\circ$) and minimum for the electron ($\pm 1^\circ$). For ions, $\delta$ is negative for $\psi = 0^\circ$ to $180^\circ$, and it is positive for $\psi = 180^\circ$ to $360^\circ$, whereas for electron, variation in $\delta$ is opposite to that of ions as they gyrate in opposite direction.

Fig. 10

The deviation in simulated magnetic mirror point latitude from their respective theoretical value as a function of initial gyro-phase for different equatorial pitch angle. a is for oxygen ion, b is for proton, and c is for electron of energy $E_k = 5 \text{MeV}$ and $L= 5$ with initial position $[x_0=L, y_0=0, z_0=0]$. The magnitude of the deviation is shown by the colour bar

The dependency of mirror point latitude on initial gyro-phase can be understood from the direction of initial force acting on the charged particle during its gyration. For example, ions (gyrating in the clockwise direction) with initial $\psi = 90^\circ$ move toward the region of a weaker magnetic field (away from the Earth) and ions with initial $\psi = 270^\circ$ move toward the stronger magnetic field (towards the Earth). It means the initial gyro-phase affects the guiding centre of the particle in the equatorial plane such that L-shell associated with their guiding centre varies with $\psi$. Hence, $\lambda _\text{m}$ varies with $\psi$. As gyration of the electron is opposite to ions, the effect of $\psi$ on $\lambda _{\text{m}}$ is opposite in case of the electron.

Discussion

The theoretical expression of magnetic mirror point latitude suggests that different charged particles with the same $\alpha _{\text{eq}}$ will bounce back from the same magnetic latitude irrespective of their mass, energy, and L-shell. However, our test particle simulations suggest a dependency of $\lambda _{\text{m}}$ on particle’s energy, L-shell, and mass. This dependence can be understood from gyro-radius ($r_{\text{L}}=\gamma m_0v_{\perp }/qB$) of the trapped particle. When the particle is at higher L-shells, it experiences the weaker ambient magnetic field resulting in larger gyro-radius. Also, as compared to the electrons, the gyro-radius of proton and oxygen is more significant due to their higher masses. Likewise, if the energy of the charged particle increases, its velocity increases, which results in more abundant gyro-radius. Therefore, higher energy, mass and L-shell of the particle results in larger gyro-radius of the particle. It may be noted that the theoretical formula given in Eq. (6) is based on the guiding centre approximation and it considers the magnetic field at one single point at the equator ($B_{\text{eq}}$) and mirror point ($B_\text{m}$) along the guiding magnetic field line. The guiding centre approximation ignores the gyration and only traces the trajectory of the guiding centre from one mirror point to another. In reality, particle gyrates around the magnetic field line, and therefore one cannot ignore the gyration associated effects. In the simulation, the trapped particles are allowed to gyrate, bounce, and drift self-consistently. For larger gyro-radius, the ambient magnetic field over a gyration is likely to vary, and the approximation used in the theoretical formulation gets violated. Hence, as demonstrated by present simulations, the actual mirror latitude of particle deviates considerably from its theoretical estimate when gyro-radius is larger.

Fig. 11

Variation of magnetic moment $\mu$ and three adiabatic invariants $J_1$, $J_2$ and $J_3$ estimated from the simulation for electron (a, d, g, j), proton (b, e, h, k), and oxygen (c, f, i, l) having energy 5 MeV, $L=6$, $\alpha _{\text{eq}}=30^{\circ }$, $\psi =90^{\circ }$. The mean value of adiabatic invariants $\langle J\rangle$ and their standard deviation $\sigma$ are mentioned in the respective panels to demonstrate their conservation

In the present study, the trapped particle performs periodic motion so that their motion along the magnetic field lines and drifts across the magnetic field lines are symmetric. In such periodic motion, the three adiabatic invariants associated with particle gyro, bounce, and drift motion are conserved. In the case of theoretical Eq. (6) as gyration is ignored, and only guiding centre trajectory is considered, the first adiabatic invariant is always conserved. $J_1$, $J_2$, and $J_3$ are the adiabatic invariants associated with gyro, bounce, and drift motion of particles, respectively. The first adiabatic invariant is associated with gyro-motion, and it is given by $J_1=2\pi m_0\langle \mu \rangle /q$, where $\mu$ is the magnetic moment, the longitudinal invariant $J_2$ is associated with the bounce motion and the flux invariant $J_3$ is associated with the azimuthal drift motion (Mukherjee and Rajaram 1981). We estimated three adiabatic invariants from the simulation and checked their conservation. These adiabatic invariants are shown as a function of time for electron, proton and oxygen in Fig. 11. In upper panels of Fig. 11a–c the variation of $\mu$ is shown for one bounce period. As particle gyrates and move from equatorial region (low magnetic field) to polar region (high magnetic field) the $\mu$ varies but average magnetic moment $\langle \mu \rangle$ over gyration is constant as shown by the black dash-dot line in upper panels of Fig 11. For electron the mean value of $\langle J_1\rangle$, $\langle J_2 \rangle$, and $\langle J_3 \rangle$ are $3.9\times 10^{-14} \pm 4.3\times 10^{-21}$, $4.2\times 10^{-12} \pm 3.2\times 10^{-18}$ and $7.1\times 10^{-11} \pm 4.6\times 10^{-14}$, for proton $3.5\times 10^{-14} \pm 4.3\times 10^{-18}$, $4.1\times 10^{-12} \pm 1.3\times 10^{-17}$ and $6.7\times 10^{-11} \pm 4.8\times 10^{-14}$, and for oxygen $1.4\times 10^{-13} \pm 3.8\times 10^{-17}$, $6.5\times 10^{-12} \pm 3.0\times 10^{-14}$ and $6.3\times 10^{-11} \pm 6.3\times 10^{-14}$, respectively. It may be noted that the standard deviation in $\langle J_1\rangle$, $\langle J_2 \rangle$, and $\langle J_3 \rangle$ are small, which implies that all three adiabatic invariants remain constant with time. As conservation of adiabatic invariants is verified, we can compare mirror point latitudes derived from the simulation with theory. One may argue that if gyro-radius is large enough to affect the underlying approximation used in (6) then the first adiabatic invariant $J_1$ may also get violated. However, all three adiabatic invariants are found to conserved in the simulation. It is because each of these invariant is obtained from the closed line integral of the canonical moment of the particle over one cycle of their respective motion. For example, the first adiabatic invariant $J_1=2\pi m_0\langle \mu \rangle /q$ is obtained by averaging it over one gyro-period, and it is conserved over a gyration.

Conclusions

In the present study, we have performed the three-dimensional test particle simulations to investigate the dynamics of charged particles trapped in the Earth’s inner magnetosphere. The magnetic field of the Earth is assumed to be dipolar and stationary with time. Since the particle is treated as a test particle, its motion does not affect the ambient magnetic field. We have simulated the trajectories of electrons, protons and oxygen ions and converted the results from Cartesian to the spherical coordinate system. We have obtained the latitudes of the magnetic mirror point through the simulation for particles of different energies, equatorial pitch angles at different L-shells and different initial gyro-phase. The chosen range of L-shell is 3–6, energy is 5 keV–5 MeV, the equatorial pitch angle between $5^{\circ }$–$85^{\circ }$ and gyro-phase between $0^{\circ }$–$360^{\circ }$. The theoretical expression of the latitude of the magnetic mirror point suggests its dependence only on the equatorial pitch angle of the charged particle (Tsurutani and Lakhina 1997; Baumjohann and Treumann 2012). However, our simulation suggests that the magnetic mirror point of the charged particle also depends on its L-shell and energy. It is because the theoretically derived magnetic mirror point latitude is based on the guiding centre approximation, which ignores the effect of gyration of particle. As compare to oxygen ions or protons ($\delta _{\max}\approx 10^{\circ }{-}16^{\circ }$), the deviation in the magnetic mirror point latitude of electrons from theory is less significant ($\delta _{\max}\approx 0.3^{\circ }$). The deviation between simulated and theoretical magnetic mirror point latitude increases with both energy and L-shell of the particle and hence, can not be ignored. For electrons, the gyro-radius is small, so the guiding centre approximation is applicable. For oxygen ions and protons, the ambient magnetic field over a gyration will not be constant due to their large gyro-radius. In such a scenario, the underlying approximation used in the theoretical formulation of the estimation of $\lambda _\text{m}$ gets violated. Thus, one has to be cautious while using theoretically derived latitudes of magnetic mirror points. The initial gyro-phase ($\psi$) affects the position of the guiding centre of the particle so that L-shell associated with their guiding centre varies with $\psi$. Hence, the magnetic mirror point latitude of the particle is found to vary with initial gyro-phase of the particle.

We have estimated the minimum vertical height ($H_{\min}$) of the trapped particles above the surface of the Earth in order to know how deeper they can enter into the ionosphere. We found that the magnetic mirror points of the trapped particles with smaller equatorial pitch angle ($\alpha _{\text{eq}} \le 10 ^{\circ }$) and smaller energy lies inside the ionosphere ($<1000 \,\text {km}$). Therefore, the particles with a smaller $\alpha _{\text{eq}}$ will get lost in the ionosphere, even in the absence of any wave–particle interaction phenomena. It is because they are entering the region where sufficient charged and neutral particles are present to absorb their energy through collisions. Such information will be useful in the atmospheric/ionospheric modelling to understand the particle precipitation in the high latitudes. The theoretical expression for the magnetic mirror point of the charged particle is often used to understand their penetration depth in the ionosphere. In this context, the simulation-based quantified information provided in the present study will be useful to the scientific community to understand the dynamics of the trapped particles in the Earth’s inner magnetosphere.