The fate of O ions observed in the plasma mantle: particle tracing modelling and Cluster observations

Ion escape is of particular interest for studying the evolution of the atmosphere on geological time scales. Previously, using Cluster-CODIF data, we investigated the oxygen ion outflow from the plasma mantle for different solar wind conditions and geomagnetic activity. We found significant correlations between solar wind parameters, geomagnetic activity (Kp index) and the O outflow. From these studies, we suggested that O ions observed in the plasma mantle and cusp have enough energy and velocity to escape the magnetosphere and be lost into the solar wind or in the distant magnetotail. Thus, this study 5 aims to investigate where do the ions observed in the plasma mantle end up. In order to answer this question, we numerically calculate the trajectories of O ions using a tracing code to further test this assumption and determine the fate of the observed ions. Our code consists of a magnetic field model (Tsyganenko T96) and an ionospheric potential model (Weimer 2001) in which particles initiated in the plasma mantle region are launched and traced forward in time. We analysed 131 observations of plasma mantle events in Cluster data between 2001 and 2007, and for each event 200 O particles were launched with an initial 10 thermal and parallel bulk velocity corresponding to the velocities observed by Cluster. After the tracing, we found that 98% of the particles are lost into the solar wind or in the distant tail. Out of these 98%, 20% escape via the dayside magnetosphere.

their particular example the ions would reach the plasma sheet around -50 Re. As the near-Earth X-line is pushed towards Earth during disturbed conditions, these ions are expected to escape in to distant tail. Models and simulations have been extensively employed to investigate polar wind and ionospheric outflow. Schunk and Sojka (1989) simulated the polar wind behaviour using a combination of a low-altitude ionosphere-atmosphere and a highaltitude hydrodynamic model in a simulated region from 120 km to 9000 km. They discovered the complexity of the polar wind 60 density structures in different altitude ranges as well as for geomagnetic variations. Polar wind behaviour during one idealised geomagnetic storm has been investigated by Schunk and Sojka (1997), who updated their model to an altitude coverage of 90 km to 9000 km for latitudes higher than 50 • . They investigated the seasonal and solar cycle variations for four idealised geomagnetic storms (winter and summer solstices and solar minimum and maximum). They found that O + upflow increases over the polar cap during the storms, while O + is the dominant ion species at all polar latitudes. These results are similar to the ones by Barakat and Schunk (2006) who studied the generalised behaviour of the polar wind, also during an idealised geomagnetic storm using a macroscopic PIC (particle-in-cell) model. Their results agreed with satellite observations. At an intermediate lower altitude of 4000 km, Horwitz et al. (1994) determined the bulk velocity and temperature profiles of O + and H + in the polar wind using a semi-kinetic outflow model. They found that centrifugal forces increase the outflowing O + flux with 2 orders of magnitude when the convection electric field is enhanced from 0 mV/m to 100 mV/m. A similar result 70 has been shown by Abudayyeh et al. (2015), who used a Monte Carlo simulation based on the Tsyganenko T96 model and included the effects of the ambipolar electric field as well as gravitational and mirror forces. Additionally, Abudayyeh et al. (2015) observed higher bulk velocities and densities (H + and O + ) in the cusp than in the polar cap.
At an altitude range of 1.2 Re to 15.2 Re, Barghouthi et al. (2016) employed the same 1-D Monte Carlo model used by Abudayyeh et al. (2015) to investigate energetic H + and O + outflows along two trajectories (from the polar cap to the cusp) 75 and compared them with Cluster data. Considering the centrifugal acceleration, the ambipolar electric field and the waveparticle interaction, they concluded that the latter was the most important mechanism especially at higher altitudes (cusp).
Finally, a statistical model of the fate of energetic ions showed that these ions are highly dependent on the magnetic field configuration. Therefore, for quiet magnetic field, more ions escape directly through the magnetopause, whereas for active magnetic field, the ions are convected towards the tail and reach the distant tail at 50 Re (Ebihara et al., 2006). Ebihara et al. the magnetosphere with 50% in the distant tail. Despite all those interesting studies, the fate of ions observed in the plasma mantle has not yet been well defined. This study aims to clarify if O + ion outflows observed in the plasma mantle will escape the magnetosphere and be lost into the solar wind as suggested previously from observations Slapak and Nilsson, 2018;Schillings et al., 2018). For a more accurate estimate of the fate of ions, the starting point should be high-altitude, so that much of the transverse heating and centrifugal acceleration are already included. In order to answer this question, we 90 traced particles in a combination of the Tsyganenko T96 (Tsyganenko, 1995) and the Weimer 2001(Weimer, 2001 models.
About 25 000 O + particles were launched from the plasma mantle with initial parameters taken from Cluster observations. This model thus incorporates the effect of the mirror force on the launched ions, centrifugal acceleration and E × B drift. It does not include any further wave particle interaction than what the ions had experienced prior to the observation point.
This paper is organised as follows: Section 2 describes the instrumentation and data set we used, followed by the method 95 and a description of our code in Section 3. Section 4 and 5 present and discuss our results respectively. Finally, Section 6, summarises our study.  (Balogh et al., 2001) with a mode sample frequency of 22.4 Hz. In our study, we use the magnetic field averaged over the spacecraft spin period (4 s).
The solar wind data were retrieved from the OMNIWeb database. This database consists of data from several satellites at diverse positions around Earth. In our simulations (see Section 3.2), we utilise the solar wind dynamic pressure and velocity, the IMF By and Bz components in high-resolution (5 min) as well as the magnetic Dst index (1 h).

Plasma mantle dataset
Our dataset consists of plasma mantle and few cusp events observed by the Cluster spacecraft 4 between 2001 and 2007. In order to only retrieve plasma mantle data, we apply several constraints on the observational data. Firstly, the CODIF O + counts are contaminated when strong proton fluxes from the magnetosheath are recorded at the same energy level as the O + ions (Nilsson et al., 2006). These false O + counts usually originate from the magnetosheath and lead to an underestimate of the O + 115 velocity moment. The technique to remove these false counts is based on the E × B drift, because the drift is neither mass nor charge dependent. Consequently, using the kinetic energy equation, the cross talk signal is seen as an O + perpendicular bulk velocity that is 1/4 of the corresponding perpendicular proton velocity, and typically the O + density is higher than 2 cm −3 (for more details, we refer the reader to Nilsson et al. (2006)). We avoid these contaminated data (and therefore magnetosheath data) using the method described by Nilsson et al. (2006). However, we slightly changed the threshold defined by Nilsson et al. 120 (2006) because over the years the quality of the Cluster data devolved and so did the threshold. The new thresholds are given To pick out only plasma mantle observations we implement different conditions for the high latitude regions. In these regions, the plasma beta β (O + and H + ) is typically higher than 0.05, whereas it is lower than 0.05 in the polar cap (Liao et al., 2010(Liao et al., , 2015Haaland et al., 2017). We use a threshold value of β > 0.1 for high latitude regions. Additionally, the perpendicular 125 temperature of the protons should be lower than 1750 eV in order to distinguish plasma sheet from plasma mantle data (Nilsson et al., 2006;Kistler et al., 2006;Slapak et al., 2017). As partly mentioned above, the O + and H + densities are restricted to n(H + ) > 10 −3 cm −3 and 10 −3 cm −3 < n(O + ) < 2 cm −3 to keep only reliable velocity estimates. In order to study the fate of ions, we take O + data with an outward flow (v > 0 or v < 0 in the southern and northern hemispheres respectively).
Finally, we use a spatial coverage restriction to remove the inner magnetosphere, which is defined by -5 Re < X GSM < 8 Re Re (see also Slapak et al. (2017)). Major geomagnetic storms data are removed to exclude other magnetospheric regions than the plasma mantle . Finally, we also excluded events associated with potential shock arrivals at Earth to avoid strong perturbations in the magnetosphere.
When all the above conditions are met, we define one event by 60 data points or more in a row. Between 2001 and 2007, our automatic routine detected 131 events that met the region criteria and the model restrictions (see Section 3.1).

Methodology
The section aims at briefly describing how our model works and its inputs and outputs.

Particle tracing simulations
We use a test particle simulation code (Gunell et al., 2019) to compute ion trajectories in the magnetic fields given by the Tsyganenko T96 model (Tsyganenko, 1995) and electric fields derived from the ionospheric potential given by the Weimer 140 2001 model (Weimer, 2001). The electric field is defined on a grid, and during the test particle trajectory calculation the electric field at the particle position is found by interpolation. Before the trajectory calculation starts, we define the electrostatic potential V on a three-dimensional grid with a cell size of 1200 km × 1200 km × 1200 km in the region −60 < X < 10 Re, |Y | < 19 Re and |Z| < 19 Re by tracing the magnetic field lines from each cell down to the ionosphere, where we retrieve the potential from the Weimer model. The electric field is then found from the relationship E = −∇V .
145 Figure 1 illustrates the magnetic field lines in light grey and the electric field grid in brown. This illustration represents the magnetosphere based on T96 and the grid for the interpolation of the electric field. Due to the limits of the Tsyganenko model, the electric field grid goes to -60 Re in the tail and 10 Re in the dayside (Tsyganenko, 1995). In the nort-south, Z, and dawn-dusk, Y , directions, the limit of the grid are at ±19 Re. Note that the illustration is not to scale. To launch a particle, its initial ion 3-D velocity is calculated using the following equation where v 0 is the bulk velocity defined by v and v E×B , the parallel and drift velocities respectively, v th is the thermal velocity of the O + ion (see Section 3.2 for more details), B the magnetic field and E the electric field. Then, the electric field at the position of the O + ion is formed by interpolation of the field saved on the grid. Finally, from the interpolated electric field, a new velocity is calculated with the Boris algorithm (e.g. Birdsall and Langdon, 1991).

155
The last step is repeated as far as the limitations of the code allows it. The tracing uses a time step based on the gyro-period, so that our time step is dt = 2π/20ω c whereas the maximum number of iterations is limited to 10000. In 99,93% of the cases, the particle stops because they reach the limit of our model (electric field grid), whereas for the 0.07% remaining the maximum number of the iterations have been completed.
The grid that defined the limit of our model is sufficient for our study because it includes the whole magnetosphere around 160 Earth (magnetopause is defined around X = 10 Re in the dayside and about |Y | = 13 Re and |Z| = 13 Re for X = 0 Re). Note that small regions near the subsolar magnetopause are not included, however, they are not relevant for our study. Further in the magnetotail, the magnetopause expands into a "cone shape" in the Y and Z direction and beyond 200 Re in the X direction. Our grid stops at X = −60 Re in the tail due to Tsyganenko model limits, moreover most of the particles reaching that distance will most likely be lost (see Sections 4 and 5 for more details). Concerning constraints, the Weimer model imposes no constraints on 165 solar wind parameters, while Tsyganenko T96 does. Therefore, when an observation meet the criteria described above it also has to match with Tsyganenko T96 constraints, which concerns the Dst index (-100 nT < Dst < 20 nT), the dynamic pressure (0.5 nPa < P dyn < 10 nPa) and IMF B z and B y (−10 nT < B z , B y < 10 nT).

Inputs and outputs of the model
The inputs of the models are (a) solar wind parameters as required by the Tsyganenko-and Weimer models, (b) the positions, the thermal and bulk parallel velocities, v and v th respectively, based on Cluster observations. The solar wind parameters (see    10 Re). We defined the escaping limit by the geocentric distance of the final positions R = X 2 fin + Y 2 fin + Z 2 fin that equals Furthermore, in Fig. 4c, we determined the minimum distance in the X direction for each trajectory as the minimum value found in the X direction of each ion trajectory (given in Re). As an example, for a trajectory of 2000 steps, the minimum X distance is the minimum value (in the X direction) within the 2000 steps, which is not necessarily the final position of the ion.
In Fig. 4c, each data point represents the minimum X distance of one trajectory. The average minimum X distance is around 225 -10 Re, which corresponds to the plasma mantle region if |Z| > 5 Re (see also on Fig. 5). This parameter is important because some particles that interact with the plasma sheet in the distant tail might come back close to Earth. However, such cases are rare because for a total of 849 trajectories having a X minimum distance beyond -50 Re only 23 trajectories finish their route close to Earth (R < 10 Re). The 826 remaining are roughly equally spread between 10 Re and 64 Re.  is observed around Earth (-3 Re < X < 3 Re) at R cyl = 23 Re. In contrast, the lower scaled flux is observed below R cyl = 10 Re and between 15 Re < R cyl < 20 Re for X lower than -20 Re.

Discussion
In our 131 events based on Cluster-CODIF observations, the parallel and thermal components of the velocities during the 245 events are taken as inputs to our forward tracing model (see Section 3.2 and Fig. 2). From these observations, we found that O + ions observed in the plasma mantle have a parallel velocity being comparable to, or lower than the thermal velocity as expected. More precisely, the ratio between the velocity components (|v |/v th ) is 0.65 ± 0.44. Since the bulk perpendicular velocity is much smaller than the parallel bulk velocity in the plasma mantle (Vaith et al., 2004), the O + ions in the mantle can be characterised by the parallel bulk velocity and the thermal velocity. The bulk perpendicular velocity is relatively small, 250 but would be important for their trajectories, in particular, storm time, and/or southward IMF (Haaland et al., 2008). Indeed, in the 13% of our events with higher convection velocity, more than half have a corresponding Dst index below -20 nT and a southward IMF.
We do not find any strong correlation between geomagnetic activity (Dst) and the final positions. For the IMF direction, we identify 48% of the events are associated with northward IMF and the final positions of these ions to be mainly spread between 255 found that for southward IMF the cold ion outflow is convected toward the plasma sheet due to strong convection, whereas for IMF directed northward convection is stagnant, so that cold ion outflow reach the far tail. Slapak et al. (2012) suggested three main routes for ion outflow; (1) cold ion that will end up mainly in the plasma sheet (Mouikis et al., 2010;Haaland et al., 2012;Liao et al., 2015), (2) energised ions from the cusp to the plasma mantle (Liao et al., 2010;Slapak et al., 2017;Schillings et al., 2019), (3) energised ions from to cusp going directly to the magnetosheath . Slapak et al. (2017); Slapak and Nilsson (2018);Schillings et al. (2019) suggested that ions observed in 265 the plasma mantle have sufficient energy and velocity to escape in the distant tail. However, our results show that very few ions reach the distant tail but instead escape directly through the magnetopause after a few minutes (∼ 17 min). These O + ions have short or middle length trajectories in our model (less than 2000 steps, see also Section 4) and represent 95% of our sample. Most (99.6%) of these O + ions reach a point where the tracing is stopped at a geocentric distance higher than 10 Re and escape the magnetosphere. For ions with trajectories longer than 2000 steps (5% of the total trajectories), 35% is 270 earthward flow due to its interaction with the plasma sheet. Most of these ions do not return to the ionosphere. Some will instead experience charge exchange, become neutral and be lost from the magnetosphere. This assumption is supported by Ebihara et al. (2006), who modelled O + trajectories and introduce a charge exchange process in their model. They estimated that 2% of the total outflow became neutral due to charge exchange with the hydrogen geocorona. Other particles will drift to the magnetopause (magnetopause shadowing) and be lost. We note that ion precipitation recorded by the DMSP spacecraft 275 (Newell et al., 2007) indicates a total precipitation of ions (H + and O + ) of the order 10 24 s −1 , which is most of the time dominated by cusp precipitation, not return flow precipitation. This is even less than the return flow estimated by Slapak and Nilsson (2018), indicating that most return flow indeed does not precipitate to the ionosphere. However, we do not study the fate of this earthward ion flow and therefore they are not considered as escaping ions in this study.
Under quiet magnetospheric conditions (Dst ≥ -20 nT), it was found that 3% of the final positions of the trajectories is 280 within a geocentric distance of 10 Re (return flow), whereas during disturbed conditions we observe only 1.6% return flow.
This result agrees with Ebihara et al. (2006), who found that under quiet time 4% to 7% of the outflowing ions return to Earth.
Under disturbed conditions, the authors estimated a smaller return of 0.6% to 0.8%.
Finally, since O + ions are launched from the plasma mantle, the particles observed by CODIF already went through transverse heating and centrifugal acceleration. Thus this model includes most of the energisation and acceleration compared to 285 other models. Moreover, the model does not include wave-particles interaction after the oxygen ion has been launched.

Summary and conclusions
Based on previous suggestions that O + ions from the plasma mantle are escaping Slapak and Nilsson, 2018;  The initial velocities and positions are determined by Cluster observations and are used as inputs for the forward tracing. Our results are summarised in the following points: 1. The O + ions observed in the plasma mantle have an initial parallel velocity lower than the thermal velocity as expected.

295
Thus, the thermal velocity dominates from the start, and through high perpendicular temperatures, the mirror force will increase the parallel velocity further downstream of the observation point.
2. 98% of the final positions (out of 26200) are located beyond a geocentric distance of 10 Re. These particles escape and are lost into the solar wind. 20% of the ions escape directly through the high-latitude dayside magnetopause.
3. 2% of the total trajectories lead back towards earth, i.e. they constitute return flow. Some of these O + ions have interacted 300 with the plasma sheet in the distant tail and eventually end up between the Earth and a geocentric distance of 10 Re.
4. Under disturbed magnetospheric conditions (Dst < -20 nT), we observe 1.6 % return flow, whereas during quiet time the return flow increases to 3%. 5. We do not find any correlation between the IMF direction, the geomagnetic disturbances and the final positions of O + in our tracing model. However, the ions ending up close to the Earth (geocentric distance smaller than 10 Re) are for 67% 305 of the time associated with southward IMF.
Code and data availability. The Cluster data can be retrieved from the Cluster Science Archive <https://csa.esac.esa.int/csa-web/>. Competing interests. The authors declare that they have no competing interests.