Multiscale hybrid modeling of the impact response of the Earth’s magnetotail to ionospheric O + outflow

Ionospheric outflow plays an important role in coupling the ionosphere with the solar wind-magnetosphere system. Previous multi-fluid MHD studies explored the global influence of oxygen ions of ionospheric origin (O +) on magnetospheric dynamics. A detailed exploration of the interaction of ionospheric ions with the magnetotail requires kinetic treatment for ions. We perform a self-consistent investigation of these processes with a three-dimensional space-time adaptive hybrid code, HYPERS, powered by an intelligent Event-driven Multi-Agent Planning System (EMAPS). By comparing simulations with and without outflow we conclude that oxygen ions, flowing from the ionosphere through the lobes into the tail or directly entering the inner magnetosphere, are able to significantly modify the magnetotail configuration and induce X-points and current sheet structures at locations where magnetic reconnection does not occur in a simulation without outflow, potentially very close to the Earth. This finding may have implications for interpreting substorms and magnetotail reconnection events observed for southward magnetic field simultaneously with significant contents of oxygen ions of ionospheric origin.


Introduction
The Earth's magnetosphere is sustained by two major plasma sources: the solar wind and the ionosphere. While the solar wind predominantly contains hydrogen (H+) and helium ions (He++), the ionospheric sources include hydrogen, helium (He+) and oxygen (O+) ions. In particular, O+ from the ionosphere is an important source for magnetospheric plasma that plays a significant role in the magnetospheric dynamics (Kistler, 2017). Since O+ originates primarily from the ionosphere, oxygen is commonly used for studying impacts of ionospheric outflows on the magnetosphere, where O+ ions are regularly observed in the lobes, the plasma sheet and the ring current (Kistler, 2017). The O+ content in all cases is higher during geomagnetically active times, with outflow rate increasing with the general planetary activity level (Kp index) (Yau and André, 1997;Cully et al., 2003a).
Ionospheric ions can be transported into the magnetosphere through two most important pathways, namely, from the nightside auroral and dayside cusp regions, which form the boundaries between open and closed magnetic field lines (e.g., Liao et al., 2010;Kistler, 2017;Toledo-Redondo et al., 2021).
Ions flowing out from the cusp ("cleft ion fountain") are heated in this region and either carried toward the magnetopause region or convected over the polar cap at high altitudes into the tail lobes along the open magnetic field lines (Liao et al., 2015). As a result, the lobe field lines above and below the plasma sheet are filled with ionospheric ions that eventually end up in the plasma sheet or become lost in the distant tail. Heavy ions from the ionosphere, accelerated in the cusp/cleft, have been identified as a direct source for the hot plasma in the plasma sheet (Liao et al., 2010). This dayside cusp outflow is a relatively consistent feature, varying with the geomagnetic conditions and the F10.7 index, a proxy of the ionizing UV radiation from the Sun . For instance, tailward streaming O+ beams are commonly observed inside the lobes even during non-storm times . The O+ outflow increases significantly during the storm main phase and drops slightly during the recovery phase, when the convection is enhanced in the magnetotail lobes (Liao et al., 2010). Cully et al. (2003b) found that under times of enhanced convection, the "in flight" ionospheric ions preferentially supply the tail region where substorm activity is thought to originate, and that a reduction in the lobe convection, as caused, for example, by a northward turning of the interplanetary magnetic field (IMF), immediately and dramatically reduces this supply, which may have consequences for substorm triggering.
The other important transport path for ionospheric ions to the magnetotail originates in the nightside auroral region. This outflow is prone to drastic intensifications during times of high geomagnetic activity. This path transports O+ directly into the near-Earth tail plasma sheet without passing through the lobe region (Liao et al., 2010). In some cases, O+ was even found to dominate the total number density of the plasma in the tail. In particular, Kistler et al. (2016) found that the plasma sheet densities closer to the Earth, at 6 − 7R E (R E is the Earth radius), were characterized with the O+/H+ ratios about a factor of 10 higher than in the 15-19R E region. Overall, in terms of mass density, the near-Earth plasma sheet is mostly dominated by the ionospheric source, owing to these O+ contributions (Toledo-Redondo et al., 2021).
At times of high solar and magnetic activity, the overall O+ fluence at the auroral latitudes increases exponentially with Kp up to ∼3 × 10 26 ions/s, with O+ ion energies reaching up to 17 KeV (Yau and André, 1997). The nightside polar cap region, where magnetic field lines are open and connected to the IMF, may also serve as an extra source of energetic O+ ions. However, the fluences of these ions ∼10 25 ions/s (Yau and André, 1997) are typically lower by an order of magnitude than those originating from the cusp and auroral region. Moreover, ions upflowing from the cusp and auroral oval typically reach higher energies than ions originating directly from the nightside polar cap. Overall outflow fluences of H+ and O+ from the high-latitude region can reach above 10 26 ions/s (Toledo-Redondo et al., 2021). These rates vary by at least one order of magnitude, increasing with stronger solar EUV and geophysical activity, especially for heavier O+ ions (Yau and André, 1997;Cully et al., 2003a).
To summarize, the outflowing ionospheric O+ is transported into the magnetosphere through the combined field-aligned (parallel) motion and magnetospheric convection. For southward IMF, magnetopause reconnection results in directing both dayside cusp and nightside auroral outflows toward the magnetotail. The cusp originating ions eventually pass over the nightside polar cap and stream into the lobes on open field lines, whereas the auroral ions may directly be injected into the magnetotail on closed field lines. In both cases, ions flowing from the high-latitude regions are also energized along their trajectories. These trajectories depend on the magnetic latitude (ML) and magnetic local time (MLT) at their initial escape surface (a sphere of radius ≥1 R E ). The points of entry of these ions into the magnetotail also depend on their initial energies and pitch angles (i.e., angles between the ion velocity and the magnetic field direction) at their local injection points.
The issue of ionospheric outflow acting as a possible substorm trigger was first raised by Baker et al. (1982), who argued that increased contents of large-gyroradius O+ ions during active times might help destabilize linear ion tearing mode in the thin plasma sheet because the growth rate is proportional to the ion gyroradius. Baker et al. (1982) also reported the asymmetries observed in substorm onset location, with onsets occurring preferentially on the duskside. Baker et al. (1982) related this finding to the observed asymmetries in the O+ content in the tail. Later observations, however, found no dawn-dusk asymmetry in the O+/H+ ratio at 7-8R E , which suggests that the overall asymmetry in substorm onset location may also have some other explanation (Kistler, 2017).
Examining the ability of the ionospheric O+ ions to significantly affect the local structure and stability of the magnetotail is critical for enabling predictions of potential outflow driven substorms and reconnection regions in the tail at ≥20R E . The effects of ionospheric outflow on sawtooth events in the magnetotail and the properties of associated substorms were investigated by Ouellette et al. (2013); Garcia et al. (2010); Brambles et al. (2010), Brambles et al. (2011. Recent in situ observations by Angelopoulos et al. (2020) have also indicated that tail reconnection may be powered by intense storms at close geocentric distances ∼10R E , where magnetic reconnection is normally expected to be suppressed by the Earth's strong dipole field. These observations also revealed a significant presence of O+ ions (up to ∼50%) in the tail, in addition to protons. What may lead to reconnection so close to Earth, however, remains an open question, together with a more general question of whether O+ may induce the triggering of reconnection in the magnetotail despite that the observations do not seem to indicate that O+ directly increases the instability of the current sheet to reconnection (Kistler, 2017). This paper derives its interest from ample evidence suggesting that many substorms may be triggered not only by changes in the solar wind but also by the ability of the ionosphere to populate the tail lobes and the near-Earth region with heavy O+ ions, producing the associated global impact on the magnetotail, as emphasized above. The outflow impact may become especially important when the IMF turns southward, even though enhanced ion outflow delivery to the plasma sheet increases gradually due to a delay in the transport time of O+ ions from their source regions, with ions with greater parallel velocities reaching further down the tail (Moore and Horwitz, 2007;Kistler et al., 2016). In this paper, for the first time, we perform three-dimensional (3D) global hybrid simulations of the Earth's magnetosphere with O+ outflows. For this purpose, we use a multiscale hybrid code, HYPERS, which treats all ions as kinetic particles and the electrons as a massless quasineutral fluid, as described by Omelchenko and Karimabadi (2012a), Omelchenko and Karimabadi, (2023); Omelchenko et al. (2021a). The purpose of these ion kinetic ("beyond MHD") simulations is to demonstrate the ability of ionospheric O+ to enter and significantly modify the Earth's magnetotail. Although the dayside magnetosphere in these simulations is not accurately resolved, it still provides a realistic global configuration for injecting O+ ions.
In Section 2, we briefly review results from the previous simulations of O+ contributions to the tail dynamics, conducted with multi-fluid MHD and local particle-in-cell (PIC) codes. In Section 3, we describe the multiscale HYPERS-Global model, which has recently become available for public use at the NASA CCMC. In Section 4, we compare results from simulations performed for two different scenarios that model the dayside cusp ("lobe") and "auroral" O+ outflows. We also perform a baseline simulation without outflow. Finally, in Section 5, we summarize the most important points of our study and discuss implications of including ionospheric sources of plasma in the future global multiscale models of the Earth's magnetosphere.

Past simulations of oxygen impact on magnetotail
Global modeling of ionospheric sources of plasma and their impact on the magnetosphere has been performed with multifluid magnetohydrodynamics (MHD) codes (e.g., Winglee, 2000;Garcia et al., 2010;Brambles et al., 2010;Brambles et al., 2011;Yu and Ridley, 2013a;Yu and Ridley, 2013b) (for more references see a review by Toledo-Redondo et al. (2021)). In particular, Brambles et al. (2011) used a multi-fluid MHD code to study the effects of auroral outflow on magnetotail dynamics. By combining simulations with quantitative analysis, they showed that the amount of hemispheric outflow flux rate may control the mode of magnetospheric convection. A polar wind outflow model (PWOM) was developed by Glocer et al. (2009) and applied to study the global magnetospheric dynamics and nightside reconnection during a magnetic storm (Ilie et al., 2015). More recently, Varney et al., 2016a,Varney et al., 2016b developed a twoway coupled ionosphere/polar wind model (IPWM) for different ionospheric outflow species.
Local full-electromagnetic Particle-in-Cell (PIC) simulations (e.g., Karimabadi et al., 2011;Liu et al., 2014;Tenfjord et al., 2019) and theoretical models  have suggested that O+ ions may significantly alter the current sheet structure and properties of magnetotail reconnection. Local multi-fluid Hall-MHD simulations (Shay and Swisdak, 2004) also demonstrated that O+ slows the outflow speed and reconnection rate. Shay and Swisdak (2004), however, argue that the increased geomagnetic activity during storms may "overpower" the local reduction in the reconnection rate due to O+. At the same time, Hesse and Birn (2004) found in their full PIC simulations that the reduction in reconnection rate from O+ is not as large as expected from the simple fluid picture.
MHD simulations with test (not self-consistently modeled) outflow particles have also been performed (e.g., Peroomian et al., 2006;Moore and Horwitz, 2007). The importance of adapting MHD-particle simulations to include self-generated fields and properly assess their effects on ion motion (especially for heavy ions) was emphasized by Moore and Horwitz (2007); Peroomian et al. (2006) after they found out that global test-particle simulations of O+ ions in dynamically evolving MHD fields produced results consistent with observations of the increasing contribution of heavy ions in the near-Earth plasma sheet. MHD simulations, however, give only averaged (Maxwellian) pictures of outflow and magnetotail plasmas since they do not include important kinetic effects observed in the magnetosphere. Peroomian et al. (2006) launched test O+ particles with an equivalent fluence of up to ∼10 26 ions/s from the cusp region in both hemispheres from just below the inner simulation boundary with a radius of 2.5R E . They argued that the dayside cusp outflow is 10-100 larger than the corresponding fluence of O+ ions from the nightside auroral zone, and therefore there was no need to include the latter in their simulation. This assumption, however, contradicts observations by Yau and André (1997), who measured net auroral O+ fluences up to 2 × 10 26 ions/s for large Kp values. An important implication of choosing the cusp source of O+ by Peroomian et al. (2006) was the generation of O+ ion streams over the nightside polar cap into the lobes, as discussed earlier. In addition, the passage of a solar wind shock and the southward turning of IMF in their simulation caused the O+ ions, injected into the lobes, to accelerate and more vigorously penetrate the plasma sheet.
Yu and Ridley (2013a), Yu and Ridley (2013b) performed multifluid MHD simulations to examine whether the ionospheric source of the O+ at the inner boundary with a radius of 2.5R E may produce a significant impact on the magnetotail. In their simulation of the cusp outflow, they chose the ionospheric plasma to flow out from a wedge-shape region at the inner boundary from 65 to 85°magnetic latitude (approximately 74.5-86.8°at the ionospheric altitude when mapped down along dipole magnetic field) and from 11 to 13 magnetic local time (MLT). The simulation of the nightside auroral outflow used a wedge-shape source region at the inner boundary that extended on the nightside from 23 to 1 MLT across the midnight, and from 55 to 75°magnetic latitude (68.7-80.6°at the ionosphere altitude). Both simulations used hemispheric O+ fluences up to 2 × 10 26 O+ ions/sec. Similarly, Garcia et al. (2010) used O+ fluences up to 1.4 × 10 26 ions/s to simulate O+ ions outflowing from an isolated patch in the nightside auroral zone. Yu and Ridley (2013b) note that the purpose of specifying such uniform and large (but still within a statistical range of observations) outflow source in their simulations aims at emphasizing the impact of O+ outflow on the magnetotail. We pursue the same goal in our hybrid simulations. Yu and Ridley (2013a) also observe that their nightside auroral outflow simulation, where O+ ions flow directly into the inner magnetosphere and plasma sheet earthward of the normal reconnection line, shows little change in the magnetotail through the entire simulation time, which implies that the mass loading in the tail produces an insignificant impact on initializing substorm onset, even though magnetic flux is constantly advected into the reconnection region.
In contrast with that simulation, the cusp outflow simulation by Yu and Ridley (2013b) indicated that the O+ ions in the tail help trigger tail reconfiguration by compressing the neutral sheet, accumulating magnetic energy in the tail, and then disrupting the tail, which was followed by another disruption. This result is consistent with results from multi-fluid simulations by Wiltberger et al. (2010), who found a second substorm in their O+ outflow model which did not occur in the no-outflow case. Yu and Ridley (2013a) finally conclude that their cusp outflow lands near the tail reconnection region and dramatically disturbs the magnetotail, while the nightside auroral outflow plasma mostly populates the near-Earth plasma sheet, "not encouraging" the release of the accumulated magnetic energy. We find these details particularly important for comparison with our hybrid simulations (see Section 4).

Hybrid equations
We model the Earth's magnetosphere with ionospheric O+ sources using a hybrid approach to approximating the full plasma-Maxwell system of equations (see, e.g., Winske et al., 2003). This approach neglects the displacement current in Maxwell's equations and treats ions as full-orbit kinetic species (macro-particles in our model), moving in self-consistent electric and magnetic fields. The plasma electrons are modeled as a massless quasi-neutral fluid with an adiabatic equation for scalar pressure. This approximation is described by the following set of equations, written in CGS units: p e = n e T e ∼ n γ e .
In Eqs. 1-5, c is the speed of light; e, T e , n e , J e are the electron absolute charge, temperature, number density, and current density, respectively; p e is the electron pressure, described by Eq. 5 with γ = 5/3; ρ i , J i are the total ion charge and current density, computed on the computational mesh by the PIC method; E is the electric field, found from the generalized Ohm's Eq 1. B, B ext , B tot are the "selfgenerated" (B| t=0 = 0), "external" (time-steady), and total magnetic fields, respectively.
This plasma approximation is implemented in a multiscale hybrid code, HYPERS (Omelchenko and Karimabadi, 2012a;Omelchenko et al., 2021a). The generalized Ohm's law Eq. 1 includes a resistive term, where the applied plasma resistivity, η is chosen to be a function of plasma density and current (Omelchenko et al., 2021a). Applied resistivity (not to confuse with numerical resistivity) in hybrid simulations serves two purposes: 1) imitate electron inertia effects in reconnection regions, 2) enable fast magnetic field propagation (diffusion) in low-density and plasma-free regions.
The latter is needed because the displacement current is excluded from the Ampere law Eq. 2 in the standard (Darwin) hybrid approximation. Failure to properly treat low-density regions in a hybrid simulation may result in producing non-physical effects (Poppe, 2019;Omelchenko et al., 2021b). In practice, the first two terms in the RHS of Eq. 1 are forced to vanish at "vacuum" cells, with n e ≤ n min , n min ≃ 0.05-0.1n 0 , where n 0 is the upstream solar wind electron number density.
The hybrid-PIC model, described by Eqs. 1-5 and coupled with the PIC method for ion macro-particles, provides important "beyond MHD" simulation capabilities, such as the Hall physics and ion kinetic/cyclotron motion effects. This model supports ion cyclotron, whistler, and kinetic-Alfvén wave modes that play a prominant role in regulating plasma turbulence in the Earth's magnetosphere and impacting its global behavior. In addition, representing ion components with Lagrangian macro-particles avoids numerical diffusion that causes meshdependent reconnection and spatial spread of ion sources in fluid simulations (e.g., Yu and Ridley, 2013a). The HYPERS model of the magnetosphere enables using arbitrary ion species for both solar wind and ionospheric outflow. Reconnection onset and dynamics are controlled by the ion kinetics, Hall physics and parameter-dependent resistivity. Ultimately, kinetic effects should be incorporated into global models of the Earth's magnetosphere, in order to better understand the influence of turbulent plasma dynamics on reconnection.
For presentation purposes, we use GSM coordinates everywhere in this paper. The solar wind plasma in this study consists of protons and electrons. We choose solar wind parameters close to those used in multi-fluid MHD studies by Yu and Ridley (2013a), Yu and Ridley (2013b). The solar wind is continuously injected from the leftmost x-boundary in the negative x-direction with the number density, n SW = 3 cm −3 and speed, v SW = 350 km/s. The initial electron and proton temperatures are chosen to be equal to 10 eV.
To allow for a magnetic flux barrier to build in the magnetosheath upstream of the dayside magnetopause (Samsonov et al., 2017), the initial IMF field, B IMF = 5 nT is chosen to be directed northward everywhere. Then, at the simulation time, t ≃ 10Ω −1 cp , where Ω cp is the proton cyclotron frequency computed with respect to the initial IMF strength, the IMF is switched southward at the injection boundary in order to enable dayside magnetopause reconnection and promote tailward convection of lobe magnetic flux, which is also consistent with favorable physical conditions for enhanced ionospheric outflow (Winglee, 2000). The O+ outflow is initiated (see Section 4) at t ≃ 44Ω −1 cp , just prior to the arrival of the southward IMF at the subsolar bow shock. All simulations were conducted for the time period ≃ 600Ω −1 cp , which is long enough to observe the outflow impact on the magnetotail. The Earth's dipole is assumed to have no tilt in this study.

Spatial scales
Normally, hybrid simulations seek to resolve spatial scales of the order of the ion inertial length, d i and ion cyclotron radius, r ci . Assuming the Earth's radius, R E ∼ 60d i , to encompass the whole magnetosphere in a global hybrid simulation, the computational domain at minimum should include the magnetopause with a typical standoff distance, R MP ≃ 10R E ∼ 600d i and the magnetotail that stretches far from the Earth to distances ∼100R E ∼ 6000d i . The CPU-intense nature of hybrid simulations currently prevents running simulations on these large scales in 3D. Therefore, in this study, we scale down the nominal magnetopause standoff distance by reducing the magnetic dipole's strength, so that R MP ≃ 100d i . This numerical dipole, however, is strong enough to resolve the most essential local kinetic features of the magnetosphere (Omidi et al., 2004). Furthermore, as was demonstrated by Tóth et al. (2017), the global magnetospheric solution is not strongly sensitive to this scaling factor, with the ion inertial length dynamics of the numerical model occurring in a self-similar manner with respect to the original plasma system.
In practice, the need to accurately account for inflow and outflow boundary conditions in the presence of the Earth's magnetic dipole requires multiplying these dimensions by a factor of 2-3. Approximating such large 3D domains with uniform meshes with cell sizes of the order of ∼1d i is prohibitively expensive in a global hybrid model. To alleviate this problem, we cover the computational domain with a logically mapped, "stretched" Cartesian mesh (Figure 1). This mesh maintains a uniform spatial resolution of Δx = Δy = Δz = 2d i in the central part ("patch") of the domain, which partially covers the lobes and near-Earth tail that are of most interest to this study. This resolution is sufficient to capture the ion kinetic effects in this region. The mesh stretching pushes physical domain boundaries away from the uniform patch, making the Earth's dipole field vanish at the inflow/outflow (x) and lateral (y and z) boundaries. This helps implement physically consistent local boundary conditions for plasma flow (Omelchenko et al., 2021a).
For all simulations in this study, we use a logical mesh of 450 × 300 × 300 cells (see Figure 1). Initially, the uniform patch is initialized with 25 solar wind macro-particles per cell. Outside the patch in the x-direction, the number of particles at each cell is increased according to the physical cell size in this direction. Along the lateral (y and z) axes, the plasma density is initialized uniform by scaling individual weights of ion macro-particles in accordance with physical cell sizes in these directions (Omelchenko et al., 2021a).
During a simulation, each particle found in the uniform patch is dynamically split (just once) into two new particles for better resolution. To maintain the initially specified solar wind plasma, protons are continuously injected into the computational domain at the positive x-boundary with a shifted Maxwellian velocity distribution. The near-Earth environment ("inner magnetosphere") in this study is not modeled in detail, being represented by a particleabsorbing spherical obstacle with a radius of R obs = 65d i = 0.65R MP . As already mentioned, this numerical obstacle size, when scaled up to realistic magnetospheric dimensions, corresponds to the physical inner domain radius, R obs = 6.5R E . Finally, to account for finite ionospheric conductivity effects, we allow the Earth's magnetic field lines to diffuse across the surface of the inner boundary by assigning a constant resistivity value to the outer shell of the obstacle.

Temporal scales
Explicit hybrid solvers need to stably integrate electromagnetic fields on the kinetic-Alfvénic and whistler timescales, in addition to the Alfvénic ("MHD") timescales. Because of the inherent particle noise, hybrid-PIC simulations typically generate whistler mode oscillations with large wave numbers, k ∼ 1/Δ and high frequencies, ω ∼ k 2 ∼ 1/Δ 2 (Δ is the local mesh size). If not accurately resolved or resistively damped, these spurious oscillations explosively grow in time and quickly abort the simulation. Therefore, predicting FIGURE 2 H+ number density and absolute value of net current density, |J| in three cross-sections in the baseline simulation (no outflow) at Ω cp t = 461: x-z ( y = 1d i ), x-y (z = 1d i ), and y-z (x = −199d i ).
optimum timesteps for computing the multiscale electromagnetic fields in hybrid simulations becomes often difficult in practice.
The problem of choosing global timesteps for the fields is further compounded by the problem of choosing proper timesteps for kinetic ions. For accuracy and stability, full-orbit particle solvers ("pushers") typically require small timesteps, Δt p , so that ΩΔt p ≪ 1, where Ω is the local particle (ion) cyclotron frequency. In the Earth's magnetosphere, where the local strength of magnetic field greatly varies in space, the need to satisfy this cyclotron constraint for all particles typically results in forcing either small timesteps for all particles or producing inaccurate trajectories for particles in strong magnetic fields. For instance, choosing particle timesteps, Ω 0 Δt p ≃ 0.05, where Ω 0 is the ion cyclotron frequency computed with respect B IMF , may result in producing large errors in the ion cyclotron orbits and drifts for particles located closer to the inner obstacle (e.g., in the cusp or even some parts of the magnetosheath), where Ω ≥ 10Ω 0 . A similar argument about timesteps can be made for non-Maxwellian particles, energized in the self-consistent electromagnetic fields in hybrid simulations. Therefore, in addition to using uniform timesteps for the fields, using the same timesteps for all particles inevitably creates another numerical bottleneck in the hybrid model.
The HYPERS code effectively overcomes this multiscale challenge by implementing an "intelligent" approach to time integration: Event-driven Multi-Agent Planning System (EMAPS), introduced by Omelchenko and Karimabadi (2006a); Omelchenko and Karimabadi (2006b); Omelchenko and Karimabadi (2007); Omelchenko and Karimabadi (2023). The time integration of cell-based fields and Lagrangian particles on a stretched mesh is performed by EMAPS in a "game of life" fashion, so that each computational quantity (cell or particle) is dynamically assigned an FIGURE 3 H+ number and absolute value of net current density, |J| in three cross-sections in the first lobe outflow simulation at Ω cp t = 461: x-z ( y = 1d i ), x-y (z = 1d i ), and y-z (x = −199d i ).
individual time increment, with cells being synchronized with other local cells on an "as needed" basis. In this process, EMAPS plays a role of a "simulation time operating system, " which self-adaptively predicts and corrects local time increments. This novel multiscale paradigm has enabled parallel, CPU-efficient 3D hybrid simulations of the Earth's magnetosphere Ng et al., 2021;Omelchenko et al., 2021a;Omelchenko et al., 2021b) and laboratory plasmas (Omelchenko and Karimabadi, 2012b;Omelchenko and Karimabadi, 2023).
Being restricted by the fast ion cyclotron time scales, 3D hybrid simulations of the Earth's magnetosphere typically run for simulation periods, t ∼ 100 − 500 Ω −1 0 . Formally, these periods span short magnetospheric times, t < 20 min, compared to hours in global MHD simulations. For convenience of comparing results to certain observations, however, some modelers prefer to multiply these simulation times by typical hybrid model scaling factors, such as the ratio of real to numerical magnetopause distance, or the ratio of simulation to real Alfvén speed (e.g., Lin et al., 2022). This upscaling procedure serves the purpose of presenting hybrid simulation results in "magnetospheric" time. We note that although this time scaling may be useful for comparing some (macroscale) simulation phenomena with appropriate observations, it is not appropriate for describing physical effects that result from ion-driven turbulence developing on the ion cyclotron scales. In what follows, we measure simulation time in the units of Ω −1 cp ≈ 2 s.

Simulations and discussion
In what follows, we present simulation data from the uniform mesh patch only, in GSM coordinates, measured in FIGURE 4 O+ temperature and O+/H+ number ratio in two cross-sections in the first lobe outflow simulation at Ω cp t = 461: x-z ( y = 1d i ) and y-z (x = −199d i ). solar wind proton inertial length (d i ) units. When converting these distances to distances expressed in terms of the Earth's radius (R E ), we assume that the nominal magnetopause standoff distance (used to calculate the strength of the magnetic dipole) is 10R E = 100d i . All simulation quantities are normalized by their appropriate upstream values: Velocities by the Alvfén speed, V A , temperatures (velocity squares) by V 2 A , number densities by n 0 , Poynting vector components by (V A /4π)B 2 IMF , and timesteps by the inverse proton cyclotron frequency, Ω −1 cp . To better assess the impact of different O+ outflows on the magnetotail structure, we first conducted a "baseline" simulation without outflow, where simulation parameters were chosen to ensure that no X-points appear in the tail (including its part that lies outside of the uniform patch shown) during the entire time period, Ω cp t ∼ 600, adopted for all simulations. Results from this baseline simulation are illustrated in Figure 2. In particular, this Figure presents a regular, well-developed proton (H+) tail, shown at a late simulation time, which is also used to present results from outflow simulations. The main focus of this paper is to demonstrate how different types of O+ outflow may produce tail configurations that are considerably different compared to the baseline case. Other simulation details, e.g., those related to the temporal dynamics of the magnetotail, are left for future discussions.

Outflow setup
Following (Yu and Ridley, 2013a;Yu and Ridley, 2013b), we initialize nightside O+ outflows in both hemispheres inside wedgeshaped regions located between 11 p.m. and 1 a.m. magnetic local time (MLT) and within specified magnetic latitudes (ML). In these simulations, O+ ions with a parallel speed of 100 km/s and thermal speed of 35 km/s are launched along the magnetic field lines from the inner spherical surface with a radius, R = 65d i , which approximately corresponds to R = 6.5R E in the full-size Earth's magnetosphere. The initial O+ energy in our simulations is consistent with observations at high altitudes (Krcelic et al., 2020) and close to O+ energy values used in the previous MHD simulations (e.g., Yu and Ridley, 2013a;Yu and Ridley, 2013b).
We consider two outflow scenarios, where O+ is injected from different nightside regions of the inner simulation boundary. In the "lobe" scenario, O+ ions are launched between 50 and 75°ML. These values respectively correspond to approximately 75 and 84°a t the ionospheric altitude (R = 1R E ), when tracing ion trajectories from the inner boundary along the unperturbed Earth's dipole lines down to the ionosphere. Physically, these injection points indirectly model the dayside cusp O+ outflow, originating from the ionospheric altitude and streaming along open field lines into the lobes, as discussed in Section 1. For this case, we perform two simulations with different total hemispheric O+ fluences: 1.9 × 10 26 FIGURE 5 H+ number density and O+/H+ number ratio in two cross-sections in the second lobe outflow simulation at Ω cp t = 314: x-z ( y = 1d i ) and y-z (x = −199d i ).
ions/s (simulation 1) and 3.8 × 10 26 ions/s (simulation 2). In the other, "auroral" scenario, we launch O+ ions between 35 and 50°ML, corresponding to 71 and 75°, respectively, at the ionospheric altitude. In this case, the total hemispheric O+ fluence is approximately 2.6 × 10 26 ions/s. The amount of outflow flux has been shown to be an important factor in determining the global state of the system . The total O+ fluences chosen in our simulations are close to the peak O+ fluence values measured in observations (Cully et al., 2003a) and maximum fluences previously used in multi-fluid MHD studies (Garcia et al., 2010;Yu and Ridley, 2013a;Yu and Ridley, 2013b).
As mentioned in Section 3, the O+ outflow is launched approximately when the southward IMF arrives at the subsolar point of the magnetosphere. Following its onset, the streaming ionospheric O+ ions mass load the plasma sheet, increase its thermal pressure, and generate turbulent electromagnetic fields that affect the initial (proton only) magnetotail configuration. Below we discuss global responses of the magnetosphere to the lobe and auroral outflows, as defined above.

Lobe outflow
In this scenario, O+ ions are injected on the nightside of the inner obstacle surface at high latitudes and arrive at the tail by flowing along open magnetic fields. Without outflow, these lines stretch into the lobes, being void of plasma under steady solar wind/IMF conditions. Figure 3 illustrates the proton tail density and absolute value of the net plasma current density, |J| at Ω cp t = 461 in the first simulation (the less intense outflow). Compared to the baseline case, (see Figure 2), this Figure clearly shows a significantly perturbed state of the magnetotail, with a turbulent, "wiggling" current sheet. Multiple X-points develop between ∼15 − 20R E during this simulation, with eventually a dominant X-point forming near ∼20R E .
In this case, the current in the y − z plane demonstrates an asymmetry of the current sheet in the dawn-dusk direction, with denser currents forming on the duskside (positive y). This correlates with the presence of more energetic O+ ions in the plasma sheet on the duskside, as shown in Figure 4 (left column). The dawn-dusk asymmetry in our outflow simulations arises due to the Hall component of electric field that tends to push largeorbit O+ towards the duskside. This effect may be responsible for asymmetric reconnection signatures in observations reported by Baker et al. (1982), as discussed in Section 1. We observe this dawn-dusk asymmetry to be more pronounced in the second lobe simulation, performed with twice as many outflowing O+ ions. The results from this simulation are illustrated in Figure 5.
As expected, this Figure shows a significantly more compressed H+ plasma (especially on the duskside), consistent with a larger current density observed in the sheet (not shown). The X-point Ion velocity (V0 stands for H+ and V1 for O+) and Poynting vector (S) components in two cross-sections in the first lobe outflow simulation at Ω cp t = 314: x-z ( y = 1d i ) and y-z (x = −199d i ).
in this case case develops at approximately the same distance of 20R E .
The stronger outflow in the second lobe simulation serves the purpose of generating a faster response of the tail to the "driver" in the absence of more sophisticated details, e.g., initial amounts of O+ in the tail. This simulation demonstrates an interesting feature, shown in the right column of Figure 5, namely, the pushing of O+ ions out of the current sheet to the outer parts if the tail. This configuration, where O+ ions mostly populate the outer regions of the current sheet, correlates with a theoretical model of multilayered current sheet . This model takes into account the contribution of non-adiabatic O+ ions as a term in the Grad-Shafranov-like system of equations, which describe the quasi-equilibrium of a multi-component (H+ and O+) tail. In this theoretical model, the bulk of the sheet current is carried by H+, while the current "wings" carried by O+.
As noted above, the hybrid model includes many physical effects not accounted for by theory or multi-fluid MHD simulations, among which ion-driven turbulence and ion energization play the most important roles. For instance, Figure 4 illustrates the transient dynamics and heating of O+ ions (left column), along with their relative population in the current sheet with respect to H+ ions (right column) in the first lobe outflow simulation. We observe that the higher temperatures of O+ ions correlate well with the regions of enhanced current density, shown in Figure 3 (right column). This property of O+ ions may serve as a useful diagnostic for future satellite observations. Kinetic ion effects account for ion-scale plasma turbulence that plays an important role in the local and global dynamics of FIGURE 7 H+ number and absolute value of net current density, |J| in three cross-sections in the auroral outflow simulation at Ω cp t = 461: x-z ( y = 1d i ), x-y (z = 1d i ), and y-z (x = −119d i ).
the dayside magnetosphere (Omelchenko et al., 2021a). Likewise, large-orbit O+ ions perturb the plasma sheet (Baker et al., 1982), enhancing the turbulent effects in the magnetotail, as evident from comparing the tail configurations with and without outflow in our simulations (see Figures 2, 3). As a result, the local inflow, redistribution, and outflow of energy in the tail plasma form a highly dynamic process, where the rate of magnetic flux conversion and concomitant energy dissipation may be enhanced or even controlled by turbulent reconnection, as previously suggested for astrophysical plasmas by Lazarian et al. (2020).  O+ temperature and number density in two cross-sections in the auroral outflow simulation at Ω cp t = 461: x-z ( y = 1d i ) and y-z (x = −119d i ).

FIGURE 9
Cell-averaged H+ and O+ particle timesteps in two cross-sections in the first lobe outflow simulation at Ω cp t = 461: x-z ( y = 1d i ) and y-z (x = −199d i ).
Frontiers in Astronomy and Space Sciences 12 frontiersin.org

Auroral outflow
In this scenario, we inject O+ ions on the nightside of the inner obstacle surface at lower latitudes, so that they directly move into the tail plasma by following closed magnetic field lines in the inner magnetosphere. In MHD simulations of this type of outflow (Yu and Ridley, 2013a), the magnetotail appeared to be steady through the entire simulation, with the X-line staying at around the same distance (∼20R E ), contrary to the dayside cusp ("lobe") outflow case, where the tail was found to experience dramatic changes. Yu and Ridley (2013a) note that the mass loading in the tail in their auroral outflow case did not have a definitive role in initializing or preventing substorm onset, even though a growing amount of magnetic flux was advected into the nominal reconnection region (∼20R E ). At the same time, their dayside cusp outflow simulation indicated that heavy O+ ions in the tail helped trigger the tail reconfiguration by compressing and disrupting the neutral sheet. Based on these results, Yu and Ridley (2013a) conclude that the degree of impact of O+ ions on the tail strongly depends on their pathways from the ionosphere into the tail.
The results from our simulation contradict their finding of the insignificant influence of auroral O+ outflow on the magnetotail. Figure 7 shows the H+ density and absolute value of net current density, |J| in the tail for this case at Ω cp t = 461. These plots clearly indicate that the O+ outflow has dramatically affected the tail stability and resulted in inducing an X-point at x ≃ −12RE, much closer to the Earth than the X-point distance in the lobe outflow case (∼20RE). Moreover, the H+ tail in this (auroral) case undergoes a dramatic disruption in the near-Earth space, with H+ ions being completely flushed out by the streaming O+ ions in some parts of the tail.
As discussed in Section 1, the recent in situ observations by Angelopoulos et al. (2020) have suggested that magnetotail reconnection may be powered by intense storms at close geocentric distances, ∼10R E , where no reconnection is normally expected due to the Earth's strong dipole field. The same observations found a significant percentage of O+ ions (up to ∼50%) in the tail in addition to protons. Although the upstream conditions in our simulation are different from those observed by Angelopoulos et al. (2020) during a storm, our simulation indicates that strong auroral outflows may serve as a plausible contributing factor to inducing reconnection in this region. This conclusion is also corroborated by results from multi-fluid simulations by Brambles et al. (2011), who demonstrated an X-line forming around ∼14R E when auroral outflow was included. Figure 8 shows temperature and density distributions of O+ ions in the auroral outflow case. By comparing the y-z density panel (O+) in this Figure with the y-z density panel (H+) in Figure 7, we conclude that plasma sheet protons are flushed away from the X-point by O+ ions. This leads to establishing a predominant population of O+ near the X-point and in the outer parts of the tail in the z-direction. This effect may also be related to the findings by Karimabadi et al. (2011), who demonstrated in their local simulations that during the non-linear stage of reconnection the initially dominating protons in the current sheet can be replaced (through a "flushing effect") by O+ ions coming from the outer parts of the sheet.

Summary and conclusion
Using a multiscale, self-adaptive hybrid code, HYPERS (Omelchenko and Karimabadi, 2012a;Omelchenko et al., 2021a), we have explored possible impacts of ionospheric O+ outflow on the structure, stability, and dynamics of magnetotail for two different scenarios: "lobe" and "auroral" outflows. In these simulations, we have treated ions of both solar wind (H+) and ionospheric (O+) origin as full-orbit particles with a self-consistent PIC approach.
Our simulations account for the important kinetic (non-Maxwellian, ion cyclotron) and Hall effects that are missing from multi-fluid MHD models. The hybrid model expands the physical richness of global modeling of the Earth's magnetosphere by adding kinetic turbulence, ion cyclotron motion, and ion energization effects on top of modeling global flow patterns that generally develop in the magnetosphere. Although this increase in physical fidelity comes at a cost of using a downsized Earth model, the physics extracted from hybrid simulations provides many important details, regularly confirmed by space observations (Omelchenko et al., 2021a) and laboratory experiments (Omelchenko and Karimabadi, 2012b;Omelchenko and Karimabadi, 2023) but not described by MHD models (Delzanno and Borovsky, 2022). As a first step towards self-consistent kinetic modeling of ionospheric outflow in the Earth's magnetosphere, in this paper we have chosen to explore O+ outflows. In the first ("lobe") scenario, we have launched O+ ions from higher magnetic latitudes. Physically, this case corresponds to the dayside cusp outflow, where O+ streams along open magnetic field lines over the polar cap on the nightside into the lobes. In the second ("auroral") scenario, we have injected O+ ions from lower magnetic latitudes into the inner magnetosphere along closed field lines. To model these outflows, we have used hemispheric outflow fluences (the total number of O+ ions injected per second) on the higher side of the spectrum of available observations, similar to those used in the previous MHD simulations.
In the case of lobe O+ ions streaming into the magnetotail, we have discovered that O+ outflow can significantly disrupt the tail structure by inducing X-points at ∼20R E distances, which is consistent with MHD results by Yu and Ridley (2013a), Yu and Ridley (2013b). Our simulations, however, emphasize the complex internal dynamics and heating of O+ ions, which are missing from the MHD simulations. In addition, by launching a stronger lobe outflow, we have also demonstrated a multi-layer tail configuration predicted by analytical theory .
In the case of auroral O+ outflow, our findings disagree with the results from the MHD simulations performed by Yu and Ridley (2013a) for a similar scenario (we note, however, that the numerical setups of outflows in our and their simulations are not identical). In particular, we have found that a strong auroral O+ outflow can induce an X-point in the near-Earth tail, which is consistent with storm observations by Angelopoulos et al. (2020) that also revealed substantial amounts of O+ near the reconnection site.
The multiscale outflow simulations presented in this paper have been powered by an event-driven "simulation time operating system", EMAPS (Omelchenko and Karimabadi, 2006a;Omelchenko and Karimabadi, 2007;Omelchenko and Karimabadi, 2023). EMAPS replaces time stepping by discrete event simulation. In HYPERS, this enables asynchronous advance of both particles and electromagnetic fields on their local timescales, in accordance with their physical rates and stretched mesh geometry. EMAPS orchestrates asynchronous computation by adaptively limiting perupdate increments of all variables to physically reasonable values. In global magnetospheric simulations, this approach prevents numerical instabilities, achieves superior accuracy, and produces significant computational gains, compared to similar time-stepping simulations, where choosing proper timesteps is always a "tug-ofwar" battle between code performance on one side and simulation accuracy and stability on the other. To illustrate this point, in Figure 9 we provide snapshots of the cell-averaged H+ and O+ particle timesteps for the time moment used to generate Figure 3 and Figure 4. Smaller particle timesteps in Figure 9 indicate the need for the simulation to increase the local temporal resolution due to the higher local particle energies and/or stronger magnetic fields. For simplicity, we do not show here the field timesteps and timestep distributions for the whole stretched mesh that also take into account the cell size variations (Omelchenko et al., 2021a).
No single computational study can answer all questions related to the observations of ionospheric O+ ions in the magnetosphere. To further explore the impact of ionospheric outflow on the magnetosphere, especially as a function of different ion sources, particle fluences, source locations, and other magnetospheric parameters, more studies need to be conducted in the future. Although there exist many different physical scenarios for exploring the effects of ionospheric ions both locally and globally, the hybrid model (kinetic ions, fluid electrons) presents a particularly useful computational venue for these studies. To improve our knowledge of the Earth's magnetosphere, results from hybrid simulations should also be compared with properly scaled MHD and full PIC simulations. By simulating the impact of O+ outflow on the magnetotail, we have made a step in this direction.

Data availability statement
The datasets used in the current study are available from the corresponding author on reasonable request, or they can be reproduced by running the HYPERS-Global model at the NASA CCMC: https://ccmc.gsfc.nasa.gov/.

Author contributions
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.