The Discrepancy Between Simulation and Observation of Electric Fields in Collisionless Shocks

Recent time series observations of electric fields within collisionless shocks have shown that the fluctuating, electrostatic fields can be in excess of one hundred times that of the quasi-static electric fields. That is, the largest amplitude electric fields occur at high frequencies, not low. In contrast, many if not most kinetic simulations show the opposite, where the quasi-static electric fields dominate, unless they are specifically tailored to examine small-scale instabilities. Further, the shock ramp thickness is often observed to fall between the electron and ion scales while many simulations tend to produce ramp thicknesses at least at or above ion scales. This raises numerous questions about the role of small-scale instabilities and about the ability to directly compare simulations with observations.


INTRODUCTION
Collisionless shock waves are an ubiquitous phenomenon in heliospheric and astrophysical plasmas. They most often manifest as a nonlinearly steepened fast magnetosonic-whistler wave that has reached a stable balance between steepening and some form of irreversible energy dissipation. If a balance is reached, a stationary shock ramp is formed. The shock ramp is the part of shock transition region between upstream and downstream with an abrupt, discontinuity-like change in number density (n S where s is the particle species), pressure 1, , quasi-static 2 magnetic field magnitude vector (B o ), and bulk flow velocity (V bulk ). The thickness of this ramp is thought to depend upon macroscopic shock parameters like the fast mode Mach number (M f ), shock normal angle, θ Bn (e.g., quasi-perpendicular shocks satisfy θ Bn ≥ 45°), and upstream averaged plasma beta (Sagdeev, 1966;Coroniti, 1970;Tidman & Krall, 1971;Galeev, 1976;Kennel et al., 1985).
The topic of interest for this study is electric fields in observations and simulations, so we will limit the discussion to wave-particle interactions and macroscopic quasi-static field effects. Further, given that the primary discrepancy between simulations and observations lies in the lack of large amplitude, high frequency electrostatic waves in the former, we will limit the discussion to high frequency electrostatic waves. Note that some PIC simulations do generate the electrostatic waves of interest but the simulations are often tailored to generate the modes (e.g., isolated simulation mimicking shock foot region) by artificially injecting known free energy sources (e.g., initialize with two counter-streaming beams). Therefore, all of the modes listed in the following discussion have been generated in PIC simulations (Dyrud and Oppenheim, 2006;Matsukiyo and Scholer, 2006;Matsukiyo and Scholer, 2012;Muschietti and Lembège, 2017;Saito et al., 2017). However, as will be shown, parameters like the wavelengths and amplitudes tend to differ from those in observations, sometimes significantly.
4 f lh f ce f cp , where f cs is the cyclotron frequency of species s ( qsBo ms where q s is the total charge, and m s is the mass of species s). 5 There are two modes radiated by the MTSI at collisionless shocks, both of which are very obliquely propagating and have real frequencies near or below f lh . The two free energy sources for the MTSIs are between incident electrons and reflected ions and incident electrons and incident ions. 1976; Lemons and Gary, 1977;Wu et al., 1983;Wu et al., 1984), electron beams (Papadopoulos and Palmadesso, 1976), and/or heat flux carrying electrons (Marsch and Chang, 1983). These modes are important for collisionless shock dynamics because they can stochastically accelerate both thermal electrons (parallel to B o ) and ions (perpendicular to B o ) to suprathermal energies (Wu et al., 1984;Cairns and McMillan, 2005).
Electrostatic IAWs have been observed in the solar wind and near collisionless shocks for over 40 years (Fredricks et al., 1968;Fredricks et al., 1970a;Gurnett and Anderson, 1977;Gurnett et al., 1979;Kurth et al., 1979). They present in spacecraft observations at frequencies, in the spacecraft frame, above the proton plasma frequency 6 (due to the Doppler effect), typically in the ∼1-10 kHz range in the solar wind near 1 AU. They are observed as linearly polarized (mostly parallel to B o but sometimes at small oblique angles), modulated sine waves with bursty wave envelopes lasting 10 s of ms (Wilson et al., 2007;Wilson et al., 2014a;Wilson et al., 2014b). They have been shown to have wavelengths on the order of a few to several Debye lengths 7, , or 10-100 s of meters near 1 AU (Fuselier and Gurnett, 1984;Breneman et al., 2013;Goodrich et al., 2018;Goodrich et al., 2019). They are thought to be driven unstable by the free energy in currents (Biskamp et al., 1972;Lemons and Gary, 1978), temperature gradients (Allan and Sanderson, 1974), electron heat flux (Dum et al., 1980;Henchen et al., 2019), or ion/ion streaming instabilities (Auer et al., 1971;Akimoto and Winske, 1985;Akimoto et al., 1985b;Goodrich et al., 2019) or they can result from a nonlinear wave-wave process (Cairns and Robinson, 1992;Dyrud and Oppenheim, 2006;Kellogg et al., 2013;Saito et al., 2017). These modes are important for collisionless shock dynamics because they can stochastically accelerate thermal electrons (parallel to B o ) generating self-similar velocity distribution functions (VDFs) or the so called "flattop" distributions (Vedenov, 1963;Sagdeev, 1966;Dum et al., 1974;Dum, 1975;Dyrud and Oppenheim, 2006). They are also capable of stochastically accelerating the high energy tail of the ion VDF (parallel to B o ) (Dum et al., 1974). Note that the generation of the flattop has recently been interpreted as evidence of inelastic collisions (Wilson et al., 2019a;Wilson et al., 2019b;Wilson et al., 2020).
ESWs present in spacecraft observations as short duration (few ms), bipolar electric field pulses parallel to B o and monopolar perpendicular (Behlke et al., 2004;Wilson et al., 2007;Wilson et al., 2014b). They tend to be on Debye scales and are thought to be BGK phase space holes (Ergun et al., 1998;Cattell et al., 2005;Franz et al., 2005;Vasko et al., 2018). ESWs can be driven unstable by electron beams (Ergun et al., 1998;Cattell et al., 2005;Franz et al., 2005), ion beams (Vasko et al., 2018), modified two-stream instability (MTSI) (Matsukiyo and Scholer, 2006), or the produce of high frequency wave decay (Singh et al., 2000). Until recently, it was thought all ESWs outside the auroral acceleration region were electron holes. However, work by (Vasko et al., 2018) and  suggest that many of the ESWs in the terrestrial bow shock are not only ion holes, they do not propagate exactly along B o as was previously thought. ESWs are important in collisionless shock dynamics because they can trap incident electrons (Dyrud and Oppenheim, 2006;Lu et al., 2008) or ions (Vasko et al., 2018;Wang et al., 2020), depending on the type of hole. They have also been shown to dramatically heat ions (Ergun et al., 1998), and/or couple to (or directly cause) the growth of IAWs (Dyrud and Oppenheim, 2006), whistler mode waves (Singh et al., 2001;Lu et al., 2008;Goldman et al., 2014), LHWs (Singh et al., 2000).
The ECDI is driven by the free energy in the relative drift between the incident electrons and shock-reflected ions (Forslund et al., 1970;Forslund et al., 1971;Lampe et al., 1972;Matsukiyo and Scholer, 2006;Muschietti and Lembège, 2013). They also range from Debye to electron thermal gyroradius scales (Breneman et al., 2013) and present in spacecraft observations as mixtures of Doppler-shifted IAWs and electron Bernstein modes. The polarization of these modes can be confusing, presenting as shaped like a tadpole or tear drop, with one part of the "tadpole" nearly parallel to B o (i.e., IAW part) and the other nearly orthogonal (i.e., the Bernstein mode part) Breneman et al., 2013;Wilson et al., 2014b;Goodrich et al., 2018). This results from the coupling between two modes that are normally orthogonal to each other in their electric field oscillations. ECDI-driven modes are important for collisionless shocks because they can resonantly interact with the bulk of the ion VDF, generate a suprathermal tail on the ion VDF, and strongly heat the electrons perpendicular to B o (Forslund et al., 1970;Forslund et al., 1972;Lampe et al., 1972;Muschietti and Lembège, 2013).
Langmuir waves have been observed upstream of collisionless shocks for decades (Gurnett and Anderson, 1977;Filbert and Kellogg, 1979;Kellogg et al., 1992;Cairns, 1994;Bale et al., 1998;Bale et al., 1999;Malaspina et al., 2009;Soucek et al., 2009;Krasnoselskikh et al., 2011). These waves have k λ e ( 1 (Soucek et al., 2009;Krasnoselskikh et al., 2011) and rest frame frequencies satisfying f rest ( f pe . Langmuir waves are driven unstable by electron beams and/or nonlinear wave decay (Pulupa et al., 2010;Kellogg et al., 2013). They tend to be linearly polarized nearly parallel to B o when electrostatic but some do exhibit circular polarization when electromagnetic (Bale et al., 1998;Malaspina and Ergun, 2008). Langmuir waves are relevant to collisionless shock dynamics in that they dissipate the free energy in reflected electron beams and can mode convert to generate free mode emissions that can serve as remote detection signatures (Cairns, 1994;Bale et al., 1999;Pulupa et al., 2010).
In summary, the most commonly observed electrostatic wave modes near collisionless shocks are IAWs, ESWs, ECDI-driven modes, and Langmuir waves. Electrostatic LHWs are less commonly observed, which may be due to instrumental effects 6 2 π f ps ns qs 2 εoms where s is the particle species. 7 λ De εo k B Te ne e 2 where n e is the electron number density.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org January 2021 | Volume 7 | Article 592634 as many electric field instruments (Bonnell et al., 2008;Bougeret et al., 2008;Cully et al., 2008) have been designed with gain rolloffs at ∼1-10 Hz (low-or high-pass filters), which happens to be the typical value of f lh in the solar wind near 1 AU. It may also be that electrostatic LHWs are just less commonly generated or damp very quickly in collisionless shocks. Langmuir waves tend to occur upstream of the shock in regions filled with shock reflected electron beams (Cairns, 1994;Bale et al., 1999;Wilson et al., 2007;Pulupa et al., 2010). Although they can be common in the upstream, they tend to be much less so in the ramp and immediate downstream region. Therefore, the remaining discussion will focus on the most commonly observed Debye-scale, electrostatic modes: IAWs, ESWs, and ECDI-driven modes. These three modes are observed in and around both quasi-parallel and quasi-perpendicular shocks. The only macroscopic shock parameters on which they appear to depend are the shock density compression ratio and M f [e.g., Wilson et al., 2007;Wilson et al., 2014a;Wilson et al., 2014b]. The ECDI-driven modes tend not to be observed for M f ( 3, since they require sufficient reflected ions to initiate the instability. Part of the reason for the lack of dependence on shock geometry is that the fluctuations in the foreshock upstream of a quasi-parallel shock, for instance, locally rotate the magnetic field to quasi-perpendicular geometries and some can even locally reflect/energize particles [e.g., Wilson et al., 2013;Wilson et al., 2016].

Spacecraft Observations
Early spacecraft electric field observers had very limited resources, compared to modern day, in memory, computational power, and spacecraft telemetry. As such, the common practice was to perform onboard computations to generate Fourier spectra for predefined frequency ranges (Fredricks et al., 1968;Fredricks et al., 1970a;Fredricks et al., 1970b;Rodriguez and Gurnett, 1975). These Fourier spectra are spectral intensity data averaged over fixed time and frequency intervals, which has been more recently shown to significantly underestimate the instantaneous wave amplitude (Tsurutani et al., 2009). The underestimation led to some confusion in multiple areas of research because the estimated wave amplitudes from the spectra were too small to noticeably impact the dynamics of the system in question. For instance, for decades the radiation belt community had relied upon such dynamic spectra and came to conclusion that the whistler mode waves (e.g., chorus and hiss) were typically in the (1 mV/m amplitude range. The advent of time series electric field data led to the discovery that some of these modes could have amplitudes in excess of ∼30 mV/m (Santolík et al., 2003). Later the STEREO spacecraft were launched and the electric field instruments were one of the first to be turned on. This led to the discovery of extremely large amplitude whistler mode waves with T200 mV/m (Cattell et al., 2008). The discovery provoked an investigation of Wind observations as it passed through the radiation belt some 60+ times early in its lifetime. The result was a series of papers using Wind and STEREO that all showed consistent observations of large amplitude whistler mode waves with T100 mV/m Breneman et al., 2011;Kellogg et al., 2011;Kersten et al., 2011;Wilson et al., 2011;Breneman et al., 2012). These results altered the design and scientific direction of NASA's Van Allen Probes mission. Similar issues arose in observations of collisionless shock waves. The early work using dynamic spectra data found electrostatic waves with spacecraft frame frequencies, f sc , greater than a few hundred hertz to have amplitudes of, at most, a few 10s of mV/m but typically smaller in the few mV/ m range (Fredricks et al., 1970b;Rodriguez and Gurnett, 1975). Numerous theoretical studies had suggested that small-scale, high frequency waves were an important dissipation mechanism to regulate the nonlinear steepening of collisionless shock waves (Sagdeev, 1966;Tidman and Krall, 1971;Papadopoulos, 1985). However, such small amplitude electric field observations raised doubts about the ability of the high frequency modes to supply sufficient dissipation to maintain a stable shock.
The first published example (of which the authors are aware) of a time series electric field component observed by a spacecraft within a collisionless shock was presented in Wygant et al. (1987) observed by the ISEE-1 probe. The observation was one of the first pieces of evidence that the dynamic spectra plots were not fully capturing the electric field dynamics because the data showed electric fields up to nearly ∼100 mV/m. Later work using the Wind spacecraft found ESWs in the terrestrial bow shock with amplitudes in excess of ∼100 mV/m (Bale et al., 1998;Bale et al., 2002). A few bow shock crossings were observed with the Polar spacecraft, which found nonlinear, electrostatic IAWs within the shock with amplitudes up to ∼80 mV/m (Hull et al., 2006). The picture starting to emerge was that high frequency, electrostatic waves were common and large amplitude in collisionless shocks. Note that the occurrence rate of electrostatic waves was already implied by studies using dynamic spectra data, but not such large amplitude. Wilson et al. (2007) examined waveform capture data of electrostatic waves above the proton cyclotron frequency, f pp , from the Wind spacecraft finding a positive correlation between peak wave amplitude and shock strength, i.e., stronger shocks had larger amplitude waves. They also observed that ion acoustic waves were the dominant electrostatic mode within the shock ramp. Shortly after, a study ) of a supercritical shock showed evidence of waves radiated by the ECDI. Since then, a series of papers using MMS (Chen et al., 2018;Goodrich et al., 2018;Goodrich et al., 2019), STEREO (Breneman et al., 2013), THEMIS (Wilson et al., 2014a;Wilson et al., 2014b), and Wind (Breneman et al., 2013) have examined these electrostatic waves in collisionless shocks.
While the discussion has almost exclusively focused on fluctuating electric fields, δE, it is critical to discuss quasi-static electric fields, E o , as well. The primary obstacle to accurate E o measurements results from the lack of a stable ground in spacecraft observations (Scime et al., 1994a;Scime et al., 1994b;Scudder et al., 2000;Pulupa et al., 2014;Lavraud and Larson, 2016) and the sheath that forms around the conducting surfaces (Ergun et al., 2010), which alters how the instrument couples to the plasma. It is beyond the scope of this study to discuss, in detail, the difficulties associated with such measurements, but some context can be gained by reviewing some recent electric field instrument papers (Wygant et al., 2013;Bale et al., 2016;Ergun et al., 2016;Lindqvist et al., 2016). In lieu of a proper E o measurement in the plasma rest frame, we can estimate the convective electric field, E c -V sw × B o (where V sw is the bulk flow solar wind velocity in the spacecraft frame) 8 . Since the parameters in most simulations are scaled or normalized, we will use the dimensionless ratio δE/E o when comparing spacecraft and simulation results. Unless otherwise specified, E o E c in these contexts.
Prior to the launch of MMS, there were several studies that attempted to measure the cross shock electric field but each suffered from inaccuracies or under resolved electric field measurements which kept the issue of its magnitude in doubt (Dimmock et al., 2011;Dimmock et al., 2012;Wilson et al., 2014a;Wilson et al., 2014b). The launch of MMS allowed researchers, for the first time, to probe E o with sufficient cadence and accuracy to properly measure the cross shock electric field in an interplanetary shock (Cohen et al., 2019). Note that the E o measured in this study peaked at (1.5 mV/m, i.e., comparable to or smaller than the magnitude of E c (which was (4 mV/m in this study). Therefore, we will assume E o as being comparable to E c in magnitude throughout and will just refer to E o instead of both. Even so, there is some discrepancy because such a measurement is extremely difficult at the terrestrial bow shock and detailed MMS observations showed that the electron dynamics seemed to be dominated by a combination of magnetosonic-whistler modes and electrostatic IAWs and ECDI waves (Chen et al., 2018).
The current picture from observations is summarized in the following. In the studies where the quasi-static electric field could be reliably measured (Cohen et al., 2019) or approximated from measurements (Wilson et al., 2014a;Wilson et al., 2014b;Goodrich et al., 2018;Goodrich et al., 2019), the findings were that δE is consistently much larger than E o , i.e., δE ≫ E o . Some of these works attempted to quantify the impact on the dynamics of the system due to δE vs. E o , finding δE dominated (Wilson et al., 2014a;Wilson et al., 2014b;Chen et al., 2018;Goodrich et al., 2018). Chen et al. (Chen et al., 2018) examined in great detail the evolution of the electron distribution through the shock finding that a magnetosonic-whistler wave accelerated the bulk of the incident distribution which rapidly became unstable to high frequency, electrostatic IAWs that scattered the electrons into the often observed flattop distribution [Wilson et al., 2019b;Wilson et al., 2019a;Wilson et al., 2020, and references therein]. This seems to somewhat contradict the results of Cohen et al. (2019) and others who argued that a quasi-static cross shock potential is dominating the shock and particle dynamics. What all of these studies do agree upon is that δE ≫ E o . Note that the purpose of comparing the fluctuating to the quasi-static field here is to help compare with simulations, which normalize the electric fields by the upstream E c value or something similar. Figure 1 shows seven waveform captures observed by the Wind spacecraft's WAVES instrument (Bougeret et al., 1995) while passing through the quasi-perpendicular terrestrial bow shock. The first column shows the x-antenna electric field (δE x ), the second the y-antenna electric field (δE y ), and the third hodograms of δE y vs. δE x . The local B o is projected onto each hodogram shown as a magenta line 9 . Each row shows a different waveform capture/snapshot that is ∼17 ms in duration. The first column contains a double-ended arrow in each panel illustrating the scale associated with 200 mV/m. The first two rows show examples of ESWs mixed with ECDI-driven waves, the third and fourth rows show ECDI-driven waves, and the fifth through seventh rows show IAWs. The distinguishing features are as follows: the ESWs have an isolated, bipolar pulse with either a linear or figure eight-like polarization and a nearly flat frequency (spacecraft frame) spectrum response in the ∼0.2-10 kHz range (not shown); the IAWs exhibit symmetric δE x and δE y about zero, are linearly polarized along B o , and have a broad frequency peak (spacecraft frame) in the ∼2-10 kHz range (not shown); and the ECDI exhibit asymmetric δE x and δE y fluctuations about zero, their polarization is not always linear along B o , and the frequency peak (spacecraft frame) is in the ∼0.5-10 kHz range with superposed cyclotron harmonics (not shown). For reference, the upstream average convective for this bow shock crossing is E c ∼ 3.5 mV/m.

Kinetic Simulations
Kinetic simulations of shocks are challenging due to the need to resolve global structure of the shock (generally associated with λ i 10 ) and the relatively long time scales associated with it simultaneously with short spatial (λ De ) and fast temporal scales (f pe ) associated with instabilities.
Early kinetic particle-in-cell (PIC) simulations were much more limited by computational constraints than those performed today. A common approach to scaling the problem in order reduce the computational load is to consider one-or twodimensional problems and to reduce ratios of the ion-to-electron mass, Mi me , and electron plasma-to-cyclotron frequency, ωpe Ωce , while keeping the plasma β and the size of the problem in units of λ i comparable to the physical system of interest. Further computational trade-offs include altering the simulation resolution (i.e., number of grid cells), the number of particles per cell for particle codes or velocity-space resolution for continuum Vlasov codes (Yang et al., 2013). Since the frequencies and the growth rates of the instabilities of interest are associated with certain characteristic time scales, such a rescaling may significantly alter the development and the role of instabilities in the simulations. For example, reducing M i me lowers the threshold for Buneman instability (Hoshino and Shimada, 2002) by reducing the difference between electron and ion thermal speeds. Values of ωpe Ωce were also expected and found to inhibit the growth of certain wave modes like Bernstein modes (Matsukiyo and Scholer, 2006;Muschietti and Lembège, 2013;Muschietti and Lembège, 2017). What's more, the M i m e ratio was shown to dramatically affect the macroscopic profile of the shock magnetic field (Scholer and Matsukiyo, 2004) and affect the growth of what are now viewed as critical instabilities like the MTSI (Umeda et al., 2012a;Umeda et al., 2012b;Umeda et al., 2014). Thus the re-scaling approach must be carefully chosen based upon its expected impact on the phenomena of interest.
Some of the first two-dimensional PIC simulations using realistic Mi me was presented by (Matsukiyo and Scholer, 2006). Since then, the community has made efforts to compromise somewhat on Mi m e in order to increase ω pe Ω ce , to more realistic values (i.e., 50-100 in solar wind near 1 AU) (Muschietti and Lembège, 2013) used ratios of M i m e 400 and ω pe Ω ce 10 to examine the higher harmonics of Bernstein modes associated with the ECDI.
More typical values for the latter fall in the ∼2-4 range for recent simulations (Umeda et al., 2014;Matsukiyo and Matsumoto, 2015;Zeković, 2019). However, much larger values have been used in cases where one can reduce the simulation to one spatial dimension (Umeda et al., 2019). Despite all of the progress made since the early full PIC simulations, there still remains two striking discrepancies between observations and many simulations: the amplitude and wavelength at which the strongest electric fields are observed and inconsistencies in the thickness of the shock ramp. The second issue is more obvious from cursory examinations of simulation results, so we will discuss it first. As previously discussed, observations consistently show that the shock ramp thickness, L sh , tends to satisfy 5 < L sh /λ e < 40 (Hobara et al., 2010;Mazelle et al., 2010). However, PIC simulations, even with realistic mass ratio, often generate FIGURE 1 | Example of six waveform captures observed by Wind while crossing the Earth's quasi-perpendicular bow shock on 1996-12-02. Each row corresponds to, by column, the X vs time, Y vs time, and Y vs X hodogram (in instrument coordinates) from the Wind/WAVES time domain sampler (TDS) instrument (Bougeret et al., 1995). Each waveform capture/snapshot is ∼17 ms in duration. The quasi-static magnetic field projection is shown as the magenta line in hodograms (Adapted from Figure 3.6 in . Frontiers in Astronomy and Space Sciences | www.frontiersin.org January 2021 | Volume 7 | Article 592634 6 shock ramps with thicknesses satisfying L sh /λ e > 43, i.e., exceeding proton inertial length (Scholer and Burgess, 2006), while some generate more realistically thin ramps (Matsukiyo and Scholer, 2012;Yang et al., 2013). Yang et al. (Yang et al., 2013) concluded that the shock ramp thickness decreased with increasing Mi me but increased with increasing ion plasma beta. Note however that (Matsukiyo & Scholer, 2012) used ∼20% finer grid resolution, twice as many particles per cell, and smaller plasma betas than (Scholer & Burgess, 2006). However, (Yang et al., 2013) used fewer particles per cell and smaller ωpe Ωce than both (Matsukiyo and Scholer, 2012) and (Scholer and Burgess, 2006). It is important to note that it's still not clear what physical or numerical parameters controls the ramp thickness in simulations or observations or even what the relevant physical scale is (e.g. λ e or λ De ).
Note that the thickness of the magnetic ramp of a collisionless shock is not significantly affected by the presence of corrugation/ripples (Johlander et al., 2016) other than the temporal dependence that can occur during reformation (Mazelle et al., 2010). The spatial extent of the entire transition region can indeed be increased by such oscillations but the magnetic gradient scale length does not appear to be significantly affected. The biggest limiation to determining the shock ramp thickness in data is time resolution.
More recent spacecraft like THEMIS (Angelopoulos, 2008) and MMS (Burch et al., 2016) have, for instance, fluxgate magnetometers that return 3-vector components 128 times every second, which is more than sufficient to resolve the bow shock ramp. The bow shock moves slower in the spacecraft frame than interplanetary shocks, so time resolution is more of a constraint for examining the shock ramp thickness of interplanetary shocks. Even so, the 128 sps of the THEMIS and MMS fluxgate magnetometers is still sufficient for most interplanetary shocks.
As previously discussed, observations consistently show δE/E o > 50 for fluctuations with wavelengths at or below a few 10s of Debye lengths (Wilson et al., 2014a;Wilson et al., 2014b;Chen et al., 2018;Goodrich et al., 2018), i.e., λ ( few 10s of λ De . Most shock simulations find values satisfying δE/E o < 10 and the scales at which the largest electric field fluctuations occur tend to satisfy kλ e < 1 (Scholer and Matsukiyo, 2004;Matsukiyo and Scholer, 2006;Scholer and Burgess, 2007;Lembège et al., 2009;Umeda et al., 2012a;Umeda et al., 2014;Matsukiyo and Matsumoto, 2015). We note that explicit fully kinetic PIC simulations tend to have spatial grid resolution of a few λ De , since such scales must be resolved for numerical stability. It has long been known that unrealistically small values of M i m e and ω pe Ω ce can lead to unrealistically large electric field amplitudes for modes with kλ e < 1 (Hoshino and Shimada, 2002;Comişel et al., 2011;Zeković, 2019). Although the three main modes discussed herein have been successfully identified in PIC simulations, they were either unrealistically small in amplitude or at different spatial scales or not excited unless the simulation was specifically tailored for that instability.
It is worth noting the severe computational costs of using fully realistic plasma parameters. The separation of spatial scales satisfies λ i /λ De Ω ce , which is again ∼ 4,000-20,000 11 . The computational cost of any given simulation scales as , where d is the number of spatial dimensions used in the simulation. Thus, one can see that increasing ω pe Ω ce from ∼10 to 100, even in a one-dimensional simulation, is at least 100 times more computationally expensive. It is also the case that simulations often use shock speeds satisfying V Tp < V Te < U shn while shocks in the solar wind tend to satisfy V Tp < U shn ≪ V Te ≪ c (see Section 3 for values and definitions). For explicit PIC codes, there are additional computational expenses since the time steps are tied to the grid cell size, which raises the order of both ωpe Ωce and Mi me by one. Therefore, it can be seen that we are approaching a computational wall and it may require new classes of simulation codes to overcome these limitations if we hope to use fully realistic plasma parameters. shows the magnetic field magnitude (black) and three normal incidence frame coordinate basis (NCB) (Scudder et al., 1986a) components vs time at 128 samples/second (sps) observed by the fluxgate magnetometer (Auster et al., 2008). Panel (B) shows the DC-coupled electric field at 128 sps observed by EFI (Bonnell et al., 2008;Cully et al., 2008). Panel (C) shows the AC-coupled electric field at 8,192 sps observed by EFI. Note that Panels (B,C) share the same vertical axis range for direct comparison (Adapted from Figure I:2 in Wilson et al. (2014a)).
11 Note that ω pe /Ω ci is larger by an additional factor of M i /m e √ .

EXAMPLE OBSERVATIONS VERSUS SIMULATIONS
In this section we will present two example observations made by the THEMIS (Angelopoulos, 2008) and MMS (Burch et al., 2016) missions to further illustrate the difference in magnitude between δE and E o . We will also present PIC simulation results with parameters representative of a wide class of simulations discussed in the literature. The purpose is to illustrate some limitations of simulations to provoke advancement in closing the gap between observations and simulations of collisionless shocks. Figure 2 shows a direct comparison between δE and E o observed by THEMIS-C during a terrestrial bow shock crossing adapted from Wilson III et al. (2014a) and Wilson III et al. (2014b). The study examined the energy dissipation rates estimated from (J · E) (i.e., from Poynting's theorem) due to the observed electric fields, E, and estimated current densities 12 , J. They expanded (j o + δj) · (E o + δE) and found that (J o · δE) was the dominant term 13 , i.e., the fluctuating fields acted to limit the low frequency currents in/around the shock. Two important things were found: the magnitude of (j o · δE) and changes in this term were much larger than (j o · E o ); the signs of the changes in (j o · δE) and (j o · E o ) are opposite. The second point was interpreted to imply that the fluctuating fields were giving energy to the particles and the quasi-static fields were gaining energy from the particles. In short, the main conclusion from (Wilson III et al., 2014a;Wilson III et al., 2014b) was to illustrate that not only are the fluctuating electric fields large, they could potentially contribute enough energy transfer to compete with quasi-static fields. Prior to this study, the view by many in the community was that these fluctuating fields were completely negligible or just a minor, secondary effect. More recent, independent studies have performed similar analyses using different spacecraft and came to similar conclusions (Chen et al., 2018;Goodrich et al., 2018;Hull et al., 2020). Figure 3 provides another example directly comparing δE and E o observed by two MMS spacecraft during a terrestrial bow shock crossing. The electric fields are shown in the De-spun, Sun-pointing, L-vector system or DSL (Angelopoulos, 2008). For each spacecraft, E o,j ( 10 mV/m was satisfied for the entire interval with most time steps satisfying E o,j ( 5 mV/m. In contrast, the peak-to-peak δE j values commonly exceed 100 mV/m in bursty, short duration, wave packet intervals. Note that even the electric field instrument on MMS has limitations in its accuracy for frequencies below ∼1 Hz (Ergun et al., 2016;Lindqvist et al., 2016). Thus, even with the significantly improved instrument technology and design of MMS, the observations consistently show δE ≫ E o .
As a practical list of reference values, we present one-variable statistics of solar wind parameters from the same data set as in Wilson et al. (2018)

and all interplanetary (IP) shocks in the Harvard
Smithsonian's Center for Astrophysics Wind shock database 14 . The following will show parameters as X 5% ( X ( X 95% ,X [units], for the entire data set, where X y% is the y th percentile andX is the median. First, the typical parameters for over 400 IP shocks are as follows:   (Russell et al., 2016). Panels (B,E) show the DC-coupled electric field components at 128 sps observed by the electric field instrument (Ergun et al., 2016;Lindqvist et al., 2016). Panels (C,F) show the AC-coupled electric field at 8,192 sps. Note that Panels (B,C,E) share the same vertical axis range for direct comparison.
12 Note that similar current densities have been found using multi-spacecraft techniques (Hull et al., 2020) supporting the results in Wilson III et al. (2014a) and Wilson III et al. (2014b). 13 Note that Q o in this context is not the quasi-static terms in quasi-linear or linear theory but that from the DC-coupled measurements. Further, δQ is the fluctuating terms from these theories but the AC-coupled measurements, thus there is no a priori requirement that 〈δQ〉 0. 14 https://www.cfa.harvard.edu/shocks/wi_data/.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org January 2021 | Volume 7 | Article 592634 8 where U shn is the upstream average flow speed in the shock rest frame and V shn is the upstream average shock speed in the spacecraft frame. Note that the values of M f , M A , and U shn will all be, on average, larger for most bow shocks in the interplanetary medium. These values are only meant to serve as a statistical baseline for reference. For example, the 11 bow shock crossings in Wilson et al. (2014a); Wilson et al. (2014b)  where f cs is the cyclotron frequency of species s, f ps is the plasma frequency of species s, V Ts is the most probable thermal speed of species s 16 , ρ cs is the thermal gyroradius of species s 17 , λ e is the inertial length of species s, and λ De is the electron Debye length. 15 Note that none of these are Doppler-shifted. 16  where 〈U shn 〉 y% is the y th percentile of U shn presented earlier in this section and c is the speed of light in vacuum. Figure 4 shows example one-dimensional cuts at three different time steps taken from a PIC simulation. The simulation parameters are as follows: θ Bn ∼ 60 deg, M A ∼ 6.5, ωpe Ωce ∼ 4, Mi me ∼ 900, Δ ∼ 1 λ De (where Δ is the grid cell size), initially 1,000 particles per cell, and λ e /λ De ∼ 8 (i.e., ∼27 times smaller than median solar wind values near 1 AU). All of the panels show data in normalized units. The electric field is normalized to the initial upstream averaged convective electric field, -V×B, i.e., the same E o referenced for spacecraft observations. Thus, in the upstream the E z component has an offset of unity. The normalization for n e and B are just the initial upstream average values of the magnitude of each. All fields are measured in the simulation frame, where the shock moves in the positive x direction with the speed of approximately 2 V A .
One can see that the largest values of |E| rarely exceed 2 (i.e., only short intervals >2 but peak only at ∼6), similar to the simulations discussed previously. Further, the spatial scales at which these fields are maximized is on λ e scales whereas observations show maximum electric fields at λ De scales. One can also see that the ramp, e.g., B in panel (D), is about L sh ∼ 28 λ e (or ∼0.9 λ i ) thick, similar to observations that typically show ramps satisfying L sh < 35 λ e (or <0.8 λ i ) (Hobara et al., 2010;Mazelle et al., 2010). The simulation does, however, generate the ubiquitous whistler precursor train upstream of the shock ramp Wilson et al., 2017). Yet it is still not clear what parameter or parameters are controlling the shock ramp thickness and electric field amplitudes at very small spatial scales in simulations.

DISCUSSION
We have presented examples illustrating that spacecraft observations of collisionless shocks consistently show δE ≫ E o where δE is due to electrostatic fluctuations satisfying k λ De ( 1 with frequencies well above f lh . In contrast, most PIC simulations of collisionless shocks show considerably smaller amplitude of electrostatic fluctuations. This is true even when the simulation uses realistic M i m e and plasma betas. Further, many simulations still generate shock ramps satisfying L sh /λ e > 43, i.e., thicker than most observations. However, much more progress has been made on this front where Yang et al. (2013) concluded that the shock ramp thickness decreased with increasing Mi me but increased with increasing ion plasma beta. There is still the question of whether ω pe Ω ce plays a role in the simulated values of L sh , though Yang et al. (2013) was able to produce realistic thicknesses despite only having ω pe Ω ce 2. Another potential issue that was not explicitly discussed in detail is that of the separation between λ e and λ De , but these are controlled by M i m e and ω pe Ω ce . As previously shown, statistical solar wind parameters satisfy λ e /λ De ∼ 215 (or λ p /λ De ∼ 9,757) while simulations often have much smaller vales of λ p /λ De ∼ 70-500 (or λ e /λ De ∼ 7-40) (Umeda et al., 2011;Umeda et al., 2012a;Umeda et al., 2014;Savoini and Lembège, 2015). It is also the case that simulations often use shock speeds satisfying V Tp < V Te < U shn < c while shocks in the solar wind tend to satisfy V Tp < U shn ≪ V Te ≪ c.
The origins of the discrepancy between the observation that δE ≫ E o for electrostatic fluctuations satisfying k λ De ( 1 remain unclear. The ratios M i m e and ω pe Ω ce are the most likely parameters since they control the separation of spatial and temporal scales between the instabilities of interest and the global shock scales in an obvious manner. A lack of spatial resolution in most simulations may also be a factor. The purpose of this work is to motivate both observational and simulation communities to bridge the gap find closure with this issue. Without an accurate reproduction of the high frequency, large amplitude waves it is not possible to determine at what scales the electric fields dominate the energy dissipation through collisionless shocks.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://cdaweb.gsfc.nasa.gov.

AUTHOR CONTRIBUTIONS
LW wrote most of the content and generated all the observational data figures presented herein. L-JC provided critical contributions to the bridge between observations and simulations. VR provided critical contributions to simulation techniques and limitations induced by the variation of different normalized parameters.

FUNDING
The work was supported by the International Space Science Institute's (ISSI) International Teams program. LW was partially supported by Wind MO&DA grants and a Heliophysics Innovation Fund (HIF) grant. Contributions of VR were supported by NASA grant 80NSSC18K1649. Computational resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.