Escaping dynamics of relativistic protons in the Earth's magnetosphere

Escaping dynamics of relativistic protons in the Earth’s magnetosphere ´ Alvaro Meseguer,1, ∗Juan C. Vallejo,2Jes´ us M. Seoane,2Francisco Marqu´ es,1and Miguel A.F. Sanju´ an2 1Departament de F´ ısica, Universidad Polit` ecnica de Catalunya, Barcelona 08034, Spain 2Nonlinear Dynamics, Chaos and Complex Systems Group, Departamento de F´ ısica, Universidad Rey Juan Carlos, Tulip´ an s/n, 28933 M´ ostoles, Madrid, Spain (Dated: April 3, 2025) 1 Abstract In this work, we study the motion of protons under the influence of the Earth’s magnetic field. We investigate the dynamics and topology of invariant spatial regions in two magnetic field models, namely, the dipolar and Luhmann tail configurations. In both cases, we analyze the motion of the protons in the relativistic regime, with kinetic energy intervals in the range 10 < E < 100 MeV. Specifically, we examine situations in which the protons undergo chaotic scattering due to their interaction with the Earth’s magnetic field. Additionally, we consider scenarios where the protons, during their interaction, fall onto the Earth’s surface. We compare the dipole and Luhmann tail cases by computing the distributions of escape and residence times, the exit basins, and the fractal dimensions of their boundaries. We present a robust symplectic numerical scheme which is suitable for analysing the impact of the nonlinear effects and the presence of strong sensitivity to initial conditions when the model goes beyond the pure dipolar field. Surprisingly, we uncover a scaling law between the coefficient of the decay law and the energy of the protons. Furthermore, our computation of the fractal dimension of the exit basin boundaries Dversus the distance between the protons and the Earth z0reveals a decrease in Das z0increases. We expect this work to be useful for a better understanding of the behavior of protons under the influence of the Earth’s magnetic field when they are in the chaotic scattering regime. I. INTRODUCTION The electrodynamics of Earth’s magnetosphere is very complex, and it has attracted the interest of many physicists and mathematicians for over a century. For a thorough review of this topic, including a comprehensive historical overview and the most relevant and up-to-date studies, we refer the reader to the standard monograph by Russell, Luhmann, and Strangeway [1], and references therein. Figure 1 shows the main characteristics of Earth’s magnetosphere. Close to the Earth, the magnetic field is very similar to a magnetic dipole, generated by the currents in the Earth’s metallic core. However, the solar wind dramatically alters the dipolar structure of the magnetic field at a few radii from the Earth. On the Sun-facing side, a bow shock is formed when the solar wind particles are deflected by Earth’s magnetic field. On the opposite side of the Earth, away from the Sun, the so-called magnetotail is formed, generated mainly by a plasma sheet near the equatorial plane. ∗alv[email protected] 2 FIG. 1. A simple schematic view of the Earth magnetosphere and its main components. The magnetic field lines in red, the Sun on the left. As the solar wind interacts with the Earth’s magnetosphere, high-energy particles are trapped and held around the Earth, forming the so-called Van Allen radiation belts, where the population of charged particles is much larger than average [2, 3]. Two belts are usually described. The inner belt mainly consists of energetic protons, with kinetic energy values ranging from 10 MeV to 100 MeV, generated through the so-called Cosmic Ray Albedo Neutron Decay (CRAND) [4, 5]. Note that, from now on, we denote by Ethe kinetic energy of the protons. This inner belt is primarily the product of cosmic ray collisions in the upper atmosphere. Some of the collision by-products can enter the atmosphere, but when the fragment is a free neutron, it can decay into a proton (which captures most of the energy), an electron, and a neutrino. Although the neutron is very fast, it can statistically decay while still close to Earth. Hence, the resulting energetic proton will remain trapped by the magnetosphere. This results in a belt with variable boundaries, typically extending from 0.2RLto 3.0RL, where RLis the Earth radius. Although the inner belt can also contain high concentrations of electrons with energies within the range of hundreds of keV, see for instance [6], it is the outer belt which consists mainly of highenergy electrons. Our work will focus on the analysis of the dynamics of the charged particles forming the inner belt as they are injected into the Earth’s magnetosphere. This analysis will use the tools of nonlinear dynamics from the perspective of an open chaotic scattering problem, aiming to get insight into the residence times for the injected particles, the topology of the invariant trapping regions where these particles remain captured, as well as the lifetimes of these regions. Scattering problems are defined as the interaction between an incident particle and a potential 3 region or massive object that scatters it. In this context, we say that the scattering is chaotic when the final state of the particle has sensitive dependence on initial conditions (see Ref. [7]). The equations of motion of the test particles that model these interactions are nonlinear, and we can find strong dependence on initial conditions and chaotic dynamics. One defines the scattering region as the domain where the particle is affected by the potential field while, outside this zone, the influence of the potential over the particle can be neglected. The scattering systems are labeled as open when the test particles may escape once they entered in the scattering region after bouncing back and forth for a while. In these open cases, chaotic scattering is regarded as a transient chaos dynamics, and, although the test particle can remain bounded, the resulting dynamics is not a simple one. Regarding the dynamics in a magnetospheric model, complex structures can be found even in the most simple case of a magnetic dipolar field. This is the so called St¨ ormer problem (see Ref. [8–10]). Here, chaos is typically present at high energies and has been observed to be independent of the initial velocity [11–13]. When considering the real Earth’s magnetosphere, the magnetic dipole component still predominates, but the dynamics of the particles and the topology of the belts are strongly dependent on the presence of the magnetotail, as shown in Fig. 2. One of the main goals of this study is therefore to analyse the effects of the presence or absence of the magnetotail regarding residence times (lifetimes) of the particles. Besides, we study how the presence or absence of the magnetotail influences the topological changes of the trapping regions by studying the fractality of the exit basins. Accordingly, the explorations reported here will consider the pure dipolar model and will compare it when adding a magnetotail. In both cases, we aim at determining how both models differ in the residence times for the injected particles and which values of the parameters drive the configuration and lifetimes of the regions where these particles remain trapped. We will also study the coexistence of escape and trapping regions of the two models, and their dependence as a function the energy and initial conditions (day-side or nightside launches). Because of the broad set of parameters modeling the magnetospheric structures, this study will focus on the dynamical study of protons. However, the mathematical formulation and numerical methodologies developed here are completely general and can be applied to other families of charged particles. When modeling particles in the magnetosphere, the guiding center approximation has been often used in the past because it reduces the required computational time, see [14] and references therein. However, here, we need to solve the full dynamics of the particles, modeling its three 4 main dynamical components, gyration, bounce and azimuthal drift, because we aim to analyse the impact of the nonlinear effects and the presence of strong sensitivity on initial conditions when the model goes beyond the pure dipolar field and both high and low energies are possible. Therefore, we will present a robust numerical scheme with reduced numerical dissipation, which is suitable for this analysis. The paper is organized as follows. In Sec. II, we introduce the dipolar and dipolar-tail magnetic models, illustrating the main features of the resulting proton dynamics. In Sec. III we provide some analytical results regarding the dipolar model, address the symmetries of the resulting dynamical system, and introduce the symplectic numerical methodologies used for the energy-preserving time integration of the proton’s motion. Section IV is devoted to characterize the topology of the space regions according to their associated residence times (i.e., the overall time spent by the protons in those regions before an eventual escape). The study of the decay law of the particles in the scattering region is carried out in Sec. V. Section VI presents both the exit basins and the basins of attractions and their fractal nature. Finally, the main conclusions and a discussion of the results are presented in Sec. VII. II. MODEL DESCRIPTION The relativistic dynamics of a charged particle of mass mand charge qin the presence of electric Eand magnetic Bfields is described by the Lorentz force md(γv) dt =q(E+v×B), γ = (1 −v2/c2)−1/2.(1) This equation can be obtained from the relativistic Lagrangian and Hamiltonian L=−mc2p1−v2/c2−q(V−v·A),(2a) H=p(p−qA)2c2+m2c4+qV, p=mγv+qA.(2b) where Vand Aare the electric and magnetic potentials, respectively. In the absence of electric field (E= 0) from (1) we see that the magnitude of the velocity is constant, and the γfactor moves outside the time derivative: mγ dv dt =qv×B(3) This equation is almost identical to the Newtonian case, except for the presence of the γfactor, which is a constant prescribed by the initial conditions of the particle. As has been observed by 5 other authors [15, 16], in the absence of electric fields we can use the modified Lagrangian and Hamiltonian given by L=1 2mγv2+qv·A,H=1 2mγ (p−qA)2,p=mγv+qA,(4) where we assume γis constant in the derivation of the equations of motion, and Ais the magnetic vector potential defined by B=∇×A. These modified Lagrangian and Hamiltonian, which resemble the Newtonian case except for the presence of the γfactor, reproduce the relativistic equations of motion (3), and we will use them because they are simpler than (2), and facilitate the discussion of symmetries and conserved quantities. In order to model the magnetic field of the magnetosphere, we will use a superposition of two different magnetic fields: the Earth dipole magnetic field, Bd, which is dominant close to the Earth, and a magnetic field modeling the magnetotail, Bt, using the Luhmann model [17]. Accordingly, the mathematical expressions of the magnetic vector potential and the magnetic fields read as follows: A=Ad+At,Ad=µˆ z×r r3, µ =µ0M 4π,At=−Btδln cosh(z/δ)ˆ y,(5) B=Bd+Bt,Bd=µ3zr r5−ˆ z r3,Bt=Bttanh(z/δ)ˆ x,(6) where Mis the Earth magnetic dipole moment pointing towards the zdirection; ˆ x,ˆ y,ˆ zare the unit vectors in the x,yand zdirections, respectively. A reasonable value for M, that we will use in the paper, is M=−8×1022 A m2, the minus sign accounting for the fact that the Geographic North Pole is approximately the magnetic south pole. This gives a value µ=µ0M/(4π) = −8×1015 T m3. Henceforth, ˆ xand ˆ yare the unit vectors in the Cartesian xand ydirections, respectively, where the xaxis points in the direction from the Sun to the Earth (so that the positive xdirection points away from the Sun). In Eq. (6), Btmeasures the strength of the magnetotail magnetic field, and the value we are using in this study is Bt=−10−8T, consistent with former models [17], whereas δis a convenient length measuring the thickness of the density current associated with Bt. The density current associated with a prescribed magnetic field is given by the Amp` ere-Maxwell equation ∇×B= µoj. For the magnetic dipole the current is located at the origin (in practice in the Earth’s core) because ∇×Bd= 0. For the magnetotail, we get a density current at the equatorial plane in the ˆ ydirection: jt=Bt µ0δcosh−2(z/δ)ˆ y≈4Bt µ0δe−2|z|/δ ˆ y.(7) 6 Because 95% of the current is concentrated in a layer of width 3δaround the equatorial plane, in this study we will use δ= 3RL. Moderate variations of the value of δdid not result in remarkable changes in the dynamics observed. (a) (b) -20 0 20 40 60 80 x -30 -15 0 15 30 z FIG. 2. Magnetic field lines of a magnetic dipole Bdincluding a magnetotail Bt. (a) Meridional plane. (b) 3D perspective view. Red dots are the points where B= 0, and the red lines are the field lines through these points. Bt=−10−8Tand δ= 3RL. Lengths are given as multiples of Earth radius. Notice that the day side is the region on the left of the Earth (corresponding to negative values of x) while the night side is the region on the right of the Earth (corresponding to positive values of x). Figure 2 shows magnetic field lines in the meridional plane y= 0 for Bt=−10−8Tand δ= 3RL. The magnetic field is usually represented by plotting the magnetic field lines, that are the solution of dr/dτ =B, where τparametrizes the magnetic lines. The value B= 0 is found at the two points marked in red in the figure, with coordinates given by r=z0(−√2,0,±1), z0= 9.442 77RL. The unstable manifold of the upper B= 0 point and the stable manifold of the lower B= 0 point define the magnetosphere and magnetotail of the Earth. They are plotted in red in Fig. 2. The magnetopause on the noonside is limited by the bow shock. The limit of the magnetopause on the noonside in the equatorial plane is very close to x0=−13.3704RL. It is convenient to render the problem dimensionless using the Earth radius RL= 6.371 × 106m, the velocity of light cand the Earth magnetic field at the equator Beq =µ/R3 L= 0.309 362 T, as units of length, velocity and magnetic field. The resulting dimensionless gov7 erning equations are dv dt =bv×B,L=1 2v2+bv·A, b =qµ mγ ,(8) A=ˆ z ×rr3−btδln cosh(z/δ)ˆ y,B=3zr r5−ˆ z r3+bttanh(z/δ)ˆ x,(9) where band btare the dimensionless versions of qµ/(mγ)and Bt/µ respectively. Parameter b depends on γand the ratio q/m, and therefore on the initial conditions and nature of the particle considered, while bt= 3.232 457 ×10−4is fixed, measuring the dimensionless strength of the magnetotail field, in good agreement with the aforementioned Luhmann model [17]. III. THEORETICAL ANALYSIS AND NUMERICAL INTEGRATION OF THE MODEL In this section, we describe the Lagrangian structure of the equations of motion, the conserved quantities derived from the Lagrangian formulation, and the numerical techniques used. To streamline the main text, the details of the analysis have been confined to the Appendix A, and we present here a summary of the relevant conclusions. Since the model does not explicitly depend on time, the Hamiltonian, that coincides with the kinetic energy of the particles, is conserved. In the dipolar model without a magnetotail, i.e., when bt= 0, there is an additional conserved quantity, related to the angular momentum of the particle ℓ=ρ2˙ θ+bρ2 r3,(10) where uis the modulus of the particle velocity (conserved), and ρ=px2+y2is the distance to the dipole axis; see equation (A2) in Appendix A. The conservation of ℓand uresults in a bounded quantity, g(ρ, z), depending only on the spatial coordinates: g(ρ, z) = a ρ−ρ r3,|g(ρ, z)| ≤ u |b|, a =ℓ/b. (11) In the case a > 0there exist bounded regions, for |g| ≤ 1/4and ρ≤2/a, where the particles remain trapped, and have a lunula shape similar to the Van Allen belts, as shown in Fig. 3. Particles trapped in these regions can penetrate along the edges of the lunulas into the atmosphere and even reach the Earth’s surface, resulting in the auroras (borealis and australis). The shape and size of the trapping regions depends on the value of a=ℓ/b, i.e. depends on the initial conditions of the particle and the nature of the particle (protons, electrons, etc.). Therefore the Van Allen belts are not fixed in space, but depend on the kind of particles coming from the 8 0 1 2 3 a; -1 -0.5 0 0.5 1 az -20 -10 0 1 FIG. 3. Contours of the function gfor a=ℓ/b > 0. The red lines correspond to g=±1/4, the green line to g= 0, and the black curves to g=±0.15. The gray region is the lunula bounded by the curves g=±0.15. solar wind, that end up populating the belts. In the tail model, there is just one conserved quantity, namely the kinetic energy of the particle. The analysis done about trapping regions in the dipole case, can not be extended in the presence of the tail. However, considering the weakness of the tail magnetic field near the Earth compared with the dipole strength, it is reasonable to assume that some of the characteristics of the trapping regions will still apply. In any case, the trajectories are very complex, often chaotic, and can only be computed numerically. Trajectories confined on the equatorial plane are exact solutions of the equations of motion, even including the magnetotail field. In fact the magnetotail field is zero on the equatorial plane, therefore both models coincide in this plane. Standard time integrators for non-stiff ODEs are generally non-symplectic and therefore do not preserve energy during the time evolution of the system. However, even symplectic integrators may fail to preserve additional conserved quantities, such as the ℓmagnitude in the dipole case. In this study, we have used and adapted the Runge-Kutta 9(8) time-stepper [23], modified with a projection method [24–26] to ensure the preservation of conserved quantities. We have chosen a high-order method with a variable time step to allow reasonable step sizes that reduce computational cost while maintaining high accuracy in the results; see details in Sec. 4 of Appendix A. Figure 4 shows the trajectory of a proton in the magnetotail field, obtained through numerical integration using the aforementioned symplectic-projection method. The trajectory appears chaotic and is confined within a lunula-shaped region. 9 V. ESCAPE TIMES AND DECAY LAWS Our system can be viewed as a scattering problem. Chaotic scattering refers to the interaction between an incident particle and a potential region or massive object, which scatters the particle [28]. Therefore, the analysis of escape times also provides insights into the dynamical properties of the system. Decay laws describe the evolution of populations of injected particles over time, as well as their likelihood of remaining bound as time progresses. In this work, we analyze the dynamical decay law of protons, which defines the time dependence of the particles’ survival probability within the scattering region. Particle motion can be either bounded or unbounded, following regular or chaotic trajectories. In conservative systems, regular motion is typically associated with Kolmogorov–Arnold–Moser (KAM) tori. In hyperbolic chaotic scattering, periodic orbits are unstable, and KAM tori do not exist in the phase space. In this case, the dynamical decay law of the particles follows an exponential form. However, in the nonhyperbolic regime, KAM tori coexist with chaotic saddles, leading to algebraic decays in the particles’ survival probability [29]. We note that we will examine how protons escape the influence of the magnetic field for energies in the range of [10,100] MeV, which corresponds to the relativistic regime. In the relativistic case [30], it has been shown that, generally, when the particle velocity is sufficiently large, the decay law of the particles becomes exponential. A. Decay laws Here, we aim to investigate whether the injected protons follow the characteristic exponential decay law typically associated with relativistic systems. To do so, we consider particles with energies in the range of [10,100] MeV, fixing the initial injection height at z0= 2.0RL. Figure 7 shows the ratio Rof bound particles for both the dipole and magnetotail models, with the latter including day-side and night-side injections. The time horizons and energies range within the intervals t∈[0,30] s and E∈[10,100] MeV, respectively. Unexpectedly, Fig. 7 shows that the magnetotail seems to have a trapping effect on the particles, albeit moderate, but clearly noticeable. Compared with the dipolar model, Fig. 7(a), Fig. 7(b) and Fig. 7(c) show that the survival ratio slightly increases when the magnetotail is present, although the particles are released from the side opposite to the most active region of the magnetotail. 16 FIG. 7. Evolution of the ratio of bound protons, R, with time, for protons injected at z0= 2RL, with the color code indicating energy levels. (a) Dipole field. (b) Magnetotail Luhman’s D-tail model (day-side injections) (c) Magnetotail Luhman’s N-tail model (night-side injections). The ratio Rof bound particles is related to the likelihood of a proton to remain bound over time. This, in turn, can be linked to the evolution of the energy spectrum, or the distribution function of the number of protons as a function of energy, which characterizes the particle population at a given time. Therefore, Fig. 8 shows the evolution of these energy spectra for both the night-side and day-side populations. From a qualitative perspective, we observe that both distributions are quite similar, although the night-side population seems to settle slightly earlier At t= 0, the population starts with a flat distribution, corresponding to all particles being injected from the night side. The figure illustrates that the most energetic protons, with E= 100 MeV, are the first to escape. As time progresses, lower-energy protons begin to escape as well. Around T≈15 s, the distribution stabilizes, with no further escape of particles occurring at a given energy. 17 FIG. 8. Evolution of the energy spectrum of the populations of protons injected at z0= 2.0RL, both from the day side, D-tail model (left) and night side, N-tail model (right). Larger energies escape earlier, but after roughly T≈15 sthe spectrum remains stable. The panels of Fig. 7 show a transient behavior for Raround T≈15 s, observed in both the dipole and magnetotail cases. However, the transient is more clearly seen in the magnetotail cases. To analyze this behavior, we now plot the evolution of the fraction Rof unbound particles for different injection heights z0, fixing the energy at E= 25 MeV and E= 100 MeV. Figure 9 shows the increase in the fraction of escaping particles, R, with energy and initial injection height, z0. It also illustrates the dependence of the asymptotic bound population on different injection heights, specifically z0= 0.5,1,2,5RL. At larger z0values, almost all injected particles escape, whereas a different behavior is observed at smaller z0. Typically, an initial fraction of particles escapes immediately after injection, while another group remains bound for some time before eventually escaping. In the asymptotic regime, only the particles that remain bound throughout the entire integration time persist. The rightmost panels of Fig. 9 provide a zoomed-in view of the evolution of this transient population. Here, two distinct slopes in the decay law curves are visible, with the slopes being clearly different for the lower energies (continuous line). At the largest z0, this asymptotically bound population disappears, as nearly all injected particles eventually escape, and no long-term bound population remains, with the decay continuing indefinitely. 18 FIG. 9. Dynamical decay law, or log −plot evolution of the fraction of particles that remain bound as the trajectory progresses. Top panels correspond to day-side injections, bottom panels to night-side injections. Rightmost panels detail the evolution at early times of the evolution up to t= 15time-units. B. Exponential behaviors at relativistic regimes In a general classical scenario, the fraction of remaining particles in the scattering region decreases linearly with energy, as shown in Ref. [31]. However, in the relativistic domain, the decay law of the particles follows an exponential behavior given by R∼e−αt, where 1/α represents the characteristic decay time due to scattering processes. We consider protons with energies in the relativistic regime, where the velocity of the particles exceeds β > 0.1. At this point, nearly all trajectories depart from the non-relativistic regime [30]. To study the potential exponential behavior of Rat different energies, we have again chosen energy values in the range E∈[10,100] MeV and fixed z0= 2RL. Figure 10 shows the trend 19 FIG. 10. Plots of the coefficient of the exponential law αversus ln Ecorresponding to protons launched from z0= 2RL. Dipole (top), magnetotail day-side injections (bottom left) and magnetotail night-side injections (bottom right). We can observe similar linear trends in all cases. The protons are escaping faster in the dipole case compared with the magnetotail day-side injection when the energy is very large, since there is not effect of a magnetotail. However, the protons escape faster in the case of the night-side magnetotail. across the three different scenarios considered. We observe that the expected exponential behavior holds in all three models. The rationale behind this linear trend is intuitively understandable. As demonstrated in our previous work [31], the fraction of trapped particles decreases with energy. Since this fraction is described by the expression e−αt, it is reasonable to expect that αscales linearly with ln E. The top panel plots the exponent of the decay law, α, against ln Efor the dipole model on the day side. The trend closely follows a straight-line fit (r2= 0.995). The middle panel of Fig. 10 shows the Luhmann tail model for day-side injections. Here, the plot exhibits a similar linear trend with r2= 0.994. However, it is noticeable that protons escape slightly faster in the 20 -2 -1 0 1 2 Log z0 -0.2 -0.15 -0.1 -0.05 0 α FIG. 11. Plots of the coefficient of the exponential law αversus log z0for the night side in which the protons have E= 25 MeV and we vary the injection height z0. A rough exponential decay can be observed where the protons escape faster insofar the distance zis increasing. dipole model. To quantify this difference, we define the fraction f=∆α ∆ ln Eas an indicator of the escape velocity. Specifically, for the dipole case, f≈0.022, while for the Luhmann tail model, f≈0.027, suggesting that the magnetotail retains particles for a shorter duration before their escape. The bottom panel depicts for the Luhmann tail model, night-side injections. The linear trend between αand ln Eis also evident, with an r2= 0.995. Here, f≈0.02, indicating that protons take longer to escape the influence of the magnetic field compared to previous cases, owing to the effect of the tail. Moreover, they take more time than on the day-side injections case. One can conclude that the influence of the magnetotail is stronger for particles injected from the night-side. This difference might arise because protons are closer to the tail on the night side. We also investigated the effect of varying the injection height on the escape time, while keeping the energy of the protons fixed. For this, we set E= 25 MeV and vary z0∈[0.25,5]RLon the night side. Figure 11 shows the exponent of the decay law, α, plotted against ln z0. An exponentiallike behavior is observed, with protons escaping more rapidly for smaller values of z0. VI. EXIT BASINS AND FRACTAL DIMENSION In the previous sections, we have analyzed the temporal behavior of the protons as they were injected. This section will numerically explore the complexity of the dynamics of the exit basins. 21 The goal is to investigate how these basins evolve with the control parameters of the model, comparing the behaviors of day-side and night-side particles, and relating the trapping structures to the temporal behaviors, residence times, and escape time distributions discussed earlier. In dissipative systems, a basin of attraction is defined as the set of points that, when taken as initial conditions, are attracted to a specific attractor. When two distinct attractors exist within a certain region of phase space, two basins are formed, separated by a basin boundary. This basin boundary may be a smooth curve or, in some cases, a fractal curve. In our study, we have observed that the injected particles can either remain bounded, fall to Earth, or escape after some time, exiting the system through a specific exit region. Therefore, in analogy to basins of attraction, we define exit basins (or escape basins) as the set of initial conditions that lead to a particular exit [7]. When two or more escape regions are possible, fractal boundaries can emerge, reflecting the system’s chaotic dynamics. In such cases, the boundary separating one basin from another is not clearly defined, leading to unpredictable behavior. Figure 12 compares the structure of the exit basins in the dipole and Luhmann tail models. The escaping particles are protons with energies of 25 MeV and 100 MeV, and their initial positions (x0, y0, z0)lie on constant z=z0planes, as indicated. The region (x0, y0)∈[−10,0]RL× [−10,10]RLis divided in 100 ×200 small cells, with particles launched from the center of each cell toward Earth. The trajectories are computed up to a time of Te= 30 s, but the time integration stops when r < 1(Earth fall, shown in red in the figure) or r > 80 (escape from Earth). When a particle escapes, we distinguish between z(Te)>0(North escape, shown in blue) and z(Te)<0 (South escape, shown in green). The remaining initial conditions, corresponding to particles that remain trapped near the Earth, are plotted in white. The time spent before escaping or falling is indicated by the grading between white and the corresponding color. The North exit, South exit, and Earth-falling basins exhibit fractal characteristics. The area of these fractal regions increases with the distance from the equatorial plane of the injection point (i.e., with increasing z0) and with with the energy of the particle. The presence of the magnetotail significantly enhances the fractality of the basins. Additionally, the time spent by the particle before falling to Earth, indicated by the intensity of the red dots within the white regions, also displays fractal behavior. In order to quantify this evolution of the fractality of the basins, we have computed the boxcounting dimension Dof the initial conditions corresponding to the confined protons, that is, the white regions of Fig. 12. Figure 13 shows the variation of the box-counting dimension Dwith 22 Dipole 25 MeV D-Tail 25 MeV N-Tail 25 MeV Dipole 100MeV D-Tail 100MeV N-Tail 100MeV z0= 0.5RL z0=RL z0= 2RL z0= 3RL FIG. 12. Exit basins for the dipole and Luhmann models. Each panel represents the initial conditions plane, x0−y0initial conditions plane in which |x0| ∈ (0,10)RLand y0∈(−10,10)RL. Green corresponds to unbounded trajectories escaping with z > 0, while blue corresponds to those escaping with z < 0. Red means particles falling to Earth. White means confined particles. From top to bottom, protons are injected at z0= 0.5,1.0,2.0,3.0RL. The first three leftmost columns correspond to 25 MeV, while the three rightmost columns to 100 MeV. 23 the height of the injection point z0. We have depicted different physical situations of interest: the night side of the tail model (denoted in blue and dark blue), the day side of the tail model (denoted in red and purple) and the dipole model (light green and green), where all of them for 25 MeV and 100 MeV, respectively. We can observe that the smaller the z0, the larger the value of the dimension of the white area, meaning it almost covers a surface. 1 2 3 4 5 z 0 0.5 1.0 1.5 2.0 D N-Tail 25 MeV N-Tail 100 MeV D-Tail 25 MeV D-Tail 100 MeV Dipole 25 MeV Dipole 100 MeV FIG. 13. Variation of the fractal dimension Dof the areas corresponding to confined particles, when launched at different z0. Here, we show different physical situations: the night side of the tail model (denoted in blue and dark blue), the day side of the tail model (denoted in red and purple), and the dipole model (denoted in light green and green), where 25 MeV and 100 MeV, respectively. A decreasing trend of Dwith z0is clearly observed both for protons with E= 25 MeV and for those with larger energy value, E= 100 MeV. The larger the energy, the smaller fraction of particles that remains bounded, hence, the smaller the box dimension. There is a steep decrease of the dimension when z0>2.0. Moreover, beyond z0≈3.5, the number of initial conditions leading to bounded trajectories is very small, and one can observe some fluctuations on the computed values of the box dimension Dat large z0. Regarding the dipole data, at the largest z0, there are no confined motions, hence no data to plot. VII. CONCLUSIONS AND DISCUSSION Our work has been focused on analyzing the dynamics of protons injected into Earth’s magnetosphere at varying heights, z0, using tools from nonlinear dynamics and chaos theory. These protons are the result of neutron decay near Earth. However, given the high velocities of the 24 neutrons, we can consider injections occurring at greater distances from Earth as well. We have first analyzed the proton dynamics analytically for the case of a pure dipolar magnetic field. However, to explore a more realistic scenario that includes a magnetic tail, a numerical approach is necessary. Our investigation revealed that the chosen integration scheme must be highly accurate to account for particles falling into the Earth’s poles, due to the strong sensitivity to initial conditions in this scenario. When analyzing residence time plots, even considering our initially limited set of horizontal launches for defining the initial conditions of the protons, we have observed that the models already produce structures resembling the inner belt. Here, the role of the tail does not seem to be very critical in creating these structures, which are formed by both bound and unbound protons. We have also analyzed the decay laws associated to those populations. In the relativistic regime, decay laws are exponential, as expected. Interestingly, the shapes of the radiation belts do not change significantly when comparing the dipole and magnetotail models. However, the magnetotail appears to have a trapping effect on the residence times of the injected particles. This effect, although moderate, is clearly noticeable. The influence of the magnetotail is stronger for particles injected from the night side, possibly because protons are closer to the tail on the night side. An exponential-like behavior can be observed, with protons escaping from Earth more rapidly for smaller values of z0. The decay laws also show that the escaping populations exhibit some transient behavior. The higher the energy of the proton, the shorter the duration of the transient. Only after some time post-injection does the energy spectrum of the populations stabilize for a given z0. Because trapped protons are more likely to interact with the Earth’s atmosphere, these low-energy populations could be the most prone to such interactions [32]. Finally, we have also analyzed the exit basins. They show the extent to which the resulting protons remain bound. The North exit, South exit, and Earth-falling basins exhibit a fractal nature, with the area of the fractal regions increasing with both the distance from the equatorial plane of the injection point (increasing z0) and the energy of the particle. The presence of the magnetotail significantly increases the fractality of the basins. A preliminary computation of the fractal dimension of the basin boundaries, D, as a function of the distance z0, shows a decrease in Dwith increasing z0. Our Luhmann model is still very simple, neglecting the tilt of the dipole, day/night side asymmetries of the magnetic lobes, and ring currents. It is also a static model, ignoring any highfrequency fluctuations that might affect the dynamics of the studied particles. However, despite 25 sonably large time steps of about 10−2, except during polar penetrations, which require a reduction of the time step by two orders of magnitude. To preserve conserved quantities, this time-stepper has been suitably modified by adding a projection of the predicted trajectory over the invariant manifold, as described in [24–26]. In all the cases explored in this work, such projection method led to phase space integration orbits that numerically preserve energy (and angular momentum, in the dipolar model) almost to machine precision. Figure 4 shows the trajectory of a proton in the magnetorail field, and first row of figure 16 shows the trajectory of a proton with the same initial conditions in the dipole field. Both trajectories look chaotic, and show penetrations along the wedges of the lunula region where the protons are confined. The differences are clearly seen in the second row of Fig. 16. This figure shows temporal series of x,θand ℓin the two cases (dipole and magnetotail). The mean angular velocity ⟨˙ θ⟩is slightly larger in the presence of the magnetotail, as shown in the θ(t)plot. Time steppers preserving the conserved quantities (energy, and ℓin the dipole case) have been used. The time series of ℓ(t)shows that in the dipole case ℓis preserved (error less than 10−14), while in the magnetotail case it shows temporal oscillations, yet very tiny. 10!10 10!5 Tol 10!2 10!1 Error h_ 3i 10!15 10!10 10!5 Tol 10!14 10!12 10!10 10!8 "` 10!15 10!10 10!5 Tol 10!15 10!10 "u FIG. 17. Figure showing error estimates as a function of resolution. We observe that a substantial reduction in error requires tolerances smaller than approximately 10−12. Error estimates and convergence tests for numerical integrations are usually conducted by comparing results obtained with different time steps. In our case, since we use variable time steps, we assess numerical accuracy based on the tolerance required in the Runge-Kutta projection method. Figure 17 presents these accuracy tests for the trajectory illustrated in Fig. 16. The second 32 and third panels show the maximum deviation of the conserved quantities along the trajectory for different tolerances, εℓ=||ℓ−ℓo||∞and εu=||u−uo||∞. We observe that for tolerances less than or equal to 10−10, the conserved quantities remain constant up to machine precision. The first panel displays the relative error of the mean azimuthal angular velocity for different tolerances. ⟨˙ θ⟩is a highly sensitive variable, as it critically depends on the accurate resolution of polar penetrations, especially for chaotic trajectories such as the one considered in Fig. 16. The error is computed by comparing with the smallest tolerance used. We observe that a substantial reduction (three orders of magnitude) in the error requires tolerances smaller than approximately 10−12. For this reason, we have used a very small tolerance of 3×10−14 in all computations. [1] C. T. Russell, J. G. Luhmann, and R. J. Strangeway. Space physics: An Introduction. (Cambridge University Press, Cambridge, 2016). [2] M. K. ¨ Ozt¨ urk. Trajectories of charged particles trapped in Earth’s magnetic field. Am. J. Phys. 80, 420 (2012). [3] M. G. Kivelson and C. T. Russell. Introduction to Space Physics (Cambridge University Press, Cambridge, 1995). [4] J. M. Albert, G. P. Ginet, and M. S. Gussenhoven. CCRRES observations of radiation belt protons: 1. Data overview and steady state radial diffusion. J. Geophys. Res. 103, 9261 (1998). [5] R. S. Selesnick, D. N. Baker, A. N. Jaynes, X. Li, S. G. Kanekal, M. K. Hudson and B. T. Kress. Observations of the inner radiation belt: CRAND and trapped solar protons. Journal of Geophysical Research: Space Physics 119(8), 6541 (2014). [6] K. Zhang, X. Li, H. Zhao, Q. Schiller, L.-Y. Khoo, Z. Xiang, R. Selesnick, M. A. Temerin, J. A. Sauvaud. Cosmic Ray Albedo Neutron Decay (CRAND) as a Source of Inner Belt Electrons: Energy Spectrum Study. Geophysical Research Letters, 46, 544 (2018). [7] J. M. Seoane, and M. A. F. Sanju´ an. New developments in classical chaotic scattering. Rep. Prog. Phys. 76, 016001 (2012) [8] R. Liu, S. Liu, F. Zhu, Q. Chen, Y. He, and C. Cai. Orbits of charged particles trapped in a dipole magnetic field. Chaos 32, 043104 (2022). [9] M. Braun. Particle motions in a magnetic field, J. Diff. Eq. 8(2), 294 (1970). [10] C. St¨ ormer. The Polar Aurora (Oxford Clarendon Press, Oxford, 1955). 33 [11] J. Chen and P. J. Palmadesso. Chaos and Nonlinear Dynamics of Single Particle-Orbits in a Magnetotaillike Magnetic Field. J. of Geophys. Res. 91, 1499 (1986). [12] J. Chen. Nonlinear Dynamics of Charged Particles in the Magnetotail. J. of Geophys. Res. 97, 15,011 (1992). [13] O. F. de Alcantara Bonfim, D. J. Griffiths, and S. Hinkley. Chaotic and hyperchaotic motion of a charged particle in a magnetic field dipole. Int. J. Bifurc. Chaos 10, 265 (2000). [14] P.K. Soni, B. Kakad, A. Kakad, Simulation study of motion of charged particles trapped in Earth’s magnetosphere. Advances in Space Research 67, 749 (2021). [15] A. J. Dragt. Trapped orbits in a magnetic dipole field. Rev. Geophys. 3, 255 (1965). [16] A. J. Dragt. Correction to Trapped orbits in a magnetic dipole field. Rev. Geophys. 4, 112 (1966). [17] J. G. Luhmann, L. Friesen. A Simple Model of the Magnetosphere. J. Geophysical Res. 84, 4405 (1979). [18] T. J. Stuchi. Symplectic integrators revisited. Brazilian J. Phys. 32, 958 (2002). [19] W. I. Newman, and A. Y. Lee. Symplectic Integration Methods and Chaos: Timestep Selection and Lyapunov Time. Bulletin of the AAS 37, 531 (2005). [20] A. Calini, N. M. Ercolani, D. W. McLaughlin ans C.M. Schober. Mel’nikov analysis of numerically induced chaos in the nonlinear Schr¨ odinger equation. Physica D 89, 227 (1996). [21] K. P. Rauch and M. Holman. Dynamical chaos in the Wisdo-Holman integrator: origins and solutions. Astron. J. 117, 1087 (1999). [22] B. M. Herbst, M. J. Ablowitz. Numerically Induced Chaos in the Nonlinear Schr¨ odinger Equation. Phys. Rev. Lett. 62, 2065 (1989). [23] J. H. Verner. Numerically optimal Runge–Kutta pairs with interpolants. Numer. Algor. 53, 383 (2010). [24] M. Calvo, D. Hern´ andez-Abreu, J. I. Montijano, and L. R´ andez. On the Preservation of Invariants by Explicit Runge–Kutta Methods. SIAM Journal on Scientific Computing 28, 868 (2006). [25] E. Hairer, C. Lubich, and G. Wanner. Geometric Numerical Integration. Structure Preserving Algorithms for Ordinary Differential Equations. (Springer, Heidelberg, 2006). [26] W. Cai, Y. Gong, and Y. Wang. An explicit and practically invariants-preserving method for conservative systems. (arXiv:2009.06877v1). [27] P. B¨ uhler. Radiation Belts. Proceedings of ESA Workshop on Space Weather, November 1998, ESA WPP-155 (1999). [28] E. Ott and T. Tel, Chaotic scattering: An introduction, Chaos 3, 4 (1993). 34 [29] A. Motter and Y.-C. Lai. Dissipative chaotic scattering. Phys. Rev. E 65, 015205(R) (2001). [30] J. D. Bernal, J. M. Seoane, and M. A. F. Sanjuan. Global relativistic effects in chaotic scattering. Phys. Rev. E 95, 032205 (2017). [31] F. Blesa, J. D. Bernal, J. M. Seoane, and M. A. F. Sanju´ an. Relativistic chaotic scattering: Unveiling scaling laws for trapped trajectories. Phys. Rev. E 109, 044204 (2024). [32] D. N. Baker. Wave–particle interaction effects in the Van Allen belts. Earth, Planets and Space 73, 1 (2021). [33] J. Naugle, D. Kniffen. The Flux and Energy Spectra of the Protons in the Inner Van Allen Belt. Phys. Rev. Lett. 7, 3 (1961). [34] J.F. Ripoll, V. Loridan, M.H. Denton, G. Cunningham, G. Reeves, G., O. Santol´ ık, J. Fennell, D.L. Turner, A.Y. Drozdov, J.S. Cervantes Villa, Y.Y. Shprits, S.A. Thaller, W.S. Kurth, C.A. Kletzing, M.G. Henderson, and A.Y. Ukhorskiy. Observations and Fokker-Planck Simulations of the L-Shell, Energy, and Pitch Angle Structure of Earth’s Electron Radiation Belts During Quiet Times. J. Geophysical Res. 124, 1125 (2019). 35