Electron Energy Partition across Interplanetary Shocks. I. Methodology and Data Product

Analyses of 15,314 electron velocity distribution functions (VDFs) within ±2 hr of 52 interplanetary (IP) shocks observed by the Wind spacecraft near 1 au are introduced. The electron VDFs are fit to the sum of three model functions for the cold dense core, hot tenuous halo, and field-aligned beam/strahl component. The best results were found by modeling the core as either a bi-kappa or a symmetric (or asymmetric) bi-self-similar VDF, while both the halo and beam/strahl components were best fit to bi-kappa VDF. This is the first statistical study to show that the core electron distribution is better fit to a self-similar VDF than a bi-Maxwellian under all conditions. The self-similar distribution deviation from a Maxwellian is a measure of inelasticity in particle scattering from waves and/or turbulence. The ranges of values defined by the lower and upper quartiles for the kappa exponents are κec ~ 5.40–10.2 for the core, κeh ~ 3.58–5.34 for the halo, and κeb ~ 3.40–5.16 for the beam/strahl. The lower-to-upper quartile range of symmetric bi-self-similar core exponents is sec ~ 2.00–2.04, and those of asymmetric bi-self-similar core exponents are pec ~ 2.20–4.00 for the parallel exponent and qec ~ 2.00–2.46 for the perpendicular exponent. The nuanced details of the fit procedure and description of resulting data product are also presented. The statistics and detailed analysis of the results are presented in Paper II and Paper III of this three-part study.


Background and Motivation
The solar wind is an ionized gas experiencing collective effects where Coulomb collisions occur, but the rates are often so low that, for instance, two constituent particle species, s′ and s, are not in thermodynamic or thermal equilibrium, i.e., (T s′ /T s ) tot ≠ 1 for s′ ≠ s, and the relevant scale lengths are orders of magnitude smaller than the collisional mean free path (e.g., Wilson et al. 2018). Therefore, for any process dependent on scales like the thermal gyroradii, ρ cs , or inertial lengths, λ s , the media is considered collisionless (see Appendix A for definitions). That the solar wind is a nonequilibrium, weakly collisional, kinetic gas results in multicomponent velocity distribution functions (VDFs) for both ions (e.g., Kasper et al. 2006Kasper et al. , 2012Kasper et al. , 2013Maruca et al. 2011;Wicks et al. 2016) and electrons (e.g., Schwartz & Marsch 1983;Maksimovic et al. 1997Maksimovic et al. , 1998Lin 1998;Pierrard et al. 1999Pierrard et al. , 2001Štverák et al. 2008, 2009Pulupa et al. 2014a).
An illustrative example, showing the three electron components typically observed in the solar wind near 1 au below ~1.2 keV, is shown in Figure 1. The component parameters are exaggerated 10 for illustrative purposes but are based on the fit results of the VDF shown in Figure 4. The core is modeled by a symmetric bi-self-similar VDF and the halo and beam/ from a bi-Maxwellian on visual inspection. The halo and beam/strahl were modeled with a bi-kappa model function for all VDFs examined since they always have a power-law tail and previous work found kappa model functions to be the best approximation (e.g., Maksimovic et al. 2005;Štverák et al. 2009).
For each IP shock, an iterative process was followed to correct for the spacecraft potential, ϕ sc (details found in Appendix B), and define fit parameter initial guess values and constraints to yield stable solutions for the most VDFs (detailed steps found in Appendix C, and list of initial guess values and constraints found in Supplemental Material ASCII files (Wilson et al. 2019c) described in Appendix F). The process of defining the initial guess values and constraints is discussed in Section 3.2, and the quantified estimates of the fit quality are discussed in Section 3.3.
A total of 15,314 electron VDFs were observed by the Wind spacecraft within ±2 hr of 52 IP shocks. Of those 15,314 VDFs, 15,210 progressed to fit analysis, and stable model function parameters were found for 14,847 (~98%) core fits, 13,871 (~91%) halo fits, and 9567 (~63%) beam/strahl fits. The reason for the large disparity in beam/strahl fits compared to the other two components will be discussed in Section 3.3 and Appendix C.

Velocity Distribution Functions
This section introduces and defines the model functions used to fit to the particle VDFs in this study with examples provided to illustrate shape and dependences on parameters.
The most common VDF used to model particle VDFs in space plasmas is the bi-Maxwellian (e.g., Feldman et al. 1979aFeldman et al. , 1979bFeldman et al. , 1983bKasper et al. 2006), given by where A M is given by where v o,j is the drift speed of the peak relative to zero along the jth component, V T, j 2 is the thermal speed given by Equation (6c), V j is the velocity ordinate of the jth component, and n o is the number density.
been around for decades (e.g., Vasyliunas 1968;Feldman et al. 1983a;Maksimovic et al. 1997;Salem et al. 2003). It is beyond the scope of this study to explain the physical interpretation/origin of this function, but there are several detailed discussions already published on the topic (e.g., Livadiotis 2015; Livadiotis et al. 2018). A generalized powerlaw particle distribution is given by a bi-kappa VDF (e.g., Mace & Sydora 2010;Livadiotis 2015), for electrons here as where A κ is given by and B κ is given by where Γ(z) is the Riemann gamma function of argument z and V Tj is again the most probable speed of a 1D Gaussian for consistency, i.e., it does not depend on κ.
The last model VDF is called a self-similar distribution, which results when a VDF evolves under the action of inelastic scattering (e.g., Dum et al. 1974;Dum 1975;Horton et al. 1976;Horton & Choi 1979;Jain & Sharma 1979;Goldman 1984) or flows through disordered porous media (e.g., Matyka et al. 2016). The symmetric form is given by where A SS is given by A slightly more general approach can be taken where the exponents are not uniform, which will be referred to as the asymmetric self-similar distribution function. The asymmetric functional form is given by where A AS is given by Again, this will reduce to a bi-Maxwellian in the limit where p → 2 and q → 2. Note that in the event that the exponents s, p, or q are not even integers, the velocity ordinates, (V ‖ − v o‖ ) and (V ⊥ − v o⊥ ), will become absolute values to avoid complex roots and negative values of f (V ‖ , V ⊥ ). Example one-dimensional cuts of these three model VDFs can be found in Figure   2 for comparison.
The self-similar exponents are mostly a new variable, since most previous work modeled the core electrons as a bi-Maxwellian (e.g., Štverák et al. 2008, 2009Bale et al. 2013;Pulupa et al. 2014b). There are a few studies that used one-dimensional self-similar functions to model a select few electron VDFs near collisionless shocks (e.g., Feldman et al. 1983aFeldman et al. , 1983b, finding values consistent with those presented in Table 2. However, these studies did not define the normalization parameter in terms of the number density and thermal speeds (see, e.g., Equation 3(a) and 4(a)), but rather found a numerical value from empirical fits, i.e., the normalization parameter was not coupled to the physical parameters of the fit function. At least one study in the solar wind did define the normalization constant, but they only considered a one-dimensional, isotropic distribution (e.g., Marsch & Livi 1985). Although several theoretical works predicted ranges of possible self-similar exponent values under various extrema scenarios (e.g., Dum et al. 1974;Dum 1975;Horton et al. 1976;Horton & Choi 1979;Jain & Sharma 1979;Goldman 1984), this is the first time the model has been used on a statistically significant set of VDFs.
The following is an illustrative example that shows how the signal-to-noise ratio of particle detectors strongly depends on the number density and thermal speed and that hot, tenuous plasmas are much more difficult to measure and accurately model. Examine the onedimensional cuts shown in Figures 2 and 4. The toy models in Figure 2 are shown to illustrate the effect of thermal speed and exponents on the model fit function peaks and shapes. Notice that increasing the thermal speed of the Maxwellian from V Te = 1500 to 5500 km s −1 drops the peak phase-space density by nearly two orders of magnitude. The cut line also passes the ±20,000 km s −1 velocity boundary (i.e., roughly the upper energy bound of the EESA Low instrument) at a phase-space density roughly one order of magnitude higher than the colder examples. That is, the change in thermal speed reduced the dynamic range of observed phase-space densities by three orders of magnitude. Suppose that one examines a more extreme example with n e = 15 cm −3 and V Te = 10,000 km s −1 . In this case, the difference between the peak and the lowest phase-space density within the ±20,000 km s −1 velocity boundary would only be a factor of ~55, i.e., slightly more than one order of magnitude. For more details about derivation and normalization constants, see the Supplemental Material (Wilson et al. 2019c).

Fit Parameter Constraints
This section involves the discussion of the constraints/limits placed on fit parameters for each electron component and justifies them based on physically significant assumptions.
As an illustrative example, Figure 3 shows the densities of the protons, alpha-particles, and three electron components (blue squares) and the associated uncertainties (red error bars) for a subcritical, quasi-perpendicular IP shock (see, e.g., Wilson et al. 2019c Note that there are two time periods after 11:00 UTC where a few fit results satisfy n eb /n eh ⩾ 1. Figure 3 is illustrative of some of the error analysis employed in the present study and the fact that the beam/strahl fit more often fails than the core or halo as evidenced by the number of points. Below the details of how the fit parameters are constrained/limited are outlined with physical arguments.
First, the present study differs from some previous studies in that the fits are performed on the two-dimensional VDF rather than separate fits on one-dimensional cuts of the twodimensional VDF (e.g., Maksimovic et al. 2005;Pulupa et al. 2014aPulupa et al. , 2014b. One of the limitations of the latter approach is that the distribution function is not necessarily a separable function, which can introduce difficulty for the physical interpretation of the results. However, the latter approach has numerous advantages, including the stability of the solutions and ease with which the solutions are found with nonlinear least-squares software, i.e., it is generally easier to fit to a one-dimensional cut than a two-dimensional distribution. The present study uses the former approach to avoid the difficulties introduced for nonseparable functions. For instance, when fitting to the parallel one-dimensional cut, the amplitude of the VDF is directly tied to the amplitude of the perpendicular cut. The amplitude of all standard model two-dimensional, gyrotropic VDFs is dependent on n s , V Ts, −1 , and V Ts, ⊥ −2 . While it is computationally possible to fix the amplitude to the observed amplitude of the data for each cut and only vary the respective thermal speeds/temperatures and exponents, the inversion to find n s can be problematic if care is not taken. For instance, the normalization constants differ for one-dimensional cuts from the two-dimensional gyrotropic VDF (see, e.g., Equation 1(a)). Although this approach involves fewer free parameters and should thus be easier to fit, it is much more restrictive in parameter space, i.e., n s only varies indirectly through the variation of the thermal speeds/temperatures and exponents.
Given that fitting to a two-dimensional gyrotropic VDF has more free parameters and orders of magnitude more degrees of freedom, a stable solution requires reasonable constraints/ limits on the variable parameters. There are some obvious boundaries determined by instrumental and physical constraints. As shown in the previous section, the difference between the highest and lowest phase-space densities is important for the signal-to-noise ratio, but it is also relevant to fitting model functions to the data. For instance, if an electron distribution had a population with V Te ⩾ 10,000 km s −1 , the weights would not provide sufficient contrast between the peak and tails to constrain a stable and reliable fit without multiple imposed constraints. In contrast, electron VDFs with thermal speeds below ~1000 km s −1 fall below the lowest energy of the detector and so would be artificially hotter if they were observed (e.g., Paschmann & Daly 1998). A similar effect is often observed by spacecraft with electrostatic analyzers designed for the magnetosphere, not the comparatively cold, fast solar wind beam (e.g., McFadden et al. 2008aMcFadden et al. , 2008bPollock et al. 2016).
Statistical studies of the solar wind have shown that the maximum range of the total electron temperature is T e,j ~ 2.29-77.2 eV or V Te,j ~ 450-2600 km s −1 (e.g., Wilson et al. 2018). Previous studies have found that the electron halo temperatures satisfy T e,j ~ 14-560 eV or V Teh,j ~ 1100-7000 km s −1 (e.g., Feldman et al. 1975Feldman et al. , 1978Feldman et al. , 1979aMaksimovic et al. 1997Maksimovic et al. , 2005Skoug et al. 2000;Tao et al. 2016aTao et al. , 2016bLazar et al. 2017). Previous studies have also found that the electron beam/strahl temperatures satisfy T eb,j ~ 20-150 eV or V Te,j ~ 1300-3600 km s −1 (e.g., Ogilvie et al. 2000;Viñas et al. 2010;Tao et al. 2016aTao et al. , 2016b. Thus, a range of allowed core thermal speeds from ~1000 to ~10,000 km s −1 can be assumed. There are similar instrumental constraints on the drift speed of the three components. The core, however, is not likely to exhibit drift speeds (in the ion rest frame) in excess of several hundred kilometers per second (e.g., Pulupa et al. 2014a). In the present work, most fit results show less than 50 km s −1 , i.e., only 1838 of 14,847, or ~12%, have drift speeds exceeding 50 km s −1 , consistent with previous work. 14 In contrast, owing to the physical interpretation of the strahl/beam component, most (8848 of 9567, or ~92%) have drift speeds in excess of 1000 km s −1 . The allowed core, halo, and beam/strahl drift speeds loosely ranged from ~1000 to ~10,000 km s −1 for most events. In some events, a lower bound was imposed to prevent unphysical fit results, e.g., beam/strahl component with near zero drift speed (see Supplemental Material ASCII files (Wilson et al. 2019c) described in Appendix F for ranges for specific events). Note that V oes,⊥ was fixed during the fitting, i.e., it was not allowed to vary. Originally this parameter was free to vary but resulted in fewer stable fits and rarely varied by more than a few kilometers per second. In some events, an explicit V oec,⊥ was set as the initial guess values determined from examination of the distributions, but this is for a small minority of events (333 of 14,847, or ~2%).
It has also been empirically found that the EESA Low detector has issues when n ce ≾ 0.5 cm −3 or n ce ≿ 50 cm −3 for typical solar wind thermal speeds. 15 This is rarely an issue, as only 41 of the 14,847 VDFs analyzed (or ~0.3%) have fit results falling outside the range ~0.5-50 cm −3 . Note that the total electron density, n e = n ec + n eh + n eb ~ n e = n p + 2n α , is constrained by the total ion density from SWE and the total electron density from the upper hybrid line observed by the WAVES radio receiver (Bougeret et al. 1995), when possible (see Appendix B for more details).
Physically, the halo and beam/strahl components are suprathermal; thus, they should not have the dominant contribution to the total phase-space density of the VDF. Therefore, it is physically consistent to assume that the fit results should satisfy n eh /n ec < 1 and n eb /n ec < 1. The solutions were constrained to satisfy n eh /n ec < 0.5 and n eb /n ec < 1 based on results found in previous studies near 1 au (e.g., Feldman et al. 1975;Maksimovic et al. 1997Maksimovic et al. , 2005Skoug et al. 2000;Štverák et al. 2009;Viñas et al. 2010;Pierrard et al. 2016;Tao et al. 2016b).
In numerous previous studies that assumed a three-component solar wind electron VDF near 1 au (e.g., Maksimovic et al. 2005;Štverák et al. 2009;Pulupa et al. 2014aPulupa et al. , 2014b, constraints were sometimes assumed such as that the fits satisfy n eb /n eh < 1. There is no restriction on this ratio 16 imposed during the fit process, and 1824 of 9313 (or ~20%) of the fits satisfy n eb /n eh ⩾ 1. In fact, it was found that imposing the constraint n eb /n eh < 1 during the fit process actually greatly reduced the number of stable solutions found for the beam/ strahl component. 17 Previous work did show that the ratio n eb /n eh decreases with increasing radial distance from the Sun, dropping below unity before 1 au, on average, but the ranges overlapped, allowing for n eb /n eh ⩾ 1 (e.g., Štverák et al. 2009).
Another constraint that is often assumed/used is that the strahl/beam component be only antisunward along B o (e.g., Maksimovic et al. 2005;Štverák et al. 2009;Pulupa et al. 2014aPulupa et al. , 2014b, though some magnetic field topologies have sunward-directed beam/strahl components (e.g., Owens et al. 2017). This constraint is imposed in this study, but it is important to note that some IP shocks examined have observable electron foreshocks. A consequence is that the halo component of the fit results effectively absorbs both the halo 15 Technically, this is an issue for nearly all electrostatic analyzers designed and flown to date. This is largely unavoidable without increasing the dynamic range of the detector significantly. 16 The number of good ratios differs from the number of beam/strahl fits because some VDFs had a stable halo or beam/strahl but not the converse. 17 Note that there was a post-fit constraint imposed limiting n eb /n eh < 3 because it was found empirically that most fits exceeding this threshold were bad/unphysical. However, not all were bad as evidenced by the example in Figure 6.
Wilson et al.

Page 12
Astrophys J Suppl Ser. Author manuscript; available in PMC 2020 July 03. and the shock-reflected electron component in the events where this is directed sunward along B o (this is very rare). If the shock-reflected electron component is directed antisunward, they will be included in the beam/strahl fit (this is much more common). The net result for the former is a smaller (T ⊥ /T ‖ ) eh and on the latter is a larger (T ⊥ /T ‖ ) eb and n eb .
The lower bound of possible κ es values is defined for mathematical/physical reasons as being ≿3/2 (e.g., Livadiotis 2015; Livadiotis et al. 2018). The upper bound is set to 100 solely because above that value the difference between a bi-Maxwellian and bi-kappa VDF is smaller than the accuracy of the measurements. Although the upper bound is allowed to extend to 100, the typical upper bound observed near 1 au is <20 (e.g., Maksimovic et al. 1997;Štverák et al. 2009;Pierrard et al. 2016;Tao et al. 2016aTao et al. , 2016bLazar et al. 2017). The range of possible values for s ec , p ec , or q ec falls between 2 and 10 for physical reasons (e.g., Dum et al. 1974;Dum 1975;Horton et al. 1976;Horton & Choi 1979;Jain & Sharma 1979;Goldman 1984).
Finally, by definition the halo and beam/strahl components represent the lowest-energy suprathermal components of the electrons. Therefore, it is natural to assume that T eh /T ec > 1.
There is no explicit restriction on this ratio imposed, and only 384 of 13,867 (or ~3%) of the fits satisfy T eh /T ec < 1, and these occur downstream of strong shocks where core heating dominates. However, there are numerous events where limits/constraints were imposed on the component temperatures individually. So the low percentage is not entirely unexpected. In contrast, there were no corresponding attempts to limit T eh /T eb in any way other than to fit to the data.

Quality Analysis
The initial approach was to use the reduced chi-squared value χ s 2 of component s (see Appendix D for definition) as a test of the quality of the fit. However, it was quickly determined that some fit lines matched well with the data but had χ s 2 > 10 while others did not fit well at all despite having χ s 2 ≲ 1. The issue is partly related to the calibration of the detector and thus the quality of the values (see Appendix B for more details). The issue is also related to fitting a gyrotropic model function to data that is not, in general, gyrotropic. A possible improvement would fold the entire VDF into a forced gyrotropy prior to fitting to improve counting statistics and the comparison between data and model functions, but that is beyond the scope of the current study. Therefore, a new quantity was defined to provide an additional definition of the quality of any given fit by direct comparison.
Let us use f (0) as the actual data and f (m) (=f (core) + f (halo) + f (beam) ) as the total model fit results. Then one can define the ratio of these two parameters as ℛ = f (0) / f (m) , which is a two-dimensional array of values. Then one calculates the median of this array, ℛ, to determine the percent deviation given by where δℛ is computed for each electron VDF. The values of δℛ were then used as uncertainties/error bars for all fit parameters for the associated VDF for all components. In general, the percent magnitude of the uncertainty in each of the six fit parameters should not be uniform as is used herein (see Appendix E for discussion of 1σ uncertainties). The uncertainty of any variable calculated using these fit parameters was propagated assuming uncorrelated errors.
Note that the δℛ value alone does not always characterize the quality of any given fit. Therefore, a combination of parameters is chosen to define a set of fit quality flags from best with a value of 10 to worst with a value of 0 (see Appendix F for definitions). In general, fits with flags of at least 2 or higher can be used, but low fit flags should be treated with caution.
Only ≾1% of all core, halo, and beam/strahl fits had flags of 1, while >95% of core, >89% of halo, and >61% of beam/strahl flags were at least 2. Figure 4 shows an example VDF that had a low χ s 2 for each component and a δℛ 3.0%, i.e., this is an example of an ideal fit. The distribution was fit using a symmetric bi-self-similar distribution for the core and a bi-kappa distribution for both the halo and beam/strahl components. The fit results are as follows: In contrast, Figure 5 shows an example VDF that had a high χ s 2 for two components yet still a small δℛ 9.4%, i.e., this is still an example of a good fit despite the bad χ s 2 values for the core and beam/strahl fits. The fit results are as follows: Further, the example VDF in Figure 5 differs from that in Figure 4 in that an asymmetric self-similar model is used for the former. The total fit lines also illustrate a weakness of the method used. Since the components are fit separately, the respective weights change with each fit to prevent the fitting software from giving too much emphasis to, for instance, the core of the distribution when fitting to the halo. 18 Thus, the resultant f (m) can exceed f (0) in some places. The software does a post-fit check for instances where either the combined or any component model fit exceeds the data by user-specified factors. 19 For most events, the threshold is set between ~2 and 4, but this varies, as some events have known issues. For instance, the known density from the upper hybrid line is 10 cm −3 , but no variation of ϕ sc yields fit results with n e ~ 10 cm −3 without the model exceeding the data at low energies.
The reason is related to known calibration issues (see Appendix B).
Finally, Figure 6 shows an example VDF that had a high χ s 2 for the core component and moderate for beam/strahl but a small δℛ 2.1%. This example VDF was chosen to illustrate a good fit even when n eb /n eh > 1. As previously discussed, there are post-fit constraints applied to the data based on statistical and physical constraints. The constraint relevant to Figure 6 is that requiring n eb /n eh < 3. This is why the fit flag value for the beam/strahl is zero and why χ tot 2 is larger than a few. The fit results are as follows: One can see from the figure that the halo component is rather weak compared to the beam/ strahl, which could be the result of an enhancement from the electron foreshock of this IP shock or the fast nature of the solar wind upstream of this IP shock. Regardless, the purpose of this example is to illustrate that stable and good fit solutions can be found that satisfy n eb /n eh > 1 even at 1 au.
After examining thousands of fit results, it was determined that δℛ with χ s 2 and χ tot 2 are consistently more reliable quantities used in combination for defining the quality of the fit than using χ s 2 alone. The value is also used as a proxy for the uncertainty of any given fit parameter, e.g., δn es = ± δℛ ⋅ n es /2 shown as the red error bars in Figure 3. Note that values of 100% correspond to fill values or bad fit results. In the following section the one-variable statistics of the χ s 2 and δℛ values are listed for reference to typical/expected values when evaluating the quality of a fit. In general, the best fits have small values for δℛ and all χ s 2 .
Further tests of consistency were also performed to validate the fit results. First, the EESA Low detector is known to saturate when the count rate exceeds ~10 7 counts s -1 (Lin et al. 1995). Examination of all VDFs found that a total of 10 energy-angle bins (from a total of 20,184,120), or ~5 × 10 −5 %, exceeded the maximum count rate. Therefore, it is not thought that saturation has a significant impact on the methodology and results of this study. Second, as illustrated in Figure 3, the total electron density satisfies n e ~ n p + 2n α for nearly all intervals. Statistically, the difference between the fit result for n e = n ec + n eh + n eb and n p + 2n α is within expectations. The median, lower quartile, and upper quartile values are 10.3%, 4.9%, and 19.0%, respectively, which is consistent with our δℛ statistics.
Finally, the total electron current, j e,tot = Σ s n es v os,‖ , in the ion rest frame should be zero to maintain a net zero current in the solar wind. The mean, median, lower quartile, and upper quartile for all data examined are ~22 km s −1 cm −3 , ~0 km s −1 cm −3 , ~−214 km s −1 cm −3 , and ~351 km s −1 cm −3 , consistent with previously published work on this data set (e.g., Bale et al. 2013;Pulupa et al. 2014a) and consistent with work in progress (C. S. Salem et al. 2019, in preparation). Normalizing j e,tot by n e times V Tec,tot yields a mean, median, lower quartile, and upper quartile for all data examined of ~0.17%, ~10 −8 %, ~−0.95%, and ~1.3%, respectively. Thus, the values are all small compared to unity. Quantitatively, ~97.5% of the j e,tot /(n e V Tec,tot ) values satisfy ≾5.5%. Figure 7 shows both j e,tot and j e,tot /(n e V Tec,tot ) versus seconds from every shock ramp center time in this study. One can see that although there are locations with significant deviation from zero (e.g., the shock ramp, which is not tremendously surprising, as that is where currents are supposed to exist), the mean (red horizontal line) and median (orange horizontal line) are small for both the raw and normalized current densities. Note that the data in Figure  7 include fit results where there may not be a solution for one or more components (see discussion of first data product ASCII file in Appendix F).
As a final note, there is the question about the validity of using a new model function to describe the thermal core. Of the 11,874 core VDFs fit with a symmetric bi-self-similar model function, there were 9559, or ~80.5%, that satisfied 2.0 ⩽ s ec ⩽ 2.05. That is, the majority of the distributions would be nearly indistinguishable from a bi-Maxwellian on visual inspection. Therefore, the use of the symmetric bi-self-similar model function is not entirely inconsistent with previous work that modeled the solar wind core with a bi-Maxwellian (e.g., Feldman et al. 1979aFeldman et al. , 1979b. In fact, these results show that most core VDFs are not far from thermal velocity distributions, consistent with results showing evidence for collisional effects on the core (e.g., Salem et al. 2003;Bale et al. 2013).

Summary of Fit Results
For the 52 IP shocks examined there were a total of 15,314 VDFs observed by Wind. Of those 15,314 VDFs, 15,210 progressed to fit analysis, and for the core only 534 (~4%) were modeled as bi-kappa VDFs, 12,095 (~80%) were modeled as symmetric bi-self-similar VDFs, and 2581 (~17%) were modeled as asymmetric bi-self-similar VDFs. All core bikappa VDFs were found in the upstream, and all downstream core VDFs used either a symmetric or asymmetric bi-self-similar model. All halo and beam/strahl components were fit to a bi-kappa model. The justifications for the use of these functions are given in Section 3 and Appendix C. Of those 15,210 that progressed to fit analysis, stable solutions were found for 14,847 (~98%) f (core) , 13,871 (~91%) f (halo) , and 9567 (~63%) f (beam) .
Recall that the fit results presented herein were performed on two-dimensional, (assumed) gyrotropic velocity distributions in the proton bulk flow rest frame. Most prior work numerically fit to one-dimensional cuts of the VDF or to one-dimensional reduced VDFs. There are benefits for either method, but here it is shown that the method employed is valid by illustrating the consistency with previous work. which are consistent with previous results near 1 au (e.g., Feldman et al. 1975Feldman et al. , 1979aFeldman et al. , 1983bMaksimovic et al. 1997;Phillips et al. 1989aPhillips et al. , 1989bNieves-Chinchilla & Viñas 2008;Skoug et al. 2000;Salem et al. 2001;Štverák et al. 2009;Pierrard et al. 2016 The purpose of listing these statistics is to provide a range of typical or expected χ s 2 and δℛ values for reference when determining the quality of any given fit. Note that the statistics for δℛ shown above were performed on arrays that excluded the lower and upper boundaries, i.e., 0.1% and 100% values. The statistical results of the model function exponent and drift speed results are presented below, and the full data product resulting from this work is described in Appendix F. Table 2 shows the one-variable statistics for the exponents from the model fits of the electron VDFs for the core (s = c), halo (s = h), and beam/strahl (s = b). The VDFs, modeled as bikappa (κ es ), symmetric bi-self-similar (s es ), and asymmetric bi-self-similar velocity distributions (p es for parallel and q es for perpendicular), are summarized for all time periods, upstream only, downstream only, low Mach number only, high Mach number only, quasiperpendicular only, and quasi-parallel only. The rows showing N/A (not available) for every entry had no fit results, i.e., the core was only modeled as a bi-kappa in the upstream and an asymmetric bi-self-similar only in the downstream, and therefore the converse had no results to examine.

Exponents and Drifts
For the VDFs fit to a bi-kappa, the core values typically lie between ~5 and 10, while the halo and beam/strahl lie in the ranges of ~3.5-5.4 and ~3.4-5.2, respectively. Only the core was fit to the bi-self-similar functions, and nearly all symmetric exponents are between ~2.00 and 2.04, while most of the asymmetric parallel and perpendicular exponents lie in the ranges of ~2.2-4.0 and ~2.0-2.5, respectively.
The κ eh and κ eb values are consistent with previous solar wind observations near 1 au (e.g., Maksimovic et al. 1997Maksimovic et al. , 2005Štverák et al. 2009;Pierrard et al. 2016;Tao et al. 2016aTao et al. , 2016bLazar et al. 2017;Horaites et al. 2018 There are several interesting things to note from Table 2. The mean, median, and lower/ upper quartile values for κ ec are slightly higher for high than for low Mach number shocks, though only the median and lower quartile values are significant. Since a bi-kappa model was only used for upstream core VDFs, this may imply that shock strength is somehow dependent on the upstream core electron distribution profiles. One possible physical interpretation would be that the sound speed depends on the polytropic index for each species, i.e., the equation of state assumed for the system. A bi-kappa core VDF could affect the estimate of the sound speed, thus altering the fast mode Mach number. However, the shape of the upstream VDFs will also affect the shock dissipation mechanisms. For instance, it is known that the existence of power-law tails improves the efficiency of shock acceleration (e.g., Trotta & Burgess 2019). Therefore, the larger κ ec associated with higher Mach number shocks may imply that lower energy particles have entered the tails, thus increasing the exponent. 20 In contrast, the asymmetric bi-self-similar exponents, only used in downstream regions, are effectively the same between low and high Mach number shocks. However, this changes when comparing quasi-parallel and quasi-perpendicular shocks. The p ec exponent has higher mean, median, and lower/upper quartile values for quasi-parallel than quasi-perpendicular shocks. The opposite is true for the q ec exponent. This is interesting, as higher p ec values are predicted to occur in the nonlinear saturation stages of ion-acoustic waves (e.g., Dum et al. 1974;Dum 1975 relative electron-ion drifts (i.e., currents) and are observed near both quasi-parallel and quasi-perpendicular shocks (e.g., Fuselier & Gurnett 1984;Wilson et al. 2007Wilson et al. , 2010Wilson et al. , 2012Wilson et al. , 2014aWilson et al. , 2014bBreneman et al. 2013), but their amplitudes increase with increasing shock strength (e.g., Wilson et al. 2007). If the largest ion-acoustic waves generate the largest values of p ec , then one would expect maximum values downstream of strong quasiperpendicular shocks, which is not the case here. This leads to the question of what fraction of energy goes to increasing p ec versus what fraction goes to increasing T ec,‖ . This would depend on the effective inelasticity of the wave-particle interactions, where larger inelasticity increases p ec and smaller inelasticity increases T ec,‖ (e.g., Dum et al. 1974;Dum 1975;Horton et al. 1976;Horton & Choi 1979;Jain & Sharma 1979;Goldman 1984). The interaction between a wave and a particle can be treated as inelastic if the particle affects the wave amplitude and kinetic energy during the interaction. Most test-particle treatments do not handle this self-consistently, and if the effect is distributed to an entire VDF, the net result can be a stochastic heating that increases p ec from 2.0 (e.g., Dum et al. 1974;Dum 1975).
Another theory predicts that flat-top electron distributions (i.e., p ec → ⩾4 and q ec → ~2-3) can result from the combined effects of a quasi-static, cross-shock electric potential and from fluctuation electric fields (e.g., Feldman et al. 1983a;Hull et al. 1998) through a process called maximal filling (e.g., Morse 1965). However, similar to the predictions for wave-driven flat tops, this theory should generate stronger flat tops (i.e., larger values of p ec ) for stronger quasi-perpendicular shocks, which we do not observe. Thus, the evolution of the electron VDFs does not seem consistent with the standard quasi-static, cross-shock electric potential, but rather in agreement with recent high-resolution observations at the bow shock (e.g., Chen et al. 2018;Goodrich et al. 2018).
Another interesting result is the difference in the κ eh values under different conditions. When the values of κ eh are larger (smaller), that implies a less (more) energized halo, i.e., softer (harder) spectra. One can see that κ eh is larger downstream than upstream and near high rather than low Mach number shocks. That is, the halo is less energized downstream of IP shocks and near strong IP shocks than the converse, which is somewhat unexpected, as strong shocks should more readily energize suprathermal particles (e.g., Malkov & Drury 2001;Treumann 2009;Caprioli & Spitkovsky 2014;Park et al. 2015;Trotta & Burgess 2019). In contrast, κ eh is slightly smaller (~10%) near quasi-parallel than quasiperpendicular shocks, which implies more energized halo electrons. Although quasi-parallel shocks are predicted (e.g., Malkov & Drury 2001;Caprioli & Spitkovsky 2014) and observed (e.g., ) to be more efficient particle accelerators, the predictions are usually specific to ions, while mildly suprathermal electrons are thought to most efficiently interact with quasi-perpendicular shocks (e.g., Wu 1984;Park et al. 2013;Trotta & Burgess 2019). Further, very recent simulation results suggest that the upstream electron suprathermal tail will become flatter (i.e., smaller kappa values) with increasing Mach number for quasi-perpendicular shocks (Trotta & Burgess 2019). This may explain why both κ eh and κ eb are smaller in the upstream than downstream. The time evolution of these kappa values will be examined in more detail in Paper III. Wilson et al. Page 20 Astrophys J Suppl Ser. Author manuscript; available in PMC 2020 July 03.
A major caveat of the above discussion is the exchange of particles between the various electron VDF components, i.e., former core electrons can be energized and move to the halo or the converse. Therefore, one needs to be careful when interpreting the change in a given component-specific parameter. This will be discussed in more detail in Paper III.
Finally, the κ eb values show a similar behavior between upstream and downstream and shock geometry as κ eh , but they differ between low and high Mach number shocks. That is, stronger shocks appear to energize the beam/strahl component more than weaker shocks. This is likely due to the electron foreshock component observed upstream of strong IP shocks (e.g., Bale et al. 1999;Pulupa & Bale 2008;Pulupa et al. 2010), combined with the usual solar wind beam/strahl component. Figure 8 shows histograms of κ es , s ec , p ec , q ec , and the drift speed magnitudes, V oes,j (s for electron components and j for parallel or perpendicular), for the three electron populations.
These histograms show distributions corresponding to the first part of One can see that, as discussed previously, the core parallel drift speeds (violet line, panel (d)) tend to fall below ~100 km s −1 , consistent with previous results (e.g., Pulupa et al. 2014a). In fact, most of the core and halo drifts are near zero, with the number of results satisfying V oec,‖ ⩽ 1 km s −1 and V oeh,‖ ⩽ 1 km s −1 being 8735 (~59%) and 7311 (~53%), respectively. Note that although there is sometimes a sizable perpendicular core drift (blue line, panel (d)) for some shock crossings, these were explicitly set after visual inspection of the VDFs during the iterative fitting process. The nonzero perpendicular drifts almost certainly result from inaccuracies in the calculation of the solar wind rest frame and a dipole correction to ϕ sc not included in the present analysis (e.g., Pulupa et al. 2014a; see Appendix B for more details).
The magnitudes of V oeh,⊥ and V oeb,⊥ never deviated from zero. 21 The magnitudes of V oeh,‖ range from ~0 to 8860 km s −1 , with a lower to upper quartile range of ~0-850 km s −1 and a mean (median) of ~580 km s −1 (~0.1 km s −1 ). The magnitudes of V oeb,‖ range from ~1000 to 9330 km s −1 , with a lower to upper quartile range of ~1750-3090 km s −1 and a mean (median) of ~2580 km s −1 (~2480 km s −1 ). As previously discussed, the lower bound for V oeb,‖ was imposed on the basis of physical arguments, while the magnitude of V oeh,‖ was allowed to go to zero. If only magnitudes satisfying V oes,‖ ⩽ 1 km s −1 are considered, the mean (median) and lower to upper quartile ranges are ~42 km s −1 (~30 km s −1 ) and ~14-52 km s −1 for V oec,‖ and ~1227 km s −1 (~903 km s −1 ) and ~362-1695 km s −1 for V oeh,‖ .
21 This was an explicit constraint imposed on all fits but would also have resulted largely from the initial guess that both V oeh,⊥ and V oeb,⊥ equal zero. That is, the fit software uses initial guesses to estimate gradient magnitudes for changes between iterations. So if the initial guess is null, the step size will be null as well.  Broiles et al. 2016). The values for s ec , p ec , and q ec are consistent with previous results as well (e.g., Feldman et al. 1983aFeldman et al. , 1983b. The interesting aspect of VDFs being well modeled by bi-self-similar functions is that such functions are used to describe the evolution of distributions for either the flow through disordered porous media (e.g., Matyka et al. 2016) or the influence of inelastic scattering (e.g., Dum et al. 1974;Dum 1975;Horton et al. 1976;Horton & Choi 1979;Jain & Sharma 1979;Goldman 1984). It is unlikely that the former applies directly, but the latter may be interpreted in the following manner. The typical approach for test-particle simulations used to examine wave-particle interactions does not include feedback from the particles on the waves. In a real plasma, the particles can alter three properties of electromagnetic waves: their amplitude (potential energy), momentum, and kinetic energy. Consider a simple scenario whereby a particle reflects off of an electromagnetic wave field along one dimension. If done self-consistently, the particle can reduce the wave amplitude in addition to affecting the field momentum and kinetic energy. In the case of a reduced wave amplitude, the resulting scattering problem can be treated as a simple inelastic collision. 22 Thus, the net result of an ensemble of particles interacting with a wave field can be stochastic (e.g., Dum et al. 1974;Dum 1975), which provides one physical justification for the use of the bi-self-similar functions. These functions are also convenient in that they reduce to bi-Maxwellians in the limit where the exponents go to 2, i.e., the deviation from a Maxwellian is a measure of inelasticity in the particles' interactions with waves and/or turbulence. 23 Further, as previously discussed, ~80.5% of the core VDFs modeled with a symmetric bi-self-similar function had exponents satisfying 2.0 ⩽ s ec ⩽ 2.05. Therefore, the majority of the core electron VDFs would be visually indistinguishable from a bi-Maxwellian, which supports previous work that used thermal distributions to model the core (e.g., Feldman et al. 1979aFeldman et al. , 1979b and work that found evidence for collisional effects in the core distribution (e.g., Salem et al. 2003;Bale et al. 2013).
The κ ec seem to correlate with 〈M f 〉 up , which may suggest a shock strength dependence on the shape of the upstream electron VDFs. In contrast with expectations from a dependence on quasi-static fields, the values of p es are higher for quasi-parallel shocks, while q es are higher for quasi-perpendicular shocks, yet neither depends on 〈M f 〉 up .
Somewhat surprisingly, the values of κ eh are larger downstream than upstream, and they increase with increasing 〈M f 〉 up . That is, the halo spectra are softer downstream and near strong shocks. Quasi-parallel shocks, however, correlate with smaller κ eh , i.e., harder halo spectra. Generally, quasi-parallel shocks are predicted to be more efficient particle accelerators for suprathermal ions and very energetic electrons 24 (e.g., Caprioli & 22 That is, the particle kinetic energy may not be preserved through the interaction even if the wave kinetic energy is conserved. 23 It is also worth noting that a finite time correlation included in wave-particle interactions, something missing from quasi-linear theory, can yield a similar VDF profile (work in progress by coauthors). 24 Suprathermal is defined here for ions in the several to tens of keV energy range, while the electrons are many tens to hundreds of keV for typical 1 au solar wind collisionless shocks.

Wilson et al. Page 23
Astrophys J Suppl Ser. Author manuscript; available in PMC 2020 July 03.
Spitkovsky 2014), but electrons in the halo energy range are predicted to be energized the most efficiently at shocks satisfying θ Bn > 80° (e.g., Park et al. 2013).
Unlike the halo, κ eb are smaller near high Mach number shocks than near low Mach number shocks. The difference is likely a twofold consequence of the combined effects from shockaccelerated foreshock electrons and the method used to fit the distributions. That is, the beam/strahl component is always fit to the antisunward, field-aligned side of the VDF, while the halo is fit to the opposite. For nearly all IP shocks at 1 au, the shock normal is antisunward in a direction that would be aligned with the nominal, ambient beam/strahl electron component. For both the halo and beam/strahl, the ratios of 〈k eh 〉 dn /〈k eh 〉 up and 〈k eb 〉 dn /〈k eb 〉 up increase with increasing 〈M f 〉 up . That is, the downstream halo and beam/ strahl spectra are softer than the upstream for stronger shocks. Again, this is likely a consequence of the foreshock electrons that are not observed upstream of weak shocks. The details of the electron component velocity moments and associated changes will be discussed further in Papers II and III.
In summary, the first part of this three-part study presented the first statistical study to find that the core electron distribution is better fit to a self-similar VDF than a bi-Maxwellian under all conditions. This is an important result for kinetic theory and solar wind evolution. This work also provides the methodology and details necessary to reproduce and qualify the results of the nonlinear least-squares fitting performed herein. In Papers II and III, the statistical and analysis results of the velocity moments will be presented in detail. These observations are relevant for comparisons with astrophysical plasmas like the intra-galaxycluster medium, and they provide a statistical baseline of electron parameters near collisionless shocks for the recent Parker Solar Probe and upcoming Solar Orbiter missions.

Appendix A: Definitions and Notation
In this appendix we define the symbols and notation used throughout. In the following, for all direction-dependent parameters we use the subscript j to represent the direction, where j = tot for the entire distribution, j = ‖ for the parallel direction, and j = ⊥ for the perpendicular direction. Note that parallel and perpendicular are with respect to the quasistatic magnetic field vector, B o (nT). The generic subscript s is used to denote the particle species (e.g., electrons, protons) 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 = eff for the effective population, and s = e for the total/ entire population. Below are the symbol/parameter definitions: s. ϕ sc ≡ the scalar, quasi-static spacecraft potential (eV) (e.g., Scime et al. 1994b;Pulupa et al. 2014a) t. E min ≡ the minimum energy bin midpoint value (eV) of an electrostatic analyzer (see, e.g., Appendices in Wilson et al. 2017Wilson et al. , 2018.
The variables that rely on multiple parameters are given in the following equations: where n e is defined as For the macroscopic shock parameters, the values are averaged over asymptotic regions away from the shock transition region.

a.
shock parameters a. subscripts up and dn ≡ denote the upstream (i.e., before the shock arrives timewise at the spacecraft for a forward shock) and downstream (i.e., the shocked region) The critical Mach numbers are phenomenologically defined as follows: for 〈M f 〉 up /M cr ⩾ 1 an ion sound wave could not phase stand within the shock ramp (e.g., Edmiston & Kennel 1984;Kennel et al. 1985), for 〈M f 〉 up /M ww ⩾ 1 a linear magnetosonic whistler cannot phase stand upstream of the shock ramp (e.g., Krasnoselskikh et al. 2002), for 〈M f 〉 up /M gr ⩾ 1 a linear magnetosonic whistler cannot group stand upstream of the shock ramp, and for 〈M f 〉 up /M nw ⩾ 1 a nonlinear magnetosonic whistler is no longer stable/stationary and will result in the shock ramp "breaking" and reforming.
These definitions are used throughout.

Appendix B: Spacecraft Potential and Detector Calibration
The electron electrostatic analyzer data suffer from several sources of uncertainty, including differences between the theoretical maximum detector efficiency and actual (e.g., Bordoni 1971;Goruganthu & Wilson 1984), unknowns regarding the detector dead time 25 (e.g., Schecker et al. 1992;Meeks & Siegel 2008), and an unknown spacecraft potential (e.g., Scime et al. 1994aScime et al. , 1994bPulupa et al. 2014a;Lavraud & Larson 2016). Although the corrections for microchannel plate (MCP) degradation, etc., have not been updated since very early in the mission, the last calibrations were performed well after the initial and most dramatic scrubbing phase that occurs when the instrument is in space (see, e.g., McFadden et al. 2008aMcFadden et al. , 2008b, for further discussions of MCP degradation over time). The currently used calibrations are those from optical geometric factor corrections, onground calibrations, and in-flight calibrations (D. Larson 2019, personal communication, 2011. Although there are expected to be corrections to these calibration values over the course of the time span examined in this work, the same data in the same time range have been presented in numerous refereed publications (including, but not limited to, Salem et al. 2001Salem et al. , 2003Wilson et al. 2009Wilson et al. , 2010Wilson et al. , 2012Wilson et al. , 2013aWilson et al. , 2013bWilson et al. , 2018Bale et al. 2013;Pulupa et al. 2014aPulupa et al. , 2014b. Updating the calibration tables is beyond the scope of this work but is actively being pursued (C. S. Salem et al. 2019, in preparation).
Although the Wind spacecraft has the capacity to measure electric fields (Bougeret et al. 1995), it does not measure the DC-coupled spacecraft potential, ϕ sc . It does, however, 25 The dead time is the time period when the detector is unable to measure incident particles owing to the channel's discharge recovery time (i.e., time to replenish electrons to wall of conductive material in the microchannel plate), preamp cycle rates, etc. 26 The cycle rate or sample rate of this preamp is listed as 2 MHz, but it is not constant.
consistently observe the upper hybrid line (also called the plasma line), which provides an unambiguous measure of the total electron density, n e . For instance, the Wind/SWE Faraday cups (FCs; Ogilvie et al. 1995) are calibrated to these measurements assuming n e = n p + 2n α . Ions are generally not significantly affected by ϕ sc , as they typically have ~1 keV of bulk kinetic energy in the solar wind.
To estimate ϕ sc , an initial guess is determined numerically from the ion density. The value of ϕ sc is then adjusted until n e = n ec + n eh + n eb from the fits roughly equals 27 n p + 2n α and/or when photoelectrons disappear from the VDF plots. 28 Once a reliable estimate of ϕ sc is determined for each VDF for each IP shock, the software is cycled through all VDFs for that event and the data are saved. This process is repeated for each IP shock event. An example time series of ϕ sc is shown in Figure 3.
Note that the values of ϕ sc determined above should not be treated as the absolute or correct spacecraft potential values. The reason is that the detector efficiency and gain calibrations suffer from the issues discussed above (e.g., Bordoni 1971;Goruganthu & Wilson 1984). Therefore, the ϕ sc values are proxies for the spacecraft potential that comprise a complicated nonlinear convolution of the real spacecraft potential and the detector dead time and efficiency. Despite this uncertainty, the ϕ sc values estimated herein are consistent with those in previously published work on the same data set within the same time span (e.g., Bale et al. 2013;Pulupa et al. 2014a). Further, the consistency checks discussed in s ec tion 3.3 provide further validation of the fit results. Table 3 provides the one-variable statistics of the ϕ sc values for all VDFs, as well as upstream and downstream only, low and high Mach number only, and quasi-parallel and quasi-perpendicular only periods. There are no dramatic differences other than that the values of ϕ sc are slightly smaller downstream than upstream, slightly higher for high than low Mach number shocks, and largest (by mean, median, and quartiles) for quasi-parallel shocks. Figure 9 shows ϕ sc versus n i as both the raw values and a renormalized version where the EESA Low detector E min is used as an offset. The data were fit to a power-law-exponential, Y = X B e CX + D, where Y = (ϕ sc + E min ) 5 (eV) and X = n i (cm −3 ). The fit parameters producing the cyan dashed line are A = 2.272 ± 0.013 (cm +3B ), B = −0.431 ± 0.019 (N/A), C = 0.00115 ± 0.00155 (cm +3 ), and D = 2.0 ± 0.0 (eV), with a reduced chi-squared value of The choice of the form of the fit line is empirical and matches the observations in trend. The typical approach is to measure the spacecraft potential and number density and then fit to a function of the spacecraft potential for the number density, i.e., n i = n i (ϕ sc ) (e.g., Scudder et al. 2000). As previously stated, Wind cannot actively measure ϕ sc , and the values shown in Figure 9 are really a proxy owing to the uncertain values for the dead time and efficiency for 27 Note that the value of n e for a constraint is taken from SWE and the upper hybrid line observed by the WAVES radio receiver (Bougeret et al. 1995), when possible. 28 When ϕ sc is too low, a discontinuous "spike" appears in the cuts of the VDF. The spike-like feature can also be seen in 1D energy spectra shown in the spacecraft frame with no adjustment for ϕ sc .

Wilson et al. Page 29
Astrophys J Suppl Ser. Author manuscript; available in PMC 2020 July 03. each detector anode. The purpose of the above approach is to find a semianalytical expression for ϕ sc that only depends on n i (or n e ) as an initial estimate. The unexpected result here is that the trend depends on E min as an offset, which is likely only reflecting a one-sided measurement boundary, preventing the detector from observing the entire VDF.
Note that similar analysis on the same data set has also found a small dipolar correction to the typical monopolar approximation used herein (e.g., Pulupa et al. 2014a). The dipole term is typically less than 1 eV, however, and only ~1.5% of all the VDFs examined in our study satisfied ϕ sc < 1.5 eV. Further, the dipole correction will only affect the odd velocity moments, i.e., the drift velocity and heat flux. We did not calculate the heat flux, but we did observe perpendicular core velocity drifts previously shown to be affected by the dipole correction (e.g., Pulupa et al. 2014a).

Appendix C: Numerical Analysis Procedure
The data are fit to a user-defined model function using the nonlinear least-squares fit algorithm called the LMA (Moré 1978). The generalized LMA software, called MPFIT (Markwardt 2009), requires at minimum the following inputs when fitting to a twodimensional array of data: The error array will be ignored if the user supplies an array of weights, . The details of the use of the software and documentation are provided by the author at https:// www.physics.wisc.edu/~raigm/idl/fitting.html and in the publication Markwardt (2009).
For the purposes of finding numerical fits to electron VDFs in the solar wind, a substantial set of wrapping routines were written for use with the MPFIT libraries and can be found at https://github.com/lynnbwilsoniii/wind_3dp_pros. The wrapping software also provides detailed documentation with extensive manual pages and numerous comments throughout.
The approach used for each electron VDF is as follows:

1.
The raw VDF data, f (0r) , are retrieved as an IDL structure with the data in units of counts. A copy is created, and the data structure tag is replaced with the square root of the number of counts, f (0cr) , i.e., Poisson statistics are assumed.

2.
A unit conversion is applied to change to units of phase-space density (i.e., cm −3 km −3 s +3 ), and then the energies are adjusted to account for the spacecraft potential (e.g., Salem et al. 2001;Wilson et al. 2014a) (details are discussed in Appendix B), giving f (0sc) and f (0csc) .

3.
Then f (0sc) and f (0csc) are transformed into the ion bulk flow rest frame (e.g., Compton & Getting 1935;Ipavich 1974) following the methods described in  using a relativistically correct Lorentz transformation. The data are then interpolated onto a regular grid using Delaunay triangulation in the plane defined by the quasi-static magnetic field, B o , along the horizontal and transverse components of the ion bulk flow velocity, V i , i.e., (B o The result is a two-dimensional gyrotropic VDF, f (0) , and the associated Poisson errors/uncertainties, f (0 c) , both as functions of the parallel, V ‖ , and perpendicular, V ⊥ , velocity with respect to B o .

4.
Numerous weighting schemes were tried, and the best results (for Wind/3DP) were achieved by defining = f (0c) −2 for the weights. 29

a.
It is important to note that the fit is not done for all components simultaneously. This was the initial approach but proved to require stringent constraints for nearly all fit parameters, and the software exited before all fit parameters were varied owing to numerical instabilities 32 (e.g., Liavas & Regalia 1999), discussed in Appendix D.

b.
Thus, the core fit, f (core) , is performed first, and then the model result is subtracted from the data to yield the first residual, f (1) .

c.
The halo fit, f (halo) , is next but only to the side of f (1) opposite to that expected for the strahl/beam, where the latter is defined as the antisunward direction along B o . The entire two-dimensional halo fit is then subtracted from f (1) to yield the second residual, f (2) , i.e., both sides are subtracted, but only one side is used for the fit.
29 Several approaches were tried for the values, but the most reliable and robust was to use Gaussian weights on Poisson errors. Reliable and robust here mean that the fitting software required the fewest number of constraints and user-imposed limits to find fit parameters that well represent the observations. 30 That is, the data are not fit to two one-dimensional cuts of a two-dimensional VDF separately, but rather both dimensions are fit simultaneously. 31 It should also be noted that initial approaches tried to fit all electron components simultaneously but failed. Later approaches tried to fit the combination of only the core and halo simultaneously, but again the analysis was too unstable. Thus, the final approach fit to each component sequentially from core to beam/strahl. 32 There is also an issue of threshold tests for convergence. The software allows the user to define the thresholds for various gradients in the Jacobian. If the gradient magnitudes fall below these thresholds, the software exits with a specific fit status parameter associated with the specific threshold. For numerous reasons, the initial approach of fitting to all three components simultaneously prevented accurate fit results owing to these thresholds being satisfied too early in the iteration process.

d.
The beam/strahl fit, f (beam) , is last and fit to only the side of f (2) that is in the antisunward direction along B o .

6.
Not all VDFs will have fit results for all three components. In fact, f (beam) is often not found either because f (halo) left too few finite elements in f (2) or for numerical instability reasons (discussed in Appendix D).
All model functions are defined with six input parameters to be varied by the LMA software in the following order: Initial guesses are defined for all elements of PARAM that are specific to each shock event determined through an iterative trial-and-error approach. For each event, a zeroth-order guess is used on a subset of all VDFs, and the PARAM arrays for each component are adjusted accordingly to maximize the number of stable fit results for all components. Note that the PARAM arrays for each component differ depending on whether the VDF is located upstream or downstream of the shock ramp. In stronger shocks, the function used also varies (i.e., use symmetric bi-self-similar upstream and asymmetric bi-self-similar downstream).
the beam or halo fits defined by small T es,j and κ es . That is, "spiky" solutions are defined as those satisfying (k es ⩽ 3) ∧((T es,‖ ⩽ 11.8) ∨ (T es,⊥ ⩽ 11.8)) for model fit parameters. As evidenced by Figures 4-6, the χ tot 2 parameter alone is not necessarily an accurate measure of the quality of the fit.
An unexpected nuance arose during the development and testing of the software. The typical phase-space density of any given element of f (0) for electrons near 1 au varies from ~10 −18 to 10 −8 cm −3 km −3 s +3 . The LMA software uses a combination of gradients by constructing a Jacobian matrix of the input model fit function. 34 This is problematic when the magnitude of the input data and output model function are much much less than unity as results in numerical instabilities (e.g., Liavas & Regalia 1999). That is, the partial derivative of a number on the order of 10 −18 with respect to a number slightly greater than unity can produce exceedingly small gradients.
While the limits of double precision are not, in general, challenged by such computations, the LMA software (Markwardt 2009) was designed such that all the inputs are near unity.
The solution was to multiply by a constant offset to increase the contrast in the Jacobian components that are used to minimize χ 2 . A consequence of this approach is that the output χ 2 , f (m) , and 1σ error estimates of the fit parameters must be renormalized by this offset factor. The more standard approach is to perform the fit in logarithmic space, which reduces the dynamic range of the data. However, as discussed in Appendix E, this does not necessarily produce better fit results.
The above approach worked well except for cases with so-called flat-top distributions (e.g., Feldman et al. 1983a;Thomsen et al. 1987), modeled using the self-similar distributions (e.g., Dum et al. 1974;Dum 1975;Horton et al. 1976;Horton & Choi 1979;Jain & Sharma 1979;Goldman 1984) given by either Equation (3(a)) or Equation (4(a)). In cases where the phase-space densities were independent of energy for the core, the use of the weights above was not sufficient to constrain the fits. In these cases, shock-specific constraints/limits were imposed on the least number of fit parameters necessary to reliably and robustly produce good results (see the Supplemental Material ASCII files in Wilson et al. 2019c, described in Appendix F for list of constraints by shock). Figure 10 shows a comparison of three different fit results to illustrate the validity of the method used herein. Given the hindsight and statistics of the results from the present analysis, more refined constraints and better initial guesses were available. The fit shown in 34 That is, the partial derivatives are with respect to the fit parameters, not the velocity coordinates.

Wilson et al. Page 33
Astrophys J Suppl Ser. Author manuscript; available in PMC 2020 July 03. panels (b) and (c), referred to as the test fit from here on, was found following the automated method used for the 52 events examined in this study, i.e., the software is given initial guesses for parameters and constraints defined by knowns like n e and then allowed to find the best fit. The test fit results shown in panels (b) and (c) were then used as initial guesses (first perturbed, of course) on the same VDF to compare the method used herein (referred to as the linear method) to the base-10 logarithm approach (referred to as the log method). A larger range of constraints were used to provide a more open parameter space. Thus, in the following a comparison between the linear and log methods is presented as an illustrative test.
Panels (d) and (e) show the fit results using the linear method with the new initial guesses and parameter constraints, while panels (f) and (g) show the log method fit results. Unexpectedly, the log method did much worse in the core fit than the linear method but did well for the halo and beam/strahl fits. The numerical fit results are as follows: Thus, one can see that the log method did not produce a better fit for this specific example, which was not the expected outcome. This is almost certainly a consequence of the large constraint ranges, and a better fit would be found for a tighter range. That is, this example is not meant to argue that the linear method is better than the log method. Rather, the example is meant to illustrate that the linear method is a viable approach.
A point should also be made about the initiation stability of the LMA software. During the course of fitting all the VDFs in the present study, it was found that the choice of initial guess parameters was critical. For instance, in the example shown in Figure  were perturbed by ~20%-30% away from these initial guesses, the log method would not initiate fit iterations owing to diverging deviates and/or diverging model results, i.e., the software could not establish an initial Jacobian. 35 Unexpectedly, the linear method was more tolerant of perturbed initial guess parameters. There are still several checks for each component fit to address this possible noninitiation error, but even so this sometimes did not fix the issue, which is one reason why not all VDFs had stable solutions.
Finally, a note about the 1σ uncertainties of every fit parameter. These values are not reported because it was found that they do not accurately or consistently reflect the quality of fit. For instance, the 1σ uncertainties of n eh and V Teh,‖ for the log method in the example VDF shown in Figure 10 (panels (f) and (g)) are ~19,988 km s −1 (i.e., ~617% error) and ~3.53 cm −3 (i.e., ~5163% error), respectively, even though χ h 2 2.51. The 1σ uncertainties for the same parameters but for the fit in panels (d) and (e) are ~110.1 km s −1 (i.e., ~3.1% error) and ~0.0047 cm −3 (i.e., ~8.5% error), and χ h 2 1.82. That is, the reduced chi-squared values differ by only ~39%, but the 1σ uncertainties differ by hundreds to thousands of percent. The 1σ uncertainties determined by the LMA software that are assigned to the output fit parameters are not representative of the actual uncertainties. The reason is related to the orthogonal basis constructed during the qr-factorization (ultimately used to minimize χ s 2 ), which is not the same basis as that for the fit parameters. The output uncertainties thus contain nonlinear convolution of 1σ uncertainties from potentially multiple fit parameters. The effect is analogous to electric field measurements from two antennas with differing noise levels. If the electric field data are rotated to a new coordinate basis from the original instrument basis, the resulting field components will have a nonlinear convolution of noise from the original components. Thus, the 1σ uncertainties were not used as errors for each parameter.
The 1σ errors are also forced to zero in the software when the fit value reaches a userdefined boundary/constraint/limit. This is reported in the fit constraints ASCII file described in Appendix F (i.e., under the heading "Peg" in the ASCII file). As previously stated, the δℛ value alone does not always characterize the quality of any given fit. Therefore, a combination of parameters were used to define fit quality flags (see Appendix F for definitions), which should be used for determining the reliability of any given fit.

Appendix F: Data Product
One of the primary purposes of this first part of this three-part study is to describe the methodology and nuances of the fit procedure to provide context and documentation for the resulting data product (Wilson et al. 2019c). This will serve as the reference document for use of the data product by the heliospheric and astrophysical communities. The nuances and details of the procedure are critical for reproducibility and quality control in the use of the data product described in this section. While Papers II and III discuss the statistics and analysis results in detail, this first part is critical for any statistical or physical interpretation of the data, and it includes analysis of the exponents and drifts.
The fit results are provided in two ASCII files. The first contains all fit parameters for the three electron components in addition to several other relevant parameters. The nonelectron data products are linearly interpolated to the midpoint time stamp of each electron VDF. The ASCII file contains a detailed header with descriptions and explanations of the parameters with associated units. The data included are as follows: UTC time of electron VDF midpoint time stamp; n p and n α measured by SWE (cm −3 ); n i measured by 3DP (cm −3 ); T p,j and T α,j measured by SWE (eV); T i,j measured by 3DP (eV); B o,j measured by MFI (nT); V p,j and V α,j measured by SWE (km s −1 ); V i,j measured by 3DP (km s −1 ); ϕ sc determined from fit process (eV); δℛ calculated from fit process (%); n es from 3DP fits (cm −3 ); T es,j from 3DP fits (eV); V oes,j from 3DP fits (km s −1 ); κ es , p es , and q es from 3DP fits (N/A); χ s 2 from 3DP fits (N/A); and the numeric fit status value for each electron component (N/A). The total reduced chi-squared values for all fits are also included in the ASCII file. The fit flags for each component fit are also included. Let Ξ ≡ Σ s χ s 2 ; then, the list is as follows:
The second ASCII file contains the fit constraints, initial guesses, whether the fit parameters reached a fit constraint boundary, the number of iterations required to reach a stable fit, the chi-squared of the fit, degrees of freedom of the inputs, and a two-letter code for the model function used.
Both ASCII files contain fit results even if they are not high-quality or reliable results, which can be determined from the combination of χ s 2 , χ tot 2 , and δℛ used to define the fit flags in the first ASCII file, as discussed previously. The entries with fill values (listed in the header) resulted because a stable fit was not found or the fit was determined to be "bad," as defined in Section 3.3 and Appendix C. When there is a significant discrepancy between n p and n i (e.g., differ by a factor exceeding ~40%), the more reliable/accurate of the two is n p . Under these circumstances, T i,j and V i,j should be subject to scrutiny as well. The model function used for the core is given in the second ASCII file.
Note that the second ASCII file will contain nonfill, fit values for the same parameters that are all fill values in the first ASCII file. Although many constraints were set as far from the expected values as possible to avoid a parameter from being limited during the fit, some were imposed after all the fits were found for a given shock crossing. These were imposed for physical reasons (see, e.g., Section 3.2) and to avoid issues during regridding and/or interpolation for comparison with other data sets (e.g., magnetic fields). These post-fit constraints are 1.5 < κ eh ⩽ 20, 1.5 < κ eb ⩽ 20, 0 ⩽ n eh /n ec ⩽ 0.75, 0 ⩽ n eb /n ec ⩽ 0.50, 0.0 ⩽ n eb /n eh ⩽ 3.0, 11.4 eV ⩽ T eh,j ⩽ 285 eV, and 11.4 eV ⩽ T eb,j ⩽ 285 eV. All statistics and fit results presented herein are with respect to the first ASCII file values, but we include all the fit results in the second ASCII file for reference. This is because some of our post-fit constraints eliminated good fits like that shown in Figure 6, which failed the n eb /n eh < 3 test.
Most of the fits that failed this specific test were clearly bad fits, but not all.
The purpose of providing the detailed inputs for the fit results is for reproducibility and for quality control/sanity checks for researchers in the heliospheric and astrophysical communities interested in future use of the data. The data product will benefit current and future missions like Parker Solar Probe, in addition to providing a statistical comparison with astrophysical shocks, which is currently not available.  Examples of one-dimensional cuts through multiple model VDF functions to illustrate the functional dependence on various parameters. The top row (panels (a) through (c)) shows the dependence on the thermal speed, denoted generically as V th here. The bottom row (panels (d) and (e)  Example IP shock crossing observed on 1996 April 2 by the Wind spacecraft. The panels are as follows from top to bottom: |B o | (nT), B o (nT, GSE); value of spacecraft potential used for fits ϕ sc (eV); n p (red line) and 100 × n α (blue line) (cm −3 , SWE); s ec (blue circles), κ eh (green circles), and κ eb (magenta circles); n ec values (blue circles) and uncertainty (red error bars) (cm −3 , 3DP fit); n eh (cm −3 , 3DP fit); and n eb (cm −3 , 3DP fit). The error bars for the four electron fit parameter panels are defined by the percent deviation discussed in Section 3.3. The error for this date satisfied 0.2% < δℛ ⩽ 54% with a median of 10.3%.    Two superposed epoch analysis plots of the total electron current density, j e,tot (Mm s −1 ) (top panel), and normalized values, j e,tot /(n e V Tec,tot ) (%) (bottom panel), vs. seconds from the shock ramp center. Shown in each panel are the lower (Q1) and upper (Q2) quartiles as magenta lines, the mean as a red line, and the median as an orange line for all data. That is, the lines are computed for the entire set of data, not at each time stamp. For reference, the axis ranges were defined as 110% of the maximum of the absolute value of X 2.5 and X 97.5 , where X 2.5 and X 97.5 are the bottom 2.5th and top 97.5th percentiles. Wilson et al. Page 48 Astrophys J Suppl Ser. Author manuscript; available in PMC 2020 July 03.

Figure 8.
Histograms of the exponents (top row) and bulk drift velocity magnitudes (bottom row) for the different electron components for all time periods as percentage of total counts. Panel (a) shows the κ es values for the core (violet), halo (blue), and beam/strahl (red) components. Panel (b) shows the s ec for the core (violet). Panel (c) shows the p ec (blue) and q ec (red) values for the core. Panels (d)-(f) show the magnitude of the parallel (violet) and perpendicular (blue) drift velocities for the core, halo, and beam/strahl components, respectively. The statistics for the exponents are listed in Table 2. Note that the tick marks are individually labeled in all panels. Wilson et al. Page 49 Astrophys J Suppl Ser. Author manuscript; available in PMC 2020 July 03.

Figure 9.
Spacecraft potential, ϕ sc , shown vs. the total ion density, n i , observed by the Wind/3DP ion electrostatic analyzer (PESA Low). The top panel shows the value of ϕ sc (eV) determined iteratively, as described in this appendix, vs. n i (cm −3 ), where the color code is defined by the IP shock data given in the lower left corner. The bottom panel shows the same data, but now ϕ sc is offset by the detector minimum energy, E min , and divided by the constant 5.0 to keep the magnitudes near unity. The E min are color-coded and date-specific, as   Astrophys J Suppl Ser. Author manuscript; available in PMC 2020 July 03.