Radial diffusion simulations of the 20 September 2007 radiation belt dropout

This is a study of a dropout of radiation belt electrons, associated with an isolated solar wind density pulse on 20 September 2007, as seen by the solid-state telescopes (SST) detectors on THEMIS (Time History of Events and Macroscale Interactions during Substorms). Omnidirectional fluxes were converted to phase space density at constant invariantsM = 700 MeV G−1 andK = 0.014RE G1/2, with the assumption of local pitch angle α ≈ 80 and using the T04 magnetic field model. The last closed drift shell, which was calculated throughout the time interval, never came within the simulation outer boundary of L = 6. It is found, using several different models for diffusion rates, that radial diffusion alone only allows the data-driven, time-dependent boundary values at Lmax = 6 andLmin = 3.7 to propagate a few tenths of anRE during the simulation; far too slow to account for the dropout observed over the broad range of L = 4–5.5. Pitch angle diffusion via resonant interactions with several types of waves (chorus, electromagnetic ion cyclotron waves, and plasmaspheric and plume hiss) also seems problematic, for several reasons which are discussed.


Introduction
Rapid, global dropouts in the radiation belts are currently of much interest, as it has become evident that they are not well understood.Somewhat parallel to the study of acceleration events, there are two dominant paradigms, both plausible: local precipitation due to pitch angle scattering by wave-particle interactions, and outward radial transport (diffusion) potentially combined with magnetopause shadowing.
The former was invoked by earlier studies of microbursts, but more recent work, based on global observations as well as unsuccessful attempts to detect precipitation, emphasizes the latter.A thorough, up-to-date review of the situation is given by Turner et al. (2012b).
This study of the 20 September 2007 dropout is based on measurements by the solid-state telescopes (SST) detectors on the five THEMIS (Time History of Events and Macroscale Interactions during Substorms) spacecraft (Angelopoulos, 2008).Because pitch angle-resolved measurements are not available for this early period in the THEMIS mission (D.Turner, personal communication, 2013), this investigation is restricted to the radial diffusion mechanism, both in its pure form and as often augmented with a simple, semi-empirical local loss term (e.g., Shprits et al., 2005).

Observations
Dropouts have been linked to sudden increases in solar wind velocity, density, or dynamic pressure P dyn (Morley et al., 2010;Shin and Lee, 2013).As seen in Fig. 1, during 20 September (263 day of year), the solar wind velocity rose moderately but the density spiked, leading to similar behavior in the dynamic pressure.This was accompanied by a moderate increase in Kp and, at its peak, a sharp increase in AE.The Dst index rose to about 38 and then decreased to about −18.In this study, the main use of these large-scale solar wind and magnetospheric quantities is in evaluating the magnetic field model and radial diffusion coefficients.Top to bottom: Dst, Kp, solar wind velocity, solar wind density, solar wind dynamic pressure, IMF B y , IMF B z , and the integrated quantities W 1 -W 6 (magenta, red, blue, cyan, green, and black, respectively) used the T04 magnetic field model.The final panel shows the last closed drift shell computed according to the T04 magnetic field model (black) and the T96 model (red).Also shown in the final panel is the plasmapause location according to several models as described in Sect.2.1.

Global setting
The underlying magnetic field model is used to calculate adiabatic invariants associated with particle flux measurements.THEMIS ephemeris data is provided in terms of the McIlwain L parameter, calculated using the T96 (Tsyganenko, 1996) model (and International Geomagnetic Reference Field -IGRF), but here the (Roederer) L * is used, computed using the T04 (Tsyganenko and Sitnov, 2005) model as implemented by the ONERA-DESP IRBEM library (Boscher et al., 2013).These field models are driven by Dst, P dyn , and IMF B y and B z ; T04 also uses the derived, timeintegrated quantities W 1 -W 6 , all shown in Fig. 1 (Qin et al., 2007).Of particular interest is the L * value of the so-called last closed drift shell (LCDS), beyond which particles are not on closed drift shells but are assumed lost to the magnetopause on a drift timescale (referred to as magnetopause shadowing).It is worth noting that the T96 and T04 models include the Shue et al. (1998) estimate of the magnetopause location (in r), which typically lies outside the LCDS.Actually, the LCDS is not unique but depends on the particle pitch angle.In fact, the situation is made much more complex by the possibility of drift-orbit bifurcation (DOB), which changes the second invariant (Öztürk and Wolf, 2007) and invalidates most calculations of L * but does not cause the same rapid loss of confinement.Since low values of equatorial pitch angle are less subject to DOB, as an expedient the LCDS is found (by iterative search) for a particle with equatorial pitch angle of 40 • at midnight.The bottom panel of Fig. 1 shows the LCDS at one hour intervals, according to the T04 model (black curve) and also according to T96 (red curve).The LCDS for this time interval never comes within L * = 6, so it should not be a direct factor in simulations with an outer boundary at L * lower than this.
Also shown in the bottom panel of Fig. 1 is the plasmapause location, based on sequences of Kp according to estimates of Carpenter and Anderson (1992) (blue curve) and O'Brien and Moldwin (2003) (green curve).Estimates of O' Brien and Moldwin (2003) based on log(AE) (cyan curve) or log(Dst) (yellow curve) are significantly higher during extended quiet intervals.

THEMIS observations
The SST flux data covers 12 energy ranges, from 26-36 keV up to 800-4000 keV, listed in Table 1.The four highest energy channels (> 300 keV) use "double detector logic" but the lower ones do not, and so are more vulnerable to proton contamination, despite extensive efforts at mitigation (D.Turner, personal communication, 2013).Also, as mentioned, only spin-averaged (i.e., not pitch angle-resolved) data are available for this period, so to proceed the flux measurements are assumed to be dominated by locally mirroring particles.
Locations of the five THEMIS spacecraft in Cartesian GSM (geocentric solar magnetospheric) coordinates are readily available at one minute intervals.Using the parameters discussed above, the IRBEM library was used to calculate values of L * at the spacecraft in the T04 field model, for approximately locally mirroring particles (pitch angle 80 • ) at intervals of 5 min.(Input parameters were linearly interpolated in time as necessary.) Figure 2 shows the spin-averaged flux vs. time at L * = 6 for a few values of energy.A dropout during 20 September (DOY 263) is quite evident for 408 and 720 keV, though weaker for 203.5 keV. Figure 3 shows the time development of flux vs. L * for fixed E = 720 keV.
Figure 4 shows zoomed-in plots of flux vs. time for a number of L and E values; the dropout occurs for all L * > 4. Also shown in each plot is a line segment whose slope was used to find the timescale of the assumed exponential decay during the dropout.Note that for E < 100 keV, fluxes increased rather than decreased.Because of the coarse time resolution, this is an overestimate of the timescale; the dropout

Radial diffusion simulation
Radial diffusion simulations evolve phase space density as a function of the adiabatic invariants, which must first be computed.The magnetic field and normalized field line integral I at the spacecraft location were obtained from IRBEM (also assuming local pitch angle 80 • ), and used to evaluate the adiabatic invariants M and J for each energy channel: where p 2 = (E/mc 2 )(2 + E/mc 2 ) and E is the center value of the channel.With B in gauss, I in R E , and mc 2 = 0.511 MeV, this gives M in MeV G −1 and J in g (cm s −1 ) R E (the units used by Brautigam and Albert, 2000).Then K, in units of R E √ G, is given by Figure 6 shows the resulting values of M, J , K, and L * for the five spacecraft.
To make the best use of the available coverage, attention is focused on the value M = 700 MeV G −1 , which dictates the value of p 2 or energy of the assumed locally mirroring particles.This also fixes the value of J , which cannot be specified independently; however, only measurements with calculated values J ≈ 4 × 10 −17 (K ≈ 0.014), within a factor of 2, are used.The flux j at the required energy is found by interpolating log j in log E. The phase space density is given by f = j/(p/mc) 2 ; the resulting values are shown in Fig. 7. (Further multiplying by the constant 6.3 × 10 −10 would give f in units of s 3 /km 6 .)Note that f (L) shows a substantial peak around L * = 5 prior to the dropout.These THEMIS values were used to generate initial and boundary values for radial diffusion simulations.A grid of 100 points was used, with L min = 3.7 and L max = 6.0.For each grid point L i , measurements with |L * −L i | < 0.1 within a time window were considered; to mitigate missing the onset of the dropout, the minimum of these f values was used.At each grid point, for initial conditions the time window was ±3 days wide, centered about t = 259.0(16 September), while for subsequent, time-dependent boundary values at L max and L min the window width was ±0.75 days.The boundary conditions were updated every 15 min; if no measurements were available within the time window, the previous values were retained.The resulting boundary values, f (t), for L max and L min are shown in Fig. 8 and show a sharp drop at L max lasting about 1.5 days.
Several different versions of radial diffusion coefficients were used, including the benchmark expressions D BA used by Brautigam and Albert (2000) and Albert et al. (2009).These represent historical estimates of substorm-associated electrostatic fluctuations and solar wind-driven electromagnetic ULF waves, both parameterized by Kp.Ozeke et al. (2014) obtained diffusion coefficients D Ozeke from much more recent, ground-based measurements and a different decomposition into electric and magnetic terms, also parameterized by Kp.Finally, a simulation of this interval (covering t = 263.0 to 265.9) was done using the Lyon-Fedder-Mobarry MHD (magnetohydrodynamics) code, and radial diffusion coefficients D MHD were constructed from the fluctuating electric and magnetic fields as in an earlier study described by Elkington et al. (2012).As shown in Fig. 9, all three expressions are quite similar over this range of L and time, at least for M = 700 MeV G −1 , although the D Ozeke values seem systematically the smallest.The D MHD values are mostly smaller than the D BA values except for a short interval around t = 264.5 (noon of 21 September).
Numerical integration of the well-known 1-D radial diffusion equation, was carried out with standard numerical methods, using the initial and time-dependent boundary values described above.
Figure 10 shows the results using D BA , while for the run of Fig. 11 D MHD was substituted when available (t = 263.0 to 265.9).In both cases, the phase space density does not decrease very much for L * < L max during the time that f (L max ) is depressed; outward radial diffusion occurs, but not nearly at the rate required to agree with the data over the range L * = 4 to 5.5 shown in Fig. 7.Even arbitrarily increasing D MHD by a factor of 3, as shown in Fig. 12, does not help much.During this same time, the simulated phase space density at L * < 4.5 actually increases.These results are not surprising.Brautigam and Albert (2000) and Albert et al. (2009) found that radial diffusion, even with data-driven boundary conditions, was unable to reproduce CRRES (Combined Release and Radiation Effects Satellite) observations of dropouts (as well as rebuilding) during the 9 October 1990 storm.Shprits and Thorne (2004) concluded more generally that "simulations with variable outer boundary conditions show that the depletion of the main phase relativistic electron fluxes at L ≤ 4 can not be explained only by variations in fluxes near geosynchronous orbit and require local lifetimes as short as 0.5 day" and that "even strong variations in the outer boundary are unable to cause the observed depletion of relativistic electron fluxes in the heart of the radiation belt (L ≤ 4) during the main phase of the storm."Su et al. (2010) also found that "combined radial diffusion and adiabatic transport contributes insignificantly to the main phase depletion of energetic electron fluxes within 5 R E ."Similar conclusions are implicit in Fig. 3a of Kim et al. (2011), Fig. 12 of Hwang et al. (2013) and Fig. S3 of Turner et al. (2012a).(Lyons and Thorne, 1973).Comparing observations to simulations with constant boundary conditions, Shprits et al. (2005) (2012) presented fits in L, E, and Kp to simple estimates of the timescale of pitch angle diffusion by chorus waves.The obvious problem with this approach is that it neglects energy diffusion, which tends to increase flux levels in the multi-keV range.Indeed, the increases seen in the THEMIS data following the dropout are most likely attributable to chorus acceleration (e.g., Li et al., 2007;Albert et al., 2009).Electromagnetic ion cyclotron (EMIC) waves can lead to very rapid pitch angle diffusion of relativistic electrons, and are expected to be proximate to the plasmapause (Summers et al., 1998;Li et al., 2007).Borovsky and Denton (2009) found from a study of geosynchronous data that dropouts, attributed to EMIC waves, coincide with the formation of a plasmaspheric drainage plume, in conjunction with a dense and hot plasma sheet.Usanova et al. (2012Usanova et al. ( , 2013) ) found from THEMIS and Cluster surveys that EMIC within plumes correlates with solar wind dynamic pressure, though the occurrence rate is still in the 5-10 % range inside of geosynchronous orbit.However, getting the minimum resonant energy below about 1 MeV requires very large values of k .Cold plasma theory provides a means for this, near ion stop bands (e.g., Albert, 2003), and some simulations have exploited this mechanism (e.g., Li et al., 2007;Su et al., 2011a, b).However, careful consideration of thermal effects (Silin et al., 2011;Chen al., 2011Chen al., , 2013) ) seems to invalidate this for sub-MeV electrons.Recent observations of EMIC waves and energy-dependent electron precipitation (Usanova et al., 2014) seem to confirm this.
Pitch angle scattering by hiss has also been considered.Lam et al. (2007) developed a statistical model of hiss from CRRES data, parameterized by Kp, encompassing both the plasmasphere and plumes; amplitudes were about 30 pT for 2 ≤ Kp < 4 and about 45 to 70 pT for Kp ≥ 4. A radial diffusion simulation found encouraging, though not decisive, agreement with observed fluxes.Li et al. (2007) provided a simple, convenient model which was originally presented for  2014) (green curves), and an MHD simulation (Elkington et al., 2012) (red curves), as described in the text.
Figure 10.Time development of phase space density from a simulation using the radial coefficient of Brautigam and Albert (2000).a 2-D study at L = 4.5, but which was quickly adopted for all L (e.g., Shprits et al., 2009).A study by Summers et al. (2008) of CRRES wave data restricted to plumes found considerably smaller wave amplitudes, but it may be argued that the plume crossings, mostly on the night side, were not favorably located.
Here, the Li et al. (2007) models of dayside chorus, nightside chorus, and whistler mode hiss were used to calculate quasi-linear pitch angle diffusion coefficients at L = 4.5 (with the additional assumption that wave normal angles had a gaussian distribution, with a characteristic width of ≤ 30 • ).For the case where L = 4.5 is outside the plasmasphere, and the hiss occurs in an extended plasmaspheric plume, Li et al. (2007) modeled these three wave populations during a storm main phase as covering 25, 25, and 15 % of a drift orbit.The case of L = 4.5 lying entirely within the plasmasphere was also modeled, by considering the hiss model only (covering the entire drift orbit).The lifetimes associated with these diffusion coefficients are shown in Fig. 13.Both sets of values, especially the hiss-only case, are of the order of magnitude of the observed rates in Fig. 5, but the modeled lifetimes increase with energy while the observed lifetimes decrease with energy.Also, compared to the empirical values of Lam et al. (2007) and Summers et al. (2008), the assumed hiss amplitude of 100 pT seems rather high.Though this event was not a major storm, a significant dropout occurred which does not seem to be reproducible by the standard diffusion-based models.As discussed, though outward radial diffusion is currently a fashionable paradigm, it does not seem sufficient to explain decreases at L * ∼ 4-5 for many events, including this one.Reasonable (though not the most recent) models of chorus and hiss, combined with quasi-linear theory, seem marginally effective enough, but it is very worrisome that the dependence on energy seems qualitatively wrong.In fact the increasing loss rates as particle energy increases is suggestive of EMIC waves, but conversely, as particle energy decreases the observed rates decrease too slowly (if current ideas about the minimum resonant energy are correct).Similar conclusions, based on different events, were recently reached by Morley et al. (2010), Turner et al. (2014), andHudson et al. (2014).
Possible resolutions include: qualitatively similar processes evaluated with better models of radial diffusion and cyclotron-resonant waves; non-diffusive radial transport (Ukhorskiy et al., 2006); non-resonant interactions with EMIC or other waves (Bortnik and Thorne, 2010); and nondiffusive wave-particle interactions (e.g., Albert, 2000;Albert et al., 2012).It has been noted in many test particle simulations that nonlinear interactions with large amplitude, resonant whistler waves lead to a preferential, rapid decrease in pitch angle, accompanied by decrease in energy (Albert, 2002;Tao et al., 2012); this could also be manifested as a loss process.

Figure 1 .
Figure1.Interplanetary and magnetospheric quantities for a time interval surrounding 20 September 2007.Top to bottom: Dst, Kp, solar wind velocity, solar wind density, solar wind dynamic pressure, IMF B y , IMF B z , and the integrated quantities W 1 -W 6 (magenta, red, blue, cyan, green, and black, respectively) used the T04 magnetic field model.The final panel shows the last closed drift shell computed according to the T04 magnetic field model (black) and the T96 model (red).Also shown in the final panel is the plasmapause location according to several models as described in Sect.2.1.

Figure 2 .
Figure 2. Particle flux, in units of #/(cm 2 s ster keV), from three energy channels of the SST detectors on the THEMIS spacecraft, at L * = 6.0 ± 0.2.The different spacecraft are identified by the color coding.

Figure 3 .
Figure 3. Observed particle fluxes vs. L * at fixed E, showing the time evolution.

Figure 4 .
Figure 4. Particle fluxes just before and after the dropout.The color coding denotes the individual spacecraft as in the previous figure.The slopes of the dotted lines were used to determine timescales.

Figure 5 .
Figure 5. Decay timescales vs. L * for several values of E, determined from the slopes indicated in the previous figure.

Figure 6 .
Figure 6.L * vs. time for the five THEMIS spacecraft, in the T04 magnetic field model.Also shown are the values of adiabatic invariants M, J , and K covered by the energy channels of the SST detector, assuming locally mirroring particles.

Figure 7 .
Figure 7. Phase space density vs. L * for fixed M = 700 MeV G −1 , calculated from the observed fluxes, showing the time development.

Figure 8 .
Figure 8. Phase space density at L = 6.0 (red) and L = 3.7 (blue) calculated from the observed fluxes, as functions of time, used as boundary conditions for radial diffusion simulations.

Figure 11 .
Figure 11.Time development of phase space density from a simulation using the MHD-based radial coefficient described in the text.

Figure 12 .
Figure 12.Time development of phase space density from a simulation using the MHD-based radial diffusion coefficient multiplied by 3.

Figure 13 .
Figure13.Decay timescales at L = 4.5 from a model of plasmaspheric hiss (red), and from a combination of chorus waves and hiss in a plasmaspheric plume (blue).

Table 1 .
Minimum, maximum, and center values of the THEMIS SST energy channels.