Early Time Evolution of Turbulence in the Space Environment by Neutral Beam Injection

The Space Measurement of A Rocket‐released Turbulence mission will demonstrate the production and evolution of turbulence in the space environment via a cascade of plasma physics processes. A sounding rocket will inject neutral barium atoms into the upper ionosphere at high speed across the magnetic field. The barium atoms photoionize, and the magnetic field confines them, leading to an ion ring distribution in velocity space that is unstable to electrostatic lower hybrid waves. Nonlinear induced scattering of lower hybrid waves produces electromagnetic magnetosonic and whistler waves that rapidly escape the source region and propagate into the radiation belts. In this paper, we present models and simulations of early time parts of this cascade until the lower hybrid waves are generated. We study the neutral atom injection, expansion, and ionization analytically and computationally. In a realistic scenario, the velocity distribution function is found to be nongyrotropic instead of a perfect ring in the plane perpendicular to the magnetic field. Using these results, we examine a representative region of the barium cloud and perform particle‐in‐cell simulations of the excitation of the lower hybrid waves. Of particular interest is the energy extracted from the kinetic energy of the barium ions, as this will ultimately affect the amplitude of the electromagnetic waves in the radiation belts. We find that for velocity distribution functions that deviate from an ideal ring distribution, the energy extracted from the fast barium ions increases.


Introduction
The SMART (Space Measurement of a Rocket-released Turbulence) mission is an upcoming sounding rocket experiment that will study the creation and evolution of plasma turbulence directly in the space environment Ganguli et al., 2019). A sounding rocket will be launched to approximately 600 km and inject (via shaped charge detonation) approximately 1.5 kg of barium (Ba) neutrals at approximately 8 km/s across the Earth's magnetic field. Barium neutrals will photoionize with a characteristic time scale of around 30 s, forming a plasma. As it ionizes, the streaming barium will form an approximate ring distribution of barium ions in velocity space that is unstable. The growth and subsequent saturation of the ion ring instability will excite lower hybrid waves , which are short wavelength and electrostatic (kc∕ p ≫ 1, where k is the wave number, c is the speed of light, and p is the plasma frequency). For large enough wave amplitudes, induced nonlinear scattering will occur and excite whistlers and magnetosonic waves, which are long wavelength and electromagnetic (kc∕ p ≲ 1) Ganguli et al., 2015;Tejero et al., 2015). These electromagnetic waves will propagate away from the source region into the magnetosphere (Crabtree et al., 2012). Figure 1 shows a diagram of the SMART mission.
The SMART experiment will trigger a cascade of plasma physics processes, and the different stages of this cascade have very different temporal and spatial scales. Table 1 lists these stages, their accompanying spatial and temporal scales, and appropriate modeling or simulation methods. The bold entries in Table 1 are the early time stages we address in this paper. The multiscale nature of the experiment means that modeling and prediction of the results are challenging endeavors. In fact, a successful prediction of this complicated process would have implications on the ability to predict space weather using physics based tools, and this is the primary motivation for the experiment. A practical way to make this problem tractable is to break up the stages and approach each one separately, using the results of the previous stage to define parameters for the next.
The barium starts as a solid, and to simulate the detonation of the shaped charge, we must take into account solid dynamics, phase change, and chemical reactions. A hydrocode such as CTH (McGlaun et al., 1990) is an appropriate tool to simulate this stage; this task is undertaken by our collaborators in the Naval Surface Warfare Center and will be reported elsewhere. The total mass of barium injected across the magnetic field at high speed is of critical importance. The barium neutrals become collisionless and expand over a few hundred kilometers, requiring a different method of simulation. Fluid and test particle simulations of the expansion will provide the barium density as well as the velocity distribution function but will not contain any of the microscopic plasma physics responsible for the waves.
To model the instability and nonlinear scattering, a kinetic plasma technique such as particle-in-cell (PIC) is necessary. A 3-D electromagnetic PIC simulation will provide the efficiency of energy extraction (via induced scattering) from the kinetic energy of the barium to the electromagnetic waves that propagate into the magnetosphere. Finally, it is necessary to track the propagation of the electromagnetic waves out of the source region to ensure that they will be detectable by a satellite in the radiation belts. A wave kinetic model is an appropriate tool for this stage. In this paper, we are focused on modeling and simulating the expansion and ionization of the barium (section 2), and the lower hybrid waves excited by the ion ring instability and the amount of energy lost by the barium ions to these waves (section 3). A 3-D electromagnetic simulation model will be necessary for addressing the remaining processes and will be discussed in a forthcoming article.

Neutral Injection
After the sounding rocket reaches the appropriate altitude (500-600 km), the shaped charge detonation will inject the neutral barium beam perpendicular to the magnetic field. The expansion and ionization of the barium cloud will take less than a minute. In this section, we model this process up until the point that the microscopic plasma physics of interest begins, namely, the generation of lower hybrid waves. We first present a fluid simulation of the release on a macroscopic scale that demonstrates that skidding, that is, creation of an electric field due to the release and ionization of barium, is not important and that structuring (fluid type instabilities) does not modify the overall structure of the cloud on either early or late time scales. Having demonstrated these facts, we are free to perform analytic and test particle calculations of the motion of the neutrals and ions.

Fluid Simulation of Injection and Expansion
We use the Naval Research Laboratory ionosphere code SAMI3/ESF (Huba et al., 2008) to model the 3-D dynamics of an ionized barium cloud in the Earth's ionosphere. The model is similar to that used to study the CRRES G-9 barium release (Huba et al., 1992). The plasma equations are modified to include the creation and injection of barium ions transverse to the geomagnetic field. First, a production term is included in the ion continuity equation for barium ions:  Ba + = Ba n Ba , where Ba = 0.0357 sec −1 is photoionization rate of barium and n Ba is the density of barium neutral atoms. Second, a source term is added to the current  given by (B∕c)Σ Ba + V Ba × ê z , where V Ba is the injection velocity of the barium atoms, Σ Ba + is the Pedersen conductance of the barium ions, and B = B ê z . The evolution of the barium neutral atoms is given by where N 0 is the number of barium atoms released, v th = 1.5 km/s is the thermal velocity of barium, = 0.05 is the initial fraction of barium ionized at t = 0, and x 0 , 0 , and z 0 is the initial position of the barium cloud. We compare and contrast two barium cloud releases. The first uses N 0 = 7.5 × 10 24 cm −3 , which is relevant to the upcoming SMART release, and the second, a very large release, uses N 0 = 1.0 × 10 27 cm −3 .
The geophysical parameters used in the simulation are for quiet time, low solar activity conditions: F10.7 = 90, F10.7A = 90, and Ap = 4. The release occurs in the ionosphere near Wallops Island, VA: latitude 37 • N, longitude 75 • W, and altitude 600 km. The day-of-year is 263 (mid-September), and the release time is 05:41 LT (i.e., predawn). The width of the longitudinal grid is 8 • . The grid resolution is n z = 401 (number of grid cells along the geomagnetic field), n = 130 (number of "magnetic field lines"), and n l = 192 (number of grid cells in longitude).
The evolution of the barium ion clouds is shown in Figure 2 as a function of distance (x) and altitude ( ). The panels on the left are for the SMART release, and on the right for the large release. The distance is in the zonal direction (i.e., longitude) and is transverse to the geomagnetic field. Color-coded contour plots of the barium ion density (log 10 cm −3 ) are shown at four times: t = 0, 20, 40, 60 s from top to bottom. At t = 0 s the barium ion clouds are initiated with a peak density at x ∼ 90 km; the barium neutral clouds are directed in the +x direction with a velocity V Ba = 10 km/s. As time progresses (from top to bottom) the barium neutral clouds ionize and the barium ions become elongated in the x direction. However, the peak barium ion density at t = 20 s is at x ∼ 108 km and it remains at this distance for the remainder of the simulation for both releases. Thus, within a photoionization time (≲28 s) the peak density of the barium ions occurs within 20 km of the release point and subsequently remains at this location, that is, the bulk of the barium ion cloud does not "skid" across the magnetic field lines. The reason for this is because the electrostatic electric field generated by the injected barium ions extends along the entire magnetic flux tube and acts to move the ionospheric plasma on this flux tube. This plasma acts as a "drag" on the barium ions because Σ i >> Σ Ba + (where Σ i is the ionospheric Pedersen conductance) and the barium ions become "frozen" to the magnetic field (Huba et al., 1992).
This point is amplified in Figure 3, which is a line plot of the barium ion density as a function of distance (x) and time (0-30 s) at the release altitude 600 km. The left panel is for the SMART release and the right panel is for the large release. At time t =5 s the peak barium ion density is at x ∼ 100 km and is associated with the transverse motion of the neutral barium atoms. However, by t = 10 s the peak ion density is at x ∼ 108 km and it remains at this location for the remainder of the simulation. However, not all of the barium atoms are ionized at t = 10 s but continue to flow in the x direction; they subsequently become ionized, and the barium ions extend to ≳400 km by 30 s (albeit at a much lower density than the peak density). The only significant difference between the two releases is that the peak density of the SMART release does not decrease very much in the first 30 s, while the peak density of the large release decreases by ∼20% after 30 s. This is because of the enhanced recombination rate associated with the higher electron density in the large release case.
Lastly, we show Figure 4, which is a color-coded contour plot of the barium ion densities as a function of latitude and altitude at times 05:41 and 05:43 LT at the release longitude. The solid lines denote the magnetic field lines. At 05:41 LT (20 s after the release) the ion clouds are roughly circular because of the uniform expansion of the barium atoms in this plane. However, at time 05:43 LT the barium ion clouds are elongated along the geomagnetic field lines because the ions can freely move parallel to B. The extent of the ion cloud along the magnetic field is larger for the large barium release because of its initial size.

Kinetic Theory of Injection and Expansion
In this section we describe the analytical solution to the evolution of the neutral and ionized barium particles, which we compare to the fully numerical solution. As seen in Koons and Pongratz (1981), the barium can be modeled as two populations: a fast beam with nonzero bulk velocity, and a slow component with zero bulk velocity that diffuses away from the release point. After a short time after the detonation the neutral barium ions have a distribution function N , and C = Δv 2 || (Δv ⟂ + Δv r ) 3∕2 , and where N is the number of barium atoms injected, is the fraction of barium atoms in the fast population, and v r is the beam velocity. The number of barium atoms in the fast distribution is and the number in the slow part of the distribution is We must also assume some spatial distribution, which we choose for convenience the following form: We assume that there are no forces acting on the neutrals over the time scale of interest; thus, the initial velocities are all constant v = v 0 and the position can be simply updated as x = x 0 +v 0 t. With these equations the evolution of neutral density may be written as where Γ ion ≈ 1∕(30 s) is the ionization frequency. This ordinary differential equation can be solved in general as N = exp(−Γ ion t) N , and then transformed back along the trajectories of neutral particles as Once the barium neutrals are photoionized, they can be kept track of using a distribution function, Ba+ , whose evolution equation is given by where the neutrals are now the source of ions and is given by equation (8). To solve this equation, we first assume that the motion of the barium ions is dictated only by the Lorentz force (we neglect any skidding electric field as the fluid simulations show is reasonable). Then the motion of the ions is described by and Then to solve equation (9), we use the evolution equations (10) and (11) to go backward in time and find the point in time that maximizes the right-hand side of equation (9). This time corresponds to the point that the particle with phase space coordinates (t, x, v) was ionized. Because the function h describing the initial density of particles is so highly localized in space, we may then expand about the maximum to arrive at a Gaussian integral and perform the integral over time analytically. We note that the right-hand side has multiple peaks which are easy to find numerically, so rather than finding the absolute maximum we simply perform the integration over all peaks. We have compared the analytical results with both numerical integration of equation (9) and with the test particle calculations and find excellent agreement between all of the methods. The analytical solution is most appropriate for investigating the fine-scale structure of the distribution function.
In Figure 5 we plot the distribution function of barium ions at 9.4 km from the release point and 2 s after the release as computed from the analytical theory. This figure can be understood as follows. Barium neutral atoms that are traveling fast arrive near the observation point first and ionize. Once ionized they begin to gyrate as slower ions are being ionized creating a continuous spiral function.
In Figure 6 we plot the distribution function at 100 km from the release point and 11 s after the release. In this case we have a very tight spiral. We also plot the distribution function as a function of v x > 5 km/s only at v = 0. One sees that each individual peak in the spiral is very cold. The Gaussian character of the initial distribution function can be discerned in this plot by comparing the height of the peaks. At this instant of time the slower particles have not yet arrived. In comparing Figures 5 and 6 one should keep in mind that the effects of velocity dispersion, that is, slower particles arriving later, are more pronounced farther from the release point; this accounts for the tightness in the spiral at 100 km as compared to 9.4 km.

Test Particle Simulation of Injection and Expansion
The fluid simulation in section 2.1 is useful for understanding any possible skidding or structuring of the cloud as well as macroscopic interactions with a realistic background ionosphere. The analysis in section 2.2 provided barium distribution functions (and thus the densities) as a function of time and space. In this section we solve the same problem computationally; the advantage of this is the inclusion of other physics such as collisions, metastable excited states of barium, chemistry, and structured electric fields, though not all of these are included in the simulations presented here.
A possible computational approach is direct simulation Monte Carlo (DSMC) (Bird, 1994), which uses computational particles to represents a large number of real atoms. DSMC is useful for rarefied gasses in which  the Knudsen number is large (Kn > 1) and the traditional Navier-Stokes equations do not apply. The Knudsen number is Kn = ∕L, where is the mean free path and L is a characteristic length scale. The barium cloud has a large Knudsen number, though it is important to note that for plasmas, as opposed to neutral fluids, a fluid simulation such as that used in section 2.1 can still adequately model macroscopic flows with a large Kn.
In DSMC, collisions randomly occur with nearby particles with an appropriate collision rate defined by the local properties of the fluid. The most relevant collisions for the barium cloud are collisions between (and among) barium atoms, the background oxygen (neutrals), and electrons, as well as photoionization of the barium atoms. Except for very early times in the release, the elastic collision time scale is longer than the photoionization time, so we neglect them and only consider random photoionization of the barium with a characteristic time scale of 30 s (Carlsten, 1975;Hoch & Hallinan, 1993); this also considerably simplifies the computation. In a model of CRRES G1 release, the beam velocity slowed due to interactions with the background by ∼15% in the first 30 s, and that release was considerably lower in altitude than that of SMART (Rairden et al., 1994). Therefore, the simulation presented in this section is not DMSC but is more appropriately called a test particle simulation. Barium atoms expand and ionize, after which they are trapped in cyclotron motion by the magnetic field.
The simulation provides the densities and distribution functions in space and time. We use this not only to understand how the cloud evolves in space but also to predict what the instrument payload will measure. Trajectories for the barium injection module (i.e., the mother) and instrument payload (i.e., the daughter) are estimated assuming free fall. Vertical acceleration is given by a(z) = −g 0 (Re∕(Re + z)) 2 , where g 0 = 9.832 m/s 2 is the acceleration due to gravity at the Earth's surface, Re = 6.378 · 10 6 m is the Earth's radius, and z is the altitude above the Earth's surface. In this case, the mother has an apogee of 600 km, a constant horizontal velocity of 400 m/s, and gravitational acceleration in the −z direction. The daughter trajectory was simulated by adding a separation velocity of 4.3 m/s at 130 s into the trajectory (at 344 km altitude) and along the direction of the velocity vector. The horizontal velocity, separation velocity, and timing for separation were taken from a previous experiment, CARE II, that launched on a similar sounding rocket (Bernhardt et al., 2017) (National Aeronautics and Space Administration designation 39.012). This resulted in the daughter being approximately 1,073 m above and 201 m downrange from the mother when the barium is injected at apogee. Figure 7 shows this trajectory and the separation between mother and daughter.
We initialize ∼17 billion computational particles according to the distribution similar to that observed by Koons and Pongratz (1981). At the beginning of the simulation, there are two populations of neutral barium atoms: a fast component with a large drift velocity and a slow component with no drift velocity. The fast component contains 11.8% of the total barium and streams at a velocity of 9.4 km/s. Each population is initially represented by a Maxwellian distribution with the following thermal velocities: v t⟂, ast = 1.717 km/s, v t||, ast = 0.8854 km/s, v t⟂,slow = 0.8567 km/s, and v t||,slow = 0.8854 km/s. Figure 8 shows the density of the barium ions 20 s after release at an altitude of 600 km. At this point in time, 48.66% of the barium has photoionized. Figure 8 shows cuts along the central planes of the barium cloud. The barium ions are free to stream parallel to the magnetic field, which results in considerably more expansion in that direction as compared to the perpendicular plane in the lower plot. In the upper plot, there are striations very close to the release point with a length scale of barium gyroradius. This purely kinetic effect, called cycloid bunching, was also seen in the CRRES releases (Bernhardt et al., 1993). Figure 9 shows the density of the barium neutrals. Figures 8 and 9 also show the trajectory and current location of the instrument payload. One can see that the payload cannot achieve much separation from the point of release relative to the scale size of the entire cloud. Therefore, there are two regions of interest from a modeling standpoint: the region near the payload for predicting and understanding data from the instruments, and a region representative of the cloud as a whole for predicting and understanding the generation of electromagnetic waves. This also shows the importance of releasing the instrument payload before the barium is injected. The fast time scales of the relevant plasma physics (see Table 1) mean that the payload must be in position before the injected barium arrives, or else the instruments would only measure a saturated state. Figure 10 is a plot of the plasma properties along the trajectory of the instrument payload. The striations near the injection point from Figure 8 are seen in the fast barium density. Unlike the vast majority of the cloud, the plasma that the payload will fly through has a higher or comparable density to the background oxygen ions for much of the flight time. Additionally, in the region close to the release point, the slow component of the barium (which has zero bulk velocity) has a higher density than the fast component. The payload will therefore see not just a ring distribution from the streaming barium, but a secondary ring with a smaller v ⟂ . The presence of the slow barium components will not change the real part of the frequency significantly; however, there will be excitation of waves with shorter wavelengths. These effects have not yet been included in the analysis of the plasma physics in section 3, primarily because the plasma physics of interest occurs within ∼1 s of the arrival of the fast ions; thus for most of the cloud (x ≿ 7km) this holds. Figure 11 shows the distribution function of barium ions in the perpendicular plane at a distance of 9.4 km away from the release point along the central axis of the cloud 2 s after release. The distribution function that forms is actually not a perfect ring but a spiral. The spiral forms because the injected beam has some thermal spread around the beam velocity. As the injected barium arrives at this point in space, a smaller number of particles traveling faster than the beam velocity form the outer part of the spiral, then a larger number of particles traveling near the beam velocity form the center of the spiral, and finally, a smaller number of particles traveling slower than the beam velocity form the inner part of the spiral. These particles execute cyclotron motion around a magnetic field line, causing the spiral distribution to appear to spin when sampled at a single point in space. This shape is still unstable and is expected to broaden the frequency range of the generated lower hybrid waves relative to an ideal cold ring distribution because there are a range of ring velocities as a function of gyrophase. Away from the central axis, the distribution is not just a spiral but a helix that wraps around the parallel axis. Note that this distribution matches very well with that predicted the analytic theory; Figure 11 can be directly compared with Figure 5. The theory is more practical for determining distribution functions at arbitrary points, whereas the simulations are more practical for finding the number densities since the moment integrals in the theory must be done numerically.

Ion Ring Instability
While section 2 models the injection and expansion of the barium cloud, it does not include any of the plasma wave generation and scattering that will be triggered by the unstable velocity distribution. In this section, we look at the stability and early time nonlinear development of an ion ring instability with theory and kinetic simulation. The sole purpose of this section is to understand the initial growth phase of the instability and to estimate the wave amplitude and amount of energy extracted from the beam in this initial phase.

Linear Theory of Partial Ring Instability
The background ions and electrons are treated as a cold fluid, since ∕k ≫ v te,i . In the electrostatic approximation, the fluid density fluctuation is given as where s is a species label and we allow for an arbitrary multicomponent cold plasma background.
We treat the beam ions with the Vlasov equation, where B is an ordering parameter to indicate that on the lower hybrid wave time scales the beam ions are essentially unmagnetized; however, on the quasi-linear evolution time scale this term could be retained to compare with magnetized PIC simulations. In the linear approximation, this leads to Using these expressions in Poisson's equation leads to the expression, are the usual Stix parameters, 2 pb = 4 e 2 n b ∕M b is the beam plasma frequency, b0 = n b0 F b , and n b0 is the heavy ion beam number density. Finally, we are interested in quasi-static beam distributions that are ordered by the magnetic field, and to connect to the traditional ring  (20)) in black, and three solutions for the partial ring (equation (21)) for three different wave vector directions that correspond to the directions in panel (a). instability, we transform to cylindrical coordinates in velocity space and real space.
where is the gyrophase and is the angle the wave vector makes with respect to the x axis. The dispersion relation can then be written as (18) where we have also set k || = 0, which maximizes the growth rate and is consistent with the two-dimensional simulations below and we have integrated over the parallel velocities. Equation (18) is the general dispersion relation of interest to this paper. To illustrate the linear instability properties with a realistic but concrete model, let us assume which describes a cold beam whose perpendicular velocity v b ( ) and density h( ) depends on the gyrophase. Assuming a full cold ring, v b = v r and h = 1∕(2 ) become constants, and in the limit of a single background ion species and pe ≫ Ω e , the dispersion relation becomes  1 where = (n r ∕n 0 )(m i ∕m r ). With the distribution function in equation (19) the dispersion relation becomes where the integral over the gyrophase can be done numerically. In the limit of the full ring equation (21) reduces to equation (20).
In Figure 12a we show the distribution function of fast barium ions 9.4 km away from the release point and 0.7 s after the release. At this point in space and time only the very fast ions have reached this point, and the faster of these ions have begun to gyrate around the magnetic field. In this figure we also draw three directions from the origin to the beam of ions corresponding to = 297 • (blue) , 332 • (orange), and 343 • (green). In panel (b) we plot the solution to the dispersion relation (equation (21)) for = 0.1 corresponding to the distribution function in panel (a) for the three different wave vector directions as a function of the wave vector magnitude. The h( ) function for this case peaks very near to = 332 • , which is where we find the maximum growth rate. At = 297 • the maximum growth rate is found at a much larger wave vector, this peak occurs roughly at a k ⟂ that is 1∕ cos( max − ) times larger than the peak at which h( ) has its maximum max . In black, we plot the solution to the full ring dispersion relation, equation (20) for = 0.1. For the full ring there is no dependence, meaning that all wave vector directions have the same growth rate, whereas for the partial ring the gyrosymmetry has been broken. Interestingly, the growth rate for the partial ring is larger by almost a factor of 2 in this case despite the beam having the same density and average velocity. This is a consequence of the gyrobunching and broken symmetry.

Kinetic Simulation of Partial Ring Instability
We use a PIC approach to simulate the growth of the ion ring instability (Birdsall & Langdon, 2004). The PIC code is 2-D electrostatic and uses a volume preserving particle mover (He et al., 2015), a high-order Figure 13. The initial distribution functions for the partial rings with respect to the gyrophase. The full ring is gyrotropic and the rest are nongyrotropic. The fraction of a ring (e.g., one third) refers to the full width at half maximum of the Gaussian. finite difference solution to Poisson's equation (Sutmann, 2007), and cubic splines for interpolation. It is impossible with current computers to simulate the entire barium cloud (which would require ∼ 10 28 cells), so we simulate a small periodic region with parameters appropriate for the rough center of the cloud in section (2.3).
It must be emphasized that while the 2-D electrostatic model used here is appropriate for assessing the instability onset and the consequent loss of beam energy that sustains the waves in the initial phase, it is not appropriate for the subsequent dynamics that in reality includes nonlinear scattering, conversion of electrostatic lower hybrid waves into electromagnetic whistlers/magnetosonic waves, and their propagation out of the source region, which are critical components of the SMART experiment. As discussed in Ganguli et al. (2010), Ganguli et al. (2015), the radiation of electromagnetic energy from the source region leads to saturation of the instability when NL (W∕nkT) = LH has been achieved as the wave amplitude grows to make W∕nkT sufficiently large. Here NL is the nonlinear scattering rate, LH is the lower hybrid wave growth rate, and W is the wave energy density. As the lower hybrid waves are generated in the ionosphere and their amplitude increases proportional to its growth rate, the nonlinear scattering rate also increases. When a critical value of the lower hybrid wave amplitude is reached, the nonlinear scattering rate becomes sufficiently large to convert all the lower hybrid waves generated into electromagnetic waves that propagate out of the source region. At this point, the rate at which the lower hybrid waves are generated equals the rate at which they are scattered into electromagnetic waves and a steady state is established, which saturates the lower hybrid wave amplitude.
Beyond this time the realistic physics is not captured in the current 2-D electrostatic simulation. Also, the nonlinear saturation occurs at a much lower amplitude compared to the quasi-linear case and results in a considerably lower heating of either the background or the beam particles since the energy escape rate is faster than the heating rate, that is, P > ei > ih , where P is the rate that whistlers propagate away from the source region, ei is the rate of electron-ion collisional damping, and ih is the ion Landau damping rate, which leads to plasma heating due to interactions between the lower hybrid waves and the oxygen plasma. Thus, in the realistic case of the SMART experiment a hierarchy of time scales exists that triggers a chain of physical processes, which are highly disparate in time and space scales and is given by It is impossible to reproduce the entire self-consistent chain of events faithfully in a simulation because of their highly disparate space and time scales. Hence, in this article we choose to understand the early time physics, which is essentially the neutral barium injection, photoionization of the neutrals to form of the ring distribution, and the consequent generation of the electrostatic lower hybrid waves. For these early time dynamics a 2-D electrostatic simulation model is sufficient. The duration of validity of this simulation to the realistic SMART experiment is until the time t = t c when W becomes sufficiently large so that the condition NL = LH is achieved. Beyond t = t c , a 3-D electromagnetic PIC model with radiative boundary condition or a wave kinetic model including parameterized nonlinear scatterings will be necessary, which will be addressed in a forthcoming article.
Since the lower hybrid wave growth time due to the ion ring instability is much smaller than the barium cyclotron period, we also expect the instability to evolve long before a barium cyclotron orbit. This means that a full cold ring, i.e. one that is gyrotropic and with a delta function distribution at the ring velocity, will never form in the actual SMART experiment. Therefore, we are also interested in how the instability changes when the ring is only partially formed (i.e., nongyrotropic). In this section, we present two types of simulations: a set of partial rings and a set of rings that form naturally via photoionization of neutrals (i.e., the barium ion density is time dependent). These cases are a more realistic representation of the distribution that will actually be generated as shown in Figure 11 than the assumption of a fully formed ring. For the partial rings, we examine five cases: a full ring, a third of a ring, a sixth of a ring, a ninth of a ring, and a twelfth of a ring. The partial rings have a Gaussian shape (at v ⟂ = v r ) with respect to the gyrophase; this is shown in Figure 13. For the photoionization simulations, we examine three cases with different photoionization times: fast photoionization pi = 3 2 cr ≈ 65, 000, medium photoionization with pi = 6 2 cr ≈ 130, 000, and slow photoionization with pi = 9 2 cr ≈ 195, 000, where cr is the ring ion cyclotron frequency and pi is the photoionization time (i.e., n r ∝ 1−exp(−t∕ pi ), where n r is the ring ion number density). Figure 14 shows the ring ion densities relative to total ion density for each photoionization simulation.
We have used parameters similar to those in Winske and Daughton (2012) for the full ring, one third ring, and one sixth ring. There are three species: electrons, background ions, and ring ions. We do not differentiate between electrons originating from photoionization and background electrons. The mass ratio is m i ∕m e = 500 and the ring ion to background ion mass ratio is m r ∕m i = 4. The charge-to-mass ratio for the electrons is −1. The fraction of ring ions to total ions is 0.3. The simulation is run in the perpendicular plane with normalized B = 1∕ √ 3. The electron plasma frequency is p,e = 1 and the electron thermal speed is v 2 t,e = 0.001. The ring speed is v r = 7v t,i and the ring is entirely cold (i.e., no spread The box length is 30 on each side for the partial ring simulations and 10 on each side for the photoionization simulations. The partial ring simulations are run until t = 15, 000, which is equivalent to about 50 lower hybrid time periods and 0.7 ring ion cyclotron periods. The photoionization simulations are run until t = 2 pi , which is approximately t = 130, 000, t = 260, 000, and t = 390, 000 respectively. The total number of particles is ∼4 billion for the partial ring simulations and ∼1 billion for the photoionization simulations. In the partial ring simulations, the ring ions are initialized to a ring distribution modified as shown in Figure 13. In the photoionization simulations, a cold beam of neutral barium particles initially streams across the periodic box with velocity v r . These neutrals will randomly photoionize, similarly to the test particle simulation in section 2.3, according to the characteristic photoionization time scale that is much longer than the other time scales of the system. Beginning with the partial rings, Figure 15 shows the ring distribution function and the electrostatic potential at the time in the simulation when the wave amplitude is largest. As the rings become more nongryotropic, the waves develop earlier and with higher amplitude. This is consistent with an increase in instability growth rate that is predicted by equation (21). For the full ring, there is no preferential direction for the wave vector, whereas for the partial rings, there is a preferential direction for the wave vector; again this is consistent with the theory in section 3.1. As the waves develop, the cold ring diffuses in velocity space until the instability saturates. As mentioned earlier, in the actual experiment, the waves will scatter into whistler/magnetosonic waves  and escape the source region, carrying turbulence energy with it. In this 2-D electrostatic model, the saturation is reached by quasi-linear relaxation, which allows large wave amplitudes that subsequently cause large diffusion of the ring particles in velocity space. In a realistic 3-D electromagnetic simulation, the saturation will occur with NL ∼ LH , which occurs much earlier than quasi-linear saturation and hence the wave amplitude will be smaller, resulting in less diffusion. Figure 16 shows the spectrum of the waves over the entire simulation for the full ring case integrated over the perpendicular plane, as well as theoretical predictions for the real and imaginary parts of the frequency. There are two relevant modes: a beam mode and a lower hybrid mode. These can be seen clearly in both the simulation spectrum and the black lines. We solve equation (20) numerically to find the predicted dispersion. The maximum growth is expected around k ⟂ v r ∕ LH = 1, which agrees with the simulation power spectrum. The lower hybrid mode has no imaginary component in the linear theory for a cold beam; a positive growth rate for this mode develops as the ring goes unstable and diffuses in velocity space. Figure 17 shows the ring ion distribution function at six snapshots in time for the slow photoionization case. There is a constant source of new particles at (v x , v ) = (v r , 0) as photoionization occurs. The most Figure 15. The barium ion distribution function (left column) and electrostatic potential (right column) for the partial ring cases at the point in time at which the wave amplitude is at its peak. Observe the increased directionality of the waves as the rings become more nongyrotropic.  (20)) is shown for the two modes in black. The theoretically predicted dispersion matches the simulation very well. Growth rate for the beam mode from the imaginary part of the solution to the dispersion relation (equation (20)) (bottom). The lower hybrid mode has no linear growth rate until the instability grows beyond the linear stage and the cold ring diffuses in velocity space. unstable part of this forming ring is near the source region. One can see in the snapshot at t LH = 150 that the instability has grown (due to the diffusion of the ring distribution function in velocity space) well before a complete ring ion cyclotron period. At this point, only 3.3% of the ring particles have photoionized. By t LH = 6000, 75% of the ring particles have photoionized. Figure 18 shows the energy in the electric field fluctuations over the course of the partial ring simulations. This field energy is the total electric field energy contained in the simulation box normalized to the total energy (field + kinetic) of the system. This confirms the observation from Figure 15 that the growth rate increases as the partial rings become more nongyrotropic. Figures 19 show the energy in the electric field fluctuations over the course of the photoionization simulations. One pane shows the whole simulation time, while the other focuses on the initial growth of the instability. Unlike the partial rings, the photoionization simulations do not have a classic growth followed by saturation and decay. New ring ions are continually being introduced, which excites waves over a much longer time period by continually driving the distribution function toward a more unstable state. Despite this, the highest amplitude waves do appear near the beginning of the simulation. The gradual decay of wave energy is the same despite differences in characteristic photoionization times.
Finally, Figure 20 shows how much of the initial kinetic energy is extracted from the ring ions. This energy flows both into the electrostatic waves and into heating the background electrons and ions. Since these are 2-D electrostatic simulations, the energy extraction predicted here will not represent what happens in the experiment with full accuracy. A third dimension adds a fast nonlinear scattering of electrostatic waves Figure 18. The energy in the electric field for the partial ring simulations. The growth rate of the instability increases as the rings become more partial. and an electromagnetic simulation will allow nonlinear scattering to whistlers/magnetosonic waves. Note that in a periodic simulation of any dimension, the whistler/magnetosonic waves will remain in the system (and interact with the plasma) instead of escaping the source region as they would in the unbounded SMART experiment. This can affect the threshold condition and hamper subsequent wave excitation, which could lower the energy extraction.
An important impact of 3-D electromagnetic simulation is that the instability will no longer saturate by quasi-linear broadening of the ring but will saturate when LH ≈ NL . This can increase the amount of energy extracted from the beam. The efficiency of the energy extraction is of paramount importance to the experiment as it determines the amplitude of the resulting whistler waves out in the radiation belts. Figure 20 demonstrates that as the ring spans a smaller angle in velocity space, the energy extraction increases significantly. Linear theory, quasi-linear theory, and other simulations predict an extraction efficiency between 2% and 10% (Kulygin et al., 1971;Scales et al., 2012;Winske & Daughton, 2012), but these have all assumed a full ring. In the photoionization case, the energy is continuously extracted as the waves are driven over the entire Figure 20. The loss of kinetic energy from the barium. The partial ring cases are shown in the left and the photoionization cases on the right; note the difference in time scale. As the ring becomes more nongyrotropic, more energy is lost to the lower hybrid waves. The photoionization case shows gradual energy extraction as more ring ions are continually created.