Thermospheric Neutral Winds as the Cause of Drift Shell Distortion in Earth’s Inner Radiation Belt

Recent analysis of energetic electron measurements from the Magnetic Electron Ion Spectrometer instruments onboard the Van Allen Probes showed a local time variation of the equatorial electron intensity in the Earth’s inner radiation belt. The local time asymmetry was interpreted as evidence of drift shell distortion by a large-scale electric field. It was also demonstrated that the inclusion of a simple dawn-to-dusk electric field model improved the agreement between observations and theoretical expectations. Yet, exactly what drives this electric field was left unexplained. We combine in-situ field and particle observations, together with a physics-based coupled model, the Rice Convection Model (RCM) Coupled Thermosphere-Ionosphere-Plasmasphere-electrodynamics (CTIPe), to revisit the local time asymmetry of the equatorial electron intensity observed in the innermost radiation belt. The study is based on the dawn-dusk difference in equatorial electron intensity measured at L = 1.30 during the first 60 days of the year 2014. Analysis of measured equatorial electron intensity in the 150–400 keV energy range, in-situ DC electric field measurements and wind dynamo modeling outputs provide consistent estimates of the order of 6–8 kV for the average dawn-to-dusk electric potential variation. This suggests that the dynamo electric fields produced by tidal motion of upper atmospheric winds flowing across Earth’s magnetic field lines - the quiet time ionospheric wind dynamo - are the main drivers of the drift shell distortion in the Earth’s inner radiation belt.


INTRODUCTION
The Van Allen Probes (Mauk et al., 2013) have provided unprecedented amounts of high quality energetic (10-100 keV) electron flux measurements near the magnetic equator in the inner belt and slot region (Reeves et al., 2016) below an equatorial altitude of about 3 Earth Radii (L < 3). These new measurements confirmed known dynamical features such as "zebra stripe" patterns in the spectrograms of energetic electrons (Imhof and Smith, 1965;Ukhorskiy et al., 2014) and energetic electron injections deep into the inner magnetosphere (Pfitzer and Winckler, 1968;Turner et al., 2015). Electron flux measurements in the 100-400 keV energy range from the Magnetic Electron Ion Spectrometer (MagEIS) instruments (Blake et al., 2013) also revealed a surprisingly persistent local time asymmetry of the equatorial electron intensity below L 1.4 during geomagnetically quiet times, with higher radiation belt electron intensities near dawn (Selesnick et al., 2016).
Radiation belt particles are magnetically trapped high energy particles whose motion exhibits three quasi-periodic types of motion occurring on three distinct timescales (e.g., Northrop and Teller, 1960). The slowest periodicity corresponds to a drift motion around the planet that is the combination of an energy-independent electric drift and a gradient-curvature magnetic drift proportional to the particle's momentum. Radiation belt particles' drift motion defines closed surfaces known as drift shells. At high enough energies (>100 keV), the role played by large-scale electric fields is usually omitted: Radiation belt particles' momentum is expected to be conserved along a drift shell and the flux is expected to be constant at a fixed L value in the absence of any significant source or loss mechanism (e.g. Roederer, 1967). The observed local time asymmetry of the equatorial electron intensity at fixed L values below L 1.4 challenged the latter, suggesting a drift shell distortion by a quasi-static electric field.
Even during geomagnetically quiet times, large-scale electric fields of the inner magnetosphere are more complex than the simple superposition of 1) a corotation electric field due to the rigid corotation of a perfectly conducting ionosphere and 2) a convection electric field set up by the coupling between the solar wind and the magnetosphere (e.g. Wolf et al., 2006). Dynamo electric fields produced by tidal motion of upper atmospheric winds flowing across the Earth's magnetic field lines-the ionospheric wind dynamo (Richmond, 1989)-are usually larger than subauroral convection electric fields below L ∼ 2 (e.g., Figure 1 in the review by Mozer (1973). Longitudinal, seasonal, and solar cycle variations of the mid-and lowlatitude ionospheric electric fields have been reported, together with a large day-to-day variability present even during geomagnetically quiet times (e.g., Fejer, 1993;Chau et al., 2010;Pfaff et al., 2010).
The objective of this study is to re-examine the time interval analyzed by Selesnick et al. (2016) in order to determine the origin of the large-scale electric field causing the observed radiation belt distortion. Observational and modeling resources are introduced in Section 2 to describe the electric field properties associated with the drift shell distortion observed. In particular, we present a method to infer electric field properties directly from an analysis of asymmetries measured in differential directional fluxes. Experimental, numerical and analytical experiments provide consistent electric potential variation estimates, of the order of 6 − 8 kV between the dawn and dusk sectors at the magnetic equator of L 1.30, i.e. at a magnetic latitude of ∼ 29°at the ionospheric field line footpoint (Section 3). The similarities between experimental data analysis and outputs from a quiet time run by the Rice Convection Model (RCM) Coupled Thermosphere-Ionosphere-Plasmasphere-electrodynamics (CTIPe) model suggests that the electric potential variation results mainly from the quiet time wind dynamo.

MATERIALS AND METHODS
Leveraging Liouville's theorem and adiabatic invariant theory, we conduct an analysis of energetic electron directional differential fluxes measured by the Van Allen Probes to infer electric field properties in the inner belt. The resulting electric field properties are compared with outputs from Van Allen Probes field measurements and from a physics-based coupled model, RCM-CTIPe.

Field and Particle Measurements Onboard the Van Allen Probes
The time interval of the study corresponds to the first 60 days of the year 2014. This time interval was also selected for the study by Selesnick et al. (2016) because 1) the two orbital legs are near dawn and dusk at L 1.3 and 2) it is geomagnetically quiet. Although some moderate (Echer et al., 2006) geomagnetic activity occurs towards the end of February, the time interval contains no significant inner belt injection. Additional details on the geomagnetic conditions associated with this time interval are provided in Section 4.
The Van Allen Probes (Fox and Burch, 2014) were twin spacecraft with similar highly elliptical orbits (perigee at ∼ 600 km, apogee at ∼ 5.8 R E ), an orbital period of about 9 h, and an inclination close to 10°. Spacecraft apogees drifted slowly so that it took a bit less than 2 years for spacecraft apogees to scan all local time sectors (precession rate of about 210°per year). Successive inbound and outbound crossings of L 1.3 occurred fast, within 30 min. During the first 60 days of 2014, inbound crossings of L 1.3 are in the dusk-premidnight region (20-22 MLT) while outbound crossings are in the dawn sector (4-6 MLT) for both Van Allen probes. The geographic longitudes of the crossings depend on time (UT).
In this study, we rely primarily on the directional differential fluxes provided by the MagEIS instrument (processed to level 3). Electric field measurements are also analyzed to compare and contrast with proposed electric field models. Experimental electric field information comes from measurements by the Electric Field and Waves (EFW) instrument (Wygant et al., 2013) and by the Electric and Magnetic Field Instrument Suite and Integrated Science (EMFISIS) (Kletzing et al., 2013). Field and particle measurements immediately following spacecraft maneuvers are omitted.

Differential Fluxes for Equatorial Electrons
The analysis focuses on equatorial electron intensity in the 150-400 keV energy range. MagEIS data below 400 keV is normally well above background. We follow the approach developed by Selesnick et al. (2016) to determine fluxes of equatorially mirroring electrons. MagEIS pitch-angle resolved measurements are extrapolated when the spacecraft is close to the magnetic equator. Specifically, directional differential fluxes are extrapolated to determine the equatorial electron intensity, j o , when the local field magnitude, B, is such that B o ≤ B < 1.08B o , where B o is the equatorial magnetic field calculated according to the International Geomagnetic Reference Field to which Olson and Pfitzer, 1977 Quiet is superimposed (Olson and Pfitzer, 1977). We define the magnetic shell parameter, L, following the definition provided by McIlwain (1961). For equatorial particles, the parameter L is such that: where B E 30, 000 nT is the magnetic equatorial field at the surface of the Earth. This definition guarantees that constant L values define curves of constant equatorial magnetic field amplitude, i.e., (B o cst.)5(L cst.), even in non-dipolar fields. When experimental information is required at a fixed L value, we interpolate as a function of L using adjacent measurements within ± 0.025 L when available. Linear interpolations of adjacent logarithmic flux measurements are also performed to estimate partial derivatives with respect to L and kinetic energy, T (with window sizes of ± 0.025 L and ± 155 keV, respectively).

Electric Drift Measurements
Electric drift measurements are pre-processed following the approach developed over the years by Lejosne and Mozer (2016a, 2016b, 2019 and briefly described below. After a slight correction applied to the orientation of the magnetometer axes (Lejosne and Mozer, 2016a), the spin-averaged ( ∼ 12 s) electric drift measurement is reformulated in a local orthonormal frame of reference (e ρ , e b , e φ ) set by the measured magnetic field direction: is the magnetic field vector at the location r. e ρ and e φ indicate the radial and azimuthal directions, respectively. The shorting factor is calibrated to a value very close to 1 (Lejosne and Mozer, 2019). Once corotation is subtracted, the electric drift measurement is projected to the magnetic equator assuming equipotential field lines. The magnetic field model is set to the International Geomagnetic Reference Field for the Earth's internal field. The choice of the external magnetic field model is unimportant when mapping at L 1.30. During data processing, the Kp-driven Tsyganenko (1989) magnetic field model was superimposed (Lejosne and Mozer, 2016b).

Numerical Model RCM-CTIPe
Electric potential and electric field values at 100 km altitude are provided by the Rice Convection Model (RCM) Coupled Thermosphere-Ionosphere-Plasmasphere-electrodynamics (CTIPe) model. This self-consistently coupled model of the magnetosphere, ionosphere and thermosphere is composed of three physical models: 1) the Coupled Thermosphere-Ionosphere-Plasmasphere with self-consistent electrodynamics (CTIPe) model Millward et al., 2001Millward et al., , 1996; 2) the Rice Convection Model-RCM (Toffoletto et al., 2003;Wolf, 1983); and 3) the global electrodynamic solver based on the National Center for Atmospheric Research Thermosphere-Ionosphere Electrodynamics General Circulation Model-NCAR-TIEGCM (Richmond and Maute, 2014). The model includes the electrodynamic coupling, interactions, and feedback between the inner magnetosphere and the thermosphere-ionosphere-magnetosphere system. A detailed description of the model coupling can be found in the article by Maruyama et al. (2011).

Theoretical Framework for the Inference of Electric Field Properties From Directional Differential Flux Analysis
The objective of this section is to show how to infer information on electric field properties from measured variations of differential fluxes of equatorially trapped particles. This theoretical framework is similar to the one developed by Lejosne and Mozer (2020) for the analysis of zebra stripe patterns. It is adapted below to the case of drift shell distortion by quasi-static electric fields.
Leveraging Liouville's theorem and adiabatic invariant theory, the variation of equatorial electron fluxes measured at different local times along the same B o cst. curve is related to trapped particle kinetic energy variation (Section 2.3.1.). Section 2.3.2. details the relationship between trapped particle kinetic energy variation and electric potential variation. Section 2.3.3. provides the equation implemented hereafter and summarizes the assumptions underlying the overall approach.

Link Between Fluctuations in Directional Differential Fluxes of Equatorially Trapped Particles and Kinetic Energy Variations
According to Liouville's theorem, the phase space density, f j o /p 2 , is constant along the continuous trajectories of the system in phase space, i.e., along the dynamical path of particles, in the absence of any source or loss. Thus, any variation of the trapped population momentum, dp/p, is related to a variation of the directional differential flux, dj o /j o along the drift shell: This equation can also be rewritten in terms of kinetic energy, T, and kinetic energy variation, dT, since: dp p where c (T + mc 2 )/mc 2 is the Lorentz factor and mc 2 511 keV is the rest mass energy for an electron. Assuming conservation of the first invariant, p 2 /2mB o cst., any variation of the momentum is related to a variation in the equatorial magnetic field along the drift shell: Noting that B o B E /L 3 (Eq. 1), the combination of Eqs 3, 4 also relates trapped particle energization, dT, and radial motion, dL: The variations discussed above occur along the same equatorial drift shell, i.e.: where (T, B o , φ 1 ) and (T + dT, B o + dB o , φ 2 ) are coordinates along the same equatorial drift shell. The phases φ 1 and φ 2 represent the magnetic local times of an L-crossing during consecutive inbound and outbound passes, respectively. In practice, the shape of the distorted drift shell is unknown, and what is measured is in fact the fluctuation in differential equatorial flux at fixed kinetic energy, T, and different magnetic local times along the same B o cst. curve. In other words, the fluctuation measured is: Assuming small distortions of the field, the fluctuation measured, δj o , is related to the fluctuation along the dynamical path of particles, dj o , in the following way: Combining Eqs 3, 4 and 8, and noting that dx/x d ln x, we obtain the relationship between the variation of equatorial electron differential fluxes measured at different local times along the same B o cst. curve, δj o , and trapped particle energy variation, dT:

Link Between Kinetic Energy Variation and Electric Potential Variation
The time rate of change of the average kinetic energy of a guiding center of charge q is provided by the energy equation (Roederer and Zhang, 2014): where V refers to the drift velocity of the guiding center and F is an external non-electromagnetic force that is omitted hereafter. The electric field, E, is described here as the sum of a corotation electric field and a subauroral electrostatic field: where Ω is the angular velocity vector of the Earth's rotation, and U is the electrostatic potential associated with the subauroral electrostatic field. Noting that zB/zt −(Ω × r) · ∇B, it follows from the combination of Eqs 10, 11 in the case of energetic electrons (q < 0) that a variation in kinetic energy, dT, along the dynamical path of particles is related to a variation in electrostatic potential, dU: where e q > 0 is the elementary charge and ς(T) corresponds to the ratio between electric and total (magnetic + electric) angular drift velocities where R E 6, 370 km is one Earth radius. Eq. 12 describes the conservation of the total energy of the guiding center along the drift shell (e.g., Whipple, 1978).

Summary of the Theoretical Framework
A change of variables from B o to L (Eq. 1) and the combination of Eqs 9, 12 provides the relationship between a measured fluctuation in the equatorial electron intensity, δj o , and an electrostatic potential variation, dU. Let us emphasize that the variations measured are between the magnetic local times of the inbound (φ 1 ) and outbound (φ 2 ) crossings of the same B o cst. curve: where , and the partial derivatives are computed outbound at (T, B o , φ 2 ). Theoretically, the electric potential variation, dU, on the left side of Eq. 14, is a quantity independent of the kinetic energy of the trapped population considered, T. Thus, the right side of Eq. 14 is expected to be independent of the variable T as well.
The assumptions underlying this equation are that there is no significant source or loss mechanism on the time scale of the equatorially trapped population drift period, no external nonelectromagnetic force, no parallel electric field, no other significant source of magnetic field time variation besides Earth's rotation, and no significant time variation of the electric potential on the timescale of the electron drift period (this includes the timescale of spacecraft motion from the inbound to the outbound crossings of L 1.3, which is of the order of 30 min).

RESULTS
Section 3.1. presents experimental information on the electric potential variation between dawn and dusk at L 1.30 (magnetic latitude of ∼ 29°at the corresponding ionospheric field line footpoint), leveraging field and particle measurements independently. Section 3.2. compares these results with estimates from analytical and numerical models. Section 3.3. Conclusion of the Experimental Analysis summarizes our findings. Figure 1 is an introduction to the approach. During the first 60 days of 2014, the Van Allen Probes crossings of L 1.30 are in the dusk-premidnight region (20-22 MLT) during inbound, φ 1 , and in the dawn (4-6 MLT) sector during outbound, φ 2 ( Figure 1A). The equatorial intensity of 226 keV electrons at L 1.30 measured by Van Allen Probes A (RBA) is ∼40% higher at dawn (in red Figure 1B) than dusk (in blue Figure 1B). The Frontiers in Astronomy and Space Sciences | www.frontiersin.org September 2021 | Volume 8 | Article 725800 corresponding equatorial intensity fluctuation at 226 keV and L 1.30, δj o /j o , is evaluated by subtracting the outbound electron intensity (the red line Figure 1B) from the inbound electron intensity (the blue line Figure 1B), and normalizing by the inbound value. This was done for 40 different consecutive sequences of near equatorial inbound/outbound crossings of L 1.3 during the 60 days considered ( Figure 1C). The partial derivatives of ln j o with respect to ln L and ln T were computed locally to provide a time series for the electric potential variation, dU, between the inbound and outbound local times of L 1.30 following Eq. 14 ( Figure 1D). The time variation of the electric potential variation, dU(t), is similar to the time variation of the intensity fluctuation, δj o /j o (t), because the coefficient by which δj o /j o is multiplied to obtain dU (Eq. 14) does not vary significantly (<10%) during the time interval considered. Analysis of the 226 keV equatorial electron intensity measured by Van Allen Probes A during the first 60 days of 2014 suggests that the average electric potential variation between dawn and dusk is 7.6 ± 2.7 kV.

Experimental Information on the Electric Potential Variation Between Dawn and Dusk Based on Measured Equatorial Electron Intensity Asymmetry
The approach is extended to both Van Allen Probes, and to all four MagEIS energy channels between 150 and 400 keV. The results are presented Figure 2. While the times series for the electric potential variations seem to depend on the measuring spacecraft (Figures 2A,B), the averages, and standard deviations of the electric potential variation, dU, appear independent of the energy channel ( Figure 2C). That the statistical characteristics of the electric potential do not depend on MagEIS energy channel is consistent with theoretical expectations (Section 2.3.3.).
The analysis suggests that the average value for the electric potential variation between dawn and dusk is dU 8 ± 4 kV during the first 60 days of 2014. In terms of average electrostatic field, since dU − E φ r o dφ, the electric potential variation corresponds to an average radial electric field between dusk (φ 1 ) and dawn (φ 2 ) of [E φ ] −0.5 ± 0.2 mV/m at the magnetic equator of L 1.3.

Experimental Information on the Electric Potential Between Dawn and Dusk Based on Electric Field Measurements
While Van Allen Probes only provide local electric field samples at L 1.3 at precisely dusk (φ 1 ) and dawn (φ 2 ) during the time interval considered, it is possible to obtain another estimate of the electric potential variation between dawn and dusk by considering Van Allen Probes electric field measurements over a longer time interval. Figure 3 provides experimental information on the characteristics of the DC electric field measured by the Van Allen Probes during the years 2013 and 2014 at L 1.3 ± 0.05 during geomagnetically quiet times (Kp < 3). The local time averages and standard deviations are computed over running windows of 1 hour size and they are represented by a black line and a shaded area, respectively. Assuming that the electrostatic electric field (−∇U) is quasi-stationary, the 24 hr-MLT average of the electric field components should be zero. Experimentally, the 24 hr-MLT averages of the radial and azimuthal electric field components measured are 0.02 mV/m and 0.20 mV/m, respectively. These small offsets could be due to small residual calibration errors (e.g. Lejosne and Mozer, 2019). They could also be due to the irregular seasonal and longitudinal samplings. These small offsets have been subtracted from the datasets before estimating the average electric field between dusk and dawn and the corresponding electric potential variation between dawn and dusk. In that context, the average azimuthal component of the electric field between dusk and dawn is −0.4 ± 0.1 mV/m , and the corresponding variation of the electric potential is 7 ± 2 kV. These are comparable to the estimates provided by the particle data analysis (Section 3.1.1.).

Model-Observation Comparison for the Electric Potential and Electric Field Components
In this section, we provide electric potential variation between dawn and dusk predicted by analytic and numerical models.

Comparison With an Analytical Expression for the Electric Potential Variation
The simple electrostatic potential proposed by Selesnick et al. (2016), U s , is: where E c 0.4 mV/m. The magnitude of E c was set in order to correctly represent the average distortion observed during the first 50 days of 2014. Under that model, the electric potential variation between dawn and dusk is of the order of 6 kV and the average azimuthal electric field component between dusk and dawn is −0.3 mV/m.

Comparison With Numerical Values During a Quiet Time Run With RCM-CTipe
24 h worth of electric field values during a quiet-time run with the RCM-CTIPe are provided in Figure 4. The quiet time run corresponds to a day in spring (set to March 17, 2013) during which the solar wind-magnetosphere dynamo is artificially forced to 0. In that context, the electric fields are driven solely by the quiet time wind dynamo. Because the magnetic local time (MLT) is a combination of universal time (UT) and longitude, the time variations observed at a given MLT are also representative of the longitudinal dependence of the electric field. Ionospheric values have been projected to the magnetic equator along dipolar equipotential field lines. This projection means that the radial (poleward) ionospheric component of the electric field has been multiplied by a factor ∼ 0.52 while the azimuthal (eastward) ionospheric component of the electric field has been multiplied by a factor ∼ 0.67 to obtain the corresponding amplitudes at the magnetic equator of L 1.3 (see also Lejosne and Mozer, 2016b, Section 2.1.3 for analytic expressions regarding electric field mapping). The modeled average electric potential variation due to the wind dynamo between dawn and dusk is 7 ± 1 kV according to RCM-CTIPe and the average azimuthal electric field component between dusk and dawn is −0.4 ± 0.1 mV. This first estimate is consistent with simulated values for the MLT distribution of the electric potential resulting solely from the quiet time motion of the neutral atmosphere of (Yamazaki and Maute, 2017, their Figure 8).

Conclusion of the Experimental Analysis
The different estimates for the dawn-dusk electric potential variation at L 1.3 are summarized in Table 1, together with the method and the time interval considered. Table 1 shows that different methods provide consistent estimates for the electric potential variation between dawn and dusk, with a magnitude of the order of 6-8 kV. The similarity is remarkable given that: 1) the neutral winds and associated wind dynamo vary from day to day and 2) every method was applied under a unique set of conditions.
An electric potential variation of 6-8 kV between dawn and dusk is ∼15 times greater than would be predicted by the usual Maynard and Chen (1975) parametrization of the Volland (1973) -Stern (1975 convection electric field model during quiet geomagnetic activity (Kp 2), suggesting that the solar windmagnetosphere dynamo is unlikely to be the main of cause of the observed drift shell distortion.
Van Allen probes electric field measurements and simulations by RCM-CTIPe are in qualitative agreement, despite differences in their time intervals. These similarities confirm that the electric potential variation results mainly from the quiet time wind dynamo. Comparing Figures 3, 4 also demonstrates that the Van Allen Probes can resolve the effects of neutral winds on global electric fields. While the modeling and experimental results provide similar estimates for the average magnitude of the dawn-dusk electric potential variation, they do not explain the high variability revealed by the particle data analysis. The time series of dawndusk electric potential variation, dU(t), has a standard deviation of the order of 7 kV ( Figure 2C). On the other hand, 1 day of quiet time run by RCM-CTIPe provides a standard deviation of the order of 1 kV only, including possible longitudinal effects. Data analyses similar to the one presented Section 3.1.1. were performed 1) using measurements by the Radiation Belt Storm Probes Ion Composition Experiment (Mitchell et al., 2013) in place of MagEIS, and/or 2) at higher L values (specifically, L 1.5, and L 1.7). While quantitative analysis of RBSPICE data is incomplete, similar trends for the time variations of the dawn-dusk electric potential were observed. These findings suggest that the time variations observed correspond to real geophysical information.
We searched for possible dependencies of the time series, dU(t), by performing correlation analyses with various indicators. We looked for possible correlations between the   1 | Estimates of the electric potential variation between dawn and dusk, and average azimuthal electric field component between the dusk and dawn regions of the magnetic equator of L 1.3 according to different experimental, numerical and analytical methods.

Method
Time interval considered Electric potential variation dU(kV) between dusk (φ 1 ) and dawn (φ 2 ) at the magnetic equator of L = 1.3 Average azimuthal electric field [E φ ] between dusk (φ 1 ) and dawn (φ 2 ) at the magnetic equator of L = 1.3 Local asymmetries in radiation belt intensity ( Newell and Gjerloev (2011) and Ohtani and Gjerloev (2020) are indicated Panel (D) in black, blue and red, respectively.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org September 2021 | Volume 8 | Article 725800 8 sampled electric potential variation, dU(t), and the spacecraft longitudinal positions at the time, t, of measurements (UT). We also looked for possible variations based on a variety of indicators (presented in Figure 5), such as F10.7, Kp, Dst, and the polar cap potential, PCP, computed according to the formula by Boyle et al. (1997). We also considered substorm onsets as defined by Forsyth et al. (2015), Newell and Gjerloev (2011) and Ohtani and Gjerloev (2020). Regardless of the parameter considered, no obvious dependency was found.
Realistic modeling of the subauroral electric field evolution is required to further investigate the origin of the variability reported here.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: The Van Allen Probes electric field data are in public access in the RBSP/EFW database: www. space.umn.edu/rbspefw-data. The MagEIS and magnetic ephemeris data are available from the ECT Science Operations and Data Center, https://rbsp-ect. newmexicoconsortium.org/rbsp_ect.php Simulation outputs can be obtained by contacting NM.

AUTHOR CONTRIBUTIONS
The first author is the lead and corresponding author. All other authors are listed in alphabetical order. We describe contributions to the paper using the CRediT (Contributor Roles Taxonomy) categories (Brand et al., 2015