The Role of Mesoscale Plasma Sheet Dynamics in Ring Current Formation

During geomagnetically active periods ions are transported from the magnetotail into the inner magnetosphere and accelerated to energies of tens to hundreds of keV. These energetic ions, of mixed composition with the most important species being H+ and O+, become the dominant source of plasma pressure in the inner magnetosphere. Ion transport and acceleration can occur at different spatial and temporal scales ranging from global quasi-steady convection to localized impulsive injection events and may depend on the ion gyroradius. In this study we ascertain the relative importance of mesoscale flow structures and the effects of ion non-adiabaticity on the produced ring current. For this we use: global magnetohydrodynamic (MHD) simulations to generate self-consistent electromagnetic fields under typical driving conditions which exhibit bursty bulk flows (BBFs); and injected test particles, initialized to match the plasma moments of the MHD simulation, and subsequently evolved according to the kinetic equations of motion. We show that the BBFs produced by our simulation reproduce thermodynamic and magnetic statistics from in situ measurements and are numerically robust. Mining the simulation data we create a data set, over a billion points, connecting particle transport to characteristics of the MHD flow. From this we show that mesoscale bubbles, localized depleted entropy regions, and particle gradient drifts are critical for ion transport. Finally we show, using identical particle ensembles with varying mass, that O+ non-adiabaticity creates qualitative differences in energization and spatial distribution while H+ non-adiabaticity has non-negligible implications for loss timescales.


INTRODUCTION
The modern understanding of transport processes, and their multiscale nature, through the magnetotail transition region has advanced considerably in the decades following the early work that first recognized their importance (e.g., Mauk and McIlwain, 1974). The definition of what constitutes the transition region varies in the literature (e.g., Ohtani and Motoba, 2017;Sergeev et al., 2018); for the purposes of this work we will take a more expansive definition of the magnetotail transition region, ∼ 8 − 20R E tailwards of Earth. It is now known that much of the plasma transport in the transition region occurs by means of transient ( ∼ 10 minute), fast ( ∼ 400 km/s) bursty bulk flows (BBFs; Baumjohann et al., 1990;Angelopoulos et al., 1992, Angelopoulos et al., 1994. These flow bursts have typical cross-tail sizes of 1-3R E (e.g., Nakamura et al., 2004;Liu et al., 2013), which we will refer to here as mesoscale to distinguish them from both global and kinetic scales in the magnetosphere. Increasingly these BBFs have been connected to localized ion injections into the ring current (Turner et al., 2017;Gabrielse et al., 2019) which can provide a substantial source of ring current pressure during stormtime (Gkioulidou et al., 2014).
A theoretical description of these flow bursts was provided by Pontius and Wolf. (1990) who identified fast flows as entropydepleted "bubbles" that would be interchange unstable within the ambient entropy profile radially-decreasing away from Earth. Confirming that flow bursts are bubbles using in situ data is complicated by both the sparsity of spacecraft in the transition region and the necessity of extrapolating flux tube entropy from point measurements of fast flows (e.g., Wolf et al., 2006). However, it has been found that inferred entropy-depletion from spacecraft measurements is a good predictor of the penetration depth of flow bursts (Dubyagin et al., 2011;Kim et al., 2012). This interpretation was supported by regional modeling, i.e., localized subdomains of the global magnetosphere, which highlighted the role of reconnection, and the resultant reduction in flux tube volume, in creating bubbles (e.g., Birn et al., 2009Birn et al., , 2011Liu et al., 2014) and the role of buoyancy (e.g., Yang et al., 2010, Yang et al., 2015Sadeghzadeh et al., 2021) in the subsequent transport of these bubbles into the inner magnetosphere. 1 Despite these advances, a full multiscale understanding of the transition region has remained elusive due to both observational and modeling limitations. Important questions remain regarding: the relative role of large-scale versus mesoscale convection; the statistical distribution of bubble lengthscales and the physical mechanisms that govern the size of bubbles; and the mechanisms of ion energization. Sparse spacecraft coverage of the transition region makes it difficult to spatially localize injections, although there has been progress using multi-mission conjunctions (Turner et al., 2017) and remote imaging (Keesee et al., 2021). Modeling multiscale transport in the transition requires: energydependent drifts, not included in global magnetohydrodynamic (MHD) models which approximate the bulk flow with the E × B drift; fast flows outside the quasi-static slow-flow approximation of typical ring current models (see e.g., Toffoletto, 2020, and references therein); global models that produce mesoscale magnetotail structures (Wiltberger et al., 2015;Cramer et al., 2017;Merkin et al., 2019), requiring high spatial resolution and low dissipation algorithms; a representation of wave-particle interactions that contribute to ion heating (e.g., Chaston et al., 2014;Cheng et al., 2020); and finite gyroradius effects, as H+ in the transition region is known to be quasi-adiabatic , while heavier ion species like O+, likely exhibit highly nonadiabatic behavior (e.g., Moebius et al., 1987;Delcourt et al., 1997;Nosé et al., 2000;Keika et al., 2013;Bingham et al., 2020). The latter point is particularly important. While the details of the O+ energization mechanisms and transport into the ring current are not yet fully understood (Ohtani et al., 2011;Keika et al., 2013;Gkioulidou et al., 2019), O+ is known to contribute a significant, if not dominant, amount of the total ring current energy density during increasingly disturbed geomagnetic conditions (Nosé et al., 2005).
Given these challenges it is perhaps not surprising that there has been little work studying bubble-mediated transport through the transition region in a global context. Ukhorskiy et al. (2018) used global MHD and test particle modeling to investigate particle interactions with a single bubble to demonstrate the importance of the mesoscale structure of the bubbles to trap the particles, whereas Cramer et al. (2017) performed a more expansive study over a range of stormtime events to demonstrate the importance of bubbles in transport. In a series of papers representing the most self-consistent work to date (Lin et al., 2017(Lin et al., , 2021Cheng et al., 2020), global hybrid simulations were used to demonstrate that reconnection creates depleted bubbles that are anisotropic, non-Maxwellian, and not in equilibrium. However, spatial rescaling is necessary for global hybrid simulations to be computationally tractable and even then the duration of these simulations remains limited.
Regional modeling has been used to reveal many aspects regarding transport through the transition region, while global modeling faces substantial intrinsic difficulties. Despite this, global modeling of the transition region is necessary to probe the complex magnetosphere-ionosphere coupling that binds the two geospace domains together. Mesoscale plasma sheet injections have global-scale consequences through their role in the ring current build-up. The partial ring current drives "Region 2" field-aligned currents (FACs; e.g., Roelof, 1989) into the ionosphere where they connect to the "Region 1" currents driven by magnetospheric processes at higher-latitudes. The mutual closure of these currents in the ionosphere occurs across the electrojets (Baumjohann, 1982), the most intense electrical currents in the ionosphere. It is in these auroral regions where strong Joule heating and momentum transfer between the ions and neutrals occur, stirring up the thermosphere. Changes in the high-latitude neutral temperature and winds are transmitted to lower latitudes through non-linear, coupled, dynamical processes, causing global-scale variations in neutral temperature, winds and composition (Li et al., 2019). This alters global ionospheric plasma densities through chemistry and plasma transport (e.g., Fuller-Rowell et al., 1994;Wang et al., 2010;Lu et al., 2016). Changes in plasma densities are also produced by energetic particle precipitation at different altitudes, depending on their energy spectra. All this results in complex ionospheric conductivity changes and electrodynamic feedback to the magnetosphere. Furthermore, Joule heating, soft electron precipitation, and plasma transport cause upwelling of ionospheric ions that enables mass coupling with the magnetosphere by providing the source for ion outflow. Therefore a systems-level understanding of geospace requires a self-consistent and cross-scale (global, meso-, and ion kinetic) treatment of the transition region.
It is the goal of this work, which extends aspects of Wiltberger et al. (2015) and Ukhorskiy et al. (2018), to better inform what aspects of the transition region transport are most critical for a global, self-consistent model to capture. To this end we utilize an approach combining MHD and test particle simulations wherein we create a test particle "mirror" of the plasma sheet, matching the spatiotemporal variations of density and temperature, whose evolution we can follow to gauge the importance of certain kinetic effects. This approach allows us to include the self-consistent formation and propagation of BBFs, within the framework of the MHD approximation, and the effects of drift physics and particle non-adiabaticity, within the framework of test particle trajectories. However, as these are test particles, there is no kinetic feedback to the global MHD electromagnetic fields. This caveat must be taken into account in the interpretation of these results. Nevertheless, our approach allows us to use very highly-resolved fluid models that produce mesoscale flow structures and to study their implications on ion transport and acceleration in the transition region.
With the appropriate caveats in mind, we will assess the importance of flow and particle characteristics to the transport of ions across the transition region. In particular we consider: flow characteristics like entropy-depletion and mesoscale structure; and particle characteristics like energy, and consequently energy-dependent drifts, and gyroradius. With the combined MHD and test particle approach we employ, we are able to mine simulation data to construct a database, exceeding a billion points, of instantaneous correlations between the inward radial transport of particle energy and properties of the ambient MHD flow. Given the importance of entropy, we find it critical to investigate the ability of our simulation to accurately reproduce the thermodynamic plasma properties of BBFs, which we verify by comparing with the THEMIS-based study of Runov et al. (2015). Next, by comparing identical test particle ensembles whose initial state is informed by MHD and are subsequently evolved with varying gyroradii, we are able to quantitatively assess the importance of non-adiabaticity and its implications to both the build-up and decay of the ring current.
This work is an extension of Ukhorskiy et al. (2018) but incorporates significant advances to the methodology. Ukhorskiy et al. (2018) seeded test particles into a single bubble and followed the evolution for ∼ 10 minutes while using a fiducial plasma sheet temperature to assess the implications to ring current buildup. Therefore, that work did not directly assess the comparative evolution of particles inside versus outside of bubbles or the result of particles interacting with multiple bubbles over longer timescales. Additionally, in this work we use flux tube entropy as a metric to explicitly study the connection between the sharp magnetic field gradients, identified as important by Ukhorskiy et al. (2018), to depleted entropy bubbles. This work is also a natural extension to the statistical data-model comparison performed by Wiltberger et al. (2015), where they used statistics derived from an idealized steady magnetospheric convection, or SMC-like period, to construct a superposed epoch analysis which they compared to Ohtani et al. (2004). The analysis of Wiltberger et al. (2015) was limited in its ability to disentangle whether the data and model differed due to differences in the properties of flow bursts or in the ambient state of the plasma sheet in their idealized simulation. Here we construct our data-model comparison to assess the relative properties of the bubbles to the ambient plasma sheet to focus on the ability of our model to reproduce how depleted bubbles are and how deeply into the inner magnetosphere they will therefore penetrate.
Having outlined the goals and limitations of our study, we now describe the plan for the remainder of the paper. In the next section, we provide an overview of the MHD (Section 2.1) and test particle (Section 2.2) simulations as well as the diagnostics (Section 2.3) we employ. Our main results are presented in Section 3 and include: a comparison of the thermodynamic properties of our modeled BBFs to observations (Section 3.1); a statistical study of the correlation between flow and particle properties to the transport and acceleration of test particle energy into the inner magnetosphere (Section 3.2); and a study of the effects of increasing non-adiabaticity relative to a guiding center ensemble (Section 3.3). Finally, we present a discussion of our results and their implications for future magnetospheric modeling in Section 4.

METHODOLOGY
Our technical approach combines two simulation methods: GAMERA (Zhang et al., 2019;Sorathia et al., 2020), a global MHD model, is used to generate electromagnetic fields; and CHIMP , a test particle code, is used to calculate particle trajectories through the time-dependent MHD FIGURE 1 | Simulation at a glance. Combined visualization of the MHD simulation and test particle data. Background color (bottom let color bar) denotes residual, i.e., non-dipolar, vertical magnetic field. Test particles are marked at the location of their field-line projection to the equatorial plane, colored by their energy (bottom right color bar), and denoted with a marker whose surface area is proportional to the logarithm of their weight. The region where test particles are seeded is marked in blue. Figure inset shows a closeup view of the test particle population in a typical BBF.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 761875 3 electromagnetic fields and weight the test particles to match a specified phase space density (PSD). In the remainder of this section we will describe in turn the details of the global magnetospheric simulation (Section 2.1), the test particle parameters and weighting (Section 2.2), and the diagnostics we will employ (Section 2.3). For convenience, however, we will provide a brief, high-level overview here. We drive the global magnetospheric simulation in an idealized configuration using pure southward interplanetary magnetic field (IMF), following a period of preconditioning. Within that global magnetospheric simulation we continuously seed new test particles in a thin arc on the nightside and weight them using the instantaneous density and temperature from the MHD simulation. In other words, we create a test particle "mirror" of the MHD flow moving through the nightside injection region. This is a mirror only within the injection region, upon moving inwards the test particles will evolve differently than the MHD flow they were created to mimic due to the kinetic equations of motion they are governed by. We seed test particles for a period of 1 h, beginning at T −1 h and ending at T 0 h, and continue to evolve those test particles that remain in the closed field region for an additional 1 h period, ending at T + 1 h. Finally, we look at the instantaneous correlation between the transport of test particle energy into the inner magnetosphere and the properties of the ambient MHD flow. In particular, we focus on the importance of the characteristic lengthscale of the ambient magnetic field, assessing the multiscale flows of the transition region, and the deviation of the flux tube entropy from an averaged background value, i.e., buoyancy, to assess the role of bubbles.
A visualization of the combined model output half-way through the simulation is illustrated in Figure 1. It depicts both the residual, or non-dipolar, magnetic field calculated by GAMERA, defined as the northward component (B Z ) with dipole contribution subtracted, and a sample of test particles with markers denoting the projection from their 3D position along the local magnetic field to the equator with the marker colored by energy and sized so that its diameter is proportional to the logarithm of the test particle's weight. The thin nightside region where test particles are continuously seeded for T < 0 is marked in blue. The inset of Figure 1 shows a zoomed-in view of a BBF having recently moved through the test particle injection region; the test particles inside the BBF are higher energy and lower weight than those outside, reflecting, as expected, the hotter, less dense plasma inside the BBF.

Global Magnetosphere
The first step in our model pipeline is the generation of electromagnetic fields and plasma moments using our global magnetosphere model GAMERA. The details of GAMERA's core MHD numerics and its validation via a standard portfolio of MHD test problems (e.g., Stone et al., 2008) were presented in detail by Zhang et al. (2019) and its first magnetospheric application by Sorathia et al. (2020). Zhang et al. (2019) also showed the critical importance of algorithmic details beyond simply the size of grid cells. In particular, the use of high-order spatial reconstruction can dramatically improve the preservation of sharp structures compared to lower-order schemes on an otherwise identical spatial grid (see Zhang et al., 2019, Section 4). For typical MHD test problems they find that lower-order reconstruction (e.g. 2 nd -order) requires 4-8× finer grid resolution as the higher-order (7 th -or 8 th -order) reconstruction used by GAMERA for the same effective resolving capability. Mapped to 3D, 4-8× the grid resolution becomes a factor of ∼ 250 − 4000 the computational cost. As we will show below, this extra resolving power is critical to capture sharp flow structures and their subsequent effect on particle evolution.
As a first study of particle transport with GAMERA, we adopt a similar approach used in previous BBF studies with GAMERA's predecessor LFM (Lyon et al., 2004). As in Wiltberger et al. (2015), we utilize an idealized magnetospheric configuration to generate an SMC-like period. Specifically, after a preconditioning phase of alternating northward and southward IMF, the IMF turns southward and remains so for the remainder of the simulation. The solar wind driving uses typical values: 5/cc density, 400 km/s velocity with only an SM-X component, and an IMF of 5 nT directed either northwards or southwards. Two hours after the final southward IMF turning, at which point steady reconnection has begun producing regular BBFs moving earthward, we begin injecting test particles whose properties are described below (Section 2.2). The inner boundary condition of the simulation is imposed at a spherical surface of 2R E via closing the field-aligned currents through a thin-shell ionosphere as described by Merkin and Lyon. (2010). In the idealized configuration we present here, we use a constant Pedersen conductance of Σ P 10 S with no Hall conductance and an ionospheric grid of 1°× 1°. It is important to note that these are ideal MHD simulations without an explicit resistivity, although there is numerical resistivity that allows reconnection. While this can have an effect on the specific details of reconnection onset (e.g., Raeder et al., 2001), we do not expect this to affect the transport and acceleration of our test particles that are seeded earthward of the reconnection line.
GAMERA, like LFM before it, uses for its magnetospheric simulations a warped spherical grid, whose spherical axis is aligned with the SM-X axis. The grid encompasses a region extending sunward to 30R E and tailward to −350R E , and to 120R E radius in the YZ-plane. The warped grid is distorted to maximize resolution in the regions where it is most important, e.g. the bow shock, the magnetopause, and the near-Earth plasma sheet. This simulation was run using 192 × 192 × 256 cells in the radial, polar, and azimuthal directions and is comparable to the highest-resolution capability of LFM, known as "OCT" resolution, which has been used to study mesoscale plasma sheet flows before (Wiltberger et al., 2015;Merkin et al., 2019). In particular, the nominal resolution in the near-Earth plasma sheet is approximately 600 km.

Test Particles and Weighting
2 h after the final southward turning, we begin injecting test particles continuously into an arc tailward of Earth (see Figure 1). The arc is of width ΔR 1R E and centered on R 18R E , extending azimuthally 8 h in MLT centered at midnight. We refer to the beginning of test particle injections as T −1 h, and inject test particles continuously into the arc for a period of 1 h. At T 0, the Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 761875 4 injections are ceased and the existing particles are evolved until T + 1 h. Over the hour of continuous injection we create 20M test particles, or ∼ 5000 per second. These test particles are created with randomly chosen pitch angles, varying between field-aligned and anti-aligned, and energies, varying between 100 eV and 100 keV.
The numerical details of how test particle trajectories are calculated and how each test particle is weighted are described in detail in  and more recent improvements to the interpolation of GAMERA electromagnetic fields by CHIMP in Sorathia et al. (2019). We, however, provide a brief overview here for convenience. For ion trajectories we advance the particle position using the Lorentz force and Boris integrator. For electrons or positrons, we use a mixed integrator alternating between a guiding center formulation, utilizing the conservation of the first invariant, and the Lorentz trajectory depending on the ratio of the particle gyroradius to local magnetic field lengthscale. When a test particle encounters a region where ϵ > 10 -2 , where ϵ is the ratio of particle gyroradius to magnetic field lengthscale (Sorathia et al., 2018, Eq. (3)), a random gyrophase is chosen and the test particle is subsequently evolved using the full Lorentz equations of motion. The test particle can be converted back to the guiding center approximation when ϵ < 10 -2 by calculating an approximate guiding center location, which we do in the manner of Snicker et al. (2010), and the first invariant at the guiding center location.
In the work we present here we consider the evolution of an ensemble of test particles with identical initial conditions, but evolved as either positrons, i. e a guiding center ring current (GCRC), H+ (PRC), or heavy ions (HRC), specifically O+. This will allow us investigate separately the adiabatic interaction of particles with MHD-produced flow structures (GCRC), as well as the effects of increasing non-adiabaticity (PRC and HRC). Note that while we use the terminology RC for ring current because our ultimate interest is the role of these particles in building the ring current, our investigation includes their behavior throughout the nightside near-Earth plasma sheet. Data Sets S1-S3 contain a random sample of test particle trajectories for GCRC, PRC, and HRC respectively.
Given a collection of test particles and their trajectories, we assign each particle a weight such that the overall collection reproduces a desired initial condition for the PSD. Here the test particle weights represent the number of real particles each test particle acts as a proxy for, i.e. macroparticles. Note, however the weighting does not alter the underlying statistics of the test particles themselves it merely assigns weights to each preexisting test particle. The flexibility of this approach is that if we have a collection of test particles that properly samples the phase space then we can easily change the choice of PSD initial condition, and thus the weighting, without recalculating any test particle trajectories.
For this study our PSD initial condition is directly informed by the MHD flow properties. Test particles are weighted using a κ-distribution with the density and temperature of the MHD flow in the location and time of the test particles injection. Here we use κ 6 consistent with ion observations for this region of the near-Earth tail (Runov et al., 2015). We note that we use the same κ value across all species so as to focus on the role of nonadiabaticity on an otherwise identical ensemble, although we do not expect sensitivity to the choice of κ .
Weighting is done on a 4D discretized phase space spanning cylindrical radius, MLT, pitch angle, and energy. The phase space is reduced by one in each of its spatial and momentum components: phase space densities are averaged over entire flux tubes, corresponding to the radius and MLT of their equatorial crossing; and over gyrophase. In this work we use a phase space, Γ(R, ϕ, α EQ , K), spanning [3, 22] R E , [0, 2π], [0, π], and [0.1, 500] keV respectively. The phase space grid is discretized linearly in the angular variables with N ϕ 72 and N α 36 and logarithmically in radius and energy with N R 30 and N K 30. Newly-seeded test particles are weighted collectively at a cadence of δT 5 s. We ensure that the resolution of the phase space and number of test particles are appropriate by verifying the approximate invariance of the results when changing the discretization or number of test particles used for weighting.
While there are a number of technical details involved in this process, the result is merely to constrain the test particles in the injection region (Figure 1, blue arc) so that they reproduce, up to the phase space discretization chosen, the distribution corresponding to a κ-distribution using the MHD moments. Effectively we create a test particle mirror of the MHD flow passing through the test particle injection region. This is similar to , which used three discrete wedge-shaped injection regions on the nightside but with much finer spatial granularity.

Diagnostics
Our goal in this study is to connect properties of the background MHD flow, e.g., deviations of integrated flux tube entropy from the background, to metrics quantifying the transport and acceleration of test particles within those flow regions, e.g., the velocity of the equatorial crossing of a test particle. Here we will define the metrics we use for the background MHD flow and the test particles. Our primary interest is the "transition region," which for our purposes we define in an expansive way to encompass the region where fast flows and particle drift physics are both non-negligible. We define a "Data Collection Region", marked in Figure 2, as −22.5 ≤ X ≤ −8 and |Y| ≤ 15 from which we will take statistics.
Interchange instability, due to flux tube entropy imbalance, has long been believed to play an important role in the earthward transport of plasma (see recent review by Toffoletto, 2020, and references therein). To identify interchange unstable regions, we use the entropy function (e.g. Birn et al., 2009): where the integral is taken along the flux tube that contains the cylindrical coorinate (R, ϕ, z 0). The quantity S is conserved in ideal MHD and related, but not identical, to the thermodynamic entropy of a flux tube (Birn et al., 2009). This distinction is not relevant for this work and we will simply refer to Eq. (1) as the flux tube entropy (FTE). Figure 2 shows a representative equatorial distribution of flux tube integrated entropy, including examples of low-entropy intruding bubbles. Supplementary Video S1 is an animation of the simulation data in the same format as Figure 2. Interchange instability is governed by deviation of the local entropy from the background (Rosenbluth and Longmire, 1957), which we can quantify using what we henceforth refer to as the "relative buoyancy", where S 0 is an average background entropy calculated using a sliding time window, namely with Δt 10 min. The timescale for averaging is chosen to lie between the timescales characterizing the propagation of individual bubbles and those governing global reconfiguration of the tail during the SMC-like event we simulate. Eq. (2), the relative buoyancy, quantifies the fractional deviation of the local entropy from the background. Negative relative buoyancy corresponds to regions that will tend to move earthwards, i.e. bubbles, and positive relative buoyancy regions will tend to move outward, i.e. blobs. An example of the relative buoyancy, δS(R, ϕ), at the same time as Figures 1, 2 is shown in Supplementary Figure S1. An example of the manner of evolution of this quantity is shown in Supplementary Figure S2, which depicts δS(ϕ, t), i.e. as a function of MLT and time, at selected fixed radii.
Previous work (Gabrielse et al., 2017;Ukhorskiy et al., 2017;Ukhorskiy et al., 2018) has identified the role of mesoscale magnetic structures within BBFs to induce "trapping" of sufficiently energetic particles within closed contours of constant magnetic field strength. An example of these "magnetic islands" can be seen in the inset of Figure 2. To identify these types of mesoscale structures in the background flow, we define the quantity where the denominator is taken to be the L 2 norm of the Jacobian, restricted to equatorial derivatives, of the magnetic field.
Turning now to the test particles, we seek to investigate the correlation of particle transport and acceleration to properties of the background flow. To this end we consider the motion of each test particle's equatorial crossing point through the nominal transition region highlighted in Figure 2. For this exercise we focus on GCRC as our interest is really in the motion of the equatorial crossing point of the particle guiding center, and using GCRC allows us to remove the additional noise that would be introduced by gyromotion. We accomplish this by projecting each test particle to the equator at a temporal cadence of ΔT 5 s, with (X EQ , Y EQ , K) i,p designating the equatorial projection of particle p at time T T 0 + iΔT and its kinetic energy. Similarly we can calculate for each particle at each time: ΔR EQ , the rate of change of the radial equatorial crossing point; ΔK, the rate of change of energy; and α EQ , the equatorial pitch angle of the projected test particle, assuming conserved first invariant. Finally, we join to each tuple of test particle characteristics the flow characteristics at the equatorial projection, i.e. L ∇B (X, Y) and δS(X, Y). In this manner we create a large collection of data points of the form, which we restrict to those for which (X EQ , Y EQ ) lie in our nominal transition region. A randomly selected subset of these is included in Data Set S4. The collection defined this way, i.e. tuples of the form given by Eq. (5) for which (X EQ , Y EQ ) is in the nominal transition region, creates a data set of over a billion points. Our data analysis benefited from the use of fast-histogram 2 as intrinsic numpy routines were often insufficient to handle data at this scale. From this data set we construct statistical populations defined by a membership criteria and weighting of each data point. These statistical populations are: • Population (P): All tuples with weights given by w p , the standard weight of test particle p. • Transport (T ): All tuples for which ΔR EQ < 0, i.e. instances of inward radial transport, with weights given by w p KΔR EQ • Acceleration (A): All tuples for which ΔK > 0, i.e. instances of particle acceleration, with weights given by ΔKw p With these definitions, P represents simply the occurrence ratio. The collection T represents all instances of inward radial transport and the weighting is such that we consider the relevance of each data point to the total keV −R E of energy transport. In other words, one 20 keV particle transported inwards 5 R E is equivalent to ten 5 keV particles transported inwards 2 R E . Finally, population A represents all instances of individual particle acceleration and the weighting quantifies the importance of each data point to the bulk energization that occurs in the nominal transition region.

RESULTS
In the sections that follow we use weighted test particles to study particle transport and acceleration due to the global and mesoscale flow, as captured by our combined MHD and test particle simulations. In our model, test particles are weighted based on the density, temperature, and flux tube volume in the MHD simulation at the location and time of their injection. Therefore, prior to our discussion of particle transport, we show that the statistical properties of the BBFs, particularly those related to the quantities used in particle weighting, produced by our model are in agreement with observations (Section 3.1).
Having established this, we next turn to an investigation of energetic, but adiabatic, particle transport through the transition region (Section 3.2). Here we show the critical role of energetic particle drifts and entropy-depleted bubbles. Given the importance of the latter, we also discuss the sensitivity of our results to numerical resolution. Finally, we go beyond the guiding center approximation and quantify the effects of increasing particle gyroradius (Section 3.3). We find that nonadiabaticity results in significant and qualitative differences in particle behavior for heavier ions, and that even for H+ there are non-negligible effects.

Bursty Bulk Flows in Global Magnetohydrodynamics and Observations
In the context of a study of BBFs using GAMERA's predecessor, the LFM, and a similar global magnetosphere simulation to ours Wiltberger et al. (2015) used simulation data to conduct a superposed epoch analysis mirroring the Geotail study of Ohtani et al. (2004). Wiltberger et al. (2015) found that the model was able to qualitatively reproduce statistical profiles of velocity, magnetic field, and density within the BBFs. Before moving on to an analysis of particle transport and acceleration, we pause here to expand upon the validation efforts of Wiltberger et al. (2015) in a similar statistical manner. Runov et al. (2015) have conducted a statistical study, using THEMIS data over a 2 year period, quantifying the relative thermodynamic and magnetic properties of intruding regions, or dipolarizing flux bundles (DFBs), to the background flow. They identify intruding regions using a criteria informed by large dB z /dt and maximum B z , indicating a sharp magnetic front, small B x to exclude events when the spacecraft are in the lobes, and decreasing number density, to help select for depleted flux tubes. These criteria select approximately 300 events which are separated based into 4 groups based on their equatorial radius: R < 9.5, R ∈ [9.5, 12], R ∈ [12, 15.5], and R ∈ [15.5, 25]. Within each of these radial groups they calculate the ratio of density, temperature, and vertical magnetic field between the intruding to background populations. Their results (see Runov et al., 2015, Table 1) are reproduced here in Table 1 under the "DATA" columns, and demonstrate clear statistical relationships that are observed in nature.
In this paper, our approach to model validation is to consider data-model comparison as a question of whether the model can reproduce observed statistical relationships. This statistical approach helps alleviate, but not eliminate, the difficulties associated with comparing in situ measurements, characterized by low spatial density over a long time duration, with highresolution models, featuring high spatial density over a very limited temporal duration. Runov et al. (2015) use 300 measurements collected over 2 years, whereas we (as described below) can generate millions, or even billions when including test particles, of measurements over only 2 hours of, in this case, synthetic southward IMF driving. Care must be taken in how we interpret comparisons of a population with sparse sampling of a large event space to dense sampling of a very limited event space. With these caveats in mind, we return to the data from our simulation.
To construct our statistics, we calculate flux tube entropy throughout the nominal transition region defined in Section 2.3 at a cadence of 5 s. We identify the intruding population using a criteria defined by the relative buoyancy and density depletion, although we return to the significance of the latter criterion shortly. The background population is defined by |δS| < 0.05 and the intruding population is defined based on the criteria δS < − 0.075 and density depletion of 10% relative to the background. We note here that in principle we can identify everything for which δS < 0 as a bubble. However, we find that for small δS the time-varying background S(R, ϕ, t) complicates reliably separating bubbles from background. Therefore we use the stricter criterion δS < − 0.075 to identify bubbles. We also note that the use of a density depletion criterion only appreciably affects the results in the innermost radial bin (R < 9.5), outside of that the ratios we get with and without the additional density depletion criterion are largely unchanged. We attribute this to the fact that there are fewer bubbles that penetrate this deeply (Supplementary Figure S2) and that the statistics of the innermost region are skewed by the induced effects of bubbles, e.g. background plasma displaced from its entropy equilibrium by braking bubbles (e.g., . This interpretation is supported by the long-wavelength striations of oscillating relative buoyancy seen inside R 8R E (Supplementary Figure S2C).
With the criteria as defined above we generate approximately 10 7 total data points, of which ≈ 10 6 are identified as intruding.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 761875 7 The fraction of intruding to background is of course higher in the model than the observational data as we are focused only on the nominal transition region during an SMC-like period. The ratios of the average density, temperature, and magnetic field between the intruding and background populations are then calculated in each of the radial groups ( Table 1) under "MODEL". The statistical relationships inferred from the data, as characterized by Table 1 can be summarized as follows: the intruding regions are depleted of mass, typically exhibiting 60% of the background mass, this relationship is largely independent of radius but exhibits large deviations; the intruding regions are hotter than the background, typically by 40%; and the intruding regions are dipolarizing, i.e. exhibit stronger magnetic fields than the background, typically twice as strong outside of R 12 R E and decreasing within. As can be seen in the comparison of MODEL and DATA shown in Table 1, our model can, for the most part, quantitatively reproduce these relationships.
The quantities we have considered in our data-model comparison are of direct relevance to how we weight test particles: the thermodynamic quantities directly inform the κ-distribution function used as an initial condition of the PSD; and the magnetic field, via the flux tube volume, is included in the phase space volume element. Beyond their importance to the initial configuration of the test particles, we expect these flow characteristics to shape their subsequent development. Hotter flows will have a higher fraction of their energy density carried by particles whose energy-dependent drifts, i.e. gradient and curvature, are comparable to the bulk E × B flow (see e.g., Sorathia et al., 2017, Figure 12). More strongly magnetized flows may have lower flux tube volumes and possibly lower flux tube entropy, S, which can result in deeper penetration of the flow (Dubyagin et al., 2011;Kim et al., 2012). Therefore our ability to reproduce these statistics in the MHD flow do not just reflect our ability to accurately create particles at the tailward edge of the transition region but to evolve those particles throughout the transition region.

Role of Mesoscale Flow Structures in Ring Current Energization and Transport
The first particle study we consider is the transport of particle energy in the adiabatic limit. To this end we use positrons, which we refer to as "GCRC" for guiding center ring current, whose gyro-and bounce-frequencies are much higher than the MHD frequencies. An animation of the combined MHD and test particle simulation is shown in Supplementary Video S2 in the same format as Figure 1. In Figure 3 we show the evolution of energy content contained within various equatorial shells. The energy content within a shell of (cylindrical) radius R is simply w p K p where the sum is taken over all test particles whose equatorial projection is within the given radius. The shaded region in Figure 3 denotes the period in which new test particles are being continuously created (T < 0), whereas for T > 0 only the existing test particles within the domain are evolved. For the outer radial shells, R > 8, there is a drop off between the maximum energy content shortly after T 0 to the energy content at the end of the simulation. This is due to the loss, without replacement, of particles on open drift orbits. Inside of R 8 the energy content increases during the period of active injection and is largely flat afterwards as particles that penetrate that far reach closed drift orbits and remain trapped. This highlights how sensitive the energy content of the inner magnetosphere is to the penetration depth of the injected particles.  Runov et al., 2015) looking at the relative density, temperature and magnetic field strength of BBFs to the background plasma. We construct analogous, although not identical, statistics (MODEL) from our simulation to compare against. Next we turn to one of our core goals in this study, namely the identification of the ambient flow properties that lead to the transport and acceleration of particle energy into the inner magnetosphere, i.e. the formation of the ring current. Using the method we describe in Section 2.3 we construct ∼ 10 9 data-points representing the instantaneous correlation between: a test particle state, through its equatorial projection, energy, and pitch angle; its transport and acceleration, through the change in its equatorial projection and energy; and the ambient MHD flow properties, through the characteristic magnetic field lengthscale (Eq. 4) and the relative buoyancy (Eq. 2).

Radial Domain [R E ]
Each of the statistical samples we defined in Section 2.3 are a collection of tuples of the form of Eq. (5), defined via a sampledependent membership criteria and weight. We define the contribution of a quantity, Q in the range Q 0 + ΔQ, as the weighted-sum of the elements that satisfy both the sample membership criteria and Q ∈ [Q 0 , Q 0 + ΔQ]. In other words, the contribution of a given relative buoyancy is: to population (P), the probability of finding a particle in a flow with that property; to transport (T ), the fraction of all particle energy transport that occurs in regions of the flow with that property; and to acceleration (A), the fraction of the total bulk energization of the particles that occurs in regions of the flow with that property.
Absent any correlation between the particle or flow characteristic to transport, the population contribution and transport contribution would be the same. Regions for which the contribution to transport/acceleration is higher than the contribution to the population represents a characteristic favorable to the delivery of energy to the inner magnetosphere, e.g., we find that particles with energies above 10 keV "punch above their weight" and contribute more to transport and acceleration than their fraction of the population. Figure 4 compares the relative contribution of particles of different energies ( Figure 4A) and in different ambient flows (Figures 4B,C) to the overall population versus transport and acceleration. Before proceeding to the more finegrained discussion of the individual characteristics, we note that contribution to transport, specifically inward radial transport from the w p KΔR EQ weighting, correlates almost identically with that of acceleration. This is perhaps unsurprising as the transition region bridges the stretched tail and dipolar inner magnetosphere and therefore adiabatic transport (for GCRC) across it would lead to both betatron and Fermi acceleration via conservation of the first and second adiabatic invariant.
In each panel of Figure 4, the shaded region corresponds to a contribution of 50% of the total transport of energy through the transition region. The overlaid percentages represent the contribution within the shaded region to the population and the 50% contribution to transport. We can summarize the results of Figure 4 as follows: ambient flows that exhibit depleted entropy or spatial mesoscale structure of the magnetic field are particularly effective at delivering particle energy into the inner magnetosphere; similarly, particles with energies above ∼ 10 keV are particularly effective at carrying that energy into the inner magnetosphere. In the remainder of this section we analyze more closely these transport critical flow structures and the characteristics of the particles they delivery into the inner magnetosphere.
From Figure 4a, particles with K > 20 keV account for half the energy transport in the statistical sample despite making up only ∼ 16.5% of the population. Previous modeling work using test particles has identified the importance of the energy-dependent gradient drift, V ∇ , whose relative importance to the energyindependent E × B drift, also the MHD bulk flow, increases with energy (Gabrielse et al., 2017(Gabrielse et al., , 2016Ukhorskiy et al., 2017. As a simple estimate we can approximate the magnitudes of the E × B and gradient drifts as V EB ∼ E/B and respectively, and define the energy for which the two velocities are equal (see Sorathia et al., 2017, Figure 12). For typical values inferred from the MHD FIGURE 4 | Efficiency of transport and acceleration in the transition region. Here we compare for various particle and flow properties their fraction of: the total population, the total inward energy transport, and the total energization as defined in Section 2.3. Regions where the contribution of a quantity to transport/acceleration is greater than population are "efficient" for ring current buildup. Shown are particle energy (A), relative buoyancy Eq. (2); (B), and magnetic lengthscale Eq. (4); (C). In each panel the shaded region is chosen to highlight approximately 50% of the total energy transport.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 761875 9 simulation, E 2.5 mV/m and L ∇B 2R E , K ∇ ∼ 25 keV. This is an estimate of the energy at which gradient drift magnitude would exceed the bulk flow, it would be non-negligible for lower energies.
Using the statistics we have derived from this simulation, we can more quantitatively assess the role of particles whose trajectories deviate appreciably from the bulk flow. We find that the fraction of transport that occurs for V ∇ /V EB > 0.5 and V ∇ /V EB > 1 is 48 and 30% respectively. In other words, we find that the majority of transport through the transition region is via particles with non-negligible energy-dependent drift speed. Ukhorskiy et al. (2018) show that for particles with K > 10 keV, comparable to where we find the contribution to transport exceeds population ( Figure 4A), magnetic gradient trapping is necessary to confine particles within azimuthally localized flow long enough to reach the inner magnetosphere. While we do not try to distinguish magnetic gradient trapping in individual cases of our billion data points, the importance of K > 10 keV particles ( Figure 4A) and L ∇B < 3R E ( Figure 4C) strongly suggests the potential importance of this mechanism.
Focusing now on Figure 4B we see the critical role of entropy in particle transport. Low-entropy "bubbles", which we define here conservatively as δS < − 0.075, account for approximately half the energy transport in the sample while making up only ∼ 15% of the population. Less conservatively, δS < 0 accounts for ∼ 50% and ∼ 65% of population and transport, respectively. The efficiency increases as we consider only deeply depleted bubbles, δS < − 0.2, which correspond to 3 and 30% of the population and energy transport respectively. This is similar to previous studies using RCM (e.g., Lemon et al., 2004;Yang et al., 2010, Yang et al., 2015 and coupled RCM-MHD (e.g., Pembroke et al., 2012;Cramer et al., 2017) which have highlighted the importance of entropy-depleted flux tubes to ring current buildup. Of note is that we find that in addition to depleted entropy regions that there is, perhaps counter-inuitively, some inward transport in regions with δS > 0. We attribute this effect, as pointed out by Yang et al. (2011), to regions directly ahead of the depleted bubbles that are initially neutrally buoyant but pushed inwards by the bubble (see also Supplementary Figure  S2). However, we can see from this analysis that while those regions do contribute to more transport than their fraction of the population, the overall effect is much smaller than the bubbles themselves.
Finally, we consider the importance of the characteristic lengthscale of the ambient magnetic field as shown in Figure  4C. Here we find that while regions with L ∇B < 3R E are efficient, the effect beyond their contribution to the population is not as significant as that of particle energy or depleted entropy. From Figure 2 (inset) we see that regions of smaller magnetic lengthscale may correspond to: the interior of the depleted entropy bubble, which is particularly efficient at transport; the region ahead of the bubble, which is somewhat efficient; or the regions azimuthally adjacent to the flow channel for which we expect, if anything, outwards flow from the vorticity at the edges of the flow channel. We can instead consider the magnetic lengthscale statistics differently and look at the characteristic lengthscales of deeply depleted bubbles, δS < − 0.2. We find that the vast majority of deeply depleted bubbles are mesoscale, between 0.5-3R E . Specifically with Pr (A|B) denoting the conditional probability of A given B.
Supplementary Figure S3 shows the probability distribution of magnetic lengthscales for all bubbles, δS < − 0.075, and deeply depleted bubbles, δS < − 0.2. For both we find a largely similar distribution of magnetic lengthscales: sharply peaked at L ∇B ≈ 1.75R E , very similar to observationally-inferred values (e.g. Nakamura et al., 2004;Liu et al., 2013), and transitioning to a rapid decline for L ∇B > 4R E . Taken together, we find: bubbles are very efficient at transport, bubbles are mesoscale, and that bubbles transport particles whose energy-dependent drifts are comparable to the ambient flow. Given the critical role of buoyancy in transporting particle energy across the transition region, and consequently building the ring current, understanding the role of numerical resolution in producing these bubbles is vital. To this end we show in Figure 5 a resolution study of the buoyancy statistics for simulations with resolutions 2× coarser (QUAD) and 2× finer (HEX) than the simulation we consider in this paper (OCT). An analogous comparison of L ∇B is shown in Supplementary  Figure S4.
From our resolution study, we find that the lowest resolution differs most notably from the two higher resolutions, and that the two higher resolution simulations are largely similar for all but the most depleted bubbles, δS < − 0.7, which account for only ≈ 1% of the total energy transport in our test particle statistics at OCT resolution. By contrast we find that the lowest resolution simulation differs appreciably for δS < − 0.3, which corresponds to ≈ 20% of the total energy transport in the OCT simulation. From this, the lower resolution simulation would clearly be FIGURE 5 | Grid resolution study of the statistics of bubble depletion during SMC-like activity. Shown are statistics calculated using a grid 2× coarser (QUAD) and 2× finer (HEX) than the grid resolution we use elsewhere in this paper (OCT).
Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 761875 missing a non-trivial amount of energy transport assuming it behaved otherwise similarly to the higher resolution simulation. We return to this point in Section 4, but expect that the lower resolution simulation misses much more due to the depressed manifestation of L ∇B < 3 R E flow structures (Supplementary Figure S4).

Sensitivity to Larmor Radius
In the preceding section we considered the transport and acceleration of positrons (GCRC) as a proxy for adiabatic ion processes. Here we broaden our experiment to assess the role of non-adiabatic effects by considering more realistic ion masses. To this end we consider the exact same ensemble of created test particles, defined by their location and time of seeding, energy, and pitch angle and evolve them with a different mass. In addition to positrons (GCRC), we also include H+ (PRC), and heavier O+ ions (HRC). While the GCRC test particles are evolved using a hybrid guiding-center and Lorentz trajectory approach, for PRC and HRC we solely use the Lorentz trajectories. The formulation of our experiment is solely to explore the role of non-adiabaticity in the evolution of an otherwise identical ensemble of particles. It does not attempt to account for differences in plasma sheet ion properties due to different source populations, e.g. ionospheric O+ versus solar wind H+. Animations of the MHD and test particle simulations for PRC and HRC are shown in Supplementary Video S2 and S3 in the same format as Figure 1.
A comparison of the energy content at different penetration depths across the test particle ensembles of varying mass is shown in Figure 6, where the shaded region denotes the period of active test particle seeding. While the bulk energization and spatial distribution of H+ (PRC) evolve relatively similarly to positrons (GCRC), there is a qualitative change in behavior when going to heavier ions (HRC). We find that HRC is more energetic at further distances and during the period of active seeding, with the deviation between HRC and GCRC/PRC diminishing closer to Earth and rapidly in time after the source of newly-injected particles is shut off.
The overlaid text on each panel of Figure 6 is approximately the ratio of energy content at peak between HRC and PRC, which we find to be: 1.75× (14R E ), 1.5× (10R E ), 1.25× (8R E ). Recall that in this comparison we have evolved the same ensemble of test particles using varying mass, so the excess energy content of HRC is a proxy for the effects of non-adiabaticity in these particles interactions with the same flow structures. While not a precise analog, we can use the O+ enhancement factor introduced by Moebius et al. (1987) and studied as a function of downtail distance using Geotail statistics by Nosé et al. (2000). The O+ enhancement factor, E O+ R O+ /R H+ , where R X is the ratio of the differential flux os species X prior to and following substorm onset. By comparing the ratio of ratios, the enhancement factor helps control for ambient differences in the properties of different ion species, e.g., due to their different sources, and makes it a good metric for comparison with our simulation that neglects these potential differences. In their study, Nosé et al. (2000) find E O+ to be, averaging over ± 1R E : E O+ 1.77 (14R E ), E O+ 1.35 (12R E ), and E O+ 1.26 (10R E ). This is very much in line with the results of our simulation, although we do find that the comparison of HRC to PRC energy content decreases more slowly towards Earth than these statistics. We suspect this is due to differences in ambient conditions; the ratio of HRC to PRC energy content in the simulation is taken after several hours of southward IMF as compared to the Geotail statistics of local dipolarizations over 3 years. In particular it is likely that our simulation has a less stretched magnetotail configuration after several hours of SMClike activity as compared to the broader statistical sample and this in turn would affect where and how strongly O+ ions are scattered.
As can be seen in Figure 6, the significant excess energy content of HRC that penetrates inside 14R E is largely absent from the region inside 8R E , i.e. the majority of this excess energy is lost prior to reaching the inner magnetosphere. Supplementary Video S2 through Supplementary Video S4 show the evolution of GCRC, PRC, and HRC respectively and Supplementary Figure S5 shows a comparison of the state of FIGURE 6 | Comparison of total energy content for identical test particle ensembles evolved with increasing gyroradius within equatorial radius R 14 (A), R 10 (B), and R 8 (C). Shaded region denotes period where test particles are being actively injected into the nightside arc at R ≈ 18. Overlaid text is the approximate ratio of energy content at peak between HRC and PRC.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 761875 GCRC, PRC, and HRC at T 0. Most evident in these comparisons is the strong duskward skew of the highly nonadiabatic HRC. To make this more quantitative we show in Figure 7 a comparison of the MLT distribution of energy inside R 12R E at T 0 for the three ensembles. There is a broad similarity across all the ensembles: a dominant peak near midnight sharply falling off towards dawn and more gradually towards dusk, and a smaller secondary peak in the post-noon sector, approximately 1400 MLT. However, we find an appreciably stronger duskward bias in HRC as compared to GCRC/PRC. The primary pressure peak of HRC is shifted several hours of MLT duskward with a broader pre-midnight extension with the secondary pressure peak extending across the dayside into the pre-noon sector. The exaggerated dawn-dusk asymmetry due to non-adiabaticity is consistent with observational statistics of plasma sheet O+ (e.g., Ohtani et al., 2011).
The duskward bias of HRC evident in Figure 7 provides a clue to the ultimate fate of the HRC energy content inside R 14R E ( Figure 6A) that is largely missing at R 8R E ( Figure 6C). Test particles are evolved until the end of the simulation or until the test particle leaves the simulation domain, in this case a spherical annulus spanning 2-30 R E . When the test particles exit the domain, they are marked out of bounds at the location they intersected the boundary. Therefore we can identify lost energy content and distinguish its manner of exit either through precipitation, for particles lost to the inner boundary, or the dusk or dawn flank based on the sign of the SM-Y coordinate. For PRC, we find the lost energy content is distributed 20%/70%/10% to precipitation, dusk, and dawn flank respectively. For HRC, we find approximately 4× more lost energy content compared to PRC and the distribution to be approximately 10%/85%/5% to precipitation, dusk, and dawn. In other words, HRC provides a great deal of non-adiabatic energy relative to PRC which reside primarily on open drift shells that lead to significant magnetopause losses on the dusk flanks. This excess energy content is rapidly lost when the source population is extinguished.
The non-adiabatic effects we find in HRC result in a dramatic difference when compared to GCRC. Additionally, despite the fact that Figure 6 shows that GCRC and PRC have very similar overall energy content and spatial distribution, we find that there are subtle but important differences in how that energy density is distributed in momentum space, i.e. pitch angle and energy spectrum.
That there should be some degree of non-adiabaticity is perhaps unsurprising as both data (e.g., Runov et al., 2017) and particle simulations (e.g., Birn et al., 2013;Ukhorskiy et al., 2018) have shown that ions at dipolarization fronts are only quasi-adiabatic. Here we find that for the energy content inside of 8R E at the end of the simulation: the energy-weighted mean equatorial pitch angle, α, goes from ≈ 45°(GCRC) to ≈ 50°(PRC); and the energy-weighted mean energy, K, goes from ≈ 50 keV (GCRC) to ≈ 60 keV (PRC). In other words we find that PRC is hotter and exhibits some perpendicular anisotropy relative to GCRC, due to the proton nonadiabaticity. The latter is potentially important as perpendicular anisotropy provides a free energy source for electromagnetic ion cyclotron (EMIC) waves known to be important in the inner magnetosphere.
For trapped ring current ions it is generally believed that charge exchange (CX) is the dominant loss mechanism for typical magnetospheric conditions (e.g., Jordanova, 2020, and references therein). The CX lifetime can be estimated as, where < n G > is the bounce-averaged geocorona density, σ CX is the CX cross-section (Lindsay and Stebbings, 2005), and v is the particle velocity. The cross-section exhibits a great deal of sensitivity to energy with σ CX ( K PRC )/σ CX ( K GCRC ) ≈ 0.6 for the mean energies calculated for PRC and GCRC above (Lindsay and Stebbings, 2005). The other critical term of τ CX is the bounce-averaged geocorona density, which will be lower for more equatorial particles whose smaller mirror latitude limits their exposure to the highest geocorona densities.
Both the equatorial anisotropy and hotter distribution will tend to increase τ CX . Here we make a simple estimate of the combined effect using the energy-weighted mean pitch angle and energy as properties of a fiducial particle at L 8 combined with: the geocorona density at L 8 from Østgaard et al. (2003), and the approximation (Smith and Bewtra, 1976), where τ EQ is the lifetime of an equatorial particle and λ M is the actual mirror latitude. This yields τ PRC /τ GCRC ≈ 1.6, an FIGURE 7 | Comparison of the magnetic longitude distribution of total energy content, within R 12 at T 0, for identical test particle ensembles evolved with increasing gyroradius.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 761875 appreciable change in the typical particle lifetime, which has important implications for longer-term ring current decay in the aftermath of active periods. This only serves to further complicate the already difficult modeling of ring current decay due to the sensitivity to different geocorona density models (Ilie et al., 2013).

DISCUSSION AND CONCLUSION
We have presented results of tracing a large number of test particles (20 M per species) in a high resolution MHD simulation of an idealized SMC-like period (Section 2.1), where particles were initialized to "mirror" the bulk properties of the MHD plasma in the plasma sheet. We continuously seed test particles in a thin (ΔR 1R E ) nightside arc at R ≈ 18 R E and weight those test particles using the density, temperature, and flux tube volume calculated from the MHD flow at that time and location (Section 2.2). Due to the importance of the MHD flow quantities to the weighting of test particles we begin with a statistical validation exercise comparing the relative thermodynamics and magnetization of BBFs in our simulation to similar statistics derived from in situ data. We use the correlations between test particles and MHD flow properties to study the transport of particle energy through the transition region, where we use an expansive definition of the transition region to encompass the domain where both fast flows and particle drift physics are both non-negligible. These correlations generate a massive database, over a billion data points connecting particle energy transport, via the motion of the equatorial projection of a particle, and local flow properties, via the magnetic lengthscale and local buoyancy (Section 2.3). Finally, we consider the importance of nonadiabaticity by comparing the evolution, using increasing mass, of identical test particle ensembles, as defined by their initial configuration and weights. Our key results are summarized as follows: • The mesoscale bursty bulk flows (BBFs) in our global MHD simulation yield thermodynamics, magnetization, and spatial scales (Table 1 and Supplementary Figure S3) that largely reproduce statistics from in situ THEMIS measurements (Nakamura et al., 2004;Runov et al., 2015). We further demonstrated that the transportcritical properties of the modelled BBFs are principally independent of the simulation resolution, as long as the simulation possesses sufficient resolving power ( Figure 5 and Supplementary Figure S4). • Mesoscale bubbles, or localized regions of depleted entropy, are critical to particle energy transport through the transition region and, for a significant fraction of the particles, the gradient and curvature drifts are nonnegligible. • The effects of ion non-adiabaticity (Section 3.3), i.e. finite gyroradius, for heavy ions (O+) are consistent with in situ measurements of relative enhancement (Nosé et al., 2000) and create qualitative differences in the bulk energization and spatial distribution (radial and MLT) of the heavy ion energy density. The overall energization and spatial distribution of H+ is largely similar to the adiabatic control run, however differences in the momentum space distribution, i.e. energy and pitch angle, have implications for recovery timescales that are non-negligible.
Although we find that MHD is able to create BBFs that reproduce observed statistical properties, this is only at sufficiently high resolution and using a low-dissipation algorithm. The importance of proper resolution to building the ring current is likely even higher than suggested in Section 3.2, where we estimate 20% of the energy transport would be lost in the lower versus higher resolution simulation by considering the energy transport for δS < − 0.3. However, Supplementary Figure S4 adds a corollary that the lower resolution simulation manifests far less structure at scales between 0.5-3 R E , which are critical for bubbles in our simulation (Supplementary Figure S3) and typically associated with in situ measurements of flow bursts (Nakamura et al., 2004). Furthermore, the larger typical lengthscale will increase K ∇ (Eq. (7)), the characteristic energy at which energy-dependent and E × B drifts are comparable; consequently, less energy density in the transition region will undergo gradient-drift dependent processes like trapping (Gabrielse et al., 2017;Ukhorskiy et al., 2017, which may result in depressed transport. In fact, we found this to be the case in our early experiments . There we demonstrated, when comparing test particle transport in simulations with the same low and high resolutions as presented in this paper, that even though the lower resolution simulation exhibited overall higher global earthward mass flux, the test particle energy density inside of R 8 R E was only half that of the higher resolution simulation. In other words, the inclusion of adiabatic particle physics without mesoscale plasma sheet dynamics will not build the same ring current. Further, these effects are not limited to the ring current. Mesoscale processes shape the wave populations of the inner magnetosphere: anisotropic ion injections provide free energy for the EMIC wave population and low energy electron injections play a similar role for VLF wave growth (e.g., Jaynes et al., 2015). The evolving plasmaspheric population correlates with the relative distribution of whistler-mode hiss and chorus waves, with important consequences to flux enhancements of energetic trapped particles in the radiation belts (e.g., Ripoll et al., 2020).
We also have found here that adiabatic particle dynamics is insufficient to capture the building, and subsequent decay, of the ring current. For H+ non-adiabaticity plays a minor role in the bulk energization and spatial distribution of particle energy in the inner magnetosphere, but even relatively small alterations to the momentum space distribution can have important implications. Non-adiabaticity leads to a more anisotropic distribution, with perpendicular bias, which can provide free energy to the wave populations of the inner magnetosphere, and a harder energy spectrum. Taken together, the more equatorial and higher energy distribution will decay more slowly due to charge exchange (Section 3.3).
More dramatic, although not unexpected (e.g., Moebius et al., 1987;Nosé et al., 2000;Delcourt et al., 1997;Keika et al., 2013; Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 761875 Bingham et al., 2020) is the non-adiabaticity exhibited by O+. Our results show that O+ non-adiabaticity creates an excess energy density, exceeding 1.5× the adiabatic energy density for R > 10 R E (Figure 6), whose contribution to the ring current is qualitatively different than H+. For O+ versus H+ we find: shallower penetration depth, greater duskward-biased asymmetry, and a rapid depletion of the excess when the source in the plasma sheet is shut off. The inclusion of O+ in ring current composition is known to affect recovery timescales via charge exchange (e.g., Hamilton et al., 1988;Ilie et al., 2013;, due to the disparate collision cross-sections and their dependence on energy. We see here that recovery timescales via magnetopause losses will also be affected by O+ composition due to nonadiabaticity. In their review, Keika et al. (2013) discuss possible mechanisms to explain the stronger O+ energization as compared to H+. These include: direct finite gyroradius particle processes, e.g. Speiser orbits or interactions with transient and localized electric fields, which will result in preferential duskward deflection; resonant interaction with inner magnetospheric waves; and differences in the underlying source population. Of these, our simulations only include the first but produce O+ enhancement factors with similar radial dependence as observed (Nosé et al., 2000), highlighting the importance of non-adiabatic interactions with the background flow. Given the importance of O+ energy density in the ring current during disturbed geomagnetic conditions (Nosé et al., 2005), modeling these non-adiabatic effects is crucial.
In the simulation presented here we did not include ring current feedback on the global magnetosphere simulation. For the relatively undisturbed magnetospheric configuration we have considered here this omission is likely acceptable. However, expanding this type of study to more disturbed periods, such as storms, will require at least coupling to an inner magnetosphere model. Even so, our work shows that models coupling global MHD magnetosphere and ring current models, while the current state of the art for simulating typical geospace timescales, are insufficient in the transition region, where both fast flows, excluded from typical inner magnetosphere models, and energy-dependent drifts, missing from MHD, are important. The results of this work stress the importance of comprehensive modeling of ion transport across the transition region which requires: producing numerically robust mesoscale bubbles with appropriate entropy depletion; resolving energetic particle interactions with the sharp boundaries and internal magnetic structure of these mesoscale bubbles; and including the effects of non-adiabatic ion particle dynamics, particularly for heavier ion species. In other words, to properly model the buildup of the H+ ring current in the inner magnetosphere, a model must account for both fast flows and energy-dependent drifts in the transition region as the transport processes operate more effectively on the higher-energy tail of the thermal distribution. Further, for heavier ions, e.g. O+, including non-adiabatic effects is necessary as they qualitatively change the character of the contribution, in the magnitude and distribution of energy density, to the ring current.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

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

FUNDING
This work was supported by NASA grants 80NSSC19K0241, 80NSSC20K1833, 80NSSC17K0679, NNX16AG73G, NNX17AI54G, and 80NSSC19K0071; the Van Allen Probes contract NNN06AA01C and by the NASA DRIVE Science Center for Geospace Storms (CGS) under grant 80NSSC20K0601.