Electron energy partition across interplanetary shocks: III. Analysis

Analysis of model fit results of 15,210 electron velocity distribution functions (VDFs), observed within $\pm$2 hours of 52 interplanetary (IP) shocks by the Wind spacecraft near 1 AU, is presented as the third and final part on electron VDFs near IP shocks. The core electrons and protons dominate in the magnitude and change in the partial-to-total thermal pressure ratio, with the core electrons often gaining as much or more than the protons. Only a moderate positive correlation is observed between the electron temperature and the kinetic energy change across the shock, while weaker, if any, correlations were found with any other macroscopic shock parameter. No VDF parameter correlated with the shock normal angle. The electron VDF evolves from a narrowly peaked core with flaring suprathermal tails in the upstream to either a slightly hotter core with steeper tails or much hotter flattop core with even steeper tails downstream of the weaker and strongest shocks, respectively. Both quasi-static and fluctuating fields are examined as possible mechanisms modifying the VDF but neither is sufficient alone. For instance, flattop VDFs can be generated by nonlinear ion acoustic wave stochastic acceleration (i.e., inelastic collisions) while other work suggested they result from the combination of quasi-static and fluctuating fields. This three-part study shows that not only are these systems not thermodynamic in nature, even kinetic models may require modification to include things like inelastic collision operators to properly model electron VDF evolution across shocks or in the solar wind.

, cometary bow shocks (e.g., Sagdeev et al. 1987), and from nonlinearly steepened electromagnetic waves radiated by instabilities (e.g., Wilson III et al. 2009. Even though the mean free path of a typical proton near 1 AU can be nearly 1 AU, i.e., orders of magnitude larger than the corresponding thermal gyroradii (ρcp) or inertial length (λp), the thickness of shock ramps -the spatial gradient scale length of the magnetic transition region -in the solar wind are often a few λe up to λp (e.g., Hobara et al. 2010;Mazelle et al. 2010). This is why shock waves in the solar wind, and in most astrophysical contexts, are called collisionless.
There are several key unresolved questions about the microphysical processes that regulate the dynamics of collisionless shock waves. One of the biggest outstanding problems in shock physics is the partition of energy between electrons and ions. A significant obstacle is that the mechanisms depositing/transferring energy are not predicted to act homogeneously, i.e., they are energy-and pitch-angle-dependent and can be speciesdependent (e.g., Artemyev et al. 2013Artemyev et al. , 2014Artemyev et al. , 2015Artemyev et al. , 2016Artemyev et al. , 2017aArtemyev et al. ,b, 2018Sagdeev 1966). Evidence supporting this prediction have been reported in some case study observations at IP shocks (e.g., Wilson III et al. 2009 and the terrestrial bow shock (e.g., Chen et al. 2018;Goodrich et al. 2018aGoodrich et al. , 2019Oka et al. 2017Oka et al. , 2019Wilson III et al. 2014a,b). Further, most collisionless shocks are subsonic to electrons, yet electrons still behave as if they have experienced a shock, even showing Mach number dependent effects (e.g., Feldman et al. 1982Feldman et al. , 1983cMasters et al. 2011;Thomsen et al. 1985Thomsen et al. , 1987Thomsen et al. , 1993Wilson III et al. 2010). The question remains, how does a collisionless shock transform the bulk flow kinetic energy into other forms like electron and/or ion heating.
To determine the energy transfer mechanism, several previous studies examined the change in average electron temperature, ∆T e,tot (see Appendix A for definitions), or the ratio of downstream-to-upstream average electron temperature,R T e,tot , across the shock. These parameters were compared with the upstream average fast Mach number, M f up, the upstream average ki-netic energy, |KE shn | up , the change in the shock normal speed, ∆Ū shn , and the square of the shock normal speed, ∆Ū shn 2 and/or the change in kinetic energy across the shock 1 . Feldman et al. (1983a) found weak, positive correlation betweenR T e,tot and M f up (technically it was with M ms up ). Feldman et al. (1983b) found weak, positive linear relationships between both ∆T e,tot and ∆T p,tot and ∆Ū shn . Thomsen et al. (1987) examined electron heating at Earth's bow shock finding a correlation between ∆T e,tot and ∆Ū shn 2 with a few to ∼10% of the kinetic energy transforming into electron heating. Later Schwartz et al. (1988) found a correlation between ∆T e,tot and ∆Ū shn 2 . They also show a strong, positive correlation for ∆T e,tot versus ∆T i,tot with a slope of ∼0.2 and y-intercept of ∼8.6 eV. Thomsen et al. (1993) examined ion and electron heating at the low Mach number, quasi-parallel Earth's bow shock finding that ∼6% of ∆KE shn was converted to electron heating. Hull et al. (2000) found a positive, linear relationship between ∆T e,tot and ∆KE shn with a slope of 0.057 +0.0041 −0.0039 , i.e., ∼6% of of ∆KE shn was converted to electron heating. Fitzenreiter et al. (2003) also found a linear relationship between ∆T e,tot and ∆Ū shn 2 . Finally, Masters et al. (2011) examined electron heating at the Kronian bow shock finding that for most crossings ∼3%-7% of KE shn up was converted to electron heating. Thus, previous work suggests that the change in bulk flow kinetic energy should be correlated with the increase in temperature. Note that previous work has either inferred (e.g., Ghavamian et al. 2007Ghavamian et al. , 2014 or showed with in situ measurements that stronger shocks heat ions more than electrons (e.g., Schwartz et al. 1988; Thomsen et al. 1993;Masters et al. 2011) with the fraction of electron heating varying with an inverse Mach number 2 . The inverse Mach number dependence was also found to be a piecewise distribution where ∼50% of the total heating goes to electrons at low Mach number and changes to 10% at higher Mach number. In this third part (Paper III) of this three-part study, the analysis of the fit results to the multi-component electron VDF analysis will be discussed in the context of the macroscopic shock parameters and some instability analysis. The results are summarized for the 52 IP shocks observed by the Wind spacecraft. The nota-tion, symbols, and data sets used herein are the same as those in Wilson III et al. (2019a) (hereafter referred to as Paper I) and Wilson III et al. (2019b) (hereafter referred to as Paper II). Paper I discussed the methodology and described the data product resulting from the application of the fit software and Paper II summarized the statistics of the fit results. This work will provide the physical analysis and interpretation of the results in the context of macroscopic shock parameters and some instability analysis.

DATA SETS AND METHODOLOGY
As in Papers I and II, all data are observed by instruments on the Wind spacecraft (Harten & Clark 1995) near 1 AU. The data used herein include quasi-static magnetic field vectors (Bo) from Wind /MFI (Lepping et al. 1995), electron and ion velocity distribution functions (VDFs) from Wind /3DP (Lin et al. 1995), and proton and alpha-particle velocity moments from the Wind /SWE Faraday Cups (Kasper et al. 2006;Ogilvie et al. 1995). The instrument details are described in Paper I. Parameters described with respect to Bo are in a field-aligned coordinate basis using a subscript j to denote the parallel (j = ), the perpendicular (j = ⊥), and total (j = tot) directions. All electron parameters are shown with a subscript s denoting the component (or sub-population) of the entire distribution where s = ec for the core, s = eh for the halo, s = eb for the beam/strahl, and s = e for the entire distribution. The combined or mixed parameters (e.g., β ef f,j ) use the subscripts s = ef f for effective and s = int for integrated parameters (see Appendix A for definitions). The analysis of the VDFs presented herein were found within ±2 hours of 52 IP shocks found in the Wind shock database from the Harvard Smithsonian Center for Astrophysics 3 between 1995-02-26 and 2000-02-20 (for full list of event dates and times, see PDF file included with additional supplemental material found at https://doi.org/10.5281/zenodo. 2875806 Wilson III et al. (2019c)). The IP shocks examined were limited to fast-forward shocks that had burst mode electron VDFs within the chosen time range about each shock.
As in Paper II, the VDF results presented herein relative to all 15,210 VDFs examined. The VDF fit results are taken from additional supplemental material in the form of two ASCII files (discussed in Paper II and described in Paper I) (Wilson III et al. 2019c). Of this total 14,418 had stable model fits (f (core) ) for the core, 13,660 stable model fits (f (halo) ) for the halo, and 3 https://www.cfa.harvard.edu/shocks/wi data/ 11,578 stable model fits (f (beam) ) for the beam/strahl. Note that all statistics presented herein are for stable fits with a fit flag for the respective component of two or higher. The selection justification of the VDFs are provided in Paper II. Similarly, this work follows Paper II with the following selection criteria: Criteria AT: All VDFs satisfying: Fit Flag {c, h, b} ≥ 2 and no violation of post-fit constraints; Criteria UP: All VDFs satisfying Criteria AT that were observed upstream of the IP shock ramp; Criteria DN: All VDFs satisfying Criteria AT that were observed downstream of the IP shock ramp; Criteria LM: All VDFs satisfying Criteria AT that were observed near IP shocks satisfying M f up < 3; Criteria HM: All VDFs satisfying Criteria AT that were observed near IP shocks satisfying M f up ≥ 3; Criteria PE: All VDFs satisfying Criteria AT that were observed near IP shocks satisfying θ Bn > 45 • ; and Criteria PA: All VDFs satisfying Criteria AT that were observed near IP shocks satisfying θ Bn ≤ 45 • . Superposed epoch analysis (SEA) of every fit parameter has been performed where the time stamps of each VDF are redefined as offsets from the associated IP shock ramp center time. Finally, to quantify the partition of energy for each component, all permutations of the difference (and ratio) between the upstream and downstream values are computed. Again, the value used for each shock is the median and the uncertainty is half the magnitude of the difference between X 5% and X 95% . The results were fit to model functions (e.g., power-law) with the dependent variable/abscissa being the different macroscopic shock parameters. All permutations are used to avoid the subjectivity introduced by selecting a user-defined time period that serves as the upstream and downstream regions. The median was used instead of an average because again the data are normally distributed and the median is less influenced by large deviations/tails. The differences or ratios are then fit to a model function using both Poisson and Gaussian weights to find the best fit results.
Finally, it should be noted that the total time range (i.e., ±2 hours) about each shock ramp is larger than the typical time range examined (i.e., ±10s of minutes). The purpose is ensure that the upstream analysis includes data that is not part of the electron foreshock. The extent of the downstream is defined purely to maintain symmetry about the shock ramp center. Although the analysis presented herein only discusses differences and analysis of the VDF solutions for the entire time range, the same analysis was performed on limited time ranges about the shock ramp for comparison (not shown). Alteration of the time range about the shock ramp did not yield significant differences for most parameters examined.

SUPERPOSED EPOCH ANALYSIS
In the following subsections multiple figures presenting superposed epoch analysis (SEA) of the various fit parameters are presented. In all sections, the subscripts s = ec is for the core, s = eh for the halo, s = eb for the beam/strahl, and s = ef f for the effective components.
The time stamps of each VDF are redefined as offsets from the associated IP shock ramp center time. The data are partitioned into 120 second time windows and then one-variable statistical analysis is performed on all data within. The time window center is shifted by 22 seconds (i.e., roughly the sample period of the 3DP instrument in survey mode) and then the analysis is repeated. In this way, the data can remain at their original sample time without affecting the magnitudes through the use of a re-gridding algorithm by averaging or interpolating data to a common set of time stamps.
For every partitioned time window,X is used as the center line (shown as a red line in Figures 1 and 2), X 5% as the lower line (shown as the lower cyan line in Figures 1 and 2), and X 95% as the upper line (shown as the upper cyan line in Figures 1 and 2). The median and percentile lines are then smoothed by selection criteria-dependent widths 4 , with Criteria AT smoothed over 10pts, Criteria LM smoothed over 20pts, Criteria HM smoothed over 30pts, Criteria PA smoothed over 30pts, and Criteria PE smoothed over 20pts. The data are also normalized by the upstream median values for each IP shock (see Tables 2-6 of additional supplemental material found in Wilson III et al. (2020) for list of values) to remove relative offsets between any two IP shocks.
The median and percentiles are used as the data are not normally distributed about a mean for each partitioned time window, as implied by the non-Gaussian histogram distributions shown in Paper II. The SEA plots provide the trend and typical ranges of the various parameters relative to the shock ramp center to illustrate changes. The parameters are also separated s ec vs.  Figure 1. Superposed epoch analysis plot of the core, sec, halo, κ eh (middle column), and beam/strahl, κ eb (right column), exponents separated into Criteria AT (first row), Criteria LM (second row), Criteria HM (third row), Criteria PA (fourth row), Criteria PE (fifth row). The red line shows the smoothedX values of the partitioned data. The lower(upper) cyan lines are the smoothed X 5% (X 95% ) values of the partitioned data. The vertical green line indicates the ramp center time. All data for a given IP shock are normalized by an upstream median value given in Table 5 found in Wilson III et al. (2020). Note that all panels share the same horizontal axis range but the first column has a different vertical range than the latter two.
by the above selection criteria to illustrate macroscopic shock parameter dependence.

Exponents
In this section, SEAs of sec/ sec up , κ eh / κ eh up , and κ eb / κ eb up are introduced and discussed in Figure 1 (see Appendix A for definitions). First, 50% of the halo and beam/strahl exponents (i.e., between the lower and upper quartiles) satisfy 3.57 κ eh 5.31 and 3.41 κ eb 5.11, which are consistent with previous solar wind observations near 1 AU (e.g., Horaites et al. 2018;Lazar et al. 2017;Maksimovic et al. 1997Maksimovic et al. , 2005Pierrard et al. 2016;Štverák et al. 2009;Tao et al. 2016a,b) (there are no studies with which to compare the sec values).
The sec exponents are relatively constant across most IP shocks, with a few exceptions for Criteria LM and Criteria PA shocks showing values up near ∼1.5. The Criteria HM are especially sparsely populated in the downstream as many of these VDFs were fit to the asym-metric self-similar model distribution (i.e., pec and qec exponents, not shown). However, the range of data between X 5% and X 95% is larger in the downstream, consistent with the interpretation of a wave-induced inelastic interaction discussed in Paper I.
The κ eh panels (middle column) are more interesting and varied. For instance, the running median clearly increases across the shock. The change is the largest for Criteria HM and Criteria PE shocks. Theory and simulation suggest that stronger and more oblique shocks are better at energizing suprathermal electrons (e.g., Caprioli & Spitkovsky 2014;Park et al. 2015;Treumann 2009;Trotta & Burgess 2019). However, if the normalization used in Figure 1 is removed (not shown), the upstream only values of κ eh between Xmin and Xmax cover a similar range for Criteria LM and Criteria HM shocks, with similar running medians as well. The downstream only κ eh values are larger at Criteria HM than Criteria LM shocks, which explains the larger one-variable statistics values of κ eh for Criteria HM than Criteria LM shocks shown in Paper I.
Note that the one-variable statistics of κ eh reported in Paper I showed a slightly larger value for Criteria PE than Criteria PA shocks in agreement with the larger change across Criteria PE than Criteria PA shocks observed in Figure 1, consistent with theory and simulations (e.g., Wu 1984;Park et al. 2013;Trotta & Burgess 2019). While the change is not in the direction expected, the magnitude of the change is certainly larger for Criteria HM than Criteria LM shocks suggesting stronger shocks have more of an impact on the halo than weaker. While the upstream kappa values have recently been predicted to increase with increasing Mach number for quasi-perpendicular shocks (Trotta & Burgess 2019), the normalization by the upstream median values for each shock removed the relative offsets to avoid this potentially misleading difference.
The increase in κ eh across the shocks could result from several effects. For instance, it could result from core electrons being energized to suprathermal energies, initially starting with a larger equivalent kappa value, thus increasing the overall suprathermal exponent value. The increase in κ eh could also result from acceleration through a quasi-static electric field (e.g., Schwartz et al. 1988Schwartz et al. , 2011Schwartz 2014;Scudder et al. 1986) for a one-dimensional VDF. The change in κ eh for a two-dimensional VDF in velocity space is complicated if the quasi-static electric field is not uniformly aligned with both component directions in velocity space, as a kappa VDF is not a simple powerlaw (e.g., Livadiotis 2015). Even so, the effect of a quasi-static electric field could result in the fit software finding a larger kappa value as the higher energy particles change velocity by a smaller fraction than the lower. However, the nonlinear relationship between the values of κs, ns, and V T s,j in affecting the shape and peak phase space density of a bi-kappa VDF (e.g., see Figure  2 of Paper I) makes it difficult to interpret the change in κ eh .
The κ eb panels (right-hand column) are less dynamic than the κ eh panels but there is important similarity. The change in the running median is largest for Criteria HM shocks, with all other selection criteria showing little-to-no change or a weak gradual change across the entire time range. The biggest difference between κ eh and κ eb is that the latter does not show a clear increase at the shock ramp for Criteria PE shocks. Similar to the discussion above for κ eh , if the normalization is removed from Figure 1 for κ eb , the difference between Criteria LM and Criteria HM shocks is mostly in the variation in the running median. For Criteria LM shocks, the running median of κ eb is very smooth and does not change much at the shock ramp (i.e., changes by less than ∼0.5) while the running median for Criteria HM shocks fluctuates with a normalized amplitude at or above unity. Again, the change in the running median does not show a strong change at the shock ramp in the unnormalized SEA (not shown).

Density Ratios
In this section, superposed epoch analyses (SEAs) of ns/n ef f are introduced and discussed, for the core (s = ec), halo (s = eh), beam/strahl (s = eb), and effective (s = ef f ). The halo-to-effective density ratios are presented in the first column of Figure 2, the halo-toeffective ratios in the second column, and the beam-toeffective ratios in the third column. Figure 2 shows the ratio of each electron component number density to the effective number density, n ef f (see Appendix A for definitions), to illustrate the relative change in the fractional composition of the electron distribution. That is, the SEA plots show whether the core, halo, and/or beam/strahl number density fraction increase or decrease across the shock for each selection criteria. It's quite obvious that the median core fraction always increases across the shock, regardless of selection criteria. The difference between selection criteria is that magnitude of the change, e.g., Criteria PA shows very little positive change in the running median and a skewness toward smaller values in the running 5 th percentile.
Unlike the core, both the halo and beam/strahl fractional densities decrease across the shock for all selection criteria. Again, the weakest change is found in the Criteria PA shocks, which is not surprising as quasi-parallel  -3600 0 +3600 -3600 0 +3600 Figure 2. Superposed epoch analysis plot of the electron component density ratios nec/n ef f (left column), n eh /n ef f (middle column), and n eb /n ef f (right column). The normalization values for each IP shock are given in Table 3 found in Wilson III et al. (2020). The format is the same as Figure  1.
shocks show smaller magnetic field compression ratios. Thus, the change in the core density is smaller which impacts the halo and beam/strahl fits. The largest change across the shock ramp occurred for Criteria LM and Criteria PE. It should also be noted that the gradient in the running median is sharpest for the beam/strahl component. The halo shows a short enhancement upstream of the shock ramp and then a decrease that continues shortly into the downstream past the ramp. The width of the upstream enhancement is larger for the Criteria HM than Criteria LM shocks, suggesting it results from some form of shock acceleration forming an electron foreshock. This is further evidenced by the stark difference in the running median profile between the Criteria PA and Criteria PE shocks, where the latter are expected to be better at energizing electrons (e.g., Wu 1984;Park et al. 2013;Trotta & Burgess 2019). The beam/strahl gradients occur almost entirely in the thin region focused on the ramp. On the time scale shown in these figures, one minor tick mark is ∼180 seconds so the halo gradient from upstream to downstream lasts upwards of ∼15 minutes. The beam/strahl gradient is much shorter at ∼6 minutes (except for Criteria PA and Criteria HM shocks which have a much gentler gradient). Recall the time windows for calculating the running median are 120 seconds and the windows β ec,// vs.  The normalization values for each IP shock are given in Table  4 found in Wilson III et al. (2020). The format is similar Figure 1 but all three columns share the same vertical axis scale.
are shifted by 22 seconds each time, thus the sharpest gradient one might expect would be a little over two minutes. These plots also show that the fraction of suprathermal particles decreases across most IP shocks, which may explain why the suprathermal exponents tend to increase across the shocks. Although the suprathermal temperatures tend toward larger values in the downstream regions (see additional SEA plots of T s,j found in Wilson III et al. (2020)), the change is very weak except for Criteria HM shocks. However, the core temperature changes across the shocks are much more dramatic than both the halo and beam/strahl. Further, the beam/strahl shows the weakest changes in T s,j for all selection criteria, i.e., the beam/strahl component's mean kinetic energy is not strongly affected by the shock.

Betas
In this section, superposed epoch analyses (SEAs) of the electron plasma betas, βs,j, for the core (s = ec), halo (s = eh), and beam/strahl (s = eb) components. The core betas are presented in Figure 3, the halo temperatures in Figure 4, and the beam/strahl temperatures in Figure 5.   The normalization values for each IP shock are given in Table  4 found in Wilson III et al. (2020). The format is the same as Figure 3.
The next thing to examine is the plasma betas, βs,j, for the core (s = ec), halo (s = eh), and beam/strahl (s = eb) components. The SEA plots in Figures 3, 4, and 5 of βs,j show the largest change across the shock compared to all other electron fit parameters examined. Of the three components, the core shows the smallest change, but all three electron components decrease across the shock. The weak decrease in βec,j, compared to β eh,j and β eb,j , is dominated by the increase in the magnitude of Bo across the shock since both T ec,j and nec increase. Although both T eh,j and T eb,j tend to increase across the shock, the increase is weak. This weak increase coupled with the stronger decrease in both the fractional n eh and n eb and increase in the magnitude of Bo makes for the corresponding decrease in suprathermal betas. Note, however, that the change in n ef f and Bo are largely governed by the Rankine-Hugoniot conservation relations while the change in any given T s,j is not. The change in T ec,j is due to the partition of free energy available due to ∆KE shn = 0, which will be discussed in Section 5.
Similar to the fractional n eh , there is evidence of a foreshock in β eh,j shown as a two-step decrease, i.e., the beta begins to drop ahead of the shock ramp, reaches a plateau or slight jump near the ramp, then continues to drop to the downstream values for all selection criteria β eb,// vs.   Table 4 found in Wilson III et al. (2020). The format is the same as Figure 3.
except Criteria HM and Criteria PA. The interesting thing is that the upstream running median is roughly near unity for all βs,j, by design of course, but the downstream running median drops to as low as ∼0.20 for β eh,j and β eb,j . That is, the running median decreases by upwards of ∼80% across the shock for the suprathermal electrons. Therefore, the ratio of the thermal-tomagnetic field energy density of the suprathermal electrons is upwards of ∼80% smaller in the downstream than upstream. Another interesting feature is that β eh,j begins to decrease upstream of the ramp center near the times when the fractional n eh increases which appears to be due to a local decrease in T eh,j . There is some bias in the steepness of the gradient due to the cluster of data near the ramp center of the IP shocks. This is due to the burst mode trigger of the 3DP instrument, which was designed to capture shock ramps in high time resolution. Thus, the large spread of values near the ramp center are dominating many of the running median calculations through that time period. In effect, this is actually reducing the variance in the running median seen in the more sparsely sampled far upstream and downstream regions. Therefore, it is worth noting to avoid over interpreting these gradients.  In this section, numerical instability thresholds will be presented. The purpose of examining various instabilities is driven, in part, by the lack of dependence on many macroscopic shock parameters of the change in many electron fit parameters (this will be discussed in Section 5) and that waves radiated by instabilities often drive VDFs away from Maxwellians (e.g., Verscharen et al. 2018Verscharen et al. , 2019, e.g., they generate power-laws and/or other non-Maxwellian features. There is also the motivation discussed in Paper I regarding the distribution of energy from a wave to the particles. The issue of inelasticity versus a more standard heating is discussed. Note that inelastic collisions have been tangentially discussed under different circumstances in previous theoretical work (e.g., Scudder & Olbert 1979). That is, inelastic scattering will increase the core exponents, sec or pec and qec, but may or may not cause a significant change in T ec,j while the standard idea of heating involves only the change in T ec,j without affecting the exponents.
Whether the core exponent changes or the temperature is critical for determining the change in the VDF profile/shape. For instance, whistler modes interacting with an initial Maxwellian VDF can generate strong deviations from Maxwellianity (e.g., Chang et al. 2013;Gary et al. 2011;Hughes et al. 2014;Saito & Gary 2007), typically resonating with the suprathermal halo and/or beam/strahl (e.g., Coroniti et al. 1982;Lengyel-Frey et al. 1994, 1996Oka et al. 2017Oka et al. , 2019Wilson III et al. 2012. Whistler waves are also of interest because there is some observational evidence that they can generate energetic electron tails (e.g., Oka et al. 2019;Wilson III et al. 2012. If the whistler is radiated by the whistler heat flux instability (WHFI) (Gary et al. 1994(Gary et al. , 1999, during the radiation of the wave fields the skewness of the VDF will reduce and the heat flux carrying electrons will scatter. The net effect will increase A eb . The wave fields can also pitch-angle scatter the halo electrons, resulting in larger A eh . If the whistler is radiated by the temperature anisotropy instability (WTAI), the act of radiating the wave fields will reduce A eh and/or A eb . The waves can propagate to another location and pitch-angle scatter the suprathermal particles which increases A eh and/or A eb much like the WHFI-driven waves.
Note that the analysis presented in the following includes a larger than typical time range (i.e., ±2 hours) about the shock ramp. Therefore, some of the results may be more indicative of the solar wind than the shock itself. Further, the following analysis uses only the fit parameters from all the electron components and both the protons and alpha-particles (Kasper et al. 2006) but does not directly include effects of non-Maxwellian features in any of the calculations. That is, the calculations use fit parameters, like density and temperature, directly in analytical formulas and numerical estimates for a given instability threshold without modification of the formula or numerical estimates. If the threshold involves only the total ion and electron populations, then the effective electron and total ion parameters are used.

Whistler Instabilities
In this section, the instability thresholds for the WHFI and WTAI are examined based on the observed properties of the electron VDFs. Specifically, this section will focus on the numerical results of Gary et al. (1994) and Gary et al. (1999). To this end, the electron heat flux is calculated which requires stable fit solutions for all three electron components, i.e., 10,983 of the total 15,210 VDFs examined satisfy this requirement (see Paper II for further details on the integration of the total electron model functions).  Figure 6 are defined as instabilities having a maximum positive growth rate, γmax, satisfying γmax > 10 −1 Ωcp. Over plotted are the electron VDF fit results from the present work. In the range of β ec, shown, there are 9362 valid VDF fit results. For reference, 90% of the data for Criteria AT satisfy the following: • 0.28 rad/s Ωcp 2.03 rad/s; • 0.49 s Ωcp −1 3.52 s; • 0.76 km ρ cef f 4.04 km; • 23.6 km ρ cef f 167 km; • 0.93 km λ cef f 3.65 km; and • 36.7 km ρ cef f 154 km. Thus, the growth times corresponding to the threshold lines in Figure 6 satisfy 4.92 s γmax −1 35.2 s or ∼0.03-0.24% of the total time examined around each IP shock. Figure 6 shows that >80% of all 10,983 VDFs are at or above the threshold for either the WHFI or WTAI. Limiting to the 9362 VDFs shown in Figure 6, ∼54% are unstable to the WHFI and ∼43% are unstable to the WTAI. That is, only ∼3% of the VDFs are stable for these criteria. Further, Gary et al. (1999) noted that in the presence of a finite electron heat flux and A eh > 1.01, the heat flux carrying electrons are always unstable to the WHFI 5 . Of the 10,983 VDFs with a calculated heat flux, ∼62% satisfied A eh > 1.01, i.e., ∼62% are at least linearly unstable to the WHFI. These rates are significantly larger than recent work using data from the ARTEMIS mission (e.g., Tong et al. 2019), which found occurrence rates of whistler waves to be 2%. However, the ARTEMIS work is limited the observations to the pristine solar wind and used A ef f instead of A eh for anisotropy threshold calculations. When they limit their occurrence rate estimates to intervals satisfying A ef f > 1, the rates jump to ∼15%.
Note that Gary et al. (1994) only used two bi-Maxwellian electron components while this work uses three non-Maxwellian electron components, which makes the comparison subject to scrutiny. For instance, the use of bi-Maxwellian instead of bi-kappa is known to cause differences in the WHFI (e.g., Lee et al. 2019;Shaaban et al. 2018) and others (e.g., Lazar et al. 2012Lazar et al. , 2013Lazar et al. , 2014Lazar et al. , 2017Lazar et al. , 2018Lazar et al. , 2019Shaaban et al. 2019). In addition, the definition of T eh ec used to create the lefthand panel in Figure 6 relied upon the halo and core components only, not a single suprathermal population like that used by Gary et al. (1994).
More than half of the VDFs are unstable to the WHFI, which is consistent with the observation that both A eb and κ eb are generally smaller than A eh and κ eh (see Paper II for values and statistics), respectively. That is, the radiation of the wave reduces A eb and κ eb and subsequent scattering can increase A eh . Though it should be noted that κ eh κ eb has been found to be true in previous solar wind studies. That is, statistically the halo exponent tends to be larger than the beam/strahl at 1 AU in some studies (e.g., Štverák et al. 2009), though others show the converse and a solar cycle dependence (e.g., Tao et al. 2016a). Therefore, the larger κ eh may be a remnant of the typical conditions in the solar wind and not directly related to local instabilities. Or the differences between the two exponents may be observed with this relationship precisely because the two populations are intrinsically coupled through whistlerlike instabilities. It may also be that quasi-static fields are affecting the core more than the halo by energizing them to suprathermal energies, thus effectively increasing the suprathermal electron exponents (e.g., Schwartz et al. 1988Schwartz et al. , 2011Schwartz 2014;Scudder et al. 1986). Kinetic simulations are required to resolve this discrepancy and are beyond the scope of this study.

Other Instabilities
In this section, the instability thresholds for short wavelength 6 , electrostatic ion acoustic waves (IAWs) will be calculated and discussed. The purpose is to calculate the probability of occurrence and determine whether such modes could be playing a role in the evolution of the electron VDFs across the examined IP shocks.
Regardless, the occurrence rates of T ec p j ≥ 3 and T ef f p j ≥ 3 were calculated for reference. For Criteria AT, T ec p j ≥ 3 is satisfied for ∼29.5%, ∼22.6%, and ∼24.4% of the VDFs for j = , ⊥, and tot, respectively. These occurrence rates jump to ∼34.8%, ∼27.2%, and ∼28.6% for T ef f p j ≥ 3. If the data are limited to Criteria UP, the occurrence rates for T ec p j ≥ 3 increase to ∼38.2%, ∼38.2%, and ∼37.7% and similarly the rates for T ef f p j ≥ 3 increase as well to ∼47.4%, ∼46.1%, and ∼42.8%. If the time range is limited to -120 s ≤ ∆t ≤ +3 s (where ∆t is the time from ramp center), the rates for T ef f p j ≥ 3 are ∼41.2%, ∼32.3%, and ∼26.4%. Thus, even if the analysis is limited to times near the shock ramps the occurrence rates of temperature ratios meeting or exceeding three are smaller than all upstream observations. The occurrence rate of IAWs has been shown to peak within collisionless shock ramps (e.g., Cohen et al. 2019;Davis et al. 2020;Fuselier & Gurnett 1984;Goodrich et al. 2018a;Wilson III et al. 2007, 2014a and this rate was found to be independent of T e p tot (e.g., Wilson III et al. 2007). That the rate of T ef f p j ≥ 3 does not peak near the ramp regions of the 52 IP shocks examined herein is consistent with the apparent lack of dependence on T e i tot found in previous work (e.g., Rodriguez & Gurnett 1975;Wilson III et al. 2007). Note that these rates are significantly higher than those estimated for the ambient solar wind in a recent study (e.g., Wilson III et al. 2018).
The critical drift speed between electrons and ions for a current-driven IAW can be analytically derived for two isotropic Maxwellians that generate a current. The critical drift is known to depend upon T e i tot (e.g., Gurnett & Bhattacharjee 2005). This critical drift threshold for current-driven IAWs, ignoring the details of the ion and electron populations, is satisfied for ∼4.78% of the 12,081 VDFs for Criteria AT that had solutions for both electron and proton data. However, as discussed in Priest & Sanderson (1972), the presence of gradients in the temperature and density reduce this idealized critical drift by factors of ∼2-8. This increases the number of VDFs satisfying the critical drift threshold to ∼5.2-10.4%. These rates require context to appreciate their magnitudes.
To provide context, some statistical calculations will be performed based upon the observations and relying upon the near ubiquity of IAWs in and around collisionless shock waves (e.g., Breneman et al. 2013;Chen et al. 2018;Fuselier & Gurnett 1984;Goodrich et al. 2018b;Gurnett et al. 1979b;Wilson III et al. 2007). The typical collisionless shock ramp thickness is anywhere from 1 λe up to 43 λe up ∼ 1 λi up (e.g., Hobara et al. 2010;Mazelle et al. 2010). For the 52 IP shocks examined herein, the following are satisfied for 90% of the events: • 1.2 km λe up 4.2 km; • 155 km/s |V shn | up 700 km/s; and • 1 km L sh 180 km.
These shock ramp thicknesses correspond to time scales of ∼0.002-1.2 seconds in the spacecraft frame, or ∼10 −5 -0.008% of the total time window examined for each IP shock. The duration of each electron VDF is ∼3 seconds 7 , which means the shock ramps at most ∼30% of the minimum cadence of the Wind 3DP instrument. Thus, any given VDF cannot parameterize only the shock ramp and the instruments cannot directly measure the shock ramp currents.
Despite the limitation of the particle data time resolution and unrealistic VDF profile assumptions in the theory, there were still nearly 600 VDFs that satisfied the critical drift threshold for current-driven IAWs. This corresponds to roughly 11 VDFs for each of the 52 IP shocks examined herein or a total duration more than 27 times that of the longest shock ramp in this study. Further, there were nearly 3500 VDFs (or ∼67 per IP shock) that satisfy T ef f p j ≥ 3, i.e., the often quoted temperature ratio threshold for IAWs. The IAWs of interest here should affect the electron VDFs on time scales much shorter than the integration time, thus the Wind 3DP instrument should only observe the post-instability form.
The Criteria DN VDFs have larger core self-similar exponents (sec, pec, and qec) than Criteria UP VDFs, with the strongest shocks consistently showing flattops (i.e., pec ≥ 4) for minutes to hours in the downstream, similar to terrestrial bow shock observations. The larger core exponents in the downstream regions combined with IAW amplitudes positively correlated with Mach number(e.g., Wilson III et al. 2007) is consistent with IAWs stochastically accelerating the core electrons (e.g., Dum et al. 1974;Dum 1975). This type of stochastic acceleration is qualitatively referred to as inelastic scattering throughout this three-part study. However, such a change in VDF profile has also been interpreted as due to the acceleration by quasi-static cross-shock electric fields (e.g., Hull et al. 2001;Schwartz et al. 1988Schwartz et al. , 2011Schwartz 2014;Scudder et al. 1986). Note that for both sec and pec the differences between Criteria LM and Criteria HM shocks are not statistically significant. That is, theory and observation suggest that the magnitude of both the quasi-static cross-shock electric fields and IAW amplitudes should increase with increasing Mach number. Therefore, if IAWs are affecting the core electrons are they increasing the 7 The instrument is actually triggered on the sun pulse from the sun sensor, which depends upon the spin rate of the spacecraft bus. The spin period for Wind has remained near ∼3 seconds for the entirety of the mission but varies by a ∼0.1 seconds depending on date and time. exponent or the temperature or both? If so, why and how? Similar questions arise regarding the quasi-static cross-shock electric field explanation. Kinetic simulations are required to resolve this discrepancy and are beyond the scope of this study.
Finally, two more instabilities will be discussed that also affect the ions since most previous instability work has focused on the ions (e.g., Kasper et al. 2013;Klein et al. 2017Klein et al. , 2018Maruca et al. 2012). The purpose is to examine the differences in the particle populations near IP shocks versus what is typically considered ambient solar wind. The threshold for both the mirror and firehose instabilities can also be calculated following the approach 8 in Chen et al. (2016). Using only VDF solutions when all five particle populations have finite velocity moments, the plasma is unstable to the firehose instability ∼1.3% of the time and mirror ∼13.5%. These rates are ∼10 and ∼20 times larger, respectively, than the rates found by Chen et al. (2016). It may not be surprising to find that VDFs are statistically more unstable near IP shocks than the ambient solar wind, since the shock itself constant is a source of multiple types of free energy. However, these rates do not significantly change if the time range of analysis is limited to ±20 minutes of the shock ramp center suggesting the enhanced rates are either due to the separation of electron components or the lack of inclusion of a secondary proton beam. Despite this the Criteria UP VDFs should be treated as generally being more unstable than time periods that intentionally avoid IP shocks.
As a final note, the instability analysis presented in this section should be interpreted with care. The rates calculated are based upon thresholds and do not directly imply anything about whether an observable wave amplitude would result. For instance, the threshold for the whistler instabilities are for growth rates at 10% of the proton cyclotron frequency, i.e., >18,000 times longer than a single electron cyclotron period. Therefore, that only ∼3% of the VDFs are stable does not imply large amplitude whistler waves should be observed for nearly all intervals examined herein 9 . That is, in the cases where the instability thresholds are barely sur-passed, the resulting fluctuation amplitudes may be so small that they are below the noise floor of many instruments. Further, the separation of the electron VDF into three populations, all non-Maxwellian, is different than the two drifting bi-Maxwellians assumed by Gary et al. (1994). Thus, the large fraction of VDFs satisfying instability thresholds presented herein should not be interpreted as most time intervals exhibiting large amplitude electromagnetic fluctuations.

Temperature Differences
In this section, the changes in the fit results across all shocks are examined. The data are presented as the median of all permutations of the difference, ∆Q, of parameter Q between the upstream and downstream values for each IP shock. The uncertainties for the fit parameters are half the magnitude of the difference between X 5% and X 95% . The uncertainties for the macroscopic shock parameters result from the standard propagation of uncertainties given by the values in the Wind shock database.
These medians with uncertainties were then fit to a model power-law function, Y = A X B + C, assuming Poisson weights, where X is one of the macroscopic shock parameters and Y is one of the medians of all permutations of the difference (or ratio) across the shock. Given that the shock must transform the change in bulk kinetic energy, the obvious shock parameter to examine is the change in kinetic energy, ∆KE shn (see Appendix A), across the shock 10 . Figure 7 shows the relationship between the change in temperature versus the change in shock kinetic energy and Table 1 shows the values of the fit parameters with goodness of fit estimates. Here, the ∆T s,j values represent the median 11 of all the permuted temperature differences between downstream and upstream for each IP shock. The ∆KE shn values are just computed from the |U shn | j values from the Wind shock database. The uncertainties are calculated using the standard propagation of uncertainties.
10 All other relevant shock parameters were examined and weaker relationships were found between M f up (not shown) and the permuted differences of β eb,j and P ec,j . Weak relationships were also observed between ∆Ū shn (not shown) and the permuted differences of the following: Aec, A ef f , T ec,j , T ef f,j , T p,j , and T i,j . 11 Note that the use of the average of all permutations of the differences did not yield significantly different results than the median.  The first thing to notice is that the relationship is not linear. In fact, several previous studies examined the change in temperature and found positive correlations between ∆T e,tot and variants of ∆KE shn (e.g., Hull et al. 2000;Schwartz et al. 1988; Thomsen et al. 1987Thomsen et al. , 1993. However, if one examines the plots of ∆T e,tot versus ∆KE shn in each of these studies, it's not clear whether the trend is linear or otherwise. For instance, examination of Figures 5 and 6 in Feldman et al. (1983b) do not appear to have a linear trend. The linear relationship between ∆T e,tot and ∆KE shn found by Hull et al. (2000) required that they ignore shocks with ∆KE shn < 100 eV in the fit. The relationship between ∆T e,tot and ∆Ū shn 2 found by Fitzenreiter et al. (2003), however, did appear to be linear. Note that many of the previous studies examined the higher Mach number terrestrial bow shock.
However, it is not unreasonable to assume that a linear relationship would exist if the energy conversion was a linear. That is, if the irreversible transformation of excess kinetic energy across the shock 12 into particle heating occurred directly with no intermediary processes 13 , 12 Some of the total energy transformation must be irreversible to initiate a shock from a nonlinearly steepening wave but once initiated, reversible processes can maintain the shock (e.g., Shu 1992  , for the effective electron temperatures (first column), core electron temperatures (second column), and total ion temperatures (third column). The first, second, and third rows show the parallel, perpendicular, and total temperature changes, respectively. The red dashed line in every panel is a power-law fit (function defined between panels a and d) to the data and the green lines show the associated uncertainty bounds. The associated fit parameters are shown in each panel. The magenta colored points were ignored during the fit process as outliers.
one may expect that a linear relationship would exist with the change in electron and ion temperature with some associated efficiency. A linear relationship also requires, in general, fewer assumptions to model. That is, the scatter plot of points in previous studies could have just as easily been fit to a nonlinear function like a power-law though it is likely these authors chose a linear relationship as that is the simplest possibility. One should also note that the exponents for ∆T ec, and ∆T ef f, are larger than all exponents of ∆T i,j . However, if the five magenta points at large ∆KE shn are included in the fits, the exponents for ∆T ec,j and ∆T ef f,j would decrease more than ∆T i,j . Note that previous work has either inferred (e.g., Ghavamian et al. 2007Ghavamian et al. , 2014 or showed with in situ measurements that stronger shocks heat ions more than electrons (e.g., Schwartz et al. 1988; Thomsen et al. 1993;Masters et al. 2011). In fact, Schwartz et al. (1988) and others have even noticed that ∆T e,tot/∆T t,tot ∝ M A up −1 . Therefore, the change in slope/trend at larger ∆KE shn may be indica-tive of differences in shock energy dissipation at stronger shocks, as suggested by trends in previous work.
Even when the five magenta points at large ∆KE shn (i.e., magenta points in upper right-hand corner of each panel of Figure 7) were included in the fit, the linear fit had much largerχ 2 or Σ f it values (not shown). Though it is also fair to argue that the fits shown in Figure 7 are not really good fits despite the lowχ 2 or Σ f it values shown in Table 1. That is, the data have large uncertainties and large relative spread for each ∆KE shn , as evidenced by the green lines on either side of the red fit lines. Therefore, it is likely that if even stronger shocks were added to this data set the power-law trend presented in Figure 7 would need to be modified by either an exponential roll-over or higher order terms or the trend would entirely fall apart. Thus, these fits should be treated with caution and/or skepticism.
To verify that the use of the median on the permutations of all differences was not causing the nonlinearity, the same temperature differences were plotted ver-sus ∆KE shn but now using ∆T s,j instead of ∆T s,j (not shown). That is, the average over each region was calculated as a single scalar prior to finding the difference between the regions. The nonlinear fit lines were still a better match to the results than a linear line. Thus, the nonlinear relationship between the temperature increase and the change in kinetic energy across the shock appears to be real at least for low ∆KE shn .
Most of the earlier work looking at the dependence of ∆T s,j on ∆KE shn focused on the Earth's bow shock, which is typically higher Mach number (i.e., M f up 3-10) than those examined herein. The Mach numbers in this study satisfy 1.01 ≤ M f up ≤ 6.4, with 90% satisfying ∼1.15-4.00, and a median of ∼1.86. Note that 45 of the 52 IP shocks examined herein satisfied M f up < 3. Further, the relationship between M f up and ∆KE shn is not linear. This begs the question of what could cause the energy transformation from ∆KE shn into ∆T s,j to be nonlinear, if the trend is real.
Regardless, the fraction of ∆KE shn distributed to each of the populations shown in Figure 7, for 90% of events (i.e., X 5% to X 95% range), satisfy the following: • 1.5% ∆T ec,j /∆KE shn 34%; • 0.8% ∆T ef f,j /∆KE shn 41%; and • 1.4% ∆T i,j /∆KE shn 63%. Although these ratios are mostly uniform for all j for both electron populations, the ∆T i, /∆KE shn range is systematically smaller than the other two components. That is, if one separates the perpendicular and total from the parallel components then 90% of the events would satisfy: • 5.8% ∆T i,{⊥,tot} /∆KE shn 63%; and • 1.4% ∆T i, /∆KE shn 40%. Similarly the median values for ∆T ec,j /∆KE shn and ∆T ef f,j /∆KE shn are in the ∼7.6-9.3% range while the ∆T i,j /∆KE shn median values satisfy ∼10.7-17.9%, with the parallel component being the smallest. In summary, the electron core and effective populations only gain ∼40-80% of what the ions do in thermal energy across the shock for these events.

Energy Density Differences
In this section, the partition of energy among the five primary constituent particle populations will be discussed. These five are the electron core (s = ec), halo (s = eh), beam/strahl (s = eb), and ion proton (s = p) and alpha-particle (s = α) populations. Similar to Section 5.1, the median of all permutations of the differences will be discussed.
The normalized pressures, ∆Πs,j and ∆ψs,j were plotted (not shown) versus θ Bn , |V shn | up , |U shn | up, M f up, M A up, M T e up, ∆Ū shn , ∆KE shn , and ∆ξ shn (i.e., every macroscopic shock parameter predicted to be of importance here). In the following, weak correlations imply there is a trend but the large scatter and multiple outliers make interpretation difficult. Moderate correlations are for clear trends but still with a significant spread in data (e.g., similar to ∆T s,j plots in Figure 7). There were no good correlations observed for any pair of parameters examined, only a few moderate correlations, and several weak correlations (not shown).
For reference, the following will show parameters as X 5% X X 95% ,X, for the 52 IP shocks examined herein: • 3.48%, ∼1.68%. Note that the Wind SWE Faraday cups have difficulty separating alpha-particles from protons and finding good nonlinear fit solutions in the immediate downstream of the stronger shocks in this study, which is likely affecting the upper bounds of both ∆Πp,j and ∆Πα,j. The reason for emphasizing this point is that previous work (e.g., Schwartz et al. 1988;Thomsen et al. 1993;Masters et al. 2011) found more energy going to the ions as the Mach number increases but the core electrons seem comparable to the protons here. A possible difference may be that previous work examined the total ion and electron VDFs, which may blur the differences in trends between the various components of each species. It is also worth noting that the ∆Πs,j values result from the absolute value of the differences, i.e., the actual change may be negative for one population.
To determine which population gained more thermal energy density across the shocks, the ratio of the elec-tron component to proton thermal energy densities are examined where the the parameters are shown as X 5% X X 95% ,X. The values the 52 IP shocks examined herein are as follows: • 391%, ∼199%. Therefore, the change in core electron thermal pressure is comparable to or larger than that for the protons even for the Criteria HM shocks, unexpectedly. For the Criteria LM shocks the core electrons receive roughly the same heating as the protons, consistent with previous results (e.g., Schwartz et al. 1988). Note that despite the use of pressure here as a measure of thermal energy density, none of these processes are occurring in a thermodynamic system. That is, there is no well defined equation of state for each of the particle populations.
In summary, none of the energy density ratios showed a clear dependence upon any macroscopic shock parameter. None of the ∆Πs,j or ∆ψs,j cared about θ Bn , i.e., the change in the ratio of thermal energy density to both the total pressure and total energy density is independent of shock geometry. In other words, the shock geometry does not appear to affect the change in the partition of energy amongst the five major particle populations. The most correlations were found with ∆Ū shn , but again none of them were good. The only good correlation was observed between ∆ξ shn and ∆ j (not shown), i.e., the change in total energy density increases with increasing change in shock kinetic energy density. This merely shows that as the available free energy increases, the total internal energy increases, as expected.
Finally, the absolute changes in normalized partial pressures were dominated by the core electrons and protons with the suprathermal electrons and alpha-particles serving as minor constituents which is again expected. However, the absolute differences discussed in this section do not inform us whether the change is positive or negative. Further, somewhat unexpectedly the core electron pressure changes were often larger than that for the protons for Criteria HM shocks, while for Criteria LM shocks the changes were closer to previous observations (e.g., Schwartz et al. 1988). Again, some of this is most likely due to the issues facing the Wind SWE Faraday cups downstream of the strongest shocks in this study. However, the larger than unity ratios of ∆P ec,j / ∆P p,j for even Criteria LM shocks were not really expected. To help address this, SEA plots are shown to illustrate the trends in ∆Πs,j versus time in the following section.

Thermal Energy Density Trends
In this section, the partition of energy among the five primary constituent particle populations will be discussed. These five are the electron core (s = ec), halo (s = eh), beam/strahl (s = eb), and ion proton (s = p) and alpha-particle (s = α) populations. Similar to Section 3, SEA plots will be presented.
Since none of the ∆Q seemed to show very good correlations with any macroscopic shock parameter predicted to be of importance, the statistical trend of the normalized energy densities directly using SEA were examined. The purpose is to see if any clear trend versus time(space) emerges that is not reflected in macroscopic differences or ratios. Note that the SEA plots involve the calculation of one-variable statistics for all points within a given time bin, while the ∆Q values are calculated on a shock-by-shock basis. Figure 8 shows the running median only from SEA of the partial thermal pressures, Πs,j, of each of the major particle populations in the solar wind, including the core electrons (blue lines), halo electrons (red lines), beam/strahl electrons (orange lines), protons (magenta lines), and alpha-particles (purple/violet lines). First, Πec,j and Πp,j dominate at all times for all selection criteria, as expected. Second, the general trend of all electron populations is for their fractional thermal energy density to reduce across the shock ramp except for Πec,j for Criteria HM shocks. Both the Π eh,j and Π eb,j electrons consistently decrease across the shock, with the weakest change across Criteria PA shocks. The Π eb,j show a continual decreasing trend across Criteria HM shocks with a distinct jump at the shock ramp in contrast to Π eh,j which basically levels off and recovers downstream. There are also several intervals where Πec,j and Πp,j seem to oscillate exactly out of phase from each other, likely owing to pressure balance features near these shocks. The interesting aspect is that they are common/strong enough to show up in a running median constructed from SEA on 52 different shocks. -3600 0 +3600 -3600 0 +3600 Figure 8. The running median (e.g., red lines in Figure  1) from the SEA of the partial thermal pressures, Π s,j (see Appendix A for parameter definitions), for each of the following particle populations: core electrons (blue lines), halo electrons (red lines), beam/strahl electrons (orange lines), protons (magenta lines), and alpha-particles (purple/violet lines). Unlike previous SEA plots, these are not normalized to an upstream median value. The vertical green line indicates the ramp center time. All panels have uniform horizontal and vertical axis ranges.
• 15% Πp,j 74%; • 1% Π eh,j 15%; • 0.2% Π eb,j 8%; and • 0.2% Πα,j 13%. To illustrate that the energy partitions shown as running medians in Figure 8 can be characteristic of most of the data, Figure 9 shows the same parameters for the same populations but as shaded regions bounded by X 5% and X 95% . What is immediately clear is that again Πec,j and Πp,j dominate at all times except for some transient excursions to large values for Π eh,j . Both Π eb,j and Πα,j have relatively large spreads in the range between percentiles, but they still follow the same trends illustrated by the running median. That is, the suprathermal electron fractional thermal energy density decreases while both ion species increase. The abrupt end in Πα,j data for Criteria PA shocks in both Figures 8 and 9 are due to low statistics owing to difficulties in finding high quality fits for the alpha-particle peak (e.g., see SWE fit flag requirements in Wilson III et al. 2018).
It is also clear that Criteria HM shocks have the largest separation between Πec,j and Πp,j for all time periods compared to other selection criteria. The downstream running median values for Criteria LM and Criteria PE shocks oscillate out of phase but this is not directly reflected in X 5% -X 95% . Criteria PA shocks exhibit the most overlap in X 5% -X 95% and Criteria HM shocks have the least. In both the runningX and X 5% -X 95% , Π eh,j shows a very clear decrease across the shock except for Criteria PA shocks where the change in X 5% and X 95% is more difficult to observe. The abrupt drop in Π eh,j and Π eb,j across shocks is likely related to the drop in n eh /n ef f and n eb /n ef f across the shocks because the associated temperatures show slight increases across the shock (see additional SEA plots of T s,j found in Wilson III et al. (2020)).

Superposed epoch analysis (SEA) of the fit results
show that the general trend of the normalized (to the upstream median) exponents for all electron components increases across the shock ramp for all selection criteria 14 . That is, the suprathermal electron tails are steeper and the core electrons at low velocities are flatter, relative to a Maxwellian, and steeper at higher energies. The normalized density, ns/n ef f , SEA plots show a general increase (relative to the upstream median) across the shock ramp for the core, but decrease for both halo and beam/strahl, for all selection criteria. That is, only the fraction of core electrons increases across the shock. Finally, the ratio of the partial-to-total pressure, Πs,j, for the protons and alpha-particles increase across the shock, while the ratio for all electron populations decrease.
An illustrative example of the possible electron VDF evolution for the weakest and strongest shocks is qualitatively shown in Figure 10. The electron VDF starts as a narrow peaked distribution with hard tails and evolves into a an almost box-like distribution with weaker, soft tails. Thus, the energy density becomes consolidated into the core population. It is clear that stronger shocks have a different downstream profile than weaker shocks. The strong shock example is indicative of the strongest shocks in this study and similar to observations downstream of the terrestrial bow shock. The larger exponents, in all three electron components, for the strong shock example produces a flatter peak in the core and steeper slopes in the halo and beam/strahl. Although the number density of the downstream weak shock ex-14 Except for Criteria UP and Criteria DN, of course ample is lower than for the strong shock example, the phase space density peak is almost an order of magnitude larger due to the smaller exponents and thermal speeds. Potential reasons for the change in profile are discussed later.
The majority of the thermal energy density is held in the core electrons with Πec,j ∼ 25-92% (i.e., fraction of the total pressure) and the protons taking up most of the rest with Πp,j ∼ 13-72% for Criteria AT. For ∼95% of the suprathermal electron and alpha-particle velocity moments, none individually satisfy Πs,j > 18%. That is, the partition of thermal energy density is completely dominated by the core electrons and protons while the suprathermal electrons and alpha-particles serve as mostly minor contributors. Further, the magnitude of the change in partial pressures, ∆Πs,j, across the 52 IP shocks examined herein (shown as Xmin X Xmax) satisfies following: • 1.7% ∆Πec,j 55.8%; • 0.6% ∆Π eh,j 35.9%; • 0.3% ∆Π eb,j 16.5%; • 1.8% ∆Π ef f,j 66.7%; • 1.8% ∆Πp,j 66.7%; and • 0.07% ∆Πα,j 10.8%. Therefore, the core electrons and protons both carry the largest Πs,j and they tend to experience the largest ∆Πs,j.
Of the three minor partial pressure populations (i.e., halo, beam/strahl, and alpha-particles), the halo electrons consistently dominate in the upstream region, especially Π eh, . In the downstream, Π eb,⊥ and Π α,⊥ can be comparable to Π eh,⊥ . Interestingly, the Πec,j and Πp,j often vary out of phase with each other suggesting a partial thermal pressure balance. This is more weakly reflected in the Πα,j variations relative to the Π eh,j and Π eb,j values.
Unexpectedly, it was found that the change in pressure of the electrons could be comparable to or larger than the protons. That is, the ratio of these changes for the 52 IP shocks examined herein (shown as Xmin X Xmax,X) satisfies following: • 16.9% ∆P ec,j / ∆P p,j 391%, ∼101%; • 1.11% ∆P eh,j / ∆P p,j 41.3%, ∼5.41%; • 0.47% ∆P eb,j / ∆P p,j 18.3%, ∼2.55%; and • 19.6% ∆P ef f,j / ∆P p,j 453%, ∼107%. Therefore, the change in core electron thermal pressure is comparable to or larger than that for the protons. This is somewhat unexpected because most IP shock observations and theory suggest that the ions should gain more thermal energy than the electrons. However, previous low Mach number bow shock observations showed that the electrons gained roughly the same fraction of thermal energy as the ions. The fractional energy density gain by the electrons increases with increasing Mach number in this study but it's unclear if the increase is influenced by difficulty of the Wind SWE Faraday cups measuring ion distributiions downstream of the strongest events.
The nonlinear trend shown in Figure 7 between ∆T s,j and ∆KE shn is apparent "by eye" but as the lower/upper uncertainty bounds on the fit lines illustrate, the trend is moderate at best. While some previous studies examined a linear relationship between these two parameters, they were primarily focused on observations at the much higher Mach number bow shock. Further, the fits were performed and shown between these two parameters because they were the best. Virtually every combination of macroscopic shock parameters was examined against nearly every ∆Q but none exhibited a good correlation except between ∆ξ shn and ∆ j (not shown). However, the good correlation between ∆ξ shn and ∆ j merely shows that as the available free energy increases, the total internal energy increases, which is not surprising.
The influence of the macroscopic shock parameters is not evident in the current dataset. In the following, a moderate correlation shows a clear trend but a significant spread in the dependent variable for any given independent variable value (e.g., similar to ∆T s,j plots in Figure 7). There is a moderate, positive correlation between ∆Ū shn and both ∆Π eb,j and ∆Πα,j. There are also moderate, positive correlations between ∆Πα,j and both ∆KE shn and M T e up. In summary, none of the energy density ratios showed a clear dependence upon any macroscopic shock parameter. Nothing was found to show any dependence upon θ Bn . That is, the change in the fractional thermal energy densities of the five major particle populations appear to be independent of the shock geometry. Therefore, microscopic instabilities were investigated because they are known to affect each population differently (e.g., Artemyev et al. 2013Artemyev et al. , 2014Artemyev et al. , 2015Artemyev et al. , 2016Artemyev et al. , 2017aArtemyev et al. ,b, 2018Chang et al. 2013;Dum 1978a,b;Hughes et al. 2014;Krall & Trivelpiece 1973;Osmane & Hamza 2012;Petkaki et al. 2003Petkaki et al. , 2006Sagdeev 1966).
Over 97% of the 9362 VDFs shown in Figure 6 are unstable to either the whistler heat flux (WHFI) or temperature anisotropy instabilities (WTAI). Roughly ∼1% and ∼14% of the 12,081 VDFs with both ion and electron data are unstable to the firehose and mirror instabilities, respectively, which is over an order of magnitude higher the ambient solar wind estimates (e.g., Chen et al. 2016). Nearly 30% of these VDFs satisfy T e p tot ≥ 3, i.e., the threshold below which ion acoustic waves (IAWs) are  predicted to experience heavy Landau damping. However, only ∼5% of the VDFs satisfied the critical drift velocity (between electrons and ions) threshold necessary to generate current-driven IAWs ignoring temperature gradients which reduce this threshold. Thus, there is sufficient statistical evidence of unstable VDFs near IP shocks, which is not surprising. There are several caveats when interpreting the above instability occurrence rates. For instance, this work used three electron, one core proton, and one alpha-particle beam populations whereas Chen et al. (2016) used one electron, two proton, and one alpha-particle beam populations. Further, the instability threshold for the WHFI and WTAI corresponds to extremely slow growth rates in excess of 18,000 electron cyclotron periods. Finally, the IAW thresholds were derived for single, isotropic Maxwellian electron and ion populations. Thus, these instability occurrence rates are only intended to serve as zeroth order proxies of a fully kinetic treatment that includes non-Maxwellian VDFs and multi-component electron and ion populations. Despite the caveats and the sometimes unexpectedly high occurrence rates, recent fully kinetic analysis of the ion VDFs in the solar wind finds that most intervals are linearly unstable (e.g., Klein et al. 2018Klein et al. , 2019. The examination of whistler and acoustic instabilities is driven by their difference in resonance energy ranges. The former tend to scatter with the suprathermal electrons affecting both As and κs (e.g., Chang et al. 2013;Gary et al. 2011;Hughes et al. 2014;Saito & Gary 2007), while the latter strongly scatters the core reducing Aec and increasing sec (or pec). There is some evidence to support the differences in A eh and A eb combined with differences in κ eh and κ eb to suggest the WHFI is present and scattering the suprathermal electrons near these IP shocks. The core electron exponents (sec, pec, and qec) are consistently larger downstream than upstream and the profile of the downstream core electrons does reach the flattop stage (i.e., pec ≥ 4) in the strongest shocks. This is consistent with the nonlinear saturation stage of IAWs interacting with electrons (e.g., Dum et al. 1974;Dum 1975;Dyrud & Oppenheim 2006;Sagdeev 1966;Vedenov 1963) but it could also be due to quasi-static cross-shock electric fields (e.g., Hull et al. 2001;Schwartz et al. 1988Schwartz et al. , 2011Schwartz 2014;Scudder et al. 1986). However, the influence of both of these effects should increase with Mach number yet the differences in core electron exponents between low and high Mach number shocks are not statistically significant. This begs the question of how much energy/momentum goes into increasing the exponent versus increasing the temperature (which does consistently and significantly change across the shocks).

CONCLUSIONS
The analysis of 15,210 VDFs observed by the Wind spacecraft within ±2 hours of 52 interplanetary (IP) shocks has been presented. Five primary constituent particle populations were included in the analysis, which are the electron core (s = ec), halo (s = eh), beam/strahl (s = eb), and proton (s = p) and alpha-particle (s = α) populations. The analysis revealed that most of the VDFs are at least linearly unstable to one or more instabilities, consistent with a similar conclusion based on different techniques applied to ion data in the ambient solar wind (e.g., Klein et al. 2018Klein et al. , 2019. For the weaker shocks, the change in core electron, effective electron, and total ion temperature is positively correlated with the change in kinetic energy across the shock in a way that appears to be nonlinear.
The evolution of the electron VDF illustrated in Figure 10 shows the qualitative trends observed in this study. The remaining question is what controls the increase in the core exponent versus temperature and why some strong shocks seem to prefer one or the other while the strongest shocks increase both. Neither instabilities or quasi-static fields alone appear to be capable of producing the observed changes in the electron VDFs. Kinetic simulations and theoretical work are required to resolve this discrepancy and are beyond the scope of this study.
Neither instabilities or quasi-static fields alone are sufficient to explain the evolution of the electron VDFs across the IP shocks. For instance, some of the stronger shocks showed significant core electron heating but little change in the core exponent (i.e., no additional flattening) while others showed significant increases in the core exponent (i.e., strong flattening) but comparatively weak heating. Theory and simulations have found that nonlinear stochastic acceleration by ion acoustic waves (i.e., referred to as inelastic collisions herein) can selfconsistently generate flattops in the core electron VDFs but it is not clear what fraction of the wave energy goes into increasing the exponent versus increasing the temperature. Further, previous studies have argued that quasi-static fields are capable of producing flattop VDFs but it's not clear what exponent value should result. Both of these explanations are plagued by the fact that the downstream core electron exponents are not positively correlated with Mach number or change in kinetic energy, quantities that have been correlated with the amplitude of both the electrostatic waves and quasistatic electric fields in collisionless shocks. That is, stronger shocks do not consistently generate flatter core electron VDFs. Addressing such fundamental questions in kinetic theory are planned in future studies.
In summary, the core electrons and protons carry the most thermal energy density and they also experience the largest changes in thermal energy density across the shocks. A moderate, positive correlation is found between ∆T s,j and ∆KE shn , for s = ec, ef f , and i. Weaker correlations are found between some VDF parameters and any other macroscopic shock parameter but nothing is correlated with the shock normal angle. Surprisingly, the change across the shock in core electron pressure relative to the proton pressure, ∆P ec,j / ∆P p,j , was found to be comparable to or larger than unity. That is, the core electron pressure change was larger than that of the protons in at least 23 of the 52 shocks, for any component. If only the parallel pressures are examined, the number of shocks increases to 28 of 52. Future work will examine the details of the kinetic physics involved in collisionless shocks to address these issues. As in Papers I and II, this appendix the symbols and notation used throughout will be defined. All directiondependent parameters we use the subscript j to represent the direction where j = tot for the entire distribution, j = for the the parallel direction, and j = ⊥ for the perpendicular direction, where parallel/perpendicular is with respect to the quasi-static magnetic field vector, Bo [nT]. The use of the generic subscript s to denote the particle species (e.g., electrons, protons, etc.) or the component of a single particle species (e.g., electron core). For the electron components, the subscript will be s = ec for the core, s = eh for the halo, s = eb for the beam/strahl, s = ef f for the effective, s = int for the integrated (see Paper II for definition), and s = e for the total/entire 15 population. Below are the symbol/parameters definitions: ) βs,j ≡ the plasma beta [N/A] of the j th component of species s (see Equations A1i and A1j) κs ≡ the kappa exponent of species s (e.g., see Wilson III et al. 2019a, for definition in model fit equation) ss ≡ the symmetric self-similar exponent of species s (e.g., see Wilson III et al. 2019a, for definition in model fit equation) ps(qs) ≡ the parallel(perpendicular) asymmetric self-similar exponent of species s (e.g., see Wilson III et al. 2019a, for definition in model fit equation) φsc ≡ the scalar, quasi-static spacecraft potential [eV] (e.g., Pulupa et al. 2014;Scime et al. 1994) (see Appendices of Paper I for more details) -Emin ≡ the minimum energy bin midpoint value [eV] of an electrostatic analyzer (e.g., see Appendices in Wilson III et al. 2017) q e, = me 2 d 3 v f e (mod) v v 2 ≡ the parallel electron heat flux [µW m −2 ] of the entire electron VDF model, f e (mod) = f (core) + f (halo) + f (beam) qeo = 3 2 me ne V T ec, 3 ≡ the free-streaming limit electron heat flux [µW m −2 ] (e.g., Gary et al. 1999) Similar to Paper I, the variables that rely upon multiple parameters are given in the following equations: T ef f,j = s ns T s,j s ns (A1a) where n ef f is defined as: For the macroscopic shock parameters, the values are averaged over asymptotic regions away from the shock transition region.
shock parameters subscripts up and dn ≡ denote the upstream (i.e., before the shock arrives time-wise at the spacecraft for a forward shock) and downstream (i.e., the shocked region) -Q j ≡ the average of parameter Q over the j th shock region, where j = up or dn -∆Q = Q dn -Q up ≡ the change in the asymptotic average of parameter Q over the j th shock region -Rns = ns dn / ns up ≡ the average shock compression ratio of species s -R Qs,j = Qs,j dn / Qs,j up ≡ the downstream-to-upstream j th component ratio of the asymptotic average of parameter Q of species s -∆Q α,β = Qα -Q β ≡ the set of all permutations of the difference between all the downstream (α) and all the upstream (β) values -∆Q ≡ the median of ∆Q α,β , i.e., the median value of all permutations of all differences across the shock of parameter Q -n sh ≡ the shock normal unit vector [N/A] θ Bn ≡ the shock normal angle 16 [deg] -|V shn | j ≡ the j th region average shock normal speed [km s −1 ] in the spacecraft frame -|U shn | j ≡ the j th region average shock normal speed [km s −1 ] in the shock rest frame (i.e., the speed of the flow relative to the shock) -|KE shn | j = 1 2 mp |U shn | j 2 ≡ the j th region average shock normal kinetic energy [eV ] in the shock rest frame -|ξ shn | j = 1 2 mp np j |U shn | j 2 ≡ the j th region average shock normal kinetic energy density [eV cm −3 ] in the shock rest frame -M A j = |U shn | j / V A j ≡ the j th region average Alfvénic Mach number [N/A] -M f j = |U shn | j / V f j ≡ the j th region average fast mode Mach number [N/A] -M T e j = |U shn | j / V T ef f,tot j ≡ the j th region average electron thermal Mach number [N/A] -M cr ≡ the first critical Mach number [N/A] For brevity the percent difference between one-variable statistics values for two different selection criteria are defined here. The percent difference between parameters satisfying Criteria DN and Criteria UP is defined as ∆Q d2u = (Q dn − Qup) /Qup × 100%, where Q is any one-variable statistic value. Similarly, the percent difference between Criteria HM and Criteria LM is defined as ∆Q h2l = (Q HM − Q LM ) /Q LM × 100% and that between Criteria PE and Criteria PA is ∆Q ⊥2 = (Q P E − Q P A ) /Q P A × 100%.
As in Paper II, integrated velocity moments refer to the velocity moments calculated by integrating over the entire model function, f e (mod) = f (core) + f (halo) + f (beam) , rather than the fit values from the components. The integrated moments are only calculated for VDFs with stable solutions for all three components using the Simpson's 1 3 Rule algorithm. The integrals are calculated in the core electron rest frame, thus the only relevant heat flux component is the parallel, q e, , because the suprathermal electrons have no finite perpendicular drift velocities (e.g., see Paper I). For further details on the integrated velocity moments, see Paper II. These definitions are used throughout.

B. EXTRA STATISTICS
In this section some extra statistics are presented in tabular form in Tables 2 and 3.  Table 2 continued   Table 3 continued Note-All values represented as percentages. The header symbols are the same as in Table 2. For symbol definitions, see Appendix A.