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.

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 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 Sorathia et al. (2018) 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 pre-existing 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 non-adiabaticity on an otherwise identical ensemble, although we do not expect sensitivity to the choice of κ (Ukhorskiy et al., 2018).

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 Sorathia et al. (2018), which used three discrete wedge-shaped injection regions on the nightside but with much finer spatial granularity.

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): S ( R , ϕ ) = ∫ p ( s ) 1 / γ d s B ( s ) , (1) 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”, δ S ( R , ϕ , t ) = S ( R , ϕ , t ) − S 0 ( R , ϕ , t ) S 0 ( R , ϕ , t ) , (2) where S 0 is an average background entropy calculated using a sliding time window, namely S 0 ( R , ϕ , t ) = 1 Δ t ∫ t − Δ t t + Δ t S ( R , ϕ , τ ) d τ , (3) 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; Sorathia et al., 2018; 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 L ∇ B = | B ⃗ | ∇ x , y B ⃗ , (4) where the denominator is taken to be the L ₂ 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 E Q , Y E Q , K ) i , p designating the equatorial projection of particle p at time T = T ₀ + 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, ( X E Q , Y E Q , K , α E Q , Δ R E Q , Δ K , L ∇ B , δ S ) , (5) 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 ² 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.

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.

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 high-resolution 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., Toffoletto et al., 2020). 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⁷ total data points, of which ≈ 1 0 6 are identified as intruding. 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.

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_pK_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.

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 ∼ 1 0 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).

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 sample-dependent membership criteria and weight. We define the contribution of a quantity, Q in the range Q ₀ + ΔQ, as the weighted-sum of the elements that satisfy both the sample membership criteria and Q ∈ [Q ₀, Q ₀ + Δ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 fine-grained discussion of the individual characteristics, we note that contribution to transport, specifically inward radial transport from the w_pKΔ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 energy-independent E × B drift, also the MHD bulk flow, increases with energy (Gabrielse et al., 2017, 2016; Ukhorskiy et al., 2017, Ukhorskiy et al., 2018). As a simple estimate we can approximate the magnitudes of the E × B and gradient drifts as V_EB ∼ E/B and V ∇ ∼ K q 1 L ∇ B 1 B , (6) respectively, and define K ∇ ∼ q E × L ∇ B , (7) the energy for which the two velocities are equal (see Sorathia et al., 2017, Figure 12). For typical values inferred from the MHD 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 Pr ( R E / 2 < L ∇ B < 3 R E ∣ δ S < − 0.2 ) = 80 % , (8) 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 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).

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 SMC-like 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 GCRC, PRC, and HRC at T = 0. Most evident in these comparisons is the strong duskward skew of the highly non-adiabatic 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 non-adiabaticity. 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, τ C X = 1 < n G > σ C X v , (9) 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 σ C X ( K ̄ P R C ) / σ C X ( K ̄ G C R C ) ≈ 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), τ C X = τ E Q ⁡ cos 3.5 ( λ M ) , (10) where τ_EQ is the lifetime of an equatorial particle and λ_M is the actual mirror latitude. This yields τ^PRC/τ^GCRC ≈ 1.6, an 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).