Electron energy partition across interplanetary shocks: I. Methodology and Data Product

TL;DR

This study analyzes electron velocity distribution functions near interplanetary shocks, demonstrating that a self-similar distribution better fits core electrons than a Maxwellian, revealing insights into particle scattering processes.

Contribution

It introduces a novel statistical methodology for fitting electron VDFs with self-similar functions, showing their superiority over Maxwellians in modeling core electrons.

Findings

Core electron distributions are better fit by self-similar functions than Maxwellians.

Kappa exponents range: 5.40–10.2 for core, 3.58–5.34 for halo, 3.40–5.16 for beam/strahl.

Detailed fit procedures and data products are described.

Abstract

Analysis of 15314 electron velocity distribution functions (VDFs) within $\pm$2 hours of 52 interplanetary (IP) shocks observed by the \emph{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 velocity distribution function, while both the halo and beam/strahl components were best fit to bi-kappa velocity distribution function. This is the first statistical study to show that the core electron distribution is better fit to a self-similar velocity distribution function 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 range of values defined by the lower and upper quartiles for the kappa exponents are $\kappa{\scriptstyle_{ec}}$ $\sim$ 5.40--10.2 for the core, $\kappa{\scriptstyle_{eh}}$ $\sim$ 3.58--5.34 for the halo, and $\kappa{\scriptstyle_{eb}}$ $\sim$ 3.40--5.16 for the beam/strahl. The lower-to-upper quartile range of symmetric bi-self-similar core exponents are $s{\scriptstyle_{ec}}$ $\sim$ 2.00--2.04, and asymmetric bi-self-similar core exponents are $p{\scriptstyle_{ec}}$ $\sim$ 2.20--4.00 for the parallel exponent, and $q{\scriptstyle_{ec}}$ $\sim$ 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.

Figures

Tables

Table 1. Table 1: Shock Parameters

Parameter	$X{}_{min}$   aaminimum	$X{}_{max}$   bbmaximum	$\overline{X}$   ccmean	$\tilde{X}$   ddmedian	$X{}_{25\%}$   eelower quartile	$X{}_{75\%}$   ffupper quartile	$\sigma {}_x$   ggstandard deviation
$⟨|𝐁{}_o{}^{}|⟩{}_{up}$ [nT]	1.04	17.4	5.96	5.59	3.99	7.10	3.01
$⟨n{}_i{}^{}⟩{}_{up}$ [$c{m}^{-3}$]	0.60	21.3	8.34	8.00	3.70	12.1	5.32
$⟨\beta {}_{tot}{}^{}⟩{}_{up}$ [N/A]	0.03	3.86	0.50	0.38	0.19	0.60	0.60
$⟨|V{}_{shn}{}^{}|⟩{}_{up}$ [km/s]	155	699	460	456	383	535	123
$⟨|U{}_{shn}{}^{}|⟩{}_{up}$ [km/s]	36.9	401	126	110	83.3	145	70.2
$\theta {}_{Bn}$ [deg]	17.1	88.6	56.8	54.6	42.7	73.3	19.5
$⟨M{}_A{}^{}⟩{}_{up}$ [N/A]	1.06	15.6	2.79	2.41	1.86	3.06	2.10
$⟨M{}_f{}^{}⟩{}_{up}$ [N/A]	1.01	6.39	2.12	1.86	1.58	2.35	0.94
$⟨M{}_f{}^{}⟩{}_{up}{}^{}/M{}_{cr}$ [N/A]	0.41	5.14	1.08	0.91	0.77	1.19	0.70
$⟨M{}_f{}^{}⟩{}_{up}{}^{}/M{}_{ww}$ [N/A]	0.06	2.49	0.36	0.18	0.11	0.32	0.51
$⟨M{}_f{}^{}⟩{}_{up}{}^{}/M{}_{gr}$ [N/A]	0.04	1.91	0.28	0.14	0.09	0.25	0.39
$⟨M{}_f{}^{}⟩{}_{up}{}^{}/M{}_{nw}$ [N/A]	0.04	1.76	0.26	0.13	0.08	0.23	0.36

Table 2. Table 2: Electron Exponent Parameters

Exponent	$X{}_{min}$   aaminimum	$X{}_{max}$   bbmaximum	$\overline{X}$   ccmean	$\tilde{X}$   ddmedian	$X{}_{25\%}$   eelower quartile	$X{}_{75\%}$   ffupper quartile
All: 15,210 VDFs
$\kappa {}_{ec}$	2.14	100.0	9.15	7.92	5.40	10.2
$s{}_{ec}$	2.00	3.00	2.03	2.00	2.00	2.04
$p{}_{ec}$	2.00	5.43	3.09	3.00	2.20	4.00
$q{}_{ec}$	2.00	3.29	2.24	2.00	2.00	2.46
$\kappa {}_{eh}$	1.51	19.7	4.62	4.38	3.58	5.34
$\kappa {}_{eb}$	1.52	20.0	4.57	4.17	3.40	5.16
Upstream Only: 6546 VDFs
$\kappa {}_{ec}$	2.14	100.0	9.15	7.92	5.40	10.2
$s{}_{ec}$	2.00	2.31	2.01	2.00	2.00	2.03
$p{}_{ec}$	N/A	N/A	N/A	N/A	N/A	N/A
$q{}_{ec}$	N/A	N/A	N/A	N/A	N/A	N/A
$\kappa {}_{eh}$	1.52	18.4	4.16	4.10	3.25	4.83
$\kappa {}_{eb}$	1.52	19.6	4.22	3.81	3.25	4.70
Downstream Only: 8664 VDFs
$\kappa {}_{ec}$	N/A	N/A	N/A	N/A	N/A	N/A
$s{}_{ec}$	2.00	3.00	2.05	2.01	2.00	2.06
$p{}_{ec}$	2.00	5.43	3.09	3.00	2.20	4.00
$q{}_{ec}$	2.00	3.29	2.24	2.00	2.00	2.46
$\kappa {}_{eh}$	1.51	19.7	4.94	4.62	3.80	5.70
$\kappa {}_{eb}$	1.53	20.0	4.82	4.45	3.61	5.44
$⟨M{}_f{}^{}⟩{}_{up}$ 3 Only: 12,988 VDFs
$\kappa {}_{ec}$	2.14	100.0	9.02	6.83	4.40	9.93
$s{}_{ec}$	2.00	3.00	2.03	2.00	2.00	2.04
$p{}_{ec}$	2.00	5.43	3.10	3.00	2.18	4.00
$q{}_{ec}$	2.00	3.14	2.26	2.01	2.00	2.49
$\kappa {}_{eh}$	1.51	19.7	4.54	4.34	3.58	5.26
$\kappa {}_{eb}$	1.52	20.0	4.62	4.20	3.46	5.19
$⟨M{}_f{}^{}⟩{}_{up}$ $\ge$ 3 Only: 2222 VDFs
$\kappa {}_{ec}$	4.32	27.2	9.30	8.60	6.89	10.4
$s{}_{ec}$	2.00	2.30	2.03	2.00	2.00	2.08
$p{}_{ec}$	2.00	5.00	3.08	2.50	2.18	4.00
$q{}_{ec}$	2.00	3.29	2.16	2.00	2.00	2.50
$\kappa {}_{eh}$	1.60	19.2	5.06	4.68	3.62	6.05
$\kappa {}_{eb}$	1.52	18.8	4.25	3.84	2.89	4.94
$\theta {}_{Bn}$ $>$ 45∘ Only: 10,940 VDFs
$\kappa {}_{ec}$	4.05	27.2	7.77	7.18	4.84	9.11
$s{}_{ec}$	2.00	2.31	2.02	2.00	2.00	2.05
$p{}_{ec}$	2.00	5.43	3.00	2.62	2.17	4.00
$q{}_{ec}$	2.00	3.29	2.28	2.04	2.00	2.56
$\kappa {}_{eh}$	1.51	19.7	4.73	4.44	3.67	5.47
$\kappa {}_{eb}$	1.52	20.0	4.67	4.20	3.33	5.33
$\theta {}_{Bn}$ $\le$ 45∘ Only: 4270 VDFs
$\kappa {}_{ec}$	2.14	100.0	16.0	11.7	10.0	14.5
$s{}_{ec}$	2.00	3.00	2.06	2.00	2.00	2.04
$p{}_{ec}$	2.00	4.28	3.29	4.00	4.00	4.28
$q{}_{ec}$	2.00	3.00	2.14	2.00	2.00	2.16
$\kappa {}_{eh}$	1.55	19.4	4.32	4.18	3.37	5.09
$\kappa {}_{eb}$	1.53	16.5	4.31	4.10	3.57	4.82

Table 3. Table 3: Spacecraft Potential Statistics

$\varphi {}_{sc}$ [eV]	$X{}_{min}$   aaHeader symbols match that of Table 2	$X{}_{max}$	$\overline{X}$	$\tilde{X}$	$X{}_{25\%}$	$X{}_{75\%}$
All: 15,144 Finite Values	1.01	26.7	7.05	6.70	5.45	7.84
Upstream Only: 6511 Finite Values	1.01	26.7	7.14	6.80	5.34	7.82
Downstream Only: 8633 Finite Values	1.92	24.8	6.43	6.45	4.00	7.37
$⟨M{}_f{}^{}⟩{}_{up}$ 3 Only: 12,932 Finite Values	1.01	26.7	6.99	6.61	5.44	7.70
$⟨M{}_f{}^{}⟩{}_{up}$ $\ge$ 3 Only: 2212 Finite Values	3.58	12.0	7.35	6.90	5.50	9.63
$\theta {}_{Bn}$ $>$ 45∘ Only: 10,894 Finite Values	1.01	26.7	6.70	6.49	5.35	7.38
$\theta {}_{Bn}$ $\le$ 45∘ Only: 4250 Finite Values	3.53	17.6	7.94	7.14	6.10	10.2

Equations

$$f\left(V{}_{\parallel },V{}_{\perp }\right)$$

$$A{}_M$$

$$f\left(V{}_{\perp },V{}_{\parallel }\right)$$

$$A{}_{\kappa }$$

$$B{}_{\kappa }$$

$$f\left(V{}_{\parallel },V{}_{\perp }\right)$$

$$A{}_{SS}$$

$$f\left(V{}_{\parallel },V{}_{\perp }\right)$$

$$A{}_{AS}$$

$$\delta \mathcal{R}=|1-\tilde{\mathcal{R}}|\cdot 100\%$$

$$T{}_{eff,j}$$

$$T{}_{s,tot}$$

$$V{}_{Ts,j}$$

$$\Omega {}_{cs}$$

$$\omega {}_{ps}$$

$$\lambda {}_{De}$$

$$\rho {}_{cs}$$

$$\lambda {}_s$$

$$\beta {}_{s,j}$$

$$\beta {}_{eff,j}$$

$$n{}_e$$

$$\chi {{}_s}^2=\sum _{i=0}^{N-1}\sum _{j=0}^{M-1}{\left(f{{}_{ij,s}}^{\left(0\right)}-f{{}_{ij,s}}^{\left(mod\right)}\right)}^2|\mathcal{W}{}_{ij,s}|$$

Full text

plasmas — shock waves — (Sun:) solar wind — Sun: coronal mass ejections (CMEs)

††journal: ApJS

1 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^{\prime}$ and $s$, are not in thermodynamic or thermal equilibrium, i.e., $\left(T{{}_{s^{\prime}}}/T{{}_{s}}\right){{}_{tot}}$ $\neq$ 1 for $s^{\prime}$ $\neq$ $s$, and the relevant scale lengths are orders of magnitude smaller than the collisional mean free path (e.g., Wilson III et al., 2018). Therefore, for any process dependent upon scales like the thermal gyroradii, $\rho{{}_{cs}}$, or inertial lengths, $\lambda{{}_{s}}$, the media is considered collisionless (see Appendix A for definitions). That the solar wind is a non-equilibrium, weakly collisional, kinetic gas results in multi-component velocity distribution functions (VDFs) for both ions (e.g., Kasper et al., 2006, 2012, 2013; Maruca et al., 2011; Maruca & Kasper, 2013; Wicks et al., 2016) and electrons (e.g., Lin, 1998; Maksimovic et al., 1997, 1998; Pierrard et al., 1999, 2001; Pulupa et al., 2014a; Schwartz & Marsch, 1983; Štverák et al., 2008, 2009).

The electron VDFs in the solar wind below $\sim$1 keV are comprised of a cold core with energies $E{{}_{ec}}$ $\lesssim$ 15 eV (e.g., Bale et al., 2013; Maksimovic et al., 1997, 1998; Pilipp et al., 1987a, b, c, 1990; Pulupa et al., 2014a), a hot, tenuous halo with $E{{}_{eh}}$ $\gtrsim$ 20 eV (e.g., Maksimovic et al., 1997, 1998; Pulupa et al., 2014a; Štverák et al., 2008, 2009), and an anti-sunward, field-aligned beam called the strahl with $E{{}_{eb}}$ $\sim$few 10s of eV (e.g., Bale et al., 2013; Crooker et al., 2003; Graham et al., 2017, 2018; Horaites et al., 2018; Štverák et al., 2009) (e.g., see Figure 1 for illustrative example). The electrons also dominate the solar wind heat flux (e.g., Bale et al., 2013; Crooker et al., 2003; Pagel et al., 2005, 2007), arising from the consistent skewness in the VDFs, specifically the halo and/or strahl components. Note that there also exists a suprathermal super halo with $E{{}_{esh}}$ $\gtrsim$ 1 keV (e.g., Lin, 1998; Wang et al., 2012, 2015), but these higher energy electrons are not examined herein.

The three electron components below $\sim$1 keV are predicted and observed to be coupled through multiple processes from wave-particle interactions (e.g., Pierrard et al., 2011, 2016; Phillips et al., 1989a, b; Saito & Gary, 2007; Saito et al., 2008; Vocks & Mann, 2003; Vocks et al., 2005; Yoon, 2014; Yoon et al., 2012, 2015, 2016) to adiabatic transport effects (e.g., Schwartz & Marsch, 1983) to collisional effects (e.g., Pilipp et al., 1987a, b, c; Schwartz & Marsch, 1983). They have also been shown to behave differently across collisionless shocks depending upon shock strength (e.g., Wilson III et al., 2009, 2010).

An illustrative example, showing the three electron components typically observed in the solar wind near 1 AU below $\sim$1.2 keV, is shown in Figure 1. The components parameters are exaggerated111The following were enhanced to increase contrast and for ease of viewing differences: parallel core temperature, perpendicular halo temperature, and parallel core drift speed. for illustrative purposes but based upon 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/strahl by a bi-kappa VDF (see Section 3.1). In this case, the self-similar exponent reduced to 2 so the VDF reduced to a bi-Maxwellian (see Section 3.1). This example is phenomenologically consistent with the majority of solar wind electron VDFs (e.g., Phillips et al., 1989a, b; Pilipp et al., 1987a, b, c; Štverák et al., 2008, 2009).

Despite its collisionless, non-equilibrium nature the solar wind can support the existence of shock waves. That the particles are in neither thermal or thermodynamic equilibrium leads to a non-homogeneous partition of energy among not only electrons and ions but also among the components of each species, e.g., the core electrons do not have the same response as the halo to collisionless shock waves. The reason for the non-homogeneous partition of energy lies in the energy-dependent mechanisms that transfer the bulk flow kinetic energy lost across the shock ramp to other forms like heat or particle acceleration (e.g., see Coroniti, 1970; Kennel et al., 1985; Sagdeev, 1966; Tidman & Krall, 1971; Treumann, 2009; Wilson III, 2016; Wilson III et al., 2017, and references therein). The mechanisms can also be dependent upon pitch-angle and species (e.g., Artemyev et al., 2013, 2014, 2015, 2016, 2017a, 2017b, 2018; Sagdeev, 1966). Most collisionless shocks are subsonic to electrons, yet electrons still respond to the shock showing even Mach number dependent effects (e.g., Feldman et al., 1982, 1983b, 1983a; Masters et al., 2011; Thomsen et al., 1985, 1987, 1993; Wilson III et al., 2010). This is all further complicated by recent observations showing that the evolution of the electron VDF through a collisionless shock is not a trivial, uniform inflation of the entire distribution, but a multi-stage process that deforms and redistributes/exchanges energy for different energies and pitch-angles at different stages (e.g., Chen et al., 2018; Goodrich et al., 2018, 2019). There is no currently known way to quantify these non-homogenous changes to capture the energy- and pitch-angle-dependent effects, therefore the next best systematic approach for a statistical study is to parameterize the electron components by their velocity moments. This is further supported by the fact that nearly all theories describing the evolution of electron VDFs rely upon either the velocity moments or a model velocity distribution function (e.g., Livadiotis, 2015, 2017; Nicolaou et al., 2018; Schunk, 1975, 1977; Schwartz & Marsch, 1983; Schwartz et al., 1988; Shizgal, 2018).

In this first part of a multi-part study we describe the methodology and numerical analysis techniques used to model the solar wind eVDFs below $\sim$1.2 keV observed by the Wind spacecraft near 1 AU around 52 interplanetary (IP) shocks. This is the first statistical study to show that the core electron distribution is better fit to a self-similar velocity distribution function than a bi-Maxwellian under all conditions. The analysis differs from numerous previous studies in its approach and the model functions used, each of which are justified herein using physically significant arguments. A benefit of the analysis is an improved, semi-analytic relationship between the spacecraft potential and ion number density. The paper also includes procedural documentation to disclose the nuances and issues associated with applying a nonlinear least squares fitting algorithm to in situ VDF data in the solar wind. This serves as a reference for use of the resulting data product described herein. In Paper II (Wilson III et al., 2019a) the statistical results of the model fits are presented with comparison to previous studies and associated discussions. In Paper III (Wilson III et al., 2019b) the analysis and interpretation of the model fit results are presented.

This paper is outlined as follows: Section 2 introduces the data sets and event selection; Section 3 introduces the methodology of the fit analysis, model functions, parameter constraints, quality control, and summary of fit results; Section 4 discusses the statistics of the fit exponents and drift velocities; and Section 5 discusses the results and interpretations with reference to further analysis in the following Papers II and III. Appendices are also included to provide additional details of the parameter definitions (Appendix A), spacecraft potential and detector calibration (Appendix B), numerical analysis procedure (Appendix C), numerical instabilities (Appendix D), direct fit method comparisons (Appendix E), and the data product produced by this effort (Appendix F).

2 Data Sets and Event Selection

In this section we introduce the instrument data sets and shock database used to examine the data observed by the Wind spacecraft (Harten & Clark, 1995) near 1 AU. The data described herein spanned from 00:55:40 UTC on 1995-02-26 to 23:04:00 UTC on 2000-02-20 (see Supplemental Material for list of dates). The symbol/parameter definitions are found in Appendix A.

Quasi-static magnetic field vectors ($\mathbf{B}{{}_{o}}$) were measured by the Wind/MFI dual, triaxial fluxgate magnetometers (Lepping et al., 1995) using the three second cadence data for each particle distribution. The components/directions of some parameters are defined with respect to $\mathbf{B}{{}_{o}}$ using the subscript $j$. That is, the parallel ($j$ $=$ $\parallel$) and the perpendicular components ($j$ $=$ $\perp$) of any vector or pseudo-tensor (e.g., temperature) are defined with respect to $\mathbf{B}{{}_{o}}$.

The electron velocity distribution functions (VDFs) were measured by the Wind/3DP low energy (i.e., few eV to $\sim$1.2 keV) electron electrostatic analyzer (Lin et al., 1995) or EESA Low. The instrument operated in both burst and survey modes for the data presented herein, which has cadences of $\sim$3 seconds and $\sim$24–78 seconds, respectively. The energy and angular resolutions are commandable but the instrument typically operates with $\Delta\ E$/$E$ $\sim$ 20% and $\Delta\ \phi$ $\sim$ 5*∘–22.5∘* depending on the poloidal anode222The ecliptic plane bins have higher angular resolution than the zenith. (e.g., see Wilson III et al., 2009, 2010, for instrument details).

The EESA Low measurements are contaminated with photoelectrons from the spacecraft, something for which must be accounted to obtain accurate velocity moments or any other results. The details of how the spacecraft potential, $\phi{{}_{sc}}$, was numerically determined for each VDF is described in Appendix B. The VDFs are transformed into the ion frame prior to any fit using relativistically correct Lorentz transformations, where the steps are as follows: (1) convert the units of the VDFs to phase space density [# cm*-3* s*+3* km*-3*]; (2) correct the energies by $\phi{{}_{sc}}$; (3) convert the energy-angle bins to velocity coordinates; and (4) transform the velocities into the ion rest frame using proper Lorentz transformations. Nothing need be done to VDFs once in units of phase space density as phase space density is a Lorentz invariant (Van Kampen, 1969) (see Appendices B and C for details).

We also examined solar wind proton and alpha-particle velocity moments determined by a nonlinear least squares fitting algorithm (e.g., Kasper et al., 2006; Maruca & Kasper, 2013) observed by the Wind/SWE Faraday Cups (Ogilvie et al., 1995). Similar quality requirements for the SWE results to that discussed in Wilson III et al. (2018) were used herein.

The VDFs examined are found within $\pm$2 hours of 52 IP shocks found in the Wind shock database from the Harvard Smithsonian Center for Astrophysics333https://www.cfa.harvard.edu/shocks/wi_data/. Of those 52 IP shocks, there were 16 quasi-parallel ($\theta{{}_{Bn}}$ $\leq$ 45*∘), 36 quasi-perpendicular ($\theta{{}_{Bn}}$ $>$ 45∘*), 45 low Mach number ($\langle M{{}_{f}}\rangle{{}_{up}}$ $<$ 3), and 7 high Mach number ($\langle M{{}_{f}}\rangle{{}_{up}}$ $\geq$ 3) shocks. The shock parameters for the 52 IP shocks examined in this three-part set of papers are shown in Table 1 (see Supplemental Material for full list of values for each shock). The IP shocks examined were selected because of burst mode 3DP availability. See Appendix A for definitions of symbols and/or parameters.

3 Fit Methodology

This section (and Appendix C) introduces and discusses the nuances of the approach and software used to numerically compute the model fit parameters for every electron VDF examined. The nuances and details are provided for reproducibility and documentation for the data product discussed in Appendix F.

The data are fit to a user defined model function using a nonlinear least squares fitting algorithm called the Levenberg-Marquardt Algorithm (LMA) (Moré, 1978). The generalized LMA software used for the present study is called MPFIT (Markwardt, 2009). The specific details for its use are outlined in Appendix C.

The components of the electron VDFs are fit to bi-Maxwellian, bi-kappa, or bi-self-similar model functions (see Section 3.1). The components can be fit separately because the solar wind is a non-equilibrium, weakly collisional, kinetic gas. That is, in the absence of a magnetic field, each electron component could, in principle, stream past the other components for nearly an astronomical unit without significant interaction. Thus, there is physical justification to fit to the sum of three model functions (see Appendix C for details).

Given that the bi-self-similar reduces to the bi-Maxwellian in the limit as the exponential argument goes to 2 and that it consistently yielded lower reduced chi-squared values, $\tilde{\chi}{{}_{s}}^{2}$, the symmetric bi-self-similar function was used as the default core model function. In the downstream of strong (i.e., $\langle M{{}_{f}}\rangle{{}_{up}}$ $\gtrsim$ 2.5) IP shocks it was found that the asymmetric bi-self-similar function produced the best results and so was the default core model function444The parallel and perpendicular profiles at low energies differ greatly in these regions and required the use of the asymmetric function to accommodate the differences. Using a symmetric function resulted in very poor fit qualities, as defined in Section 3.3.. Note that of all the core VDFs fit to a symmetric bi-self-similar function, $\sim$80.5% that satisfied 2.0 $\leq$ $s{{}_{ec}}$ $\leq$ 2.05. That is, the majority of the distributions would be nearly indistinguishable 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, $\phi{{}_{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 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 is discussed in Section 3.3.

A total of 15,314 electron VDFs were observed by the Wind spacecraft within $\pm$2 hours 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($\sim$98%) core fits, 13,871($\sim$91%) halo fits, and 9567($\sim$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.

3.1 Velocity Distribution Functions

This section introduces and defines the model functions used to fit to the particle velocity distribution functions (VDFs) in this study with examples provided to illustrate shape and dependences on parameters.

The most common velocity distribution function (VDF) used to model particle VDFs in space plasmas is the bi-Maxwellian (e.g., Feldman et al., 1979b, a, 1983a; Kasper et al., 2006), given by:

[TABLE]

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 A1c, $V{{}_{j}}$ is the velocity ordinate of the jth component, and $n{{}_{o}}$ is the number density.

The second most popular model VDF is the kappa distribution. The kappa velocity distribution has gained popularity in recent years owing to improvements in particle detectors and the ubiquitous non-Maxwellian tails observed for both ions and electrons (e.g., Lazar et al., 2015b, a, 2016, 2017, 2018; Livadiotis, 2015; Livadiotis et al., 2018; Mace & Sydora, 2010; Pulupa et al., 2014b; Saeed et al., 2017; Shaaban et al., 2018), but references to and use of kappa or kappa-like (e.g., modified Lorentzian) distributions have been around for decades (e.g., Feldman et al., 1983b; Maksimovic et al., 1997; Salem et al., 2003; Vasyliunas, 1968). 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 power-law particle distribution is given by a bi-kappa VDF (e.g., Livadiotis, 2015; Mace & Sydora, 2010), for electrons here as:

[TABLE]

where $\Gamma\left(z\right)$ 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 upon $\kappa$.

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; Goldman, 1984; Horton et al., 1976; Horton & Choi, 1979; Jain & Sharma, 1979) or flows through disordered porous media (e.g., Matyka et al., 2016). The symmetric form is given by:

[TABLE]

Note that $V{{}_{Tj}}$ is again the most probable speed of a 1D Gaussian for consistency, i.e., it does not depend upon $s$. Further, one can see that Equation 3a reduces to Equation 1a in the limit where $s\rightarrow 2$. The function in Equation 3a will be referred to as the symmetric self-similar distribution function.

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:

[TABLE]

Again, this will reduce to a bi-Maxwellian in the limit where $p\rightarrow 2$ and $q\rightarrow 2$. Note that in the event that the the exponents $s$, $p$, or $q$ are not even integers, the velocity ordinates, $\left(V{{}_{\parallel}}-v{{}_{o\parallel}}\right)$ and $\left(V{{}_{\perp}}-v{{}_{o\perp}}\right)$, will become absolute values to avoid complex roots and negative values of $f\left(V{{}_{\parallel}},V{{}_{\perp}}\right)$. 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., Bale et al., 2013; Pulupa et al., 2014b; Štverák et al., 2008, 2009). 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., 1983b, a) 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 (e.g., see Equations 3a and 4a) 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; Goldman, 1984; Horton et al., 1976; Horton & Choi, 1979; Jain & Sharma, 1979), 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 upon the number density and thermal speed and that hot, tenuous plasmas are much more difficult to measure and accurately model. Examine the one-dimensional 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 km/s to 5500 km/s drops the peak phase space density by nearly two orders of magnitude. The cut line also passes the $\pm$20,000 km/s 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 one examines a more extreme example with $n{{}_{e}}$ $=$ 15 $cm^{-3}$ and $V{{}_{Te}}$ $=$ 10,000 km/s. In this case, the difference between the peak and the lowest phase space density within the $\pm$20,000 km/s velocity boundary would only be a factor of $\sim$55, i.e., slightly more than one order of magnitude.

For reference, the list of potential free parameters are as follows (see Appendix A for symbol definitions):

• Core

• –

$n{{}_{ec}}$

• –

$V{{}_{Tec,j}}$ or $T{{}_{ec,j}}$

• –

$v{{}_{oec,j}}$

• –

$s{{}_{ec}}$

• –

$p{{}_{ec}}$

• –

$q{{}_{ec}}$

• –

$\kappa{{}_{ec}}$

• Halo

• –

$n{{}_{eh}}$

• –

$V{{}_{Teh,j}}$ or $T{{}_{eh,j}}$

• –

$v{{}_{oeh,j}}$

• –

$\kappa{{}_{eh}}$

• Beam/Strahl

• –

$n{{}_{eb}}$

• –

$V{{}_{Teb,j}}$ or $T{{}_{eb,j}}$

• –

$v{{}_{oeb,j}}$

• –

$\kappa{{}_{eb}}$

For more details about derivation and normalization constants, see the Supplemental Material.

3.2 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 Supplemental Materials for shock parameters) observed by Wind on 1996-04-02 at 10:07:57.525 UTC. For this event, the plasma parameters are listed below in the following form Min–Max(Mean)[Median]

• Upstream

• –

$\lvert\mathbf{B}{{}_{o}}\rvert$ $\sim$ 0.53–3.14(1.96)[1.53] nT;

• –

$n{{}_{p}}$ $\sim$ 11.3–15.8(12.0)[11.9] $cm^{-3}$;

• –

$n{{}_{\alpha}}$ $\sim$ 0.06–0.18(0.10)[0.11] $cm^{-3}$;

• –

$s{{}_{ec}}$ $\sim$ 2.00–2.09(2.00)[2.00] N/A;

• –

$\kappa{{}_{eh}}$ $\sim$ 2.83–12.2(4.46)[4.40] N/A;

• –

$\kappa{{}_{eb}}$ $\sim$ 1.67–12.6(4.85)[5.10] N/A;

• –

$n{{}_{ec}}$ $\sim$ 10.7–13.0(11.7)[11.5] $cm^{-3}$;

• –

$n{{}_{eh}}$ $\sim$ 0.06–1.44(0.69)[0.54] $cm^{-3}$;

• –

$n{{}_{eb}}$ $\sim$ 0.02–0.17(0.09)[0.09] $cm^{-3}$;

• Downstream

• –

$\lvert\mathbf{B}{{}_{o}}\rvert$ $\sim$ 3.45–5.99(4.85)[5.19] nT;

• –

$n{{}_{p}}$ $\sim$ 14.9–19.7(18.0)[18.1] $cm^{-3}$;

• –

$n{{}_{\alpha}}$ $\sim$ 0.14–0.27(0.19)[0.19] $cm^{-3}$;

• –

$s{{}_{ec}}$ $\sim$ 2.00–2.07(2.01)[2.01] N/A;

• –

$\kappa{{}_{eh}}$ $\sim$ 2.72–6.96(4.39)[4.29] N/A;

• –

$\kappa{{}_{eb}}$ $\sim$ 2.74–7.27(4.45)[4.50] N/A;

• –

$n{{}_{ec}}$ $\sim$ 13.6–18.4(16.7)[16.8] $cm^{-3}$;

• –

$n{{}_{eh}}$ $\sim$ 0.02–2.53(0.56)[0.44] $cm^{-3}$;

• –

$n{{}_{eb}}$ $\sim$ 0.01–0.29(0.12)[0.11] $cm^{-3}$;

Note that there are two time periods after 11:00 UTC where a few fit results satisfy $n{{}_{eb}}/n{{}_{eh}}$ $\geq$ 1. Figure 3 is illustrative of some of the error analysis employed in the present study and 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 two-dimensional VDF (e.g., Maksimovic et al., 2005; Pulupa et al., 2014a, b). 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 non-separable 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 upon $n{{}_{s}}$, $V{{}_{Ts,\parallel}}^{-1}$, and $V{{}_{Ts,\perp}}^{-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 (e.g., see Equation 1a). 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}}$ $\geq$ 10,000 km/s 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 $\sim$1000 km/s 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., 2008b, a; Pollock et al., 2016).

Statistical studies of the solar wind have shown that the maximum range of the total electron temperature is $T{{}_{e,j}}$ $\sim$ 2.29–77.2 eV or $V{{}_{Te,j}}$ $\sim$ 450–2600 km/s (e.g., Wilson III et al., 2018). Previous studies have found that the electron halo temperatures satisfy $T{{}_{e,j}}$ $\sim$ 14–560 eV or $V{{}_{Teh,j}}$ $\sim$ 1100–7000 km/s (e.g., Feldman et al., 1975, 1978, 1979a; Lazar et al., 2017; Maksimovic et al., 1997, 2005; Skoug et al., 2000; Tao et al., 2016a, b). Previous studies have also found that the electron beam/strahl temperatures satisfy $T{{}_{eb,j}}$ $\sim$ 20–150 eV or $V{{}_{Te,j}}$ $\sim$ 1300–3600 km/s (e.g., Ogilvie et al., 2000; Tao et al., 2016a, b; Viñas et al., 2010). Thus, a range of allowed core thermal speeds from $\sim$1000 km/s to $\sim$10,000 km/s 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 km/s (e.g., Pulupa et al., 2014a). In the present work, most fit results show less than 50 km/s, i.e., only 1838 of 14847 or $\sim$12% have drift speeds exceeding 50 km/s, consistent with previous work555Note that in the present work the dipole correction to $\phi{{}_{sc}}$ was not applied, which affects the drift velocity and heat flux velocity moments. Thus, the core drift velocities in our work suffer the greatest from this correction.. In contrast, owing the physical interpretation of the strahl/beam component most (8848 of 9567 or $\sim$92%) have drift speeds in excess of 1000 km/s. The range of allowed core, halo, and beam/strahl drift speeds loosely ranged from $\sim$1000 km/s to $\sim$10,000 km/s 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 described in Appendix F for ranges for specific events). Note that $V{{}_{oes,\perp}}$ 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 more than few km/s. In some events, an explicit $V{{}_{oec,\perp}}$ was set as the initial guess values determined from examination of the distributions, but this is for a small minority of events (333 of 14847 or $\sim$2%).

It has also been empirically found that the EESA Low detector has issues when $n{{}_{ce}}$ $\lesssim$ 0.5 $cm^{-3}$ or $n{{}_{ce}}$ $\gtrsim$ 50 $cm^{-3}$ for typical solar wind thermal speeds666Technically, 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.. This is rarely an issue as only 41 of the 14847 or $\sim$0.3% VDFs analyzed have fit results falling outside the range $\sim$0.5–50 $cm^{-3}$. Note that the total electron density, $n{{}_{e}}$ $=$ $n{{}_{ec}}$ $+$ $n{{}_{eh}}$ $+$ $n{{}_{eb}}$ $\sim$ $n{{}_{e}}$ $=$ $n{{}_{p}}$ $+$ 2$n{{}_{\alpha}}$, 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 upon results found in previous studies near 1 AU (e.g., Feldman et al., 1975; Maksimovic et al., 1997, 2005; Pierrard et al., 2016; Skoug et al., 2000; Štverák et al., 2009; Tao et al., 2016a; Viñas et al., 2010).

In numerous previous studies that assumed a three component solar wind electron VDF near 1 AU (e.g., Maksimovic et al., 2005; Pulupa et al., 2014a, b; Štverák et al., 2009), constraints were sometimes assumed such as that the fits satisfy $n{{}_{eb}}/n{{}_{eh}}$ $<$ 1. There is no restriction on this ratio777The 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. imposed during the fit process and 1824 of 9313 or $\sim$20% of the fits satisfy $n{{}_{eb}}/n{{}_{eh}}$ $\geq$ 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 component888Note 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.. 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}}$ $\geq$ 1 (e.g., Štverák et al., 2009).

Another constraint that is often assumed/used is that the strahl/beam component be only anti-sunward along $\mathbf{B}{{}_{o}}$ (e.g., Maksimovic et al., 2005; Pulupa et al., 2014a, b; Štverák et al., 2009), 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 and the shock-reflected electron component in the events where this is directed sunward along $\mathbf{B}{{}_{o}}$ (this is very rare). If the shock-reflected electron component is directed anti-sunward they will be included in the beam/strahl fit (this is much more common). The net result for the former is a smaller $\left(T{{}_{\perp}}/T{{}_{\parallel}}\right){{}_{eh}}$ and on the latter a larger $\left(T{{}_{\perp}}/T{{}_{\parallel}}\right){{}_{eb}}$ and $n{{}_{eb}}$.

The lower bound of possible $\kappa{{}_{es}}$ values is defined for mathematical/physical reasons as being $\gtrsim$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., Lazar et al., 2017; Maksimovic et al., 1997; Pierrard et al., 2016; Štverák et al., 2009; Tao et al., 2016a, b). 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; Goldman, 1984; Horton et al., 1976; Horton & Choi, 1979; Jain & Sharma, 1979).

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 13867 or $\sim$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.

3.3 Quality Analysis

The initial approach was to use the reduced chi-squared value $\tilde{\chi}{{}_{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 $\tilde{\chi}{{}_{s}}^{2}$ $>$ 10 while others did not fit well at all despite having $\tilde{\chi}{{}_{s}}^{2}$ $\lesssim$ 1. The issue is partly related to the calibration of the detector and thus the quality of the $\mathcal{W}$ 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^{\left(0\right)}$ as the actual data and $f^{\left(m\right)}$ ($=$ $f^{\left(core\right)}$ $+$ $f^{\left(halo\right)}$ $+$ $f^{\left(beam\right)}$) as the total model fit results. Then one can define the ratio of these two parameters as $\mathcal{R}$ $=$ $f^{\left(0\right)}/f^{\left(m\right)}$, which is a two-dimensional array of values. Then one calculates the median of this array, $\tilde{\mathcal{R}}$, to determine the percent deviation given by:

[TABLE]

where $\delta\mathcal{R}$ is computed for each electron VDF. The values of $\delta\mathcal{R}$ 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 one-sigma uncertainties). The uncertainty of any variable calculated using these fit parameters was propagated assuming uncorrelated errors.

Note that the $\delta\mathcal{R}$ 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 at least 2 or higher can be used but low fit flags should be treated with caution. Only $\lesssim$1% of all core, halo, and beam/strahl fits had flags of 1 while $>$95% of core, $>$89% of halo, and $>$61% of beam/strah flags were at least 2.

Figure 4 shows an example VDF that had a low $\tilde{\chi}{{}_{s}}^{2}$ for each component and a $\delta\mathcal{R}$ $\sim$ 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 for both the halo and beam/strahl component. The fit results are as follows:

• •

$n{{}_{e\{c,h,b\}}}$ $=$ {15.43, 2.01, 0.056} $cm^{-3}$;

• •

$V{{}_{Te\{c,h,b\},\parallel}}$ $=$ {1959.6, 2500.0, 3964.7} $km\ s^{-1}$;

• •

$V{{}_{Te\{c,h,b\},\perp}}$ $=$ {1937.9, 2575.5, 4516.2} $km\ s^{-1}$;

• •

$V{{}_{oe\{c,h,b\},\parallel}}$ $=$ {$+$44.58, -0.00, -3898.7} $km\ s^{-1}$;

• •

$V{{}_{oe\{c,h,b\},\perp}}$ $=$ {-0.00, -0.00, -0.00} $km\ s^{-1}$;

• •

$\{s{{}_{ec}},\kappa{{}_{eh}},\kappa{{}_{eb}}\}$ $=$ {2.00, 4.58, 2.57}, where $s{{}_{ec}}$ is the self-similar exponent and $\kappa{{}_{es}}$ is the kappa value;

• •

$\tilde{\chi}{{}_{e\{c,h,b\}}}^{2}$ $=$ {1.07, 1.36, 0.41 };

• •

$\tilde{\chi}{{}_{tot}}^{2}$ $=$ 6.14; and

• •

Fit Flag {c,h,b} $=$ {10, 10, 10}.

In contrast, Figure 5 shows an example VDF that had a high $\tilde{\chi}{{}_{s}}^{2}$ for two components yet still a small $\delta\mathcal{R}$ $\sim$ 9.4%, i.e., this is still an example of a good fit despite the bad $\tilde{\chi}{{}_{s}}^{2}$ values for the core and beam/strahl fits. The fit results are as follows:

• •

$n{{}_{e\{c,h,b\}}}$ $=$ {4.41, 0.57, 0.32} $cm^{-3}$;

• •

$V{{}_{Te\{c,h,b\},\parallel}}$ $=$ {3882.6, 2624.5, 4574.5} $km\ s^{-1}$;

• •

$V{{}_{Te\{c,h,b\},\perp}}$ $=$ {2728.2, 2986.3, 2387.6} $km\ s^{-1}$;

• •

$V{{}_{oe\{c,h,b\},\parallel}}$ $=$ {-0.00, -594.9, $+$2000.0} $km\ s^{-1}$;

• •

$V{{}_{oe\{c,h,b\},\perp}}$ $=$ {-0.00, -0.00, -0.00} $km\ s^{-1}$;

• •

$\{p{{}_{ec}},q{{}_{ec}},\kappa{{}_{eh}},\kappa{{}_{eb}}\}$ $=$ {4.00, 2.00, 2.27, 4.61}, where $p{{}_{ec}}$($q{{}_{ec}}$) is the parallel(perpendicular) self-similar exponent and $\kappa{{}_{es}}$ is the kappa value;

• •

$\tilde{\chi}{{}_{e\{c,h,b\}}}^{2}$ $=$ {28.5, 0.55, 14.4};

• •

$\tilde{\chi}{{}_{tot}}^{2}$ $=$ 14.40; and

• •

Fit Flag {c,h,b} $=$ {4, 6, 5}.

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 halo999That is, the weights for the halo and beam/strahl fits are modified to force the software to examine only one-side of the velocity distribution at a time. The weights also remove elements from the core fit to avoid including the core in the fit.. Thus, the resultant $f^{\left(m\right)}$ can exceed $f^{\left(0\right)}$ 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 factors101010For instance, below $\sim$1000 km/s in Figure 5 the magnitude of $f^{\left(m\right)}/f^{\left(0\right)}$ stays below $\sim$1.7 and exceed 2.0 on the anti-parallel side above $\sim$10,000 km/s. The latter was not flagged by the software because it resulted from the beam/strahl fit and that is only fit to the parallel side for this VDF.. For most events, the threshold is set between $\sim$2–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 $\phi{{}_{sc}}$ yields fit results with $n{{}_{e}}$ $\sim$ 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 $\tilde{\chi}{{}_{s}}^{2}$ for the core component and moderate for beam/strahl but a small $\delta\mathcal{R}$ $\sim$ 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 upon 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 $\tilde{\chi}{{}_{tot}}^{2}$ is larger than a few. The fit results are as follows:

• •

$n{{}_{e\{c,h,b\}}}$ $=$ {3.37, 0.03, 0.14} $cm^{-3}$;

• •

$V{{}_{Te\{c,h,b\},\parallel}}$ $=$ {2609.8, 5293.2, 4686.9} $km\ s^{-1}$;

• •

$V{{}_{Te\{c,h,b\},\perp}}$ $=$ {2286.9, 5494.9, 2516.2} $km\ s^{-1}$;

• •

$V{{}_{oe\{c,h,b\},\parallel}}$ $=$ {-0.00, -222.8, $+$3273.0} $km\ s^{-1}$;

• •

$V{{}_{oe\{c,h,b\},\perp}}$ $=$ {-0.00, -0.00, -0.00} $km\ s^{-1}$;

• •

$\{s{{}_{ec}},\kappa{{}_{eh}},\kappa{{}_{eb}}\}$ $=$ {2.00, 3.83, 3.53};

• •

$\tilde{\chi}{{}_{e\{c,h,b\}}}^{2}$ $=$ {17.84, 0.17, 5.14 };

• •

$\tilde{\chi}{{}_{tot}}^{2}$ $=$ 13.17; and

• •

Fit Flag {c,h,b} $=$ {6, 6, 0}.

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 the combination of $\delta\mathcal{R}$ with $\tilde{\chi}{{}_{s}}^{2}$ and $\tilde{\chi}{{}_{tot}}^{2}$ are consistently more reliable quantities used in combination for defining the quality of the fit than using $\tilde{\chi}{{}_{s}}^{2}$ alone. The value is also used as a proxy for the uncertainty of any given fit parameter, e.g., $\delta n{{}_{es}}$ $=$ $\pm\ \delta\mathcal{R}\cdot 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 $\tilde{\chi}{{}_{s}}^{2}$ and $\delta\mathcal{R}$ 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 $\delta\mathcal{R}$ and all $\tilde{\chi}{{}_{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 $\sim$107 counts/second (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 $\sim 5\times 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}}$ $\sim$ $n{{}_{p}}$ $+$ 2$n{{}_{\alpha}}$ for nearly all intervals. Statistically, the difference between the fit result for $n{{}_{e}}$ $=$ $n{{}_{ec}}$ $+$ $n{{}_{eh}}$ $+$ $n{{}_{eb}}$ and $n{{}_{p}}$ $+$ 2$n{{}_{\alpha}}$ are within expectations. The median(lower quartile)[upper quartile] values are 10.3%(4.9%)[19.0%], which is consistent with our $\delta\mathcal{R}$ statistics.

Finally, the total electron current, $j{{}_{e,tot}}$ $=$ $\sum_{s}\ n{{}_{es}}\ v{{}_{os,\parallel}}$, 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 $\sim$22 $km/s\ cm^{-3}$, $\sim$0 $km/s\ cm^{-3}$, $\sim$-214 $km/s\ cm^{-3}$, and $\sim$351 $km/s\ cm^{-3}$, consistent with previously published work on this dataset (e.g., Bale et al., 2013; Pulupa et al., 2014a) and consistent with work in progress [Salem et al., 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 are $\sim$0.17%, $\sim$10*-8*%, $\sim$-0.95%, and $\sim$1.3%, respectively. Thus, the values are all small compared to unity. Quantitatively, $\sim$97.5% of the $j{{}_{e,tot}}$/($n{{}_{e}}$ $V{{}_{Tec,tot}}$) values satisfy $\lesssim$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 includes 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 $\sim$80.5% that satisfied 2.0 $\leq$ $s{{}_{ec}}$ $\leq$ 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., 1979b, a). 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., Bale et al., 2013; Salem et al., 2003).

3.4 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($\sim$4%) were modeled as bi-kappa VDFs, 12,095($\sim$80%) were modeled as symmetric bi-self-similar VDFs, and 2581($\sim$17%) were modeled as asymmetric bi-self-similar VDFs. All core bi-kappa 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 is given in Section 3 and Appendix C. Of those 15,210 that progressed to fit analysis stable solutions were found for 14,847($\sim$98%) $f^{\left(core\right)}$, 13,871($\sim$91%) $f^{\left(halo\right)}$, and 9567($\sim$63%) $f^{\left(beam\right)}$.

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. The statistical results of the densities are summarized below in the following form lower quartile–upper quartile(Mean)[Median]

• All

• –

$n{{}_{ec}}$ $\sim$ 6.44–19.5(13.7)[11.3] $cm^{-3}$;

• –

$n{{}_{eh}}$ $\sim$ 0.21–0.63(0.52)[0.36] $cm^{-3}$;

• –

$n{{}_{eb}}$ $\sim$ 0.09–0.27(0.21)[0.16] $cm^{-3}$;

• Upstream

• –

$n{{}_{ec}}$ $\sim$ 4.06–12.5(8.90)[8.09] $cm^{-3}$;

• –

$n{{}_{eh}}$ $\sim$ 0.17–0.49(0.42)[0.27] $cm^{-3}$;

• –

$n{{}_{eb}}$ $\sim$ 0.09–0.26(0.22)[0.16] $cm^{-3}$;

• Downstream

• –

$n{{}_{ec}}$ $\sim$ 8.44–24.2(17.3)[16.6] $cm^{-3}$;

• –

$n{{}_{eh}}$ $\sim$ 0.26–0.70(0.59)[0.44] $cm^{-3}$;

• –

$n{{}_{eb}}$ $\sim$ 0.09–0.28(0.21)[0.17] $cm^{-3}$.

which are consistent with previous results near 1 AU (e.g., Feldman et al., 1975, 1979a, 1983a; Maksimovic et al., 1997; Nieves-Chinchilla & Viñas, 2008; Phillips et al., 1989a, b; Pierrard et al., 2016; Salem et al., 2001; Skoug et al., 2000; Štverák et al., 2009). The full statistical results and associated histograms are presented in Paper II.

The statistical results of the quality analysis are listed below in the following form lower quartile–upper quartile(mean)[median]

• All

• –

$\delta\mathcal{R}$ $\sim$ 5.4%–15.0%(11.5%)[9.1%];

• –

$\tilde{\chi}{{}_{c}}^{2}$ $\sim$ 0.89–4.28(6.67)[1.95];

• –

$\tilde{\chi}{{}_{h}}^{2}$ $\sim$ 0.41–1.61(2.17)[0.72];

• –

$\tilde{\chi}{{}_{b}}^{2}$ $\sim$ 0.37–1.31(1.56)[0.66];

• –

$\tilde{\chi}{{}_{tot}}^{2}$ $\sim$ 2.82–9.40(694)[4.90];

• Upstream

• –

$\delta\mathcal{R}$ $\sim$ 7.4%–16.7%(13.3%)[11.3%];

• –

$\tilde{\chi}{{}_{c}}^{2}$ $\sim$ 0.62–2.43(2.05)[1.38];

• –

$\tilde{\chi}{{}_{h}}^{2}$ $\sim$ 0.31–1.01(2.03)[0.50];

• –

$\tilde{\chi}{{}_{b}}^{2}$ $\sim$ 0.32–0.98(1.24)[0.55];

• –

$\tilde{\chi}{{}_{tot}}^{2}$ $\sim$ 2.08–5.51(1557)[3.17];

• Downstream

• –

$\delta\mathcal{R}$ $\sim$ 6.3%–15.7%(12.1%)[10.0%];

• –

$\tilde{\chi}{{}_{c}}^{2}$ $\sim$ 1.22–8.55(10.1)[2.67];

• –

$\tilde{\chi}{{}_{h}}^{2}$ $\sim$ 0.52–1.95(2.26)[0.91];

• –

$\tilde{\chi}{{}_{b}}^{2}$ $\sim$ 0.41–1.55(1.78)[0.79];

• –

$\tilde{\chi}{{}_{tot}}^{2}$ $\sim$ 3.81–12.8(54.1)[6.98].

The purpose of listing these statistics is to provide a range of typical or expected $\tilde{\chi}{{}_{s}}^{2}$ and $\delta\mathcal{R}$ values for reference when determining the quality of any given fit. Note that the statistics for $\delta\mathcal{R}$ 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.

4 Exponents and Drifts

Table 2 shows the one-variable statistics for the exponents from the model fits of the electron VDFs are introduced and discussed, for the core ($s$ $=$ $c$), halo ($s$ $=$ $h$), and beam/strahl ($s$ $=$ $b$). The VDFs, modeled as bi-kappa ($\kappa{{}_{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, quasi-perpendicular 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 therefore the converse had no results to examine..

For the VDFs fit to a bi-kappa, the core values typically lie between $\sim$5–10 while the halo and beam/strahl lie between $\sim$3.5–5.4 and $\sim$3.4–5.2, respectively. Only the core was fit to the bi-self-similar functions and nearly all symmetric exponents are between $\sim$2.00–2.04 while most of the asymmetric parallel and perpendicular exponents lie $\sim$2.2–4.0 and $\sim$2.0–2.5, respectively.

\startlongtable

The $\kappa{{}_{eh}}$ and $\kappa{{}_{eb}}$ values are consistent with previous solar wind observations near 1 AU (e.g., Horaites et al., 2018; Lazar et al., 2017; Maksimovic et al., 1997, 2005; Pierrard et al., 2016; Štverák et al., 2009; Tao et al., 2016a, b). The $\kappa{{}_{ec}}$ values are also consistent with previous solar wind observations (e.g., Broiles et al., 2016; Nieves-Chinchilla & Viñas, 2008).

There are several interesting things to note from Table 2. The mean, median, and lower/upper quartile values for $\kappa{{}_{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 upon the upstream core electron distribution profiles. One possible physical interpretation would be that the sound speed depends upon the polytropic index for each species, i.e., the equation of state assumed for the system. A bi-kappa core VDF could effect 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 $\kappa{{}_{ec}}$ associated with higher Mach number shocks may imply that lower energy particles have entered the tails thus increasing the exponent111111Recall that $\kappa{{}_{ec}}$ values only exist for upstream VDF fits, so the dependence on Mach number is not about thermalization..

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). Such waves are driven by relative electron-ion drifts (i.e., currents) and are observed near both quasi-parallel and quasi-perpendicular shocks (e.g., Breneman et al., 2013; Fuselier & Gurnett, 1984; Wilson III et al., 2007, 2010, 2012, 2014a, 2014b) but their amplitudes increase with increasing shock strength (e.g., Wilson III 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 quasi-perpendicular 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,\parallel}}$. This would depend upon the effective inelasticity of the wave-particle interactions, where larger inelasticity increases $p{{}_{ec}}$ and smaller increases $T{{}_{ec,\parallel}}$ (e.g., Dum et al., 1974; Dum, 1975; Goldman, 1984; Horton et al., 1976; Horton & Choi, 1979; Jain & Sharma, 1979). 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 flattop electron distributions (i.e., $p{{}_{ec}}$ $\rightarrow$ $\geq$4 and $q{{}_{ec}}$ $\rightarrow$ $\sim$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., 1983b; Hull et al., 1998) through a process called maximal filling (e.g., Morse, 1965). However, similar to the predictions for wave-driven flattops this theory should generate stronger flattops (i.e., larger values of $p{{}_{ec}}$) for stronger quasi-perpendicular shocks, which we do not observe. Thus, the evolution of the electron VDFs do 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 $\kappa{{}_{eh}}$ values under different conditions. When the values of $\kappa{{}_{eh}}$ are larger(smaller), that implies a less(more) energized halo, i.e., softer(harder) spectra. One can see that $\kappa{{}_{eh}}$ is larger downstream than upstream and near high 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., Caprioli & Spitkovsky, 2014; Malkov & Drury, 2001; Park et al., 2015; Treumann, 2009; Trotta & Burgess, 2019). In contrast, $\kappa{{}_{eh}}$ is slightly smaller ($\sim$10%) near quasi-parallel than quasi-perpendicular shocks, which implies more energized halo electrons. Although quasi-parallel shocks are predicted (e.g., Caprioli & Spitkovsky, 2014; Malkov & Drury, 2001) and observed (e.g., Wilson III et al., 2016) 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 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 $\kappa{{}_{eh}}$ and $\kappa{{}_{eb}}$ are smaller in the upstream than downstream. The time-evolution of these kappa values will be examined in more detail in Paper III.

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 $\kappa{{}_{eb}}$ values show a similar behavior between upstream and downstream and shock geometry as $\kappa{{}_{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 $\kappa{{}_{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 Table 2, i.e., all VDF solutions. In many of the panels there are isolated, dominant peaks, nearly all of which result from constraints imposed for specific events, not necessarily an underlying physical reason. For instance, the peaks for $p{{}_{ec}}$ $=$ 3 and $=$ 4 in panel c are for strong shocks exhibiting flattop VDFs in the downstream where the fit routines were not finding stable solutions without imposing constraints on both the exponents and the minimum number density for the core distribution.

One can see that, as discussed previously, the core parallel drift speeds (violet line, panel d) tend to fall below $\sim$100 km/s, 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,\parallel}}$ $\leq$ 1 km/s and $V{{}_{oeh,\parallel}}$ $\leq$ 1 km/s are 8735($\sim$59%) and 7311($\sim$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 non-zero perpendicular drifts almost certainly result from inaccuracies in the calculation of the solar wind rest frame and a dipole correction to $\phi{{}_{sc}}$ not included in the present analysis (e.g., Pulupa et al., 2014a) (see Appendix B for more details).

The magnitudes of $V{{}_{oeh,\perp}}$ and $V{{}_{oeb,\perp}}$ never deviated from zero121212This was an explicit constraint imposed on all fits but would also have resulted largely from the initial guess that both $V{{}_{oeh,\perp}}$ and $V{{}_{oeb,\perp}}$ 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.. The magnitudes of $V{{}_{oeh,\parallel}}$ range from $\sim$0–8860 km/s with a lower to upper quartile range of $\sim$0–850 km/s and a mean(median) of $\sim$580 km/s($\sim$0.1 km/s). The magnitudes of $V{{}_{oeb,\parallel}}$ range from $\sim$1000–9330 km/s with a lower to upper quartile range of $\sim$1750–3090 km/s and a mean(median) of $\sim$2580 km/s($\sim$2480 km/s). As previously discussed, the lower bound for $V{{}_{oeb,\parallel}}$ was imposed on the basis of physical arguments while the magnitude of $V{{}_{oeh,\parallel}}$ was allowed to go to zero. If only magnitudes satisfying $V{{}_{oes,\parallel}}$ $>$ 1 km/s are considered, the mean(median) and lower to upper quartile ranges are $\sim$42 km/s($\sim$30 km/s) and $\sim$14–52 km/s for $V{{}_{oec,\parallel}}$ and $\sim$1227 km/s($\sim$903 km/s) and $\sim$362–1695 km/s for $V{{}_{oeh,\parallel}}$.

5 Discussion

A total of 15,314 electron VDFs were observed by the Wind spacecraft within $\pm$2 hours of 52 IP shocks of which 15,210 had a stable solution for at least one component. Stable model function parameters were found for 14,847($\sim$98%) core fits, 13,871($\sim$91%) halo fits, and 9567($\sim$63%) beam/strahl fits. The fit parameters are consistent with previous studies and will be discussed in detail in the following two parts of this study. Of the 15,210 VDFs examined herein, the core was modeled as a bi-kappa for 534($\sim$4%) VDFs, as a symmetric bi-self-similar for 12,095($\sim$80%) VDFs, and as an asymmetric bi-self-similar for 2581($\sim$17%) VDFs. This is the first statistical study to find that the core electron distribution is better fit to a self-similar velocity distribution function than a Maxwellian under all conditions.

The exponents are summarized below in the following form lower quartile–upper quartile(Mean)[Median]

• All

• –

$s{{}_{ec}}$ $\sim$ 2.00–2.04(2.03)[2.00];

• –

$p{{}_{ec}}$ $\sim$ 2.20–4.00(3.09)[3.00];

• –

$q{{}_{ec}}$ $\sim$ 2.00–2.46(2.24)[2.00];

• –

$\kappa{{}_{ec}}$ $\sim$ 5.40–10.2(9.15)[7.92];

• –

$\kappa{{}_{eh}}$ $\sim$ 3.58–5.34(4.62)[4.38];

• –

$\kappa{{}_{eb}}$ $\sim$ 3.40–5.16(4.57)[4.17];

• Upstream

• –

$s{{}_{ec}}$ $\sim$ 2.00–2.03(2.01)[2.00];

• –

$p{{}_{ec}}$ $\sim$ N/A;

• –

$q{{}_{ec}}$ $\sim$ N/A;

• –

$\kappa{{}_{ec}}$ $\sim$ 5.40–10.2(9.15)[7.92];

• –

$\kappa{{}_{eh}}$ $\sim$ 3.25–4.83(4.16)[4.10];

• –

$\kappa{{}_{eb}}$ $\sim$ 3.25–4.70(4.22)[3.81];

• Downstream

• –

$s{{}_{ec}}$ $\sim$ 2.00–2.06(2.05)[2.01];

• –

$p{{}_{ec}}$ $\sim$ 2.20–4.00(3.09)[3.00];

• –

$q{{}_{ec}}$ $\sim$ 2.00–2.46(2.24)[2.00];

• –

$\kappa{{}_{ec}}$ $\sim$ N/A;

• –

$\kappa{{}_{eh}}$ $\sim$ 3.80–5.70(4.94)[4.62];

• –

$\kappa{{}_{eb}}$ $\sim$ 3.61–5.44(4.82)[4.45];

Overall the $\kappa{{}_{eh}}$ and $\kappa{{}_{eb}}$ values are consistent with previous solar wind observations near 1 AU (e.g., Horaites et al., 2018; Lazar et al., 2017; Pierrard et al., 2016; Štverák et al., 2009). The $\kappa{{}_{ec}}$ values are also consistent with previous solar wind observations (e.g., Broiles et al., 2016; Nieves-Chinchilla & Viñas, 2008). The values for $s{{}_{ec}}$, $p{{}_{ec}}$, and $q{{}_{ec}}$ are consistent with previous results as well (e.g., Feldman et al., 1983b, a).

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 either for 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; Goldman, 1984; Horton et al., 1976; Horton & Choi, 1979; Jain & Sharma, 1979). 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 collision131313That is, the particle kinetic energy may not be preserved through the interaction even if the wave kinetic energy is conserved.. 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 two, i.e., the deviation from a Maxwellian is a measure of inelasticity in the particles interactions with waves and/or turbulence141414It 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].. Further, as previously discussed, $\sim$80.5% of the core VDFs modeled with a symmetric bi-self-similar function had exponents satisfying 2.0 $\leq$ $s{{}_{ec}}$ $\leq$ 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., 1979b, a) and work that found evidence for collisional effects in the core distribution (e.g., Bale et al., 2013; Salem et al., 2003).

The $\kappa{{}_{ec}}$ seem to correlate with $\langle M{{}_{f}}\rangle{{}_{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 upon $\langle M{{}_{f}}\rangle{{}_{up}}$.

Somewhat surprisingly the values of $\kappa{{}_{eh}}$ are larger downstream than upstream and they increase with increasing $\langle M{{}_{f}}\rangle{{}_{up}}$. That is, the halo spectra are softer downstream and near strong shocks. Quasi-parallel shocks, however, correlate with smaller $\kappa{{}_{eh}}$, i.e., harder halo spectra. Generally, quasi-parallel shocks are predicted to be more efficient particle accelerators for suprathermal ions and very energetic electrons151515Suprathermal is defined here for ions in the several to 10s of keV energy range while the electrons are many 10s to 100s of keV for typical 1 AU solar wind collisionless shocks. (e.g., Caprioli & Spitkovsky, 2014) but electrons in the halo energy range are predicted to be energized the most efficiently at shocks satisfying $\theta{{}_{Bn}}$ $>$ 80*∘* (e.g., Park et al., 2013).

Unlike the halo, $\kappa{{}_{eb}}$ are smaller near high Mach number shocks than near low Mach number shocks. The difference is likely a two-fold consequence of the combined effects from shock-accelerated foreshock electrons and the method used to fit the distributions. That is, the beam/strahl component is always fit to the anti-sunward, field-aligned side of the VDF while the halo to the opposite. For nearly all IP shocks at 1 AU, the shock normal is anti-sunward 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 $\langle\kappa{{}_{eh}}\rangle{{}_{dn}}/\langle\kappa{{}_{eh}}\rangle{{}_{up}}$ and $\langle\kappa{{}_{eb}}\rangle{{}_{dn}}/\langle\kappa{{}_{eb}}\rangle{{}_{up}}$ increase with increasing $\langle M{{}_{f}}\rangle{{}_{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 velocity distribution function than a bi-Maxwellian under all conditions. This is an important result for kinetic theory and solar wind evolution. This work aslo 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-galaxy-cluster medium and they provide a statistical baseline of electron parameters near collisionless shocks for the recent Parker Solar Probe and upcoming Solar Orbiter missions.

The authors thank A.F.- Viñas and D.A. Roberts for useful discussions of basic plasma physics and C. Markwardt for helpful feedback on the usage nuances of his MPFIT software. The work was supported by the International Space Science Institute’s (ISSI) International Teams programme. L.B.W. was partially supported by Wind MO&DA grants and a Heliophysics Innovation Fund (HIF) grant. L.-J.C. and S.W. were partially supported by the MMS mission in addition to NASA grants 80NSSC18K1369 and 80NSSC17K0012, NSF grants AGS-1619584 and AGS-1552142, and DOE grant DESC0016278. D.L.T. was partially supported by NASA grant NNX16AQ50G. M.L.S. was partially supported by grants NNX14AT26G and NNX13AI75G. J.C.K. was partially supported by NASA grants NNX14AR78G and 80NSSC18K0986. D.C. was partially supported by grants NNX17AG30G, GO8-19110A, 80NSSC18K1726, 80NSSC18K1218, and NSF grant 1714658. S.J.S. was partially supported by the MMS/FIELDS investigation. C.S.S. was partially supported by NASA grant NNX16AI59G and NSF SHINE grant 1622498. S.D.B. and C.S.S. were partially supported by NASA grant NNX16AP95G. M.P.P. and K.A.G. were supported by Parker Solar Probe instrument funds.

Appendix A Definitions and Notation

In this appendix we define the symbols and notation used throughout. In the following, all direction-dependent parameters we use the subscript $j$ to represent the direction where $j$ $=$ $tot$ for the entire distribution, $j$ $=$ $\parallel$ for the the parallel direction, and $j$ $=$ $\perp$ for the perpendicular direction. Note that parallel and perpendicular are with respect to the quasi-static magnetic field vector, $\mathbf{B}{{}_{o}}$ [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, and $s$ $=$ $eff$ for the effective, and $s$ $=$ $e$ for the total/entire population. Below are the symbol/parameters definitions:

• one-variable statistics

• –

$X{{}_{min}}$ $\equiv$ minimum

• –

$X{{}_{max}}$ $\equiv$ maximum

• –

$\bar{X}$ $\equiv$ mean

• –

$\tilde{X}$ $\equiv$ median

• –

$X{{}_{25\%}}$ $\equiv$ lower quartile

• –

$X{{}_{75\%}}$ $\equiv$ upper quartile

• fundamental parameters

• –

$\varepsilon{{}_{o}}$ $\equiv$ permittivity of free space

• –

$\mu{{}_{o}}$ $\equiv$ permeability of free space

• –

$c$ $\equiv$ speed of light in vacuum [$km\ s^{-1}$] $=$ $\left(\varepsilon{{}_{o}}\ \mu{{}_{o}}\right)^{-1/2}$

• –

$k{{}_{B}}$ $\equiv$ the Boltzmann constant [$J\ K^{-1}$]

• –

$e$ $\equiv$ the fundamental charge [$C$]

• plasma parameters

• –

$n{{}_{s}}$ $\equiv$ the number density [$cm^{-3}$] of species $s$

• –

$m{{}_{s}}$ $\equiv$ the mass [$kg$] of species $s$

• –

$Z{{}_{s}}$ $\equiv$ the charge state of species $s$

• –

$q{{}_{s}}$ $\equiv$ the charge [$C$] of species $s$ $=$ $Z{{}_{s}}\ e$

• –

$T{{}_{s,j}}$ $\equiv$ the scalar temperature [$eV$] of the jth component of species $s$

• –

$\left(T{{}_{s^{\prime}}}/T{{}_{s}}\right){{}_{j}}$ $\equiv$ the temperature ratio [N/A] of species $s$ and $s^{\prime}$ of the jth component

• –

$\left(T{{}_{\perp}}/T{{}_{\parallel}}\right){{}_{s}}$ $\equiv$ the temperature anisotropy [N/A] of species $s$

• –

$V{{}_{Ts,j}}$ $\equiv$ the most probable thermal speed [$km\ s^{-1}$] of a one-dimensional velocity distribution (see Equation A1c)

• –

$\mathbf{v}{{}_{os}}$ $\equiv$ the drift velocity [$km\ s^{-1}$] of species $s$ in the plasma bulk flow rest frame

• –

$C{{}_{s}}$ $\equiv$ the sound or ion-acoustic sound speed [$km\ s^{-1}$] (see Supplemental Material for definitions)

• –

$V{{}_{A}}$ $\equiv$ the Alfvén speed [$km\ s^{-1}$] (see Supplemental Material for definitions)

• –

$V{{}_{f}}$ $\equiv$ the fast mode speed [$km\ s^{-1}$] (see Supplemental Material for definitions)

• –

$\Omega{{}_{cs}}$ $\equiv$ the angular cyclotron frequency [$rad\ s^{-1}$] (see Equation A1d)

• –

$\omega{{}_{ps}}$ $\equiv$ the angular plasma frequency [$rad\ s^{-1}$] (see Equation A1e)

• –

$\lambda{{}_{De}}$ $\equiv$ the electron Debye length [$m$] (see Equation A1f)

• –

$\rho{{}_{cs}}$ $\equiv$ the thermal gyroradius [$km$] (see Equation A1g)

• –

$\lambda{{}_{s}}$ $\equiv$ the inertial length [$km$] (see Equation A1h)

• –

$\beta{{}_{s,j}}$ $\equiv$ the plasma beta [N/A] of the jth component of species $s$ (see Equations A1i and A1j)

• –

$\phi{{}_{sc}}$ $\equiv$ the scalar, quasi-static spacecraft potential [eV] (e.g., Pulupa et al., 2014a; Scime et al., 1994b)

• –

$E{{}_{min}}$ $\equiv$ the minimum energy bin midpoint value [eV] of an electrostatic analyzer (e.g., see Appendices in Wilson III et al., 2017, 2018)

The variables that rely upon multiple parameters are given in the following equations:

[TABLE]

For the macroscopic shock parameters, the values are averaged over asymptotic regions away from the shock transition region.

• shock parameters

• –

subscripts $up$ and $dn$ $\equiv$ 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)

• –

$\langle Q\rangle{{}_{j}}$ $\equiv$ the average of parameter $Q$ over the $j^{th}$ shock region, where $j$ $=$ $up$ or $dn$

• –

$\mathbf{n}{{}_{sh}}$ $\equiv$ the shock normal unit vector [N/A]

• –

$\theta{{}_{Bn}}$ $\equiv$ the shock normal angle [deg], defined as the acute reference angle between $\langle\mathbf{B}{{}_{o}}\rangle{{}_{up}}$ and $\mathbf{n}{{}_{sh}}$

• –

$\langle\lvert V{{}_{shn}}\rvert\rangle{{}_{j}}$ $\equiv$ the $j^{th}$ region average shock normal speed [$km\ s^{-1}$] in the spacecraft frame

• –

$\langle\lvert U{{}_{shn}}\rvert\rangle{{}_{j}}$ $\equiv$ 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)

• –

$\langle M{{}_{A}}\rangle{{}_{j}}$ $\equiv$ the $j^{th}$ region average Alfvénic Mach number [N/A] $=$ $\langle\lvert U{{}_{shn}}\rvert\rangle{{}_{j}}/\langle V{{}_{A}}\rangle{{}_{j}}$

• –

$\langle M{{}_{f}}\rangle{{}_{j}}$ $\equiv$ the $j^{th}$ region average fast mode Mach number [N/A]

• –

$M{{}_{cr}}$ $\equiv$ the first critical Mach number [N/A]

• –

$M{{}_{ww}}$ $\equiv$ the linear whistler (phase) Mach number [N/A]

• –

$M{{}_{gr}}$ $\equiv$ the linear whistler (group) Mach number [N/A]

• –

$M{{}_{nw}}$ $\equiv$ the nonlinear whistler Mach number

The critical Mach numbers are phenomenologically defined as follows: for $\langle M{{}_{f}}\rangle{{}_{up}}/M{{}_{cr}}$ $\geq$ 1 an ion sound wave could not phase stand within the shock ramp (e.g., Edmiston & Kennel, 1984; Kennel et al., 1985); for $\langle M{{}_{f}}\rangle{{}_{up}}/M{{}_{ww}}$ $\geq$ 1 a linear magnetosonic-whistler cannot phase stand upstream of the shock ramp (e.g., Krasnoselskikh et al., 2002); for $\langle M{{}_{f}}\rangle{{}_{up}}/M{{}_{gr}}$ $\geq$ 1 a linear magnetosonic-whistler cannot group stand upstream of the shock ramp; and for $\langle M{{}_{f}}\rangle{{}_{up}}/M{{}_{nw}}$ $\geq$ 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 deadtime161616The deadtime is the time period when the detector is unable to measure incident particles due 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. (e.g., Meeks & Siegel, 2008; Schecker et al., 1992), and an unknown spacecraft potential (e.g., Lavraud & Larson, 2016; Pulupa et al., 2014a; Scime et al., 1994b, a). Significant advances in understanding the response and calibration of electrostatic analyzers have been made in recent years with the development and launch of the Magnetospheric Multiscale (MMS) mission (e.g., Gershman et al., 2016, 2017; Pollock et al., 2016). However, the improvements resulted from an exhaustive ground calibration campaign that most other missions, including Wind, have not had. Further, the electronic deadtime171717The cycle rate or sample rate of this preamp is listed as 2 MHz but it is not constant. of the EESA Low preamp (i.e., AMPTEK A111) depends upon the pulse height distribution of the previous pulse [J.P. McFaddon, Personal Communication, July 18, 2011].

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 (e.g., see McFadden et al., 2008b, a, for further discussions of MCP degradation over time). The currently used calibrations are those from optical geometric factor corrections, on-ground calibrations, and in-flight calibrations [D. Larson, Personal Communication, July 18, 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 has been presented in numerous refereed publications including but not limited to Bale et al. (2013); Pulupa et al. (2014a, b); Salem et al. (2001, 2003); Wilson III et al. (2009, 2010, 2012, 2013a, 2013b, 2018). Updating the calibration tables is beyond the scope of this work but is actively being pursued [Salem et al., 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, $\phi{{}_{sc}}$. It does, however, 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}}$ $+$ 2$n{{}_{\alpha}}$. Ions are generally not significantly affected by $\phi{{}_{sc}}$ as they typically have $\sim$1 keV of bulk kinetic energy in the solar wind.

To estimate $\phi{{}_{sc}}$ an initial guess is determined numerically from the ion density. The value of $\phi{{}_{sc}}$ is then adjusted until $n{{}_{e}}$ $=$ $n{{}_{ec}}$ $+$ $n{{}_{eh}}$ $+$ $n{{}_{eb}}$ from the fits roughly equals181818Note 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. $n{{}_{p}}$ $+$ 2$n{{}_{\alpha}}$ and/or when photoelectrons disappear from the VDF plots191919When $\phi{{}_{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 $\phi{{}_{sc}}$.. Once a reliable estimate of $\phi{{}_{sc}}$ 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 $\phi{{}_{sc}}$ is shown in Figure 3.

Note that the values of $\phi{{}_{sc}}$ determined above should not be treated as the absolute or correct spacecraft potential values. The reason being that the detector efficiency and gain calibrations suffer from the issues discussed above. (e.g., Bordoni, 1971; Goruganthu & Wilson, 1984). Therefore, the $\phi{{}_{sc}}$ values are proxies for the spacecraft potential that comprise a complicated nonlinear convolution of the real spacecraft potential and the detector deadtime and efficiency. Despite this uncertainty, the $\phi{{}_{sc}}$ values estimated herein are consistent with those in previously published work on the same dataset within the same time span (e.g., Bale et al., 2013; Pulupa et al., 2014a). Further, the consistency checks discussed in Section 3.3 provide further validation of the fit results.

Table 3 provides the one-variable statistics of the $\phi{{}_{sc}}$ values for all VDFs, 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 $\phi{{}_{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 $\phi{{}_{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^{C\ X}+D$, where $Y$ $=$ $\left(\phi{{}_{sc}}+E{{}_{min}}\right)/5$ [eV] and $X$ $=$ $n{{}_{i}}$ [$cm^{-3}$]. The fit parameters producing the cyan dashed line are $A$ $=$ 2.272$\pm$0.013 [$cm^{+3\ B}$], $B$ $=$ -0.431$\pm$0.019 [N/A], $C$ $=$ 0.00115$\pm$0.00155 [$cm^{+3}$], and $D$ $=$ 2.0$\pm$0.0 [eV], with a reduced chi-squared value of $\tilde{\chi}^{2}$ $\sim$ 0.144.

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 then fit to a function of the spacecraft potential for the number density, i.e., $n{{}_{i}}$ $=$ $n{{}_{i}}\left(\phi{{}_{sc}}\right)$ (e.g., Scudder et al., 2000). As previously stated, Wind cannot actively measure $\phi{{}_{sc}}$ and the values shown in Figure 9 are really a proxy due to the uncertain values for the deadtime and efficiency for each detector anode. The purpose of the above approach is to find a semi-analytical expression for $\phi{{}_{sc}}$ that only depends upon $n{{}_{i}}$ (or $n{{}_{e}}$) as an initial estimate. The unexpected result here is that the trend depends upon $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 dataset 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 $\sim$1.5% of all the VDFs examined in our study satisfied $\phi{{}_{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 Levenberg-Marquardt Algorithm (LMA) (Moré, 1978). The generalized LMA software, called MPFIT (Markwardt, 2009), requires at minimum the following inputs when fitting to a two-dimensional array of data:

FUNC

A scalar [string] defining the model function routine file name;

X(Y)

N(M)-element [numeric] array defining the first(second) dimension coordinate abscissa values;

Z

NxM-element [numeric] array defining the dependent data associated with X and Y abscissa values;

ERR

NxM-element [numeric] array defining the error associated with each element of Z; and

PARAM

K-element [numeric] array defining the initial guesses for the fit parameters supplied to the model function routine FUNC.

The error array will be ignored if the user supplies an array of weights, $\mathcal{W}$. The details of the use of the software and documentation are provided by the author at

https://www.physics.wisc.edu/c̃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:

• •

The raw VDF data, $f^{\left(0r\right)}$, is 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^{\left(0cr\right)}$, i.e., Poisson statistics are assumed.

• •

A unit conversion is applied to change to units of phase space density [i.e., cm*-3* km*-3* s*+3*] then the energies are adjusted to account for the spacecraft potential (e.g., Salem et al., 2001; Wilson III et al., 2014a, 2016) (details are discussed in Appendix B) giving $f^{\left(0sc\right)}$ and $f^{\left(0csc\right)}$.

• •

Then $f^{\left(0sc\right)}$ and $f^{\left(0csc\right)}$ are transformed into the ion bulk flow rest frame (e.g., Compton & Getting, 1935; Ipavich, 1974) following the methods described in Wilson III et al. (2016) (see also the associated Supplemental Material) 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, $\mathbf{B}{{}_{o}}$, along the horizontal and the transverse component of the ion bulk flow velocity, $\mathbf{V}{{}_{i}}$, i.e., $\left(\mathbf{B}{{}_{o}}\times\mathbf{V}{{}_{i}}\right)\times\mathbf{B}{{}_{o}}$. The result is a two-dimensional gyrotropic VDF, $f^{\left(0\right)}$, and the associated Poisson errors/uncertainties, $f^{\left(0c\right)}$, both as functions of the parallel, $V{{}_{\parallel}}$, and perpendicular, $V{{}_{\perp}}$, velocity with respect to $\mathbf{B}{{}_{o}}$.

• •

Numerous weighting schemes were tried and the best results (for Wind/3DP) were achieved by defining $\mathcal{W}$ $=$ $\left(f^{\left(0c\right)}\right)^{-2}$ for the weights202020Several approaches were tried for the $\mathcal{W}$ 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..

• •

Every $f^{\left(0\right)}$ is fit to the sum of three model functions in two-dimensions212121That is, the data are not fit to two one-dimensional cuts of a two-dimensional VDF separately but rather both dimensions are fit simultaneously. for the core, halo, and beam/strahl components. Again, the components can be fit separately because the solar wind is a non-equilibrium, weakly collisional, kinetic gas222222It 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.. The allowed model functions (defined in Section 3.1) and are bi-Maxwellian (e.g., Kasper et al., 2006), bi-kappa (e.g., Livadiotis, 2015; Mace & Sydora, 2010; Vasyliunas, 1968), symmetric bi-self-similar (e.g., Dum et al., 1974; Dum, 1975), and asymmetric bi-self-similar (defined in Section 3.1).

• –

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 due to numerical instabilities232323There 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 due to these thresholds being satisfied too early in the iteration process. (e.g., Liavas & Regalia, 1999), discussed in Appendix D.

• –

Thus, the core fit, $f^{\left(core\right)}$, is performed first and then the model result subtracted from the data to yield the first residual, $f^{\left(1\right)}$.

• –

The halo fit, $f^{\left(halo\right)}$, is next but only to the side of $f^{\left(1\right)}$ opposite to that expected for the strahl/beam, where the latter is defined as the anti-sunward direction along $\mathbf{B}{{}_{o}}$. The entire two-dimensional halo fit is then subtracted from $f^{\left(1\right)}$ to yield the second residual, $f^{\left(2\right)}$, i.e., both sides are subtracted but only one side is used for the fit.

• –

The beam/strahl fit, $f^{\left(beam\right)}$, is last and fit to only the side of $f^{\left(2\right)}$ that is in the anti-sunward direction along $\mathbf{B}{{}_{o}}$.

• •

Not all VDFs will have fit results for all three components. In fact, $f^{\left(beam\right)}$ is often not found either because $f^{\left(halo\right)}$ left too few finite elements in $f^{\left(2\right)}$ or 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: PARAM[0] is the number density, $n{{}_{s}}$ [$cm^{-3}$]; PARAM[1] and PARAM[2] are the parallel and perpendicular thermal speeds, $V{{}_{Ts,j}}$ [$km\ s^{-1}$]; PARAM[3] and PARAM[4] are the parallel and perpendicular drift speeds, $V{{}_{os,j}}$ [$km\ s^{-1}$]; and PARAM[5] is the function exponent. The exponent input is ignored for the bi-Maxwellian routine as it is always 2.0 but can vary in the other routines. For the asymmetric bi-self-similar routine PARAM[4] is the parallel exponent and PARAM[5] is the perpendicular exponent (see Section 3.1 for functional form).

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).

Appendix D Numerical Instability

The Levenberg-Marquardt Algorithm (LMA) software works by minimizing the the chi-squared value given by:

[TABLE]

where $f{{}_{s}}^{\left(mod\right)}$ is the model fit function of component $s$ returned by the model function routine FUNC (see Section C), $\chi{{}_{s}}^{2}$ is the chi-squared value of the fit of component $s$, and the $i$ and $j$ subscripts correspond to the indices of the parallel and perpendicular velocity space coordinates, respectively.

A total reduced chi-squared, $\tilde{\chi}{{}_{tot}}^{2}$, value was also calculated for all VDFs analyzed herein. The difference in calculation is that the weights were not offset and the model function and distribution function are for the entire VDF, not the components. Further, unlike the components, the $\tilde{\chi}{{}_{tot}}^{2}$ values used all data points in $f^{\left(0\right)}$ and $\mathcal{W}$ even if they were excluded during the fit process242424Specific energy-angle bins were excluded for various physical reasons in some VDFs including, for instance, energy and/or pitch-angle range constraints to avoid “contamination” by other components as is done to examine the halo-only and beam/strahl-only parts of the VDF.. However, the $\tilde{\chi}{{}_{tot}}^{2}$ calculation excluded data below the nine-count level to avoid non-Gaussian weights in low-count values removed “spiky” solutions in the beam or halo fits defined by small $T{{}_{es,j}}$ and $\kappa{{}_{es}}$. That is, “spiky” solutions are defined as those satisfying $\left(\kappa{{}_{es}}\leq 3\right)$ $\wedge$ $\left(\left(T{{}_{es,\parallel}}\leq 11.8\right)\vee\left(T{{}_{es,\perp}}\leq 11.8\right)\right)$ for model fit parameters. As evidenced by Figures 4–6, the $\tilde{\chi}{{}_{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^{\left(0\right)}$ for electrons near 1 AU varies from $\sim$ 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 function252525i.e., the partial derivatives are with respect to the fit parameters, not the velocity coordinates. This is problematic when the magnitude of the input data and output model function are much much less than unity as it 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 be near unity. The solution was to multiply $\mathcal{W}$ by a constant offset to increase the contrast in the Jacobian components that are used to minimize $\chi^{2}$. A consequence of this approach is that the output $\chi^{2}$, $f^{\left(m\right)}$, and one-sigma error estimates of the fit parameters must be re-normalized 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 flattop distributions (e.g., Feldman et al., 1983b; Thomsen et al., 1987), modeled using the self-similar distributions (e.g., Dum et al., 1974; Dum, 1975; Goldman, 1984; Horton et al., 1976; Horton & Choi, 1979; Jain & Sharma, 1979) given by either Equation 3a or 4a. 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 Supplemental Material ASCII files described in Appendix F for list of constraints by shock).

Appendix E Numerical Method Comparisons

As stated in Appendix D, the standard approach to avoiding numerical instabilities due to the small magnitude of $f^{\left(0\right)}$ usually involves fitting to the logarithm of $f^{\left(0\right)}$ (e.g., Štverák et al., 2009). To illustrate the validity of the method used herein, an example VDF was chosen from a different study [Farrugia et al., in preparation] that examines a single shock-magnetic-cloud system.

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 Panels b and c, referred to as the test fit from hereon, 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}}$ 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 linear method) to the base-10 logarithm approach (referred to as 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:

• Test Fit (Panels b and c)

• –

$n{{}_{ec(h)[b]}}$ $\sim$ 1.407(0.054)[0.060] $cm^{-3}$;

• –

$V{{}_{Tec(h)[b],\parallel}}$ $\sim$ 2028.0(3621.7)[4183.7] $km\ s^{-1}$;

• –

$V{{}_{Tec(h)[b],\perp}}$ $\sim$ 1927.2(3486.6)[2833.1] $km\ s^{-1}$;

• –

$V{{}_{oec(h)[b],\parallel}}$ $\sim$ $+$50.4(0.0)[-3752.6] $km\ s^{-1}$;

• –

$V{{}_{oec(h)[b],\perp}}$ $\sim$ 0.0(0.0)[0.0] $km\ s^{-1}$;

• –

$s{{}_{ec}}$ $\sim$ 2.002;

• –

$\kappa{{}_{eh}}$ $\sim$ 1.908;

• –

$\kappa{{}_{eb}}$ $\sim$ 5.151;

• –

$\delta\mathcal{R}$ $\sim$ 14.1%;

• –

$\tilde{\chi}{{}_{c(h)[b]}}^{2}$ $\sim$ 4.52(1.82)[2.98];

• –

$\tilde{\chi}{{}_{tot}}^{2}$ $\sim$ 1.30;

• –

Fit Flag {c,h,b} $=$ {8, 8, 8}.

• Linear Method Fit (Panels d and e)

• –

$n{{}_{ec(h)[b]}}$ $\sim$ 1.122(0.051)[0.055] $cm^{-3}$;

• –

$V{{}_{Tec(h)[b],\parallel}}$ $\sim$ 2183.7(3694.7)[4154.0] $km\ s^{-1}$;

• –

$V{{}_{Tec(h)[b],\perp}}$ $\sim$ 1947.4(3557.0)[2863.1] $km\ s^{-1}$;

• –

$V{{}_{oec(h)[b],\parallel}}$ $\sim$ 0.0(0.0)[-3960.1] $km\ s^{-1}$;

• –

$V{{}_{oec(h)[b],\perp}}$ $\sim$ 0.0(0.0)[0.0] $km\ s^{-1}$;

• –

$s{{}_{ec}}$ $\sim$ 2.000;

• –

$\kappa{{}_{eh}}$ $\sim$ 1.901;

• –

$\kappa{{}_{eb}}$ $\sim$ 5.073; and

• –

$\delta\mathcal{R}$ $\sim$ 12.5%;

• –

$\tilde{\chi}{{}_{c(h)[b]}}^{2}$ $\sim$ 3.73(1.82)[2.98];

• –

$\tilde{\chi}{{}_{tot}}^{2}$ $\sim$ 1.01;

• –

Fit Flag {c,h,b} $=$ {8, 8, 8}.

• Log Method Fit (Panels f and g)

• –

$n{{}_{ec(h)[b]}}$ $\sim$ 1.086(0.089)[0.062] $cm^{-3}$;

• –

$V{{}_{Tec(h)[b],\parallel}}$ $\sim$ 3248.5(2938.9)[3762.2] $km\ s^{-1}$;

• –

$V{{}_{Tec(h)[b],\perp}}$ $\sim$ 2086.3(3043.6)[2652.4] $km\ s^{-1}$;

• –

$V{{}_{oec(h)[b],\parallel}}$ $\sim$ 0.0(0.0)[-4206.0] $km\ s^{-1}$;

• –

$V{{}_{oec(h)[b],\perp}}$ $\sim$ 0.0(0.0)[0.0] $km\ s^{-1}$;

• –

$s{{}_{ec}}$ $\sim$ 2.000;

• –

$\kappa{{}_{eh}}$ $\sim$ 1.852;

• –

$\kappa{{}_{eb}}$ $\sim$ 4.016;

• –

$\delta\mathcal{R}$ $\sim$ 20.0%;

• –

$\tilde{\chi}{{}_{c(h)[b]}}^{2}$ $\sim$ 82.7(2.51)[2.15];

• –

$\tilde{\chi}{{}_{tot}}^{2}$ $\sim$ 0.71;

• –

Fit Flag {c,h,b} $=$ {2, 6, 6}.

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 10, the initial guess values used for the core fit were $n{{}_{ec}}$ $\sim$ 2.0 $cm^{-3}$, $V{{}_{Tec,\parallel[\perp]}}$ $\sim$ 2297[2297] km/s (i.e., 15 eV temperatures), $V{{}_{oec,\parallel[\perp]}}$ $\sim$ $+$10.0[0.0] km/s, and $s{{}_{ec}}$ $\sim$ 2.0. If any of the parameters were perturbed by $\sim$20–30% away from these initial guesses, the log method would not initiate fit iterations due to diverging deviates and/or diverging model results, i.e., the software could not establish an initial Jacobian262626This is associated with a fit status code of -16 as reported in the fit results ASCII file discussed in Appendix F.. 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 non-initiation 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 one-sigma uncertainties of every fit parameter. These values are not reported because it was found they do not accurately or consistently reflect the quality of fit. For instance, the one-sigma uncertainties of $n{{}_{eh}}$ and $V{{}_{Teh,\parallel}}$ for the log method in the example VDF shown in Figure 10 (Panels f and g) are $\sim$19,988 km/s (i.e., $\sim$617% error) and $\sim$3.53 $cm^{-3}$ (i.e., $\sim$5163% error), respectively, even though $\tilde{\chi}{{}_{h}}^{2}$ $\sim$ 2.51. The one-sigma uncertainties for the same parameters but for the fit in Panels d and e are $\sim$110.1 km/s (i.e., $\sim$3.1% error) and $\sim$0.0047 $cm^{-3}$ (i.e., $\sim$8.5% error) and $\tilde{\chi}{{}_{h}}^{2}$ $\sim$ 1.82. That is, the reduced chi-squared values differ by only $\sim$39% but the one-sigma uncertainties by 100s to 1000s percent. The one-sigma 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 $\tilde{\chi}{{}_{s}}^{2}$) is not the same basis as that for the fit parameters. The output uncertainties thus contain nonlinear convolution of one-sigma uncertainties from potentially multiple fit parameters. The effect is analogous to electric field measurements from two antenna 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 one-sigma uncertainties were not used as errors for each parameter.

The one-sigma errors are also forced to zero in the software when the fit value reaches a user-defined 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 $\delta\mathcal{R}$ 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. 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 non-electron 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{{}_{\alpha}}$ measured by SWE [$cm^{-3}$]; $n{{}_{i}}$ measured by 3DP [$cm^{-3}$]; $T{{}_{p,j}}$ and $T{{}_{\alpha,j}}$ measured by SWE [$eV$]; $T{{}_{i,j}}$ measured by 3DP [$eV$]; $B{{}_{o,j}}$ measured by MFI [$nT$]; $V{{}_{p,j}}$ and $V{{}_{\alpha,j}}$ measured by SWE [$km\ s^{-1}$]; $V{{}_{i,j}}$ measured by 3DP [$km\ s^{-1}$]; $\phi{{}_{sc}}$ determined from fit process [$eV$]; $\delta\mathcal{R}$ 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}$]; $\kappa{{}_{es}}$, $p{{}_{es}}$, and $q{{}_{es}}$ from 3DP fits [N/A]; $\tilde{\chi}{{}_{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 $\Xi$ $\equiv$ $\sum_{s}\ \tilde{\chi}{{}_{s}}^{2}$, then the list is as follows:

• •

Fit Flag $\{c,h,b\}$ $=$ 0 : $\left(100\%\leq\delta\mathcal{R}\right)$ $\vee$ non-finite for any of the following: $\Xi$, $\tilde{\chi}{{}_{s}}^{2}$, $\delta\mathcal{R}$

• •

Fit Flag $\{c,h,b\}$ $=$ 1 : $\left(100\leq\tilde{\chi}{{}_{tot}}^{2}<10^{30}\right)$ $\wedge$ $\left(\left(\Xi<200\right)\vee\left(\tilde{\chi}{{}_{s}}^{2}\leq 200\right)\right)$ $\wedge$ $\left(\delta\mathcal{R}<95\%\right)$

• •

Fit Flag $\{c,h,b\}$ $=$ 2 : $\left(0\leq\tilde{\chi}{{}_{tot}}^{2}<100\right)$ $\wedge$ $\left(\left(\Xi<100\right)\vee\left(\tilde{\chi}{{}_{s}}^{2}\leq 100\right)\right)$ $\wedge$ $\left(\delta\mathcal{R}<75\%\right)$

• •

Fit Flag $\{c,h,b\}$ $=$ 3 : $\left(0\leq\tilde{\chi}{{}_{tot}}^{2}<100\right)$ $\wedge$ $\left(\left(\Xi<50\right)\vee\left(\tilde{\chi}{{}_{s}}^{2}\leq 40\right)\right)$ $\wedge$ $\left(\delta\mathcal{R}<55\%\right)$

• •

Fit Flag $\{c,h,b\}$ $=$ 4 : $\left(0\leq\tilde{\chi}{{}_{tot}}^{2}<100\right)$ $\wedge$ $\left(\left(\Xi<40\right)\vee\left(\tilde{\chi}{{}_{s}}^{2}\leq 30\right)\right)$ $\wedge$ $\left(\delta\mathcal{R}<50\%\right)$

• •

Fit Flag $\{c,h,b\}$ $=$ 5 : $\left(0\leq\tilde{\chi}{{}_{tot}}^{2}<100\right)$ $\wedge$ $\left(\left(\Xi<30\right)\vee\left(\tilde{\chi}{{}_{s}}^{2}\leq 20\right)\right)$ $\wedge$ $\left(\delta\mathcal{R}<45\%\right)$

• •

Fit Flag $\{c,h,b\}$ $=$ 6 : $\left(0\leq\tilde{\chi}{{}_{tot}}^{2}<100\right)$ $\wedge$ $\left(\left(\Xi<20\right)\vee\left(\tilde{\chi}{{}_{s}}^{2}\leq 10\right)\right)$ $\wedge$ $\left(\delta\mathcal{R}<40\%\right)$

• •

Fit Flag $\{c,h,b\}$ $=$ 7 : $\left(0\leq\tilde{\chi}{{}_{tot}}^{2}<30\right)$ $\wedge$ $\left(\left(\Xi<15\right)\vee\left(\tilde{\chi}{{}_{s}}^{2}\leq 9\right)\right)$ $\wedge$ $\left(\delta\mathcal{R}<30\%\right)$

• •

Fit Flag $\{c,h,b\}$ $=$ 8 : $\left(0\leq\tilde{\chi}{{}_{tot}}^{2}<30\right)$ $\wedge$ $\left(\left(\Xi<10\right)\vee\left(\tilde{\chi}{{}_{s}}^{2}\leq 7\right)\right)$ $\wedge$ $\left(\delta\mathcal{R}<20\%\right)$

• •

Fit Flag $\{c,h,b\}$ $=$ 9 : $\left(0\leq\tilde{\chi}{{}_{tot}}^{2}<15\right)$ $\wedge$ $\left(\left(\Xi<7\right)\vee\left(\tilde{\chi}{{}_{s}}^{2}\leq 5\right)\right)$ $\wedge$ $\left(\delta\mathcal{R}<15\%\right)$

• •

Fit Flag $\{c,h,b\}$ $=$ 10 : $\left(0\leq\tilde{\chi}{{}_{tot}}^{2}<7\right)$ $\wedge$ $\left(\left(\Xi<5\right)\vee\left(\tilde{\chi}{{}_{s}}^{2}\leq 3\right)\right)$ $\wedge$ $\left(\delta\mathcal{R}\leq 10\%\right)$

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, the 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 $\tilde{\chi}{{}_{s}}^{2}$, $\tilde{\chi}{{}_{tot}}^{2}$, and $\delta\mathcal{R}$ 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 $\sim$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 non-fill, 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 (e.g., see Section 3.2) and to avoid issues during regridding and/or interpolation for comparison with other datasets (e.g., magnetic fields). These post-fit constraints are 1.5 $<$ $\kappa{{}_{eh}}$ $\leq$ 20, 1.5 $<$ $\kappa{{}_{eb}}$ $\leq$ 20, 0 $\leq$ $n{{}_{eh}}/n{{}_{ec}}$ $\leq$ 0.75, 0 $\leq$ $n{{}_{eb}}/n{{}_{ec}}$ $\leq$ 0.50, 0.0 $\leq$ $n{{}_{eb}}/n{{}_{eh}}$ $\leq$ 3.0, 11.4 eV $\leq$ $T{{}_{eh,j}}$ $\leq$ 285 eV, and 11.4 eV $\leq$ $T{{}_{eb,j}}$ $\leq$ 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 users interested in future use by the heliospheric and astrophysical communities. The data product will benefit current and future missions like Parker Solar Probe in addition to providing a statistical comparison with astrophysical shocks, which was currently not available.