Modeling of ion dynamics in the inner geospace during enhanced magnetospheric activity

We investigate the effect of magnetic disturbances on the ring current buildup and the dynamics of the current systems in the inner geospace by means of numerical simulations of ion orbits during enhanced magnetospheric activity. For this purpose, we developed a particle-tracing model that solves for the ion motion in a dynamic geomagnetic field and an electric field due to convection, corotation and Faraday induction and which mimics reconfigurations typical to such events. The kinematic data of the test particles is used for analyzing the dependence of the system on the initial conditions, as well as for mapping the different ion species to the magnetospheric currents. Furthermore, an estimation of Dst is given in terms of the ensemble-averaged ring and tail currents. The presented model may serve as a tool in a Sunto-Earth modeling chain of major solar eruptions, providing an estimation of the inner geospace response.


Introduction
During each solar cycle, sequences of eruptive flares are followed by coronal mass ejections and interplanetary shocks, some of which arrive near Earth.At times when the solar wind enters into Earth's magnetosphere, these solar eruptions modify the dynamic conditions in geospace and trigger space weather effects like geomagnetic storms and magnetospheric substorms (Daglis, 2004;Schwenn, 2006;Pulkkinen, 2007).Geomagnetic storms occur when the energy transfer from the Sun to geospace intensifies, as a result of the occurrence of magnetic reconnection at the dayside magnetopause during periods when the interplanetary magnetic field (IMF) has a strong and prolonged southward component (Akasofu, 1981).Magnetospheric substorms are caused by the variability in the north-south orientation of the IMF and evolve as energy loading-dissipation cycles (Baker et al., 1996).This kind of activity brings up a configuration change in the magnetosphere, including the ionosphere, and enhances the ring current and corresponding current systems flowing on the magnetopause, along the magnetotail and in the Birkeland regions (see in Pulkkinen et al., 2005).The associated dynamic processes evolve in a variety of timescales, from days for geomagnetic storms and hours for magnetospheric substorms down to minutes, or even seconds, for local plasma instabilities.
Solar eruptions are characterized as geoeffective when the magnetospheric response amounts to large electromagnetic perturbations, with severe consequences for the performance of ground-based power and communication networks as well as for spacecraft and weather satellites (Daglis, 2004).A major goal in space weather research is to predict the dynamic state of the geospace from measured solar wind and IMF data, so as to timely distinguish those events that are harmful.In this respect, the simulation of physical processes dominating extreme space weather conditions, such as magnetic reconnection, convective plasma transport and charged particle acceleration, is required (Daglis et al., 2009).For numerous events, magnetospheric activity can be described by means of a few geomagnetic indices, like Kp and Dst, which can in principle be derived from solar wind and IMF values.However, these indices include systematic and/or statistical errors which limit the capability to establish consistent (causal) correlations (Rostoker, 1972).Therefore, more detailed, large-scale numerical solvers of the coupled solar wind-magnetosphere system may be employed; such models have advanced with the increased availability of computer resources.
There are cases where global fluid and magnetohydrodynamic (MHD) simulations reproduce the observed changes in the magnetic topology to quite good accuracy (Moretto et al., 2006;Honkonen et al., 2013).However, their use has limitations due to missing physics for the description of non-collisional processes in a multi-species plasma.Kinetic solvers and test-particle simulations, with a description of the plasma motions in adjustable physics detail, increase the reliability at smaller scales by properly addressing effects like thermal instabilities and anomalous transport (Buneman et al., 1992;Khazanov, 2010).The drawback of microscopic models for global simulations is the large demand on computer resources; to cope with this, it is customary to separately model each source playing an important role in the dynamics: the ring current, the near-Earth tail currents, the radiation belts and the magnetopause.In this frame, the origin and transport of ring current and radiation belt particles during storm time, their interaction with the tail current, the escape of high-energy particles and the dynamic connection with the substorm phases have been addressed through the analysis of observations and dedicated numerical simulations (Chen and Schulz, 1996).
For the description of the ring current dynamics, the plasma current distributions at the near-Earth region have been modeled in terms of the bounce-averaged, drift-kinetic equation.Fok et al. (2001) described the particle drifts in the storm-time field in terms of the initial and boundary particle distributions, with the coefficients in the kinetic equation calculated from the Hamiltonian description of motion.Jordanova et al. (2010) modeled the radially diffusive plasma dynamics in self-consistence with the fields by coupling the kinetic equation with a 3-D, force-balanced magnetic equilibrium code and a MHD solver for the convection electric field.In another self-consistent treatment, Lemon et al. (2004) employed a collisionless kinetic code together with a model for the electrostatic potential, taking into account the current closure with the ionosphere.The specific model has been coupled to the code of Fok et al., where it serves as the solver for the electric field.
The method we adopt in this work is to directly follow the 3-D particle trajectories under the effect of the electric and magnetic forces during the dynamic phases of the disturbance (see, for example, Delcourt et al., 1990;Ganushkina and Pulkkinen, 2002;Ebihara et al., 2003).An advantage of studying the individual particle motions is the physics insight gained, as well as the statistics built from ensembles of particles.In such models, the Lorentz equation of motion is solved, either in its full form or reduced in terms of the guiding-center approximation, and the driving forces are (as above) the dynamic magnetic field coming from the superposition of the Earth's terrestrial magnet with the fields gen-erated by the magnetospheric current sources (Tsyganenko, 2013), and the electric field due to large-scale plasma convection and corotation with the Earth (Volland, 1973).An important factor, however, is the modeling of the electric field component induced by the time variation in the magnetic field.The specific field is involved in the strong acceleration of charged particles which is observed during geomagnetic disturbances; however, there are a relatively low number of test-particle-based studies which have been performed in this direction (like, for example, Delcourt, 2002).
It becomes apparent that the modeling of the near-Earth plasma response to geoeffective solar events is of major importance for the improvement of space weather prediction.The model requirements are a consistent description of the geomagnetic and electric fields, the computation of the Sun-driven plasma dynamics and the assessment of the numerical data for the estimation of parameters related to space weather, including benchmarks against ground-based and satellite observations.In this paper, we present results from the simulation of the electric and magnetic fields and of the energetic particles in the inner magnetosphere, focusing on the ring current buildup and decay when disturbances are occurring.The physics of our model for the forces driving the plasma dynamics are cast in a form suitable for use with 3-D test-particle codes.Provided that there are suitable simulation data, a statistical evaluation for the dynamics of the different ion types is performed over the initial conditions, and an ensemble-averaged estimation of the Dst index stemming from the ring and tail current populations is given.
The structure of the paper is as follows: in Sect.2, the physics model for the geomagnetic and electric fields is explained, accompanied by field-line tracing and equipotential contour simulations, and, following that, we describe the main aspects of the particle-tracing model.In Sect. 3 we present the numerical results: the different types of ion motion found in the disturbed magnetosphere, the statistical analysis of the particle dynamics and the estimation of the Dst index.Finally, in the concluding section, the merits of this work are summarized, the limitations of our model are discussed and further studies are proposed.

Geomagnetic field
The magnetic field in geospace is expressed as the sum of two contributions: the first one is from the Earth's terrestrial field, whereas the second comes from the external field generated by the electric currents flowing inside the magnetosphere (including the magnetopause).The Earth's magnetic field is well approximated as the one of a tilted dipole magnet with inverse polarity (Parks, 1991).In geocentric solar magnetospheric (GSM) Cartesian coordinates, the expression of the dipole field is where R E = 6378 km is the Earth radius, B E = 31 000 nT is the value of the magnetic field on the surface, r = r/r is the direction vector and θ t is the tilt angle.We note here that, since B ter varies very slowly (through B E and θ t ) in comparison to the solar activity and its geomagnetic response, it may be considered time-independent in the context of our computations.
The second component of the geomagnetic field, denoted by B ext , is generated by the electric currents which result from the interaction of the magnetospheric plasma with the solar wind (Parks, 1991).The most important components are (i) The magnetopause current, which is controlled by the solar wind's dynamic pressure P dyn , (ii) the magnetotail currents, which extend from 10 R E to well beyond 100 R E , and (iii) the ring current, flowing around the Earth inside a toroidal band approximately within [3,9] R E .The external field's spatial dependence is defined by the distribution of the current sources, and its time dependence by the evolution of these sources during quiet time and the events.
The mainstream of models for the static part of B, due to the Tsyganenko algorithms T89, T96 and TS05 (see Tsyganenko, 2013, and references therein), follows a data-based approach towards correlation with parameters of geomagnetic activity like P dyn , the IMF vector B, the planetary index Kp and the disturbance storm-time index Dst (Rostoker, 1972).In T89, a physics-based description of the magnetospheric currents and the corresponding vector potential was introduced; the T96 model improved T89 in the description of the magnetopause geometry and the equatorial tail physics, whereas TS05 upgraded T96 with the inclusion of storm and substorm dynamics.All models require specific parameter values at input: in T89, Kp and θ t are to be given; in T96, apart from θ t , it is P dyn , Dst and B y , B z , whereas TS05 requires the input of T96 plus six additional parameters, S 1 , S 2 , . .., S 6 , related to the storm-time effects.
For the visualization of the magnetic field, one employs the standard field-line tracing technique (Parks, 1991).In 2-D, the field-line map is a clear picture only on the x − z planes, as a result from the existing symmetries of the magnetosphere's geometry in the GSM system: the x axis is the line connecting Sun and Earth, whereas the z axis may always be placed on the magnetic dipole axis.In Fig. 1 we show the map of the total magnetic field on the x − z plane defined by the meridian y = 0, as computed with the T89 model.We present two cases with different values of the Kp index, one relevant to quiet time (Kp = 1) and one reflecting storm-time conditions (Kp = 5), for the typical inclination of the Earth's dipole (θ t = 11.5 • ).The typical properties of the geomagnetic field appear in the results; for example, in the second case, where Kp is larger, the field lines are more dense close to the Equator and towards Earth due to the rise in the convection intensity.
For proper application of the Tsyganenko models, the role of the differences between the models and the properties of the computed physics, especially in strongly disturbed cases, has to be investigated.The benchmark of these models against observations is an issue that has been addressed by a number of authors.Woodfield et al. (2007) performed a comparison of T89 and T96 with magnetic field data from the Cluster mission, and the results have shown noticeable deviations only in the outer ring current region on the nightside and near the cusp.Boschini et al. (2013) utilized T96 and TS05 in a geomagnetic backtracing code and benchmarked against AMS-02 data, finding significant differences near Earth only for storm conditions (Kp > 5 or P dyn > 3 nPa).Also, McCollough et al. (2008) performed a detailed statistical comparison of all established models and, based on the results, the use of models that include magnetospheric asymmetry is encouraged when Kp > 4 in regions including the dayside and the dawn-dusk neighborhood.
A comparison of T89, T96 and TS05 is presented in Fig. 2. The reference case involves a strongly perturbed, non-tilted dipole, and the exact input for each model is given in Table 1.The differences in the computation of B by the different models are quantified in terms of the relative de- Within the limits set by the differences in the input of T89 with respect to the other models, the quantitative comparison does not exhibit very large deviations in most of the region of interest (approximately within [−25R E , 10R E ] along x and [−15R E , 15R E ] along z, and always inside the magnetopause boundary).Noticeable deviations, ranging from 50 % to 150 %, appear in the outer region on the nightside, near the cusp on the dayside and in the far-Earth magnetopause, which is in agreement with the benchmarks presented above.Consequently, the specific choice of field model is not expected to play a crucial role in the test-particle results.
The dynamic part of the magnetic field is determined by the modification of the geomagnetic parameters in time.With introduction of the vector G = [G j ], with components the input parameters required for each Tsyganenko model (e.g., G = [θ t ,Kp] for T89), the partial derivative of where the functions G j (t) may be specified analytically, in terms of an approximation by continuous functions, or directly as a time series of observations (the derivatives then being computed as discrete-time finite differences).In principle, the terms ∂B ext /∂G j are not available in analytic form; these could be discretized and computed by repetitive usage of the numerical field model for the parameter values of interest, but, in this fashion, the computing cost heavily increases.
In order to simplify the computation, these terms are approximated by the variation in B ext within the start and stop times of the event t i and t f , and by a normalized profile function F B (t).In this frame, the total field is expressed as In our model, an event starts at time t i , stops at t f and, during this interval, it evolves in phases described by the function F B and the values of G at t i , t f .F B (t) is defined on the basis of the properties of the magnetic field, as observed in measurements.Here, we refer to events which have an initial "growth" period where the field strength is increasing to high values, followed by a (shorter) "relaxation" phase where B returns to its previous levels (Metallinou, 2008).To this end, F B is chosen to be In the above, ) is a product of Heaviside step functions and t g is the time stamp of the growth phase, whereas w j are fitting coefficients.In Fig. 3 we illustrate F B (t) for a substorm event with time stamps t i = 0 min, t g = 30 min, and t f = t r = 35 min, and five fitting parameters, w 1 = w 2 = 0, w 3 = 10, w 4 = −15, and w 5 = 6.

Electric field
The electric field is divided into three components (Russell, 2000): the first one is due to plasma convection in the magnetosphere, the second stems from near-Earth plasma corotation, and the third one is generated by the dynamic variation in the geomagnetic field during the events.The slow timescale of the convection and corotation processes in comparison to the overall plasma dynamics allows for their consideration as electrostatic.A variety of models have been developed for the calculation of the electrostatic potential cc that generates the convection-corotation field: (i) the Volland-Stern-Maynard-Chen (VSMC) model (Volland, 1973;Stern, 1975), based on an empirical dawn-dusk potential distribution with Kp dependence and magnetopause shielding; (ii) the E5D model, a transpolar, Kp-driven analytical approximation of the convection potential (McIlwain, 1986); (iii) the Boyle-Reiff-Hairston (BRH) model, which describes the convection field with a polar-cap potential function driven by the solar wind and the IMF (Boyle et al., 1997); and (iv) the Weimer (WM) model, which is derived from a combination of low-altitude measurements of the convection velocities at high latitudes (Weimer, 2005).
The effect of the model differences to the computed dynamics in the transition to stormy conditions is again put under question.In Khazanov et al. (2004), against the background of plasma kinetic simulations, the BRH and VSMC models were compared and the results did not yield measurable differences, except from regions near the magnetopause and the distant tail.In the same manner, in the context of MHD plasmapause simulations (Pierrard et al., 2008), the comparisons between the VSMC, E5D and WM models con-cluded in a similar picture for the near-Earth convection.An indirect benchmark of the VSMC and E5D models was performed by using these, together with the Tsyganenko models, as input to gyro-particle simulations (Woelffle et al., 2011).It was shown that the differences in the magnetic field do not influence the computation as much as the ones in the electric field, which was highlighted by significant variations in the particle trajectory shape and the energy variation during transport.One concludes that the choice of model should be made according to the performance under conditions implied by the event under study; for example, VSMC offers a good global description of transport in the plasma sheet, whereas E5D predicts the magnetopause position better.
From the aforementioned tools we choose to employ the VSMC model, which combines sufficient accuracy in the physics description with simplicity in the computer imple- In Eq. ( 5), ω E = 2π/24 rad h −1 is Earth's rotation frequency; γ is the magnetopause shielding factor; and υ 1 , υ 2 , and υ 3 are constant parameters, which are calculated in terms of data fitting over magnetic field measurements in the inner tail region.On the right-hand side of Eq. ( 5), the leftmost term represents the potential for the convection field, in which the fraction involving Kp determines the field intensity, whereas the rightmost term is the potential generating the corotation field.
Vector fields coming from a scalar potential are represented in terms of their equipotential (contour) surfaces.The contour surfaces of cc are calculated by solving Eq. ( 5) with respect to the coordinates on a certain potential level, i.e. cc (x, y, z) = l .In Fig. 4 we perform a 2-D visualization of the contour lines on the equatorial plane (z = 0), for geospace-related parameter values γ = 2, υ 1 = 0.045 kV m −2 , υ 2 = 0.0093 and υ 3 = −0.159, in two cases of solar activity level with different intensity: (a) quiet time (Kp = 1) and (b) disturbed (Kp = 5).The main physics properties of the convection field are well reproduced by the model, like, for example, the global increase in the field values as Kp increases.
The role of the electric field component induced by the dynamic variation in the magnetic field in properly modeling the solar-driven perturbations is very important.This is due to the fact that it has a short space/timescale, which is effective in accelerating ions to very high energies (as observed during storms and substorms), whereas the convection process forms a distribution of plasma currents of comparatively low energy.In this context, the total electric field is expressed in terms of the potentials According to Eq. ( 6), the calculation requires knowledge of the vector potential A ext , which is the generating function of It is known, however, that, given arbitrary magnetic field, an analytic solution for the vector potential is, in most cases, not possible.T89 involves simplifications in the description of the plasma current sources which allow the analytic calculation of A ext , whereas the later models are based on a more complicated formulation, including spherical harmonic expansion and integrals of special functions, and thus cannot fall in this category.

Particle tracing
The test-particle model computes the near-Earth ion dynamics during the geomagnetic disturbance by following the 3-D trajectories under the effect of the associated electric and magnetic fields.The particle trajectory is traced by solving numerically the Lorentz equation including the gravitational force In Eq. ( 7), g = g ER 2 E r/r 2 is the gravitational acceleration (g E = 9.81 m s −2 , its value on Earth's surface) and m and q are the particle mass and electric charge.For electrons it is m e = 9.31 × 10 −31 kg and q e = −1.6 × 10 −19 Cb, while for an ion of atomic mass A i and ionization rate s i it is m i = 1837A i m e and q i = s i |q e |.
The particle motions may also be evaluated in terms of a reduction to the full model, depending on the validity of the guiding-center (GC) approximation over the simulated region.The GC trajectory describes the overall motion well in cases where the electric/magnetic field variations remain sufficiently small over each revolution (Parks, 1991).This translates to relations of the Larmor radius ρ L and the rotation frequency f L with the spatiotemporal scales of E, B Equation ( 8) suggests that the GC approach is invalid when the field-line curvature is comparable to the Larmor radius, as well as for heavy ions that exhibit large periods of gyration.
If the approximation is valid, the GC equation is employed in the following form (Northrop, 1961) with v gc the velocity of the GC (the symbols || and ⊥ refer to the parallel and perpendicular components) and µ gc = v 2 gc,⊥ /(2B) the particle's magnetic moment, which here is an adiabatic invariant.In Eq. ( 9), the terms on the right-hand side refer to the effect of the electric field, the gravitational force, the magnetic field gradient and the magnetic curvature on the GC drift motion.
The particle-tracing scheme combines the models presented so far: T89 is employed for the static part of B and the function of Eq. ( 4) for its time variation, VSMC is used for the electric field due to convection and the induced part is computed on the guidelines described near (Eq.6), and the particle motion is followed by solving the Lorentz equation or by adopting the GC model, with the ability to interplay between the two orbit solvers.The computation is interrupted if the particle leaves far from the inner magnetosphere, either by crashing onto Earth (r ≤ R E ), crossing the magnetopause or reaching a tailward distance further than 70R E , with different stop codes so that each case is distinguished.
The orbits may be traced with the Lorentz equation, with no simplification adopted at any stage of the computation.The GC model, in the regions where it is valid according to the conditions (Eq.8), is an efficient method to provide a simpler trajectory calculation.In such a scheme, in principle, the GC conditions of validity should be checked at every time step and, depending on the outcome, the physics model to be applied should be chosen.However, since this tactic radically decreases the code speed, in practice the orbit solver is interchanged whenever the radial position of the particle becomes less that an empirically set limit R fm .In this frame, an issue which should be investigated is the differences in the orbits with respect to the computation using the full model, particularly in conditions of amplified disturbances.
In the literature, the benchmarking between the different magnetospheric particle solvers in the presence of intense electric and magnetic fields is not sufficiently extensive and the results are contradictory.Cladis and Francis (1989) performed a comparison of the GC and the Lorentz solutions for heavy ions under the effect of a geomagnetic field given by the T89 model, and the results showed acceptable deviations in the particle flux rates and the drift paths.However, Shibahara and Nose (2009) computed energetic ion motions using the TS05 and VSMC models for the fields, and found measurable differences in the occurrence of large ion gyroradii and pitch angle values close to π/2.In such cases, some of the adiabatic invariants are broken and the validity of the GC approximation becomes questionable.
In order to clarify this issue, we compute a specific trajectory for several values of the threshold distance R fm , and compare the results in Fig. 5.We evaluate the orbit with R fm ranging from R E (full orbit) up to 60R E , and for a plain GC simulation we set R fm = 100R E .In Fig. 5a and 5b, the position r and the kinetic energy E k are plotted vs. t for several values of R fm .One observes that the deviations between the Lorentz, GC and mixed computations appear after the event has ended and are measurable for E k and less important for r.The deviation of the results for E k from intermediate values of R fm exhibits an irregular behavior; indicatively, using R fm = 18R E one is still close to the full model, whereas for R fm = 10R E the deviation is larger and for R fm = 8R E the particle follows a completely different orbit.This picture is verified by Figure 5c, where the maximum relative error from all orbit quantities is computed as a function of the threshold radius.Due the sensitivity of the results on the interchanging procedure, one should be cautious with the choice of R fm ; for this, we choose to use only the Lorentz model in order to provide the most reliable approach.

Numerical results
In this section, the results from test-particle simulations are shown and analyzed.The space weather scenario under study involves the occurrence of a single magnetospheric disturbance.The growth phase of the event starts at t i = 0 with a quiet magnetosphere, indexed with Kp (t i ) = 1, and completes after t g − t i = 30 min by reaching a disturbed state with index Kp (t g ) = 5.Then, the relaxation phase follows immediately and completes after t r − t g = 5 min, during which Kp returns to its initial value, i.e.Kp (t r ) = Kp(t i ).The particle starts its flight at the time stamp t 0 , which may be before (t 0 − t i < 0) or after (t 0 − t i > 0) the onset of the growth phase, interacts with the disturbance until t r and continues moving under the effect of the restored fields until the time stamp t 1 .
In the disturbed magnetosphere, three primary types of ion trajectories are met: (i) orbits which become trapped inside the ring current, (ii) orbits that precipitate into Earth's atmosphere, and (iii) orbits escaping tailward or by crossing the magnetopause.We have computed these types by sampling various initial conditions for the ion position and energy, and the results are shown in GSM coordinates in Figs. 6 and 7.In Fig. 6, we have the planar projections of the orbit of an O + ion that eventually integrates to the ring current.The motion initiates at t 0 = −8 min, before the event, with initial radius r(t 0 ) = 20R E , magnetic local time (MLT) φ(t 0 ) = 24 h, latitude θ (t 0 ) = 25 • , pitch angle α(t 0 ) = π/2 and kinetic energy E k (t 0 ) = 4 keV, and is followed until t 1 = 180 min (∼ 2.5 h after the event).In Fig. 7, we represent the other types of motion in 3-D space for different ion species with the same input as before except for E k : the precipitating orbit is of an H + ion with E k (t 0 ) = 0.5 keV, whereas the escaping orbit is of an O + ion with E k (t 0 ) = 7 keV.
In Fig. 6, the O + ion is launched from the plasma sheet, driven towards Earth by the disturbed electric fields, and finally gets trapped in the ring current.A careful examination of the numerical data yields that, in its course to the ring current region, the ion is considerably accelerated from the energy exchange with the electric field, whereas its pitch angle has a random behavior before the entrance to the ring current and afterwards varies periodically.In the case of the H + ion that crashes onto Earth, the orbit of which is shown in Fig. 7a, the particle begins with a low initial energy and is intensely accelerated, and probably due to the relation of its pitch angle with the loss cone it ends up on the terrestrial atmosphere at t = 33 min (well before t 1 ), with a final energy as large as in the previous case.Finally, in Fig. 7b, the O + ion starts with a relatively high value of energy and escapes from the inner magnetosphere, along the meridian at 22:00 MLT, before t 1 (at t = 40 min) with a velocity gain.The different behavior of the oxygen ions for different initial energy is an indicator of the sensitivity of the ion dynamics to the initial conditions.
Regarding the particle acceleration, in Fig. 8 we examine the kinetic energy and the pitch angle for the motions in Figs. 6 and 7.In Fig. 8a, the energy of the trapped O + ion appears to have a gain of about 1.3 orders of magnitude.The largest part of the increase occurs during the relaxation phase (t g < t < t r ), where the magnetic field exhibits a steep decrease and, consequently, the induced electric field attains large values and accelerates the ions.Another incidence of energy gain occurs a little after t = 2 h well inside the ring current region.This is connected to an intense pitch angle variation, as seen in Fig. 8b, where α(t) is displayed, which is induced by the structure (i.e. the gradients over time and space) of the local fields at the specific time.The kinetic energy of the H + ion that crashes onto Earth also appears to have a sizeable gain (almost 2.6 orders of magnitude) at the time of reaching the atmosphere, after nearly 45 min of flight, whereas the energy of the escaping O + ion appears a gain of nearly 1.6 orders of magnitude at the time it crosses the magnetopause, close to the end of the relaxation phase.
For the statistical analysis of the motions, numerical data have been produced over the trajectories of the ion species keV, moving until t 1 = 3 h during a disturbance with t g = 30 min, t r = 35 min, Kp (t i ) = 1, and Kp (t g ) = 5.
relevant to each territory of the magnetosphere, in loops of different initial conditions for r, φ, θ , α and E k where, each time, only one or more of these quantities were varied.In Fig. 9 we present results for the final kinetic energy and the pitch angle from different simulations with an ensemble of N ens = 1000 O + ions, where the initial conditions varied are r(t 0 ) and E k (t 0 ).In the first computation, r(t 0 ) took values in a loop from 2 R E to 30 R E , whereas, in the other case, E k (t 0 ) ranged from 0.5 to 20 keV, and all the rest of the input quantities were equal to the values already defined: t 0 = −8 min, φ(t 0 ) = 24 h, θ (t 0 ) = 25 • , and α(t 0 ) = π/2.The general picture (also implied from the above results) is that the dependence of the particle dynamics on the initial conditions is very sensitive, which is imprinted in the wide regions where E k (t 1 ) and α(t 1 ) exhibit an irregular, non-smooth variation over the values at t = t 0 (see especially Fig. 9a).However, in all cases one identifies consecutive regions where the ions either get accelerated or remain at low energy.In Fig. 9a, there is a spatial region from 14 R E to almost 17 R E where all injected ions gain significant amounts of energy, as well as one within 21 R E and 25 R E where nearly all particles do not appear to have a net energization.In Fig. 9b the probability of acceleration appears to be larger for O + ions with low energy at the event start (E k (t 0 ) < 6 keV) than for initially energetic ions (having, for example, E k (t 0 ) > 15 keV).The specific result reveals the role of the plasma sheet as a reservoir of oxygen ions which, on the course of storms/substorms, get accelerated and enhance the ring current (Metallinou, 2008, and references therein).Finally, Fig. 9c, there are regions where α(t 0 ) varies rapidly, suggesting rotational behavior, as well as regions where the variation is slow, denoting motions close to ballistic.Most of the former ions, as implied by the inbound direction of motion driven by the large parallel velocities of the specific pitch angles, are candidates of joining the ring current.
We also analyze the kinematics of H + ions using an ensemble of 1000 particles with varying r(t 0 ) and E k (t 0 ) for the same input as above.In Fig. 10, we plot the final kinetic energy and pitch angle as a function of the initial values.The overall behavior resembles that of the (heavier) O + ions; notice, for example, the regions of quasiperiodic and quasiballistic motion in Fig. 10c, similar to the ones in Fig. 9c.Nevertheless, the effect of acceleration, as imprinted in Fig. 10a, is found to be much weaker.This is connected to the difference in charge / mass ratio of the different ion species, and verifies the known storm-time composition for the energy density of the ring current, which is dominated by the O + ions coming from the plasma sheet, in contrast to the situation in quiet time where H + is the majority species (see, for example, Korth et al., 2002).In Fig. 10c one observes that the regions of periodic-like pitch angle behavior are now more narrow, which is in accordance with the fact that hydrogen ions launched from the plasma sheet are not effective in assimilating into the ring current region.Going one step further, we estimate the statistical weight of each one of the populations formed by distributing the ions launched from the plasma sheet to the types of orbits described in the above (ring current, near-Earth tail, precipitating and escaping), and the outcome is shown in Table 2.The simulations involved two different ensembles of oxygen and hydrogen ions with N ens = 10 000 particles each, which were injected from the plasma sheet with random initial conditions for r(t 0 ) within [18R E , 22R E ], for E k (t 0 ) within [2, 6] keV and for θ (t 0 ) within [20,30] • , and all the other input being the same as above.One should notice that the conclusions drawn from Figs. 9 and 10, on the basis of single-particle dynamics, are verified.We highlight that, according to the computations, 46 % of the O + ions of the ensemble are incorporated into the ring current, in contrast to nearly 4 % of the H + ions, whereas a little more than 20 % in both species occupy the near-tail region; however, 25 % of the O + particles and 39 % of the H + ones escape the inner geospace.
The magnetic perturbation and the connection of Dst to the ring curren, as well as the contribution of each current source during the event phases, are under debate.In many cases, Dst is assumed to be correlated with the ring current energy from storm maximum well into recovery, on the basis that ring current ions provide the primary contribution to the storm-time Dst depression (Greenspan and Hamilton, 2000).However, it is suggested that Dst is also related to other sources, the effect of which may become important during disturbances.Based on ground measurements, Arykov and Maltsev (1993) indicate circumstances where the tail currents dominate the Dst development during storms.Turner et al. (2000) assess the effect of the tail currents on Dst by introducing a correction to the total current density, in terms of subtracting the magnetic curl in the tail regions as calculated by T89 and T96.The tail current was found to be most dominant in the end of the growth phase, and the accretion to Dst was estimated to scale up to 25 %.
A straightforward approach to calculate the Dst index from test particles involves the computation of the electric current densities from the particle velocities and the derivation of the generated magnetic fields (according to Ampere's law).However, the increased difficulty in the computation of surface current densities from particle orbits and the complexity of calculating the magnetic field from the electric currents, as well as the requirement of including the real positions of the ground-based sensors, imply a poor modeling performance.For studies related to the inner magnetosphere, the connection of Dst with the energy of the ring current has been described in terms of the Dessler-Parker-Sckopke (DPS) relation (Dessler and Parker, 1959;Sckopke, 1966).The advantage here is that, at input, the kinetic energy of the local plasma is required, which is a scalar quantity and simple to deduct from the test-particle results.
The original DPS relation, which connects the energy E pp stored in a specific plasma population of the magnetosphere with the associated magnetic field perturbation, takes into account only those energetic particles that gyrate around the magnetic field lines and drift longitudinally due to the field gradient.In this context, the DPS formalism provides a sufficient estimation of the perturbations due to the ring current dynamics, as well as a well-balanced one of the near-Earth tail current contribution.With the introduction of a correction term for the magnetopause current in the original relation (O'Brien and McPherron, 2000), one obtains a modified equation for the Dst index, where b dps is associated to the magnetopause correction and c dps to the quiet-time energy level.Equation ( 10) yields that, in order to estimate Dst, with the values of P dyn , b dps and c dps for a specific event or scenario, one only requires the computation of E pp for the ring and tail plasma populations.
In the simulations, the ring current particles are assumed to be confined inside a torus with radii R rc = 6R E and r rc = 3R E (i.e.extending from 3 R E to 9 R E ), whereas the near-Earth tail region is defined as the remaining area in the simulation box ranging within [r rc , r rc +R ntl ] along the Sun-Earth axis and [−R ntl , R ntl ] in the other two directions, with R ntl = 20R E .The energies E rc and E ntl of the ring and near-tail current particles are described over the average energy of H + and O + test ions that belong to these currents.This is quantified by E pp = V pp i n pp,i E k pp,i , where n pp,i is the plasma density of each ion species in each population and V pp is the volumes of the regions occupied by the plasma populations, accordingly given by V rc = 2π 2 R rc r 2 rc and V ntl = R 3 ntl − V rc .In the formula for E pp , the (ensemble) average value of the kinetic energy for the O + and H + ions in each current is computed, at each time step, over the particles contained inside the specific region at that time.
The results of the Dst computation using the scheme described above are presented in Fig. 11.The event scenario explored is again the one introduced in Sect.3, i.e. a single disturbance that begins at t = 0 from a quiet state with Kp = 1, reaches its maximum level Kp = 5 at t = 30 min and returns to the quiet state until t = 35 min.Two thousand test particles were used for the computation, divided in two different ensembles: one of 1000 oxygen ions launched from the plasma sheet, and one of 1000 hydrogen ions started in the ring current.The initial conditions for the MLT and the pitch angle were the same for both species, φ(t 0 ) = 24 h and α(t 0 ) = 90 • , whereas the initial radii, latitudes and kinetic energies were assigned randomly within different ranges for each species: for O + , r(t 0 ) ranged in [18R E , 22R E ], θ (t 0 ) in [20,30] • and E k (t 0 ) in [2, 6] keV, while for H + the corresponding intervals were [5R E , 7R E ], [0, 5] • and [1, 3] keV.The test ions were traced from 8 min before the beginning of the event until t = 120 min, and at each time step the Dst index was computed from Eq. ( 10) and associated relations, where it was assumed that P dyn = 4 nPa, b dps = 7.26 nT/nPa 1/2 , c dps = −11 nT, n rc,O = n rc,H = 10 cm −3 and n ntl,O = n ntl,H = 1 cm −3 .
In Fig. 11a we plot the dynamic evolution of Dst, and the contributions from the ring current and the near-Earth tail plasma distributions to its value are given for comparisons.A qualitative agreement with the usual evolution of Dst during a substorm is seen: during the growth phase, Dst decreases continuously, with the most rapid variation being around the interchange from growth to relaxation phase, and values of Dst indicating magnetic perturbation persist for some time after the event termination.Overall, the contribution of the ring current to Dst is larger than the one coming from the energetic particles in the near-Earth tail region.The contribution of the tail current is measurable up to t = 1.2 h, with a maximum near the end of the event growth (t = 0.55 h), and from there on the Dst is essentially determined only by the ring current; this is in agreement with the behavior stated in Arykov and Maltsev (1993) and Greenspan and Hamilton (2000).The contribution to Dst by the tail current is found equal to 30 % on the average.This is a little larger than the reported figure of 25 % in the literature; however, such deviations are justified considering the adopted assumptions in these models.
A comparison of Kp and Dst during geomagnetic events is necessary for assessing their differences in response to different storm-time current systems.In cases where the dynamic pressure and the IMF both refer to the same category in storm magnitude, the minimum Dst is expected to decrease as a function of Kp.This has been verified in terms of an additional computation, where the maximum value of Kp during the event, attained right at the end of the growth phase, was modified from 1 to 7 (these are the minimum and maximum disturbance levels allowed by the T89 model) and, in each case, the minimum Dst value was recorded.
The correlation of the maximum values of |Dst| and Kp is shown in Fig. 11b and, as expected, has an increasing monotony.In the same figure, our result is compared to the linear regression curves derived from the statistical evaluation of data from substorm events during 1987-1996(Rostocker, 2000) ) and storms in the period 1996-1999 (Huttunen et al., 2002).The comparison with the results of Rostocker shows an agreement only in the range of values 4 < Kp < 5, and with Huttunen et al.only for Kp > 5+, which correspond to moderate and intense events.The main sources of disagreement in the other ranges are probably connected to the difference of the reference values of the geomagnetic disturbance (dynamic pressure, IMF, plasma density) in the analyzed data with respect to the input given to the code, as well as to the differences with the corresponding values in the data set employed by the T89 model.

Discussion and conclusions
In this paper, we employ a collection of models for the electric and magnetic field in the inner magnetosphere for the investigation of the dynamic evolution of the ring current and the near-Earth tail ion population during the occurrence of magnetospheric disturbances.Within this research framework, we have developed an orbit-solving code which computes the test-particle motion due to convection, corotation and Faraday induction in the dynamic magnetic and electric fields of the magnetosphere.We have used the code to study the ion dynamics, and in particular the dependence of ion acceleration on the initial conditions.Furthermore, we performed a numerical estimation of the Dst index based on the test-particle energies.The results of all computations have been found to be in qualitative agreement with previous studies on the ring current evolution during magnetospheric activity.
The ion motions have been traced by solving the nonrelativistic Lorentz equation, without adopting simplifications at any stage of the computation.In practice, one usually shifts to the GC equations when the particle reaches a distance smaller than a threshold radius, from where on the GC approximation is empirically assumed to be valid.In this respect, the choice of retaining the full-orbit description prevents inaccuracies from occurring in cases where some of the adiabatic invariants are broken.During intense disturbances, such cases have the potential to occur locally in space/time, and we have verified this situation by finding major deviations in the computation of a specific trajectory for several values of the threshold distance.
The analysis of test-particle orbits reveals fragments of the ion dynamics during the disturbance.We have identified three main types of ion orbits: orbits getting trapped around Earth, orbits precipitating in the Earth's atmosphere, and others escaping from the inner geospace.During the event, a percentage of oxygen ions launched from the plasma sheet are found to be accelerated and become trapped in the ring current.However, hydrogen ions (which are known to populate the ring current during quiet times), mainly escape from the inner geospace when launched from the plasma sheet.The largest part of the O + acceleration occurs during the relaxation phase, where the magnetic field exhibits a steep decrease and, consequently, the induced electric field attains large values.The addition of this component to the convection field provides a mechanism for the observed energization levels of ions which drift towards the ring current region, contrary to an electric field purely due to plasma convection (a similar result was found in Fok et al., 1999).
Further analysis of the ion motions reveals a sensitive dependence of the particle dynamics on the initial conditions.We have found regions in geospace, including the plasma sheet, from where injected oxygen ions get preferentially accelerated, while ions starting from other regions may or may not appear a net energization depending on the initial energy.For O + launched from the plasma sheet, the possibility for acceleration is found to be larger for ions having low energy at the beginning of the event.Consequently, the composition of the ring current may be modified by oxygen ions, the majority of which are initially in specific phase-space regions, which get accelerated and drift towards Earth.These findings are consistent with the results of previous studies on the role of substorms on the ring current dynamics, and have been verified here by an additional simulation.
For the effect of each current source to the Dst index during the event phases, we have concluded that one should, in principle, also account for the dependence of Dst on other effective sources apart from the ring current energy.Our computation of the Dst, in terms of the Dessler-Parker-Sckopke relation and test-particle results, indicates a measurable contribution from the near-Earth tail current of 30 % on the average, and yields a fair agreement with other estimations indicated in the literature (∼ 25 %).In the course of the event, the largest contribution of the tail current occurs during the growth phase, and persists for some time past its maximum.Thereafter, the effect of the tail currents gradually fades away, and the value of Dst is driven only by the ring current.Dst retains small values (related to meaningful disturbances) for long times after the event termination.A more accurate estimation of Dst may be achieved with the inclusion of the physics of loss mechanisms (collisions, cyclotron emission) and wave-particle interactions.
The present work may serve as the final link in a Sunto-Earth modeling chain of major solar eruptions, providing an estimation of the inner geospace response once the solar burst reaches Earth.In this frame, a comparison of our results with other available models (e.g., Fok et al., 2001;Jordanova et al., 2010), as well as with data from observations, may act as further validation.In a relevant work by Ganushkina et al. (2012), a benchmark of models was conducted, and the results showed that the computed ring current, for moderate and intense disturbances, depends on the field models.In our work, the benchmark of the T89 model has shown few differences within the simulated region in comparison to the later models, the facilitation of which may, however, increase the accuracy in the magnetic field.Also, the Kp index is used as a parameter to describe the geomagnetic field during the disturbances, which is partially in contrast to the global character of the specific index (see, for example, Rostoker, 1972).In order to assess our results, the correlation of Kp and Dst was followed and a connection was found with previous studies.

Figure 1 .
Figure 1.Geomagnetic field map on the GSM x −z plane, as calculated with the T89 model, for θ t = 11.5 • in a case where the magnetosphere is (a) quiet (Kp = 1) and (b) disturbed (Kp = 5).

Table 1 .
Input to the Tsyganenko models T89, T96 and TS05 for the benchmark case computations presented in Fig.2.Model Kp DstP dyn B y , B z S 1 , S 2 , . . .,

Figure 2 .
Figure 2. Computation of the relative deviation between the computations of the geomagnetic field, for zero tilt angle in a disturbed magnetosphere (Kp = 5 and Dst = −70 nT), using the models (a) T89 and T96 and (b) T89 and TS05.

Figure 5 .
Figure 5. Application of the model interchanging technique in testparticle computations for different values of the threshold radius: (a) r t, (b) E k vs. t, and(c) maximum relative error vs. R fm .

Figure 7 .
Figure 7. Three-dimensional plot of ion orbits which conclude outside the inner magnetosphere, with initial conditions same as in Fig. 6 apart from (a) E k (t 0 ) = 0.5 keV for H + and (b) E k (t 0 ) = 7 keV for O + .

Figure 8 .
Figure 8. Dynamic evolution of the (a) kinetic energy for the ion orbits analyzed in Figs. 6 and 7 and (b) pitch angle for the trapped O + orbit of Fig. 6.

Figure 10 .
Figure 10.Final kinetic energy of N ens = 1000 H + ions as a function of the initial (a) radial coordinate, (b) kinetic energy, and (c) final pitch angle as a function of initial kinetic energy, for varying initial conditions and all the remaining quantities same as in Fig. 9.

Figure 11 .
Figure 11.Analysis of Dst based on an ensemble of 2000 test ions, 1000 O + launched from the plasma sheet and 1000 H + in the ring current, during the event introduced in Sect.2: (a) dynamic evolution of Dst and of its ring/tail current contributions and (b) correlation of the maxima of |Dst| and Kp, computed by varying only Kp in (a), as compared to formerly derived results.

Table 2 .
Distribution to the ring current (N rc ), near-tail (N ntl ), precipitating (N pr ) and escaping (N esc ) O + and H + populations of test ions injected from the plasma sheet, with initial conditions for r in [18R E , 22R E ], E k in[2, 6]keV and θ in[20, 30]• .