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

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

Background and Motivation

The solar wind is an ionized gas experiencing collective effects where Coulomb collisions occur, but the rates are often so low that, for instance, two constituent particle species, s′ and s, are not in thermodynamic or thermal equilibrium, i.e., (T/T)tot ≠ 1 for s′ ≠ s, and the relevant scale lengths are orders of magnitude smaller than the collisional mean free path (e.g., Wilson et al. 2018). Therefore, for any process dependent on scales like the thermal gyroradii, ρ, or inertial lengths, λ, the media is considered collisionless (see Appendix A for definitions). That the solar wind is a nonequilibrium, weakly collisional, kinetic gas results in multicomponent velocity distribution functions (VDFs) for both ions (e.g., Kasper et al. 2006, 2012, 2013; Maruca et al. 2011; Maruca & Kasper 2013; Wicks et al. 2016) and electrons (e.g., Schwartz & Marsch 1983; Maksimovic et al. 1997, 1998; Lin 1998; Pierrard et al. 1999, 2001; Štverák et al. 2008, 2009; Pulupa et al. 2014a).

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

Figure 1.

Illustrative example VDFs of the core, halo, and beam/strahl components of the electron VDFs in the solar wind near 1 au. The top row (panels (a) through (c)) shows contours of constant phase-space density (cm−3 km−3 s+3) of a two-dimensional cut through a three-dimensional VDF. The plane and coordinate basis are defined by the quasi-static magnetic field, , and the ion bulk flow velocity, . The vertical axis is defined by the unit vector ( × ) ´ and the horizontal axis by . The bottom row (panels (d) through (f)) shows one-dimensional cuts of the VDF along the horizontal (solid red line) and along the vertical (solid blue line). The locations of these cuts are defined by the color-coded crosshairs in the top row of panels. The VDF is shown in the ion bulk flow rest frame.

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

An illustrative example, showing the three electron components typically observed in the solar wind near 1 au below ~1.2 keV, is shown in Figure 1. The component parameters are exaggerated[10] for illustrative purposes but are based on the fit results of the VDF shown in Figure 4. The core is modeled by a symmetric bi-self-similar VDF and the halo and beam/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 qbi-Maxwellian (see Section 3.1). This example is phenomenologically consistent with the majority of solar wind electron VDFs (e.g., Pilipp et al. 1987a, 1987b, 1987c; Phillips et al. 1989a, 1989b; Štverák et al. 2008, 2009).

Figure 4.

Example VDF observed at 02:55:41.008 UTC on 1999 August 4 by the Wind/3DP EESA Low detector. Panel (a) shows a 2D cut through the 3D VDF as contours of constant phase-space density, where the cut plane is defined by the unit vectors ( × ) × on the vertical and on the horizontal, where = (+6.41, −7.64, −8.48) (nT, GSE). The origin in velocity space is defined by = (−388.38, +3.13, −32.63) (km s−1, GSE). The value of ϕsc for this VDF is 6.35 eV. Projected onto panel (a) are the following vectors: ion bulk flow velocity or bulk (purple arrow), (cyan arrow), shock normal vector sh (green arrow), and the Sun direction (magenta arrow). The small cyan circles show the location of actual measurements prior to regularized gridding with Delaunay triangulation. Panels (b) and (c) show the 1D parallel cuts along the horizontal (solid red line is data in both panels), and panels (d) and (e) show the 1D perpendicular cuts along the vertical (solid blue line is data in both panels). Panels (b) and (d) show the individual electron component fit results, while panels (c) and (e) show the sum of the fit results all as dashed lines and with color-coded labels. Panel (c) shows the one-count level for reference.

Despite its collisionless, nonequilibrium nature, the solar wind can support the existence of shock waves. That the particles are in neither thermal nor thermodynamic equilibrium leads to a nonhomogeneous partition of energy not only among 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 nonhomogeneous 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 (see, e.g., Coroniti 1970; Tidman & Krall 1971; Sagdeev 1966; Kennel et al. 1985; Treumann 2009; Wilson 2016; Wilson et al. 2017, and references therein). The mechanisms can also be dependent on pitch angle and species (e.g., Sagdeev 1966; Artemyev et al. 2013, 2014, 2015, 2016, 2017a, 2017b, 2018). 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, 1983a, 1983b; Thomsen et al. 1985, 1987, 1993; Wilson et al. 2010; Masters et al. 2011). 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 multistage 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 nonhomogenous 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 on either the velocity moments or a model VDF (e.g., Schunk 1975, 1977; Schwartz & Marsch 1983; Schwartz et al. 1988; Livadiotis 2015, 2017; Nicolaou et al. 2018; Shizgal 2018).

In this first part of a multipart study we describe the methodology and numerical analysis techniques used to model the solar wind eVDFs below ~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 VDF 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 is justified herein using physically significant arguments. A benefit of the analysis is an improved, semianalytic 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 et al. 2019a) the statistical results of the model fits are presented with comparison to previous studies and associated discussions. In Paper III (Wilson 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).

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 February 26 to 23:04:00 UTC on 2000 February 20. Additional supplemental material, including a PDF file containing the list of interplanetary shock event dates with associated parameters, shock parameter definitions, shock normal technique definitions, additional statistics in the form of histograms, and additional information about the model VDFs used herein, can be found at [doi:10.5281/zenodo.2875806] (Wilson et al. 2019c). The supplemental material also includes two ASCII files of fit results described in Appendix F. The symbol/parameter definitions are found in Appendix A.

Quasi-static magnetic field vectors () were measured by the Wind/MFI dual, triaxial fluxgate magnetometers (Lepping et al. 1995) using the 3 s cadence data for each particle distribution. The components/directions of some parameters are defined with respect to using the subscript j. That is, the parallel (j = ‖) and the perpendicular components (j = ⊥) of any vector or pseudo-tensor (e.g., temperature) are defined with respect to .

The electron VDFs were measured by the Wind/3DP low-energy (i.e., few eV to ~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 have cadences of ~3 s and ~24–78 s, respectively. The energy and angular resolutions are commandable, but the instrument typically operates with ΔE/E ~ 20% and Δϕ ~ 5°−22°.5 depending on the poloidal anode[11] (see, e.g., Wilson et al. 2009, 2010, for instrument details).

The EESA Low measurements are contaminated with photoelectrons from the spacecraft, something that must be accounted for to obtain accurate velocity moments or any other results. The details of how the spacecraft potential, ϕsc, was numerically determined for each VDF are 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 ϕ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 those discussed in Wilson et al. (2018) were used herein.

The VDFs examined are found within ±2 hr of 52 IP shocks found in the Wind shock database from the Harvard Smithsonian Center for Astrophysics.[12] Of those 52 IP shocks, there were 16 quasi-parallel (θ ⩽ 45°), 36 quasi-perpendicular (θ > 45°), 45 low Mach number (〈M〉up < 3), and 7 high Mach number (〈M〉up ⩾ 3) shocks. The shock parameters for the 52 IP shocks examined in this three-part set of papers are shown in Table 1 (see, e.g., Wilson et al. 2019c, for a 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.

Table 1

Shock Parameters

Parameter	Xmin[a]	Xmax[b]	X¯[c]	X˜[d]	X25%[e]	X75%[f]	σx[g]
〈|Bo|〉up (nT)	1.04	17.4	5.96	5.59	3.99	7.10	3.01
〈ni〉up (cm−3)	0.60	21.3	8.34	8.00	3.70	12.1	5.32
〈βtot〉up (N/A)	0.03	3.86	0.50	0.38	0.19	0.60	0.60
〈|Vshn|〉up (km s−1)	155	699	460	456	383	535	123
〈|Ushn|〉up (km s−1)	36.9	401	126	110	83.3	145	70.2
θBn (deg)	17.1	88.6	56.8	54.6	42.7	73.3	19.5
〈MA〉up (N/A)	1.06	15.6	2.79	2.41	1.86	3.06	2.10
〈Mf〉up (N/A)	1.01	6.39	2.12	1.86	1.58	2.35	0.94
〈Mf〉up/Mcr (N/A)	0.41	5.14	1.08	0.91	0.77	1.19	0.70
〈Mf〉up/Mww (N/A)	0.06	2.49	0.36	0.18	0.11	0.32	0.51
〈Mf〉up/Mgr (N/A)	0.04	1.91	0.28	0.14	0.09	0.25	0.39
〈Mf〉up/Mnw (N/A)	0.04	1.76	0.26	0.13	0.08	0.23	0.36

Notes. For symbol definitions, see Appendix A.

Minimum.

Maximum.

Mean.

Median.

Lower quartile.

Upper quartile.

Standard deviation.

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 nonequilibrium, 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, , the symmetric bi-self-similar function was used as the default core model function. In the downstream of strong (i.e., 〈M〉up ≿ 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 function.[13] Note that of all the core VDFs fit to a symmetric bi-self-similar function, ~80.5% satisfied 2.0 ⩽ s ⩽ 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, ϕsc (details found in Appendix B), and define fit parameter initial guess values and constraints to yield stable solutions for the most VDFs (detailed steps found in Appendix C, and list of initial guess values and constraints found in Supplemental Material ASCII files (Wilson et al. 2019c) described in Appendix F). The process of defining the initial guess values and constraints is discussed in Section 3.2, and the quantified estimates of the fit quality are discussed in Section 3.3.

A total of 15,314 electron VDFs were observed by the Wind spacecraft within ±2 hr of 52 IP shocks. Of those 15,314 VDFs, 15,210 progressed to fit analysis, and stable model function parameters were found for 14,847 (~98%) core fits, 13,871 (~91%) halo fits, and 9567 (~63%) beam/strahl fits. The reason for the large disparity in beam/strahl fits compared to the other two components will be discussed in Section 3.3 and Appendix C.

Velocity Distribution Functions

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

The most common VDF used to model particle VDFs in space plasmas is the bi-Maxwellian (e.g., Feldman et al. 1979a, 1979b, 1983b; Kasper et al. 2006), given by where A is given by where v is the drift speed of the peak relative to zero along the jth component, is the thermal speed given by Equation (6c), V is the velocity ordinate of the jth component, and n 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., Mace & Sydora 2010; Pulupa et al. 2014b; Lazar et al. 2015a, 2015b, 2016, 2017, 2018; Livadiotis 2015; Livadiotis et al. 2018; 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., Vasyliunas 1968; Feldman et al. 1983a; Maksimovic et al. 1997; Salem et al. 2003). It is beyond the scope of this study to explain the physical interpretation/origin of this function, but there are several detailed discussions already published on the topic (e.g., Livadiotis 2015; Livadiotis et al. 2018). A generalized power-law particle distribution is given by a bi-kappa VDF (e.g., Mace & Sydora 2010; Livadiotis 2015), for electrons here as where A is given by and B is given by where Γ(z) is the Riemann gamma function of argument z and V is again the most probable speed of a 1D Gaussian for consistency, i.e., it does not depend on κ.

The last model VDF is called a self-similar distribution, which results when a VDF evolves under the action of inelastic scattering (e.g., Dum et al. 1974; Dum 1975; Horton et al. 1976; Horton & Choi 1979; Jain & Sharma 1979; Goldman 1984) or flows through disordered porous media (e.g., Matyka et al. 2016). The symmetric form is given by where A is given by Note that V is again the most probable speed of a 1D Gaussian for consistency, i.e., it does not depend on s. Further, one can see that Equation 3(a) reduces to Equation 1(a) in the limit where s → 2. The function in Equation 3(a) 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 where A is given by Again, this will reduce to a bi-Maxwellian in the limit where p → 2 and q → 2. Note that in the event that the exponents s, p, or q are not even integers, the velocity ordinates, (V‖ − v) and (V⊥ − v), will become absolute values to avoid complex roots and negative values of f (V‖, V⊥). Example one-dimensional cuts of these three model VDFs can be found in Figure 2 for comparison.

Figure 2.

Examples of one-dimensional cuts through multiple model VDF functions to illustrate the functional dependence on various parameters. The top row (panels (a) through (c)) shows the dependence on the thermal speed, denoted generically as Vth here. The bottom row (panels (d) and (e)) shows the exponent dependencies. Panel (a) shows bi-Maxwellian VDFs (Equation 1(a)), panels (b) and (d) show bi-kappa VDFs (Equation 2(a)), and panels (c) and (e) show bi-self-similar VDFs (Equation 3(a) and 4(a)). All examples shown have the same number density of 5 cm−3, denoted generically as n here.

The self-similar exponents are mostly a new variable, since most previous work modeled the core electrons as a bi-Maxwellian (e.g., Štverák et al. 2008, 2009; Bale et al. 2013; Pulupa et al. 2014b). There are a few studies that used one-dimensional self-similar functions to model a select few electron VDFs near collisionless shocks (e.g., Feldman et al. 1983a, 1983b), finding values consistent with those presented in Table 2. However, these studies did not define the normalization parameter in terms of the number density and thermal speeds (see, e.g., Equation 3(a) and 4(a)), but rather found a numerical value from empirical fits, i.e., the normalization parameter was not coupled to the physical parameters of the fit function. At least one study in the solar wind did define the normalization constant, but they only considered a one-dimensional, isotropic distribution (e.g., Marsch & Livi 1985). Although several theoretical works predicted ranges of possible self-similar exponent values under various extrema scenarios (e.g., Dum et al. 1974; Dum 1975; Horton et al. 1976; Horton & Choi 1979; Jain & Sharma 1979; Goldman 1984), this is the first time the model has been used on a statistically significant set of VDFs.

Table 2

Electron Exponent Parameters

Exponent	Xmin[a]	Xmax[b]	X¯[c]	X˜[d]	X25%[e]	X75%[f]
All: 15,210 VDFs
κec	2.14	100.0	9.15	7.92	5.40	10.2
sec	2.00	3.00	2.03	2.00	2.00	2.04
pec	2.00	5.43	3.09	3.00	2.20	4.00
qec	2.00	3.29	2.24	2.00	2.00	2.46
κeh	1.51	19.7	4.62	4.38	3.58	5.34
κeb	1.52	20.0	4.57	4.17	3.40	5.16
Upstream Only: 6546 VDFs
κec	2.14	100.0	9.15	7.92	5.40	10.2
sec	2.00	2.31	2.01	2.00	2.00	2.03
pec	N/A	N/A	N/A	N/A	N/A	N/A
qec	N/A	N/A	N/A	N/A	N/A	N/A
κeh	1.52	18.4	4.16	4.10	3.25	4.83
κeb	1.52	19.6	4.22	3.81	3.25	4.70
Downstream Only: 8664 VDFs
κec	N/A	N/A	N/A	N/A	N/A	N/A
sec	2.00	3.00	2.05	2.01	2.00	2.06
pec	2.00	5.43	3.09	3.00	2.20	4.00
qec	2.00	3.29	2.24	2.00	2.00	2.46
κeh	1.51	19.7	4.94	4.62	3.80	5.70
κeb	1.53	20.0	4.82	4.45	3.61	5.44
〈Mf〉up < 3 Only: 12,988 VDFs
κec	2.14	100.0	9.02	6.83	4.40	9.93
sec	2.00	3.00	2.03	2.00	2.00	2.04
pec	2.00	5.43	3.10	3.00	2.18	4.00
qec	2.00	3.14	2.26	2.01	2.00	2.49
κeh	1.51	19.7	4.54	4.34	3.58	5.26
κeb	1.52	20.0	4.62	4.20	3.46	5.19
〈Mf〉up ⩾ 3 Only: 2222 VDFs
κec	4.32	27.2	9.30	8.60	6.89	10.4
sec	2.00	2.30	2.03	2.00	2.00	2.08
pec	2.00	5.00	3.08	2.50	2.18	4.00
qec	2.00	3.29	2.16	2.00	2.00	2.50
κeh	1.60	19.2	5.06	4.68	3.62	6.05
κeb	1.52	18.8	4.25	3.84	2.89	4.94
θBn > 45° Only: 10,940 VDFs
κec	4.05	27.2	7.77	7.18	4.84	9.11
sec	2.00	2.31	2.02	2.00	2.00	2.05
pec	2.00	5.43	3.00	2.62	2.17	4.00
qec	2.00	3.29	2.28	2.04	2.00	2.56
κeh	1.51	19.7	4.73	4.44	3.67	5.47
κeb	1.52	20.0	4.67	4.20	3.33	5.33
θBn ⩽ 45° Only: 4270 VDFs
κec	2.14	100.0	16.0	11.7	10.0	14.5
sec	2.00	3.00	2.06	2.00	2.00	2.04
pec	2.00	4.28	3.29	4.00	4.00	4.28
qec	2.00	3.00	2.14	2.00	2.00	2.16
κeh	1.55	19.4	4.32	4.18	3.37	5.09
κeb	1.53	16.5	4.31	4.10	3.57	4.82

Notes. For symbol definitions, see Appendix A.

Minimum.

Maximum.

Mean.

Median.

Lower quartile.

Upper quartile.

The following is an illustrative example that shows how the signal-to-noise ratio of particle detectors strongly depends on the number density and thermal speed and that hot, tenuous plasmas are much more difficult to measure and accurately model. Examine the 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 = 1500 to 5500 km s−1 drops the peak phase-space density by nearly two orders of magnitude. The cut line also passes the ±20,000 km s−1 velocity boundary (i.e., roughly the upper energy bound of the EESA Low instrument) at a phase-space density roughly one order of magnitude higher than the colder examples. That is, the change in thermal speed reduced the dynamic range of observed phase-space densities by three orders of magnitude. Suppose that one examines a more extreme example with n = 15 cm−3 and V = 10,000 km s−1. In this case, the difference between the peak and the lowest phase-space density within the ±20,000 km s−1 velocity boundary would only be a factor of ~55, i.e., slightly more than one order of magnitude.

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

Core

n

V or T

v

s

p

q

κ

Halo

n

V or T

v

κ

Beam/Strahl

n

V or T

v

κ.

For more details about derivation and normalization constants, see the Supplemental Material (Wilson et al. 2019c).

Fit Parameter Constraints

This section involves the discussion of the constraints/limits placed on fit parameters for each electron component and justifies them based on physically significant assumptions.

As an illustrative example, Figure 3 shows the densities of the protons, alpha-particles, and three electron components (blue squares) and the associated uncertainties (red error bars) for a subcritical, quasi-perpendicular IP shock (see, e.g., Wilson et al. 2019c, for shock parameters) observed by Wind on 1996 April 2 at 10:07:57.525 UTC. For this event, the plasma parameters are listed below in the form Min–Max (Mean)[Median]:

Figure 3.

Example IP shock crossing observed on 1996 April 2 by the Wind spacecraft. The panels are as follows from top to bottom: || (nT), (nT, GSE); value of spacecraft potential used for fits ϕsc (eV); n (red line) and 100 × n (blue line) (cm−3, SWE); s (blue circles), κ (green circles), and κ (magenta circles); n values (blue circles) and uncertainty (red error bars) (cm−3, 3DP fit); n (cm−3, 3DP fit); and n (cm−3, 3DP fit). The error bars for the four electron fit parameter panels are defined by the percent deviation discussed in Section 3.3. The error for this date satisfied with a median of 10.3%.

Upstream

|| ~ 0.53–3.14(1.96)[1.53] nT;

n ~ 11.3–15.8(12.0)[11.9] cm−3;

n ~ 0.06–0.18(0.10)[0.11] cm−3;

s ~ 2.00–2.09(2.00)[2.00] N/A;

κ ~ 2.83–12.2(4.46)[4.40] N/A;

κ ~ 1.67–12.6(4.85)[5.10] N/A;

n ~ 10.7–13.0(11.7)[11.5] cm−3;

n ~ 0.06–1.44(0.69)[0.54] cm−3;

n ~ 0.02–0.17(0.09)[0.09] cm−3;

Downstream

|| ~ 3.45–5.99(4.85)[5.19] nT;

n ~ 14.9–19.7(18.0)[18.1] cm−3;

n ~ 0.14–0.27(0.19)[0.19] cm−3;

s ~ 2.00–2.07(2.01)[2.01] N/A;

κ ~ 2.72–6.96(4.39)[4.29] N/A;

κ ~ 2.74–7.27(4.45)[4.50] N/A;

n ~ 13.6–18.4(16.7)[16.8] cm−3;

n ~ 0.02–2.53(0.56)[0.44] cm−3;

n ~ 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/n ⩾ 1. Figure 3 is illustrative of some of the error analysis employed in the present study and the fact that the beam/strahl fit more often fails than the core or halo as evidenced by the number of points. Below the details of how the fit parameters are constrained/limited are outlined with physical arguments.

First, the present study differs from some previous studies in that the fits are performed on the two-dimensional VDF rather than separate fits on one-dimensional cuts of the two-dimensional VDF (e.g., Maksimovic et al. 2005; Pulupa et al. 2014a, 2014b). One of the limitations of the latter approach is that the distribution function is not necessarily a separable function, which can introduce difficulty for the physical interpretation of the results. However, the latter approach has numerous advantages, including the stability of the solutions and ease with which the solutions are found with nonlinear least-squares software, i.e., it is generally easier to fit to a one-dimensional cut than a two-dimensional distribution.

The present study uses the former approach to avoid the difficulties introduced for nonseparable functions. For instance, when fitting to the parallel one-dimensional cut, the amplitude of the VDF is directly tied to the amplitude of the perpendicular cut. The amplitude of all standard model two-dimensional, gyrotropic VDFs is dependent on n, , and . 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 can be problematic if care is not taken. For instance, the normalization constants differ for one-dimensional cuts from the two-dimensional gyrotropic VDF (see, e.g., Equation 1(a)). Although this approach involves fewer free parameters and should thus be easier to fit, it is much more restrictive in parameter space, i.e., n 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 ⩾ 10,000 km s−1, the weights would not provide sufficient contrast between the peak and tails to constrain a stable and reliable fit without multiple imposed constraints. In contrast, electron VDFs with thermal speeds below ~1000 km s−1 fall below the lowest energy of the detector and so would be artificially hotter if they were observed (e.g., Paschmann & Daly 1998). A similar effect is often observed by spacecraft with electrostatic analyzers designed for the magnetosphere, not the comparatively cold, fast solar wind beam (e.g., McFadden et al. 2008a, 2008b; Pollock et al. 2016).

Statistical studies of the solar wind have shown that the maximum range of the total electron temperature is T ~ 2.29–77.2 eV or V ~ 450–2600 km s−1 (e.g., Wilson et al. 2018). Previous studies have found that the electron halo temperatures satisfy T ~ 14–560 eV or V ~ 1100–7000 km s−1 (e.g., Feldman et al. 1975, 1978, 1979a; Maksimovic et al. 1997, 2005; Skoug et al. 2000; Tao et al. 2016a, 2016b; Lazar et al. 2017). Previous studies have also found that the electron beam/strahl temperatures satisfy T ~ 20–150 eV or V ~ 1300–3600 km s−1 (e.g., Ogilvie et al. 2000; Viñas et al. 2010; Tao et al. 2016a, 2016b). Thus, a range of allowed core thermal speeds from ~1000 to ~10,000 km s−1 can be assumed.

There are similar instrumental constraints on the drift speed of the three components. The core, however, is not likely to exhibit drift speeds (in the ion rest frame) in excess of several hundred kilometers per second (e.g., Pulupa et al. 2014a). In the present work, most fit results show less than 50 km s−1, i.e., only 1838 of 14,847, or ~12%, have drift speeds exceeding 50 km s−1, consistent with previous work.[14] In contrast, owing to the physical interpretation of the strahl/beam component, most (8848 of 9567, or ~92%) have drift speeds in excess of 1000 km s−1. The allowed core, halo, and beam/strahl drift speeds loosely ranged from ~1000 to ~10,000 km s−1 for most events. In some events, a lower bound was imposed to prevent unphysical fit results, e.g., beam/strahl component with near zero drift speed (see Supplemental Material ASCII files (Wilson et al. 2019c) described in Appendix F for ranges for specific events). Note that V was fixed during the fitting, i.e., it was not allowed to vary. Originally this parameter was free to vary but resulted in fewer stable fits and rarely varied by more than a few kilometers per second. In some events, an explicit V was set as the initial guess values determined from examination of the distributions, but this is for a small minority of events (333 of 14,847, or ~2%).

It has also been empirically found that the EESA Low detector has issues when n ≾ 0.5 cm−3 or n ≿ 50 cm−3 for typical solar wind thermal speeds.[15] This is rarely an issue, as only 41 of the 14,847 VDFs analyzed (or ~0.3%) have fit results falling outside the range ~0.5–50 cm−3. Note that the total electron density, n = n + n + n ~ n = n + 2n, is constrained by the total ion density from SWE and the total electron density from the upper hybrid line observed by the WAVES radio receiver (Bougeret et al. 1995), when possible (see Appendix B for more details).

Physically, the halo and beam/strahl components are suprathermal; thus, they should not have the dominant contribution to the total phase-space density of the VDF. Therefore, it is physically consistent to assume that the fit results should satisfy n/n < 1 and n/n < 1. The solutions were constrained to satisfy n/n < 0.5 and n/n < 1 based on results found in previous studies near 1 au (e.g., Feldman et al. 1975; Maksimovic et al. 1997, 2005; Skoug et al. 2000; Štverák et al. 2009; Viñas et al. 2010; Pierrard et al. 2016; Tao et al. 2016b).

In numerous previous studies that assumed a three-component solar wind electron VDF near 1 au (e.g., Maksimovic et al. 2005; Štverák et al. 2009; Pulupa et al. 2014a, 2014b), constraints were sometimes assumed such as that the fits satisfy n/n < 1. There is no restriction on this ratio[16] imposed during the fit process, and 1824 of 9313 (or ~20%) of the fits satisfy n/n ⩾ 1. In fact, it was found that imposing the constraint n/n < 1 during the fit process actually greatly reduced the number of stable solutions found for the beam/strahl component.[17] Previous work did show that the ratio n/n decreases with increasing radial distance from the Sun, dropping below unity before 1 au, on average, but the ranges overlapped, allowing for n/n ⩾ 1 (e.g., Štverák et al. 2009).

Another constraint that is often assumed/used is that the strahl/beam component be only antisunward along (e.g., Maksimovic et al. 2005; Štverák et al. 2009; Pulupa et al. 2014a, 2014b), though some magnetic field topologies have sunward-directed beam/strahl components (e.g., Owens et al. 2017). This constraint is imposed in this study, but it is important to note that some IP shocks examined have observable electron foreshocks. A consequence is that the halo component of the fit results effectively absorbs both the halo and the shock-reflected electron component in the events where this is directed sunward along (this is very rare). If the shock-reflected electron component is directed antisunward, they will be included in the beam/strahl fit (this is much more common). The net result for the former is a smaller (T⊥/T‖) and on the latter is a larger (T⊥/T‖) and n.

The lower bound of possible κ values is defined for mathematical/physical reasons as being ≿3/2 (e.g., Livadiotis 2015; Livadiotis et al. 2018). The upper bound is set to 100 solely because above that value the difference between a bi-Maxwellian and bi-kappa VDF is smaller than the accuracy of the measurements. Although the upper bound is allowed to extend to 100, the typical upper bound observed near 1 au is <20 (e.g., Maksimovic et al. 1997; Štverák et al. 2009; Pierrard et al. 2016; Tao et al. 2016a, 2016b; Lazar et al. 2017). The range of possible values for s, p, or q falls between 2 and 10 for physical reasons (e.g., Dum et al. 1974; Dum 1975; Horton et al. 1976; Horton & Choi 1979; Jain & Sharma 1979; Goldman 1984).

Finally, by definition the halo and beam/strahl components represent the lowest-energy suprathermal components of the electrons. Therefore, it is natural to assume that T/T > 1. There is no explicit restriction on this ratio imposed, and only 384 of 13,867 (or ~3%) of the fits satisfy T/T < 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/T in any way other than to fit to the data.

Quality Analysis

The initial approach was to use the reduced chi-squared value 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 while others did not fit well at all despite having . The issue is partly related to the calibration of the detector and thus the quality of the values (see Appendix B for more details). The issue is also related to fitting a gyrotropic model function to data that is not, in general, gyrotropic. A possible improvement would fold the entire VDF into a forced gyrotropy prior to fitting to improve counting statistics and the comparison between data and model functions, but that is beyond the scope of the current study. Therefore, a new quantity was defined to provide an additional definition of the quality of any given fit by direct comparison.

Let us use f (0) as the actual data and f ( (=f (core) + f (halo) + f (beam)) as the total model fit results. Then one can define the ratio of these two parameters as , which is a two-dimensional array of values. Then one calculates the median of this array, , to determine the percent deviation given by where is computed for each electron VDF. The values of were then used as uncertainties/error bars for all fit parameters for the associated VDF for all components. In general, the percent magnitude of the uncertainty in each of the six fit parameters should not be uniform as is used herein (see Appendix E for discussion of 1σ uncertainties). The uncertainty of any variable calculated using these fit parameters was propagated assuming uncorrelated errors.

Note that the value alone does not always characterize the quality of any given fit. Therefore, a combination of parameters is chosen to define a set of fit quality flags from best with a value of 10 to worst with a value of 0 (see Appendix F for definitions). In general, fits with flags of at least 2 or higher can be used, but low fit flags should be treated with caution. Only ≾1% of all core, halo, and beam/strahl fits had flags of 1, while >95% of core, >89% of halo, and >61% of beam/strahl flags were at least 2.

Figure 4 shows an example VDF that had a low for each component and a , i.e., this is an example of an ideal fit. The distribution was fit using a symmetric bi-self-similar distribution for the core and a bi-kappa distribution for both the halo and beam/strahl components. The fit results are as follows:

n = {15.43, 2.01, 0.056} cm−3;

V = {1959.6, 2500.0, 3964.7} km s−1;

V = {1937.9, 2575.5, 4516.2} km s−1;

V = {+44.58, −0.00, −3898.7} km s−1;

V = {−0.00, −0.00, −0.00} km s−1;

{s, κ, κ} = {2.00, 4.58, 2.57}, where s is the self-similar exponent and κ is the kappa value;

;

; and

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

In contrast, Figure 5 shows an example VDF that had a high for two components yet still a small , i.e., this is still an example of a good fit despite the bad values for the core and beam/strahl fits. The fit results are as follows:

Figure 5.

Another example VDF observed at 20:22:43.490 UTC on 1999 January 22 by the Wind/3DP EESA Low detector in burst mode. The format is the same as in Figure 4, where this VDF has = (−6.95, +9.78, −8.77) (nT, GSE), = (−619.12, +26.66, +21.19) (km s−1, GSE), and ϕsc = 9.45 eV.

n = {4.41, 0.57, 0.32} cm−3;

V = {3882.6, 2624.5, 4574.5} km s−1;

V = {2728.2, 2986.3, 2387.6} km s−1;

V = {−0.00, −594.9, +2000.0} km s−1;

V = {−0.00, −0.00, −0.00} km s−1;

{p, q, κ, κ} = {4.00, 2.00, 2.27, 4.61}, where p(q) is the parallel (perpendicular) self-similar exponent and κ is the kappa value;

;

; 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 halo.[18] Thus, the resultant f ( can exceed f (0) in some places. The software does a post-fit check for instances where either the combined or any component model fit exceeds the data by user-specified factors.[19] For most events, the threshold is set between ~2 and 4, but this varies, as some events have known issues. For instance, the known density from the upper hybrid line is 10 cm−3, but no variation of ϕsc yields fit results with n ~ 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 for the core component and moderate for beam/strahl but a small . This example VDF was chosen to illustrate a good fit even when n/n > 1. As previously discussed, there are post-fit constraints applied to the data based on statistical and physical constraints. The constraint relevant to Figure 6 is that requiring n/n < 3. This is why the fit flag value for the beam/strahl is zero and why is larger than a few. The fit results are as follows:

Figure 6.

Another example VDF observed at 18:23:06.116 UTC on 1999 January 22 by the Wind/3DP EESA Low detector in burst mode. The format is the same as in Figures 4 and 5, where this VDF has = (−0.89, −0.32, −10.57) (nT, GSE), = (−626.59, +93.06, +76.13) (km s−1, GSE), and ϕsc = 10.67 eV.

n = {3.37, 0.03, 0.14} cm−3;

V = {2609.8, 5293.2, 4686.9} km s−1;

V = {2286.9, 5494.9, 2516.2} km s−1;

V = {−0.00, −222.8, +3273.0} km s−1;

V = {−0.00, −0.00, −0.00} km s−1;

{s, κ, κ} = {2.00, 3.83, 3.53};

;

; 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/n > 1 even at 1 au.

After examining thousands of fit results, it was determined that with and are consistently more reliable quantities used in combination for defining the quality of the fit than using alone. The value is also used as a proxy for the uncertainty of any given fit parameter, e.g., 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 and values are listed for reference to typical/expected values when evaluating the quality of a fit. In general, the best fits have small values for and all .

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 ~107 counts s–1 (Lin et al. 1995). Examination of all VDFs found that a total of 10 energy-angle bins (from a total of 20,184,120), or ~5 × 10−5%, exceeded the maximum count rate. Therefore, it is not thought that saturation has a significant impact on the methodology and results of this study. Second, as illustrated in Figure 3, the total electron density satisfies n ~ n + 2n for nearly all intervals. Statistically, the difference between the fit result for n = n + n + n and n + 2n is within expectations. The median, lower quartile, and upper quartile values are 10.3%, 4.9%, and 19.0%, respectively, which is consistent with our statistics.

Finally, the total electron current, j = Σ n v, in the ion rest frame should be zero to maintain a net zero current in the solar wind. The mean, median, lower quartile, and upper quartile for all data examined are ~22 km s−1 cm−3, ~0 km s−1 cm−3, ~−214 km s−1 cm−3, and ~351 km s−1 cm−3, consistent with previously published work on this data set (e.g., Bale et al. 2013; Pulupa et al. 2014a) and consistent with work in progress (C. S. Salem et al. 2019, in preparation). Normalizing j by n times V yields a mean, median, lower quartile, and upper quartile for all data examined of ~0.17%, ~10−8%, ~−0.95%, and ~1.3%, respectively. Thus, the values are all small compared to unity. Quantitatively, ~97.5% of the j/(nV) values satisfy ≾5.5%.

Figure 7 shows both j and j/(nV) versus seconds from every shock ramp center time in this study. One can see that although there are locations with significant deviation from zero (e.g., the shock ramp, which is not tremendously surprising, as that is where currents are supposed to exist), the mean (red horizontal line) and median (orange horizontal line) are small for both the raw and normalized current densities. Note that the data in Figure 7 include fit results where there may not be a solution for one or more components (see discussion of first data product ASCII file in Appendix F).

Figure 7.

Two superposed epoch analysis plots of the total electron current density, j (Mm s−1) (top panel), and normalized values, j/(nV) (%) (bottom panel), vs. seconds from the shock ramp center. Shown in each panel are the lower (Q1) and upper (Q2) quartiles as magenta lines, the mean as a red line, and the median as an orange line for all data. That is, the lines are computed for the entire set of data, not at each time stamp. For reference, the axis ranges were defined as 110% of the maximum of the absolute value of X2.5 and X97.5, where X2.5 and X97.5 are the bottom 2.5th and top 97.5th percentiles.

As a final note, there is the question about the validity of using a new model function to describe the thermal core. Of the 11,874 core VDFs fit with a symmetric bi-self-similar model function, there were 9559, or ~80.5%, that satisfied 2.0 ⩽ s ⩽ 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. 1979a, 1979b). In fact, these results show that most core VDFs are not far from thermal velocity distributions, consistent with results showing evidence for collisional effects on the core (e.g., Salem et al. 2003; Bale et al. 2013).

Summary of Fit Results

For the 52 IP shocks examined there were a total of 15,314 VDFs observed by Wind. Of those 15,314 VDFs, 15,210 progressed to fit analysis, and for the core only 534 (~4%) were modeled as bi-kappa VDFs, 12,095 (~80%) were modeled as symmetric bi-self-similar VDFs, and 2581 (~17%) were modeled as asymmetric bi-self-similar VDFs. All core 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 are given in Section 3 and Appendix C. Of those 15,210 that progressed to fit analysis, stable solutions were found for 14,847 (~98%) f (core), 13,871 (~91%) f (halo), and 9567 (~63%) f (beam).

Recall that the fit results presented herein were performed on two-dimensional, (assumed) gyrotropic velocity distributions in the proton bulk flow rest frame. Most prior work numerically fit to one-dimensional cuts of the VDF or to one-dimensional reduced VDFs. There are benefits for either method, but here it is shown that the method employed is valid by illustrating the consistency with previous work. The statistical results of the densities are summarized below in the form lower quartile–upper quartile(Mean)[Median]: which are consistent with previous results near 1 au (e.g., Feldman et al. 1975, 1979a, 1983b; Maksimovic et al. 1997; Phillips et al. 1989a, 1989b; Nieves-Chinchilla & Viñas 2008; Skoug et al. 2000; Salem et al. 2001; Štverák et al. 2009; Pierrard et al. 2016). The full statistical results and associated histograms are presented in Paper II.

All

n ~ 6.44–19.5(13.7)[11.3] cm−3;

n ~ 0.21–0.63(0.52)[0.36] cm−3;

n ~ 0.09–0.27(0.21)[0.16] cm−3;

Upstream

n ~ 4.06–12.5(8.90)[8.09] cm−3;

n ~ 0.17–0.49(0.42)[0.27] cm−3;

n ~ 0.09–0.26(0.22)[0.16] cm−3;

Downstream

n ~ 8.44–24.2(17.3)[16.6] cm−3;

n ~ 0.26–0.70(0.59)[0.44] cm−3;

n ~ 0.09–0.28(0.21)[0.17] cm−3;

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

All

;

;

;

;

;

Upstream

;

;

;

;

;

Downstream

;

;

;

;

.

The purpose of listing these statistics is to provide a range of typical or expected and values for reference when determining the quality of any given fit. Note that the statistics for shown above were performed on arrays that excluded the lower and upper boundaries, i.e., 0.1% and 100% values. The statistical results of the model function exponent and drift speed results are presented below, and the full data product resulting from this work is described in Appendix F.

Exponents and Drifts

Table 2 shows the one-variable statistics for the exponents from the model fits of the electron VDFs for the core (s = c), halo (s = h), and beam/strahl (s = b). The VDFs, modeled as bi-kappa (κ), symmetric bi-self-similar (s), and asymmetric bi-self-similar velocity distributions (p for parallel and q 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, and therefore the converse had no results to examine.

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

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

There are several interesting things to note from Table 2. The mean, median, and lower/upper quartile values for κ are slightly higher for high than for low Mach number shocks, though only the median and lower quartile values are significant. Since a bi-kappa model was only used for upstream core VDFs, this may imply that shock strength is somehow dependent on the upstream core electron distribution profiles. One possible physical interpretation would be that the sound speed depends on the polytropic index for each species, i.e., the equation of state assumed for the system. A bi-kappa core VDF could affect the estimate of the sound speed, thus altering the fast mode Mach number. However, the shape of the upstream VDFs will also affect the shock dissipation mechanisms. For instance, it is known that the existence of power-law tails improves the efficiency of shock acceleration (e.g., Trotta & Burgess 2019). Therefore, the larger κ associated with higher Mach number shocks may imply that lower energy particles have entered the tails, thus increasing the exponent.[20]

In contrast, the asymmetric bi-self-similar exponents, only used in downstream regions, are effectively the same between low and high Mach number shocks. However, this changes when comparing quasi-parallel and quasi-perpendicular shocks. The p exponent has higher mean, median, and lower/upper quartile values for quasi-parallel than quasi-perpendicular shocks. The opposite is true for the q exponent.

This is interesting, as higher p 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., Fuselier & Gurnett 1984; Wilson et al. 2007, 2010, 2012, 2014a, 2014b; Breneman et al. 2013), but their amplitudes increase with increasing shock strength (e.g., Wilson et al. 2007). If the largest ion-acoustic waves generate the largest values of p, 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 versus what fraction goes to increasing T. This would depend on the effective inelasticity of the wave–particle interactions, where larger inelasticity increases p and smaller inelasticity increases T (e.g., Dum et al. 1974; Dum 1975; Horton et al. 1976; Horton & Choi 1979; Jain & Sharma 1979; Goldman 1984). The interaction between a wave and a particle can be treated as inelastic if the particle affects the wave amplitude and kinetic energy during the interaction. Most test-particle treatments do not handle this self-consistently, and if the effect is distributed to an entire VDF, the net result can be a stochastic heating that increases p from 2.0 (e.g., Dum et al. 1974; Dum 1975).

Another theory predicts that flat-top electron distributions (i.e., p → ⩾4 and q → ~2–3) can result from the combined effects of a quasi-static, cross-shock electric potential and from fluctuation electric fields (e.g., Feldman et al. 1983a; Hull et al. 1998) through a process called maximal filling (e.g., Morse 1965). However, similar to the predictions for wave-driven flat tops, this theory should generate stronger flat tops (i.e., larger values of p) for stronger quasi-perpendicular shocks, which we do not observe. Thus, the evolution of the electron VDFs does not seem consistent with the standard quasi-static, cross-shock electric potential, but rather in agreement with recent high-resolution observations at the bow shock (e.g., Chen et al. 2018; Goodrich et al. 2018).

Another interesting result is the difference in the κ values under different conditions. When the values of κ are larger (smaller), that implies a less (more) energized halo, i.e., softer (harder) spectra. One can see that κ is larger downstream than upstream and near high rather than low Mach number shocks. That is, the halo is less energized downstream of IP shocks and near strong IP shocks than the converse, which is somewhat unexpected, as strong shocks should more readily energize suprathermal particles (e.g., Malkov & Drury 2001; Treumann 2009; Caprioli & Spitkovsky 2014; Park et al. 2015; Trotta & Burgess 2019). In contrast, κ is slightly smaller (~10%) near quasi-parallel than quasi-perpendicular shocks, which implies more energized halo electrons. Although quasi-parallel shocks are predicted (e.g., Malkov & Drury 2001; Caprioli & Spitkovsky 2014) and observed (e.g., Wilson 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 that the upstream electron suprathermal tail will become flatter (i.e., smaller kappa values) with increasing Mach number for quasi-perpendicular shocks (Trotta & Burgess 2019). This may explain why both κ and κ 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 κ values show a similar behavior between upstream and downstream and shock geometry as κ, 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 κ, s, p, q, and the drift speed magnitudes, V (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 = 3 and 4 in panel (c) are for strong shocks exhibiting flat-top 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.

Figure 8.

Histograms of the exponents (top row) and bulk drift velocity magnitudes (bottom row) for the different electron components for all time periods as percentage of total counts. Panel (a) shows the κ values for the core (violet), halo (blue), and beam/strahl (red) components. Panel (b) shows the s for the core (violet). Panel (c) shows the p (blue) and q (red) values for the core. Panels (d)–(f) show the magnitude of the parallel (violet) and perpendicular (blue) drift velocities for the core, halo, and beam/strahl components, respectively. The statistics for the exponents are listed in Table 2. Note that the tick marks are individually labeled in all panels.

One can see that, as discussed previously, the core parallel drift speeds (violet line, panel (d)) tend to fall below ~100 km s−1, consistent with previous results (e.g., Pulupa et al. 2014a). In fact, most of the core and halo drifts are near zero, with the number of results satisfying V ⩽ 1 km s−1 and V ⩽ 1 km s−1 being 8735 (~59%) and 7311 (~53%), respectively. Note that although there is sometimes a sizable perpendicular core drift (blue line, panel (d)) for some shock crossings, these were explicitly set after visual inspection of the VDFs during the iterative fitting process. The nonzero perpendicular drifts almost certainly result from inaccuracies in the calculation of the solar wind rest frame and a dipole correction to ϕsc not included in the present analysis (e.g., Pulupa et al. 2014a; see Appendix B for more details).

The magnitudes of V and V never deviated from zero.[21] The magnitudes of V range from ~0 to 8860 km s−1, with a lower to upper quartile range of ~0–850 km s−1 and a mean (median) of ~580 km s−1 (~0.1 km s−1). The magnitudes of Voeb,‖ range from ~1000 to 9330 km s−1, with a lower to upper quartile range of ~1750–3090 km s−1 and a mean (median) of ~2580 km s−1 (~2480 km s−1). As previously discussed, the lower bound for V was imposed on the basis of physical arguments, while the magnitude of V was allowed to go to zero. If only magnitudes satisfying V ⩽ 1 km s−1 are considered, the mean (median) and lower to upper quartile ranges are ~42 km s−1 (~30 km s−1) and ~14–52 km s−1 for V and ~1227 km s−1 (~903 km s−1) and ~362–1695 km s−1 for Voeh,‖.

Discussion

A total of 15,314 electron VDFs were observed by the Wind spacecraft within ±2 hr 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 (~98%) core fits, 13,871 (~91%) halo fits, and 9567 (~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 (~4%) VDFs, as a symmetric bi-self-similar for 12,095 (~80%) VDFs, and as an asymmetric bi-self-similar for 2581 (~17%) VDFs. This is the first statistical study to find that the core electron distribution is better fit to a self-similar VDF than a Maxwellian under all conditions.

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

All

s ~ 2.00–2.04(2.03)[2.00];

p ~ 2.20–4.00(3.09)[3.00];

q ~ 2.00–2.46(2.24)[2.00];

κ ~ 5.40–10.2(9.15)[7.92];

κ ~ 3.58–5.34(4.62)[4.38];

κ ~ 3.40–5.16(4.57)[4.17];

Upstream

s ~ 2.00–2.03(2.01)[2.00];

p ~ N/A;

q ~ N/A;

κ ~ 5.40–10.2(9.15)[7.92];

κ ~ 3.25–4.83(4.16)[4.10];

κ ~ 3.25–4.70(4.22)[3.81];

Downstream

s ~ 2.00–2.06(2.05)[2.01];

p ~ 2.20–4.00(3.09)[3.00];

q ~ 2.00–2.46(2.24)[2.00];

κ ~ N/A;

κ ~ 3.80–5.70(4.94)[4.62];

κ ~ 3.61–5.44(4.82)[4.45].

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

The interesting aspect of VDFs being well modeled by bi-self-similar functions is that such functions are used to describe the evolution of distributions for either the flow through disordered porous media (e.g., Matyka et al. 2016) or the influence of inelastic scattering (e.g., Dum et al. 1974; Dum 1975; Horton et al. 1976; Horton & Choi 1979; Jain & Sharma 1979; Goldman 1984). It is unlikely that the former applies directly, but the latter may be interpreted in the following manner. The typical approach for test-particle simulations used to examine wave–particle interactions does not include feedback from the particles on the waves. In a real plasma, the particles can alter three properties of electromagnetic waves: their amplitude (potential energy), momentum, and kinetic energy. Consider a simple scenario whereby a particle reflects off of an electromagnetic wave field along one dimension. If done self-consistently, the particle can reduce the wave amplitude in addition to affecting the field momentum and kinetic energy. In the case of a reduced wave amplitude, the resulting scattering problem can be treated as a simple inelastic collision.[22] Thus, the net result of an ensemble of particles interacting with a wave field can be stochastic (e.g., Dum et al. 1974; Dum 1975), which provides one physical justification for the use of the bi-self-similar functions. These functions are also convenient in that they reduce to bi-Maxwellians in the limit where the exponents go to 2, i.e., the deviation from a Maxwellian is a measure of inelasticity in the particles’ interactions with waves and/or turbulence.[23] Further, as previously discussed, ~80.5% of the core VDFs modeled with a symmetric bi-self-similar function had exponents satisfying 2.0 ⩽ s ⩽ 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. 1979a, 1979b) and work that found evidence for collisional effects in the core distribution (e.g., Salem et al. 2003; Bale et al. 2013).

The κ seem to correlate with 〈M〉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 are higher for quasi-parallel shocks, while q are higher for quasi-perpendicular shocks, yet neither depends on 〈M〉up.

Somewhat surprisingly, the values of κ are larger downstream than upstream, and they increase with increasing 〈M〉up. That is, the halo spectra are softer downstream and near strong shocks. Quasi-parallel shocks, however, correlate with smaller κ, i.e., harder halo spectra. Generally, quasi-parallel shocks are predicted to be more efficient particle accelerators for suprathermal ions and very energetic electrons[24] (e.g., Caprioli & Spitkovsky 2014), but electrons in the halo energy range are predicted to be energized the most efficiently at shocks satisfying θ > 80° (e.g., Park et al. 2013).

Unlike the halo, κ are smaller near high Mach number shocks than near low Mach number shocks. The difference is likely a twofold consequence of the combined effects from shock-accelerated foreshock electrons and the method used to fit the distributions. That is, the beam/strahl component is always fit to the antisunward, field-aligned side of the VDF, while the halo is fit to the opposite. For nearly all IP shocks at 1 au, the shock normal is antisunward in a direction that would be aligned with the nominal, ambient beam/strahl electron component. For both the halo and beam/strahl, the ratios of 〈k〉/〈k〉up and 〈k〉/〈k〉up increase with increasing 〈M〉up. That is, the downstream halo and beam/strahl spectra are softer than the upstream for stronger shocks. Again, this is likely a consequence of the foreshock electrons that are not observed upstream of weak shocks. The details of the electron component velocity moments and associated changes will be discussed further in Papers II and III.

In summary, the first part of this three-part study presented the first statistical study to find that the core electron distribution is better fit to a self-similar VDF than a bi-Maxwellian under all conditions. This is an important result for kinetic theory and solar wind evolution. This work also provides the methodology and details necessary to reproduce and qualify the results of the nonlinear least-squares fitting performed herein. In Papers II and III, the statistical and analysis results of the velocity moments will be presented in detail. These observations are relevant for comparisons with astrophysical plasmas like the intra-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.

Table 3

Spacecraft Potential Statistics

ϕsc(eV)	Xmin[a]	Xmax	X¯	X˜	X25%	X75%
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
〈Mf〉up < 3 only: 12,932 finite values	1.01	26.7	6.99	6.61	5.44	7.70
〈Mf〉up ⩾ 3 only: 2212 Finite Values	3.58	12.0	7.35	6.90	5.50	9.63
θBn > 45° only: 10,894 finite values	1.01	26.7	6.70	6.49	5.35	7.38
θBn ⩽ 45° only: 4250 finite values	3.53	17.6	7.94	7.14	6.10	10.2

Notes. For symbol definitions, see Appendix A.

Header symbols match those of Table 2.