Electron phase space holes (EHs) associated with electron trapping are commonly observed as bipolar electric field signatures in both space and laboratory plasma. Until recently, it has not been possible to resolve EHs in electron measurements. We report observations of EHs in the plasma sheet boundary layer, here identified as the separatrix region of magnetic reconnection in the magnetotail. The intense EHs are observed together with an electron beam moving toward the X line, showing signs of thermalization. Using the electron drift instrument onboard the satellites of the Magnetospheric Multiscale mission, we make direct millisecond measurements of the electron particle flux associated with individual electron phase space holes. The electron flux is measured at a millisecond cadence in a narrow parallel speed range within that of the trapped electrons. The flux modulations are of order unity and are direct evidence of the strong nonlinear wave–electron interaction that may effectively thermalize beams and contribute to transforming directed drift energy to thermal energy.

## I. INTRODUCTION

Wave–particle interactions are one of the fundamental ways energy is transferred and redistributed in fully or partially collisionless plasmas. Electron phase space holes (EHs) are one example of strong wave–particle interactions. They are nonlinear plasma structures formed as a result of electron trapping by electric fields.^{1,2} EHs are typically considered as time-stationary solutions to the Vlasov–Poisson equations and are often referred to as Bernstein–Green–Kruskal (BGK) modes.^{1} They are common in both space^{3–11} and laboratory plasmas,^{12,13} where they are typically identified by their divergent bipolar electric field parallel to the ambient magnetic field or positive electrostatic potential. The positive potential structure is due to a local depletion of electrons with respect to ions at the center of the EHs. EHs are a type of electrostatic solitary waves (ESWs) and typically form during the nonlinear stage of streaming instabilities, such as electron two-stream instabilities,^{14} Buneman-type instabilities,^{15,16} and ion-acoustic instabilities.^{17} Their strong wave–particle interaction can lead to electron heating, beam thermalization,^{14,18–20} electron scattering,^{21} and electron–ion coupling.^{17,18,22} As they are often formed in regions of strong plasma flows, they are commonly observed during magnetic reconnection,^{23} perhaps predominantly along the separatrices,^{20,24,25} but also inside the exhaust.^{7,26}

Although ESWs are commonly observed, their identification as EHs in space plasmas has historically been indirect through their electric field signature. While a single EH is typically observed over a few milliseconds, particle instruments have accumulation times of tens of ms to seconds, making time-resolved measurements of EHs difficult. Not until recently, using the Magnetospheric Multiscale (MMS) mission was the actual depletion in electron phase space observed. By using sub sweep measurements of the dual electron spectrometer (DES) of the Fast Plasma Investigation^{28} (FPI), Mozer *et al.*^{27} made a superposed epoch analysis of a group of several EHs. That is, the parallel electron phase space was recreated by combining measurements at varying energy sampled at different times. In this paper, using the Electric Drift Instrument^{29} (EDI) onboard the MMS spacecraft, we employ fast, millisecond sampling of the electron distribution function at a single but constant parallel speed where electrons efficiently interact with EHs. Therefore, we are able to directly observe the continuous spatiotemporal variations in the electron distribution function and quantify the strong EH–electron interaction.

In a broader context, we know that wave–particle interactions of varying strengths are ubiquitous in plasmas; for example, they occur in turbulent environments in planetary magnetospheres and in interplanetary, stellar, and interstellar space. Due to the often small time and length scales involved, studying the physics of these wave–particle interactions is only possible *in situ* within the Heliosphere and requires that instrumentation dedicated to wave–particle interaction (see, e.g., Dombrowski *et al.* and references therein^{30}) of varying strength and character be a part of future space-plasma missions.

This manuscript is organized in the following manner: in Sec. II, we present an overview of the magnetotail event and detailed measurements of the ESW properties that identify them as EHs; in Sec. III, we use the measured wave properties to reconstruct a 1D model of electron phase space and compare it to the observed electron flux; and in Sec. IV, we discuss our results.

## II. OBSERVATIONS

We present burst mode data from July 6, 2017 when MMS crossed the plasma sheet boundary layer (PSBL) in the Earth's magnetotail, at about $[\u221225,1,4]$ Earth radii in geocentric solar ecliptic (GSE) coordinates. The electric field, sampled at 8192 Hz, is from the electric field double probes (EDP) comprising the spin-plane double probes (SDP)^{31} and the axial double probes (ADP).^{32} The magnetic field, sampled at 128 Hz, is from the FluxGate Magnetometer (FGM).^{33} The plasma distributions and moments are from the FPI,^{28} where the ion and electron data are sampled at 6.7 and 33 Hz, respectively. The electron flux from the EDI^{29} is sampled at 1000 Hz. The data interval examined is 13:53:44–13:54:15. All times are given in coordinated universal time (UTC). While all measurements are made in the frame of the spacecraft, the spacecraft moves at a speed of <1 km/s and, thus, its speed can be neglected.

In the plasma sheet, prior to entering the southern PSBL, tailward (*v _{x}* < 0) ion and electron flows are observed [Fig. 1(b)], indicating the presence of a magnetic reconnection X line Earthward of MMS.

^{34}Inside the PSBL, approximately between 13:54:00 and 13:54:15, the reduced electron distribution along $B$, $fe(t,v\u2225)=\u222b\u222b\u2212\u221e\u221efe(t,v\u2225,v\u22a5,1,v\u22a5,2)dv\u22a5,1dv\u22a5,2$, where $v\u2225$ and $v\u22a5,1,2$ are the speeds parallel and perpendicular to the ambient magnetic field

**B**, respectively, shows counterstreaming electron populations [Fig. 1(c)]. A lower energy population ($\u221220\xd7103<v\u2225<5\xd7103$ km/s) streams toward the X line, while a higher energy electron population ($5\xd7103<v\u2225<80\xd7103$ km/s) streams away from the X line. The largest amplitude parallel electric fields [Fig. 1(e)] are observed at around 13:54:06, indicated by the yellow shaded region, inside a density cavity [Fig. 1(d), $n\u22480.04$ cm

^{−3}]. This environment is similar to previous observations from the PSBL and magnetic reconnection separatrices

^{20,35,36}and is consistent with the separatrix environment in simulations of symmetric magnetic reconnection.

^{37,38}The counterstreaming populations at magnetic reconnection separatrices are often unstable to streaming instabilities and the generation of plasma waves.

^{20,22,25,39,40}

### A. Properties of electrostatic solitary waves

The parallel electric fields inside the density cavity are largely composed of solitary bipolar structures, often referred to as ESWs [Fig. 2(a)]. For many of the ESWs, we also observe a finite perpendicular electric field component, with magnitudes typically below, but comparable to the parallel component [Fig. 2(b)]. This indicates that the parallel length scales are comparable to the perpendicular length scales of the structures.^{41–44}

During the event, MMS is in a tetrahedron formation with an average inter-spacecraft separation of 16 km. For comparison, the local Debye length is $\lambda De\u22481$ km. Thirteen of the ESWs shown in Fig. 2 are observed by all four spacecraft. Performing multi-spacecraft timing,^{45} we can unambiguously determine the propagation velocity of these ESWs. The propagation speeds *v _{ph}* range between −10 500 and −7300 km/s, with an average of $vph,av\u2248\u22128500$ km/s. The obtained velocities are all within $19\xb0$ of $\u2212B$.

Figure 3 shows the range of *v _{ph}* overlaid onto the reduced one- and two-dimensional electron distributions $fe(v||)$ and $fe(v||,v\u22a5,1,2)$ sampled during 30 ms around the time 13:54:05.562. The unit vectors of the parallel and perpendicular velocity components are given by $v||=[\u22120.99,\u20090.08,\u22120.12]$ $v\u22a5,1=[\u22120.07,\u22121.00,0.00]$ and $v\u22a5,2=[\u22120.12,0.00,0.99]$. The ESWs all propagate in the direction of a warm component of the electron distribution inside the PSBL [Figs. 3(a) and 3(b), inner dashed circle], slightly above its bulk speed. Together with the warm population, we also observe a hot population moving away from the X line (outer dashed circle). In comparison, the ion bulk speed is $|vix|\u227250$ km/s, and the ion thermal speed is $vti,\u2225\u223c300$ km/s. The relatively large

*v*s suggest that the ESWs were formed as a consequence of a streaming instability related to the warm streaming electron population.

_{ph}^{20}However, the measured population is not unstable to wave generation in the range of measured

*v*s. It is, therefore, likely that the electrons have already been affected by the waves' presence and heated. In $fe(v||)$ [Fig. 3(c)], we can see that the warm population observed at the same time as the ESWs forms a plateau and is heated relative to the cold population observed closer to the lobe. Although $fewarm(v||)$ is clearly not Maxwellian, we estimate characteristic thermal energies to be $Te,||warm\u223c400$ eV. $fewarm(v\u22a5),\u2009felobe(v||)$, and $felobe(v\u22a5)$ (not shown) are closer to Maxwellian with weak thermal tails or kappa distributions. Their characteristic energies are $Te,\u22a5warm\u223c200$ eV and $Te,||lobe\u223cTe,\u22a5lobe\u223c70$ eV. These observations also suggest that the electron populations have been affected by the waves. We note that both the parallel and perpendicular temperatures of the warm population exceed those of the lobe population, indicating the electrons might have been heated not only in the direction parallel to the magnetic field, but also perpendicular to it. Such scattering could have been made possible by the finite perpendicular electric field of the ESWs.

_{ph}^{21}

Using the measured propagation velocity, we can determine length scales and electrostatic potentials, along **B**, of the ESWs: $\varphi =\u2212\u222bE\u2225dl\u2225=\u222bE\u2225vphdt$. We note that both the peak-to-peak lengths and the maximum potentials for a given ESW vary between the spacecraft. This suggests they are either temporally evolving on timescales comparable to the delay between the spacecraft ($\u223c1$ ms) or that they have an irregular spatial structure. However, we find no correlation between peak-to-peak lengths and time delay, which you would expect for an expanding or shrinking structure. The mean peak-to-peak lengths for each ESW vary between *L _{pp}* = 5.4 and 10.5 km with an average of 7.0 km. The maximum potentials $\varphi max$ for each ESW vary between 140 and 410 V with an average of 240 V. To calculate the electrostatic potential shown in Fig. 2(c), we use $vph,av$, and first high pass-filter $E\u2225$ at 30 Hz. Thereafter, for the benefit of the phase space modeling in Sec. III B, $\varphi $ is detrended such that $\varphi =0$ at semi-regular intervals [shown by * at the bottom of Fig. 6(b) for MMS 1]. This further removes the effects of low-frequency components of $E||$. This does neither mean that such effects are not present nor that they are unimportant.

^{7,46}In fact, a weak parallel electric field that aids the acceleration of electrons toward the X line is expected in the magnetic reconnection separatrix region.

^{47}However, due to the low amplitude of such a field, it is hard to measure accurately.

While the ESWs do show multi-dimensional properties, estimating their 3D structure is not straightforward. The ESWs are embedded in a boundary layer with finite perpendicular electric fields as well as additional wave modes. Using a method similar to Tong *et al.*,^{42} we find that the two perpendicular length scales generally differ and are typically at least 2–3 times larger than the parallel length scale. These aspect ratios are toward the more oblate end but comparable to other studies from a similar region.^{42,44} However, we note that the uncertainties are large, and these numbers are only to be considered as rough estimates.

From **E**, we can also calculate the density perturbation (i.e., the deviation from charge-neutrality) associated with the ESWs, $\delta n=ni\u2212ne=(\epsilon 0/e)\u2207\xb7E$, using single- and multi-spacecraft techniques. With four spacecraft, we can calculate the divergence of **E**. In the centers of the structures, $\delta n>0$ [Fig. 2(d)], corresponding to a relative excess of positive charge or absence of negative charge. The relative density perturbation is $\delta n/n0\u223c1%$, where $n0=0.04$ cm^{−3} is the time- and spacecraft-averaged density inside the density cavity [Fig. 1(d)]. The positive potentials and charge perturbations are consistent with the electron depletion of EHs. For some of the ESWs, we can also see a negative charge perturbation at the edges of the structures, indicating the excess of negative charge acting as a charge-neutralizing cloud.^{42} Since the spacecraft separation is comparable to, but slightly larger than the ESW length scales, *δn* calculated from four spacecraft is likely underestimated. We, therefore, also estimate *δn* using the separate spacecraft individually. Since the relative trajectories of the spacecraft through the ESWs are roughly parallel to **B**, we use the parallel component of the divergence of the electric field, i.e., the first term in $\delta n=(\epsilon 0/e)(\u2202\u2225E\u2225+\u2207\u22a5\xb7E\u22a5)=(\epsilon 0/e)(\u2212vph\u22121\u2202tE\u2225+\u2207\u22a5\xb7E\u22a5)$. The results, also shown in Fig. 2(d), with $\delta n/n0\u223c4%$, are larger than the method using all four spacecraft. The negative charge perturbations at the edges of the ESWs are also more visible when using the single spacecraft technique. Since we cannot directly access the perpendicular electric field derivatives using a single spacecraft, we cannot directly evaluate the contribution of the $\u2207\u22a5\xb7E\u22a5$ term to *δn*. However, since we know that $E\u22a5=0$ for many of the ESWs, it is likely that this term is also important. The relative importance of the $\u2207\u22a5\xb7E\u22a5$ term depends on the 3D structure of the ESWs.^{48–52} For example, considering Poisson's equation for a potential structure in the form of a triple Gaussian $\varphi =\varphi 0\u2009exp(\u2212x2/2l\u22a512\u2212y2/2l\u22a522\u2212z2/l||2)$, the density difference at the center of the structure is $\delta n=(\epsilon 0/e)\varphi 0(l\u22a51\u22122+l\u22a52\u22122+l||\u22122)$. Here, $l||=Lpp/2$ corresponds to half the peak-to-peak length scale. Given a moderately oblate structure with $l\u22a51=2l\u22a52=4l||$, as opposed to a 1D structure with $l\u22a51=l\u22a52\u2192\u221e$, the difference in *δn* is a factor 1.3125. We note that with such a modification, the density perturbation would still remain only a fraction of the background density.

Another proxy for the electron density, or rather the electron flux, is the spacecraft potential measured using the SDP probe potentials.^{53,54} In a sunlit plasma, photons impact the spacecraft, and as a result, electrons are ejected, charging the spacecraft positively. Ambient electrons impacting the spacecraft instead charge the spacecraft negatively. Figure 2(e) shows the spacecraft potential $\delta Vsc$ high-pass filtered at 100 Hz. There are rapid variations closely anticorrelated with *δn*. At the center of the ESWs, $\delta Vsc>0$ as a result of less electrons impacting the spacecraft. The decreases in $\delta Vsc$ seen at the edges of several ESWs instead indicate a relative increase in electrons impacting the spacecraft. Due to the low plasma density, the charging time of the spacecraft ($\u223c10$ ms) is estimated to be larger than the observation timescale of the ESWs ($\u223c1\u20132$ ms). Therefore, the changes in the spacecraft potential can only be used as a qualitative estimate of the density variations.

### B. Electron flux observed by the electric drift instrument

EDI consists of two detectors, one on either side of the spacecraft. Each detector has 32 channels with $11.25\xb0$ width in the instrument azimuth plane out of which four adjacent channels can be selected. Electrostatic optics allow the selection of the polar look angle with respect to the instrument symmetry axis in $0.35\xb0$ step sizes and can be configured to look parallel or perpendicular to the magnetic field in a variety of ways. On June 7, 2017, EDI was configured to observe parallel and antiparallel electrons by centering channel 1 on (against) the magnetic field and having channels 2–4 extend toward higher (lower) pitch angles. Configured in such a way, EDI was collecting electrons at a range of pitch angles $\theta =[0,45]\xb0$ and $\theta =[135,180]\xb0$, divided into four sectors each. The energy range of EDI was set to $EEDI\xb112\Delta EEDI$, where *E _{EDI}* = 500 eV and $\Delta EEDI=50$ eV.

EDI counts are calibrated to fluxes by correcting for instrument dead time, flat-fielding, and then applying an absolute scaling factor. Dead time correction takes into account the known refresh rate of the electronics between samples. Flat fielding ensures each channel produces the same measurement under constant conditions; as the spacecraft spin changes the detector orientation with respect to the magnetic field, the active channels change to maintain a constant look direction with respect to the magnetic field. In quiet conditions when the ambient electron flux is constant, a correction to constant counts can be determined as a function of azimuth (channel) and polar angle. The absolute scaling factor transforms counts to fluxes. Because of the complicated optical properties of the instrument, there is no single geometric factor that can be applied to the data. Instead, EDI counts are trended against FPI DES parallel (0°–60° pitch angle) and anti-parallel (120°–180° pitch angle) electron fluxes from the nearest energy channel (typically 485 eV) to produce a single scale factor. Errors associated with the scale factor are not quantified but are estimated to be on the order of 20%.

Correcting for additional electron acceleration by the positively charged spacecraft of *V _{sc}* = 36 V, the EDI energy interval corresponds to field aligned electron speeds $v\u2225,EDI\xb1=vEDI\u2009cos\u2009\theta $, where $vEDI\u2212=12\u2009430$ km/s and $vEDI+=13\u2009120$ km/s ($vEDI+\u2212vEDI\u2212=690$ km/s). The collecting interval is 1 ms, providing electron particle fluxes $jeEDI$ at 1000 samples/s. In comparison, the ESWs peak-to-peak timescales are roughly 1 ms.

In Sec. II A, we noted that the variations in spacecraft potential associated with the ESWs were potentially underestimated due to the relatively slow spacecraft charging time compared to the ESW observation time. Even if these changes are underestimated by an order of magnitude, say $\delta Vsc*\u223c2V$ instead of the measured $\delta Vsc\u223c0.2$ V, they are small in comparison to both the sampling energy of EDI, $\delta Vsc*/EEDI\u223c0.004$, and the average spacecraft potential $\delta Vsc*/Vsc\u223c0.05$. The changes in spacecraft potential due to the ESWs will, therefore, have a negligible effect on the electron measurements by EDI.

Figure 2(f) shows the antiparallel electron flux $jeEDI$ in the pitch angle range $\theta =[168.75,180]\xb0$. The electron flux is characterized by large amplitude variations, $\delta jeEDI/jeEDI\u223c1$, on all four spacecraft. The flux approaches zero at regular intervals. The variations in $jeEDI$ are anti-correlated with the variations seen in $\varphi $ and *δn*. To highlight the anti-correlation, a side-by-side comparison of $\varphi $ and $jeEDI$ for MMS 1 is shown in Fig. 4. This anti-correlation supports the interpretation of the ESWs as EHs, since the electron depletion in the center of EHs should be associated with decreased electron flux. The parallel electron flux measured by EDI (not shown) is very low in comparison ($\u22720.5\xd7106s\u22121cm\u22122sr\u22121$) for all four spacecraft and shows no correlation with the EHs.

To better illustrate the inter-spacecraft differences, and compare to FPI, Fig. 2(g) shows $jeEDI$ resampled to the FPI cadence. MMS 1 and MMS 2 almost consistently measure the lowest and highest fluxes, respectively, with $je,MMS2EDI/je,MMS1EDI\u22722$. Figure 2(h) shows the antiparallel electron flux $jeFPI$ as measured by FPI in the energy and pitch angle ranges 485 ± 65 eV and $\theta =[157.50,180]\xb0$, respectively. This FPI energy channel is the closest one to that of EDI. We choose to plot an increased pitch angle range (with respect to $\Delta \varphi EDI=11.25\xb0$) to reduce some of the noise related to low counts caused by the relatively small geometric factor of the FPI detectors. The time-averaged ratios for the displayed interval are $\u27e8jeEDI/jeFPI\u27e9=[1.2,\u20092.3,\u20091.0,\u20092.3]$ for MMS 1–4, respectively. Based on the comparison between $jeEDI$ and $jeFPI$, and the small spacecraft separation, it is possible that some of the inter-spacecraft differences in $jeEDI$ may be due to instrumental effects. Other effects may be due to the 3D structure of the EHs and the beam structure, which are outside the scope of this present paper. For this time interval, $jeEDI$ from MMS1 agrees well with the flux measured by FPI, and we, therefore, use it in the following analysis.

## III. MODELLING OF ELECTRON FLUX

To investigate if the flux oscillations observed by EDI are consistent with the variation in flux expected for EHs with the observed $\varphi $ and *v _{ph}*, we model the expected phase space density under the simplifying assumption of a 1D EH structure. This approach has been applied for environments where the electrons are strongly magnetized, such that only motion parallel to the ambient magnetic field needs to be considered.

^{48,49,55}For this event ($B\u223c20$ nT), the thermal gyroradii for the warm ($Te,\u22a5\u223c200$ eV) electron population is $rewarm=2.4$ km. Assuming perpendicular length scales $l\u22a5=3l\u2225=3Lpp/2=10.5$ km, where we have used the average observed

*L*= 7.0 km, we obtain that $l\u22a5=4.4rewarm$. The electrons are, thus, moderately magnetized, and it is possible that 3D effects such as pitch angle scattering

_{pp}^{21}and temporary trapping

^{52}may be present. However, including 3D effects complicates the modeling and is only meaningful if we have good estimates of the 3D structure of the ESWs, which we lack. Possible 3D effects are further discussed in Sec. IV.

While the literature regarding the modeling of phase space structures is extensive, here, we connect them to the direct observations of phase space. The general solution to the unmagnetized, one-dimensional time-stationary Vlasov equation can be written in terms of the constants of motion $U=me2(v\u2212vph)2\u2212e\varphi $ and $\sigma =sign(v\u2212vph)$: $f=f(U,\sigma )$. The electrons can be divided into two classes: trapped and passing. If *U *>* *0, the kinetic energy of the electrons exceeds the potential energy of the structure, and the electrons are free to pass the electrostatic potential structure. If *U *<* *0, the potential energy exceeds the kinetic energy, and the electrons are trapped in the electrostatic potential structure.

### A. Modeling a single electron phase space hole

Trapping of electrons by an ESW potential structure is illustrated in Fig. 5 for one of the ESWs observed by MMS1 during the time 13:54:05.5837–05.5873. The phase speed of the ESW is $\u223c8600$ km/s and $\varphi max\u223c310$ V [Fig. 5(a)]. The trapped (passing) region of phase space ($x\u2212v\u2225$, where *x* and $v\u2225$ are the position and speed parallel to **B**, respectively) are indicated by blue (red) in Fig. 5(b). To obtain the phase space density of the passing electrons *f _{ep}*, we use Liouville mapping, where we use the fact that $df/dt=0$ along the phase space trajectories defined by the constant

*U*. That is, $fep(x,v)=f0(x0,v0)$, where

*f*

_{0}is an unperturbed phase space distribution at a point

*x*

_{0}where the potential $\varphi (x0)=0$. The velocity $v0=vph\xb1((v\u2212vph)2\u22122e\varphi (x)/me)1/2$ is obtained from the conservation of the total energy $U(x,v)=U(x0,v0)$. We have chosen an

*f*

_{0}that consists of three Maxwellian distributions with temperatures $Te=[50,100,\u2009and\u20091500]$ eV, drift speeds $vd=[\u221210\u2009000,\u22121000,and\u20096000]$ km/s, and densities $n=[0.016,\u20090.010,\u2009and\u20090.014]$ cm

^{−3}, respectively, such that the total density is $n0=0.04$ cm

^{−3}(Fig. 7, yellow line). While

*f*

_{0}, in principle, could have a large number of degrees of freedom, here, we have limited it to a relatively simple model that fairly well reproduces the flux observed by EDI, and the averaged flux observed over a train of ESWs by FPI. Some effects related to the choice of

*f*

_{0}will be discussed later.

The resulting *f _{ep}* is shown in Fig. 5(c) outside the black line (

*U*>

*0). The density of the passing electrons*

*n*decreases locally as the electrons are first accelerated and thereafter decelerated in the potential structure [Fig. 5(d)]. The trapped density is subsequently given by

_{ep}and is shown in Fig. 5(d). Note that $(\epsilon 0/e)\u22072\varphi =\u2212(\epsilon 0/e)vph\u22121\u2202E||/\u2202t$ is the density perturbation *δn* calculated using a single spacecraft and shown in Fig. 2(d). As mentioned before, since we do not take into account perpendicular derivatives, this density perturbation is to be considered as a lower estimate. In the frame of the waves, the ion kinetic energy largely exceeds the potential energy of the ions due the wave structure $mi2(vi,||+vti,||\u2212vph,av)2/e\varphi \u2273103$, where $|vi,|||\u227250$ km/s is the ion bulk drift speed and $vti,||\u2248300$ km/s is the thermal speed of the ions. The density perturbation of the ions is, therefore, assumed to be negligible, and we, thus, treat it as a constant $ni=n0=0.04$ cm^{−3}.

To model the phase space density of the trapped electrons, *f _{et}*, we employ the integral approach following Bernstein

*et al.*

^{1}The trapped phase space density

*f*is related to

_{et}*n*as

_{et}where we have used $\varphi (x)=0$ for $x\u2192\u221e$. As suggested by Hutchinson,^{57} we do not solve for $fet(U)$ directly, but for the difference $fet(U)\u2212fe(U=0)$, where $fe(U=0)$ is the phase space density at the separator between free and trapped electrons. $netflat(\varphi )$ is, therefore, the difference between *n _{et}* and the density corresponding to a distribution with constant phase space density $fe(U=0)$: $netflat(\varphi )=net(\varphi )\u2212fe(U=0)8e\varphi /me$. This integral equation is solved numerically for each $\varphi (x)$ and the range of

*U*corresponding to the trapped trajectories. Since $\varphi $ and, therefore,

*n*[as seen in Figs. 5(a) and 5(d)] are obtained at a relatively low cadence, using finite differences to calculate $dnetflat/d\varphi $ gives a noisy result. We, therefore, make an 8th order polynomial fit to $netflat(\varphi )$ and take the analytical derivative of this function to obtain $dnetflat/d\varphi $. The accuracy to which the derivative can be estimated can affect the structure of the phase space hole. For example, while changing to a second-order polynomial has next to no effect on the flux in the EDI energy range, it affects the trapped density significantly. In particular, the high order polynomial is necessary to accurately fit values at low $\varphi $. The modeled phase space density

_{et}is shown in Fig. 5(c), and the total modeled density $nemod$ is shown in Fig. 5(d) (red). Good agreement between $nemod$ and $neobs$ (derived from $E||$) serves as a proof of method or a consistency check, i.e., it indicates that the phase space density was well modeled, and that $dnetflat(x)/d\varphi (x)$ was adequately estimated. The phase space density at the center of the EH is 30% of that at the trapped/passing boundary.

The differential particle flux $vfemod$ is shown in Fig. 5(e). The flux within a given velocity interval is subsequently calculated as $jemod=\u222bv\u2225\u2212v\u2225+vfemoddv$. For the antiparallel and parallel directions $v\u2225\xb1=\u2212vEDI\xb1$ and $vEDI\u2213$, respectively [shown as two pairs of horizontal black lines in Fig. 5(e)]. The obtained flux (red) is shown together with $jeEDI$ (blue) in Fig. 5(f). The model reproduces the large amplitude variations observed by EDI. In addition, the model shows flux enhancements at the trapped/passing boundary that are not captured by EDI. This demonstrates the qualitative agreement between the model and observations. For a quantitative comparison between the model, which is intrinsically 1D, and EDI, which measures the field-of-view flux in a limited solid angle and energy range, we require knowledge of the velocity distribution in the entire perpendicular plane $(v\u22a5,1,v\u22a5,2)$. We can only obtain this at a lower cadence through FPI that measures the full three-dimensional distribution. This is further investigated in Sec. III B, where we extend our model to a larger time interval.

### B. Modeling a time series of multiple electron phase space holes

We have performed the same modeling of the phase space density for the entire train of ESWs, using $vph,av=\u22128500$ km/s to calculate $\varphi $ [Figs. 6(a)–6(c)]. In order to calculate $dnetflat/d\varphi $ more precisely, we have divided the time series into sub-intervals [marked by * in Fig. 6(b)] for which we calculate $dnetflat/d\varphi $ separately. This sub-division is done by detrending the calculated potential such that $\varphi =0$ at semi-regular intervals. Forcing $\varphi $ to zero ensures that the trapped populations of individual EHs are separate, and that phase space is continuous. This affects the amplitudes of $\varphi $ by 40–110 V.

As can be seen, some EHs are modeled independently, while some are modeled in pairs or triplets. Because all ESWs are slightly different, $dnetflat/d\varphi $ used to construct $femod$ is not as well determined for the pairs and triplets. This can be seen in the discrepancies in the density perturbation *δn* [Fig. 6(d)]. The largest deviations between $nemod$ and $neobs$ are seen at the edges of the EHs where $neobs<0$. The reason why we do not separate all EHs into individual intervals is that it requires a higher level of detrending, which affects the amplitudes of $\varphi $ to an even greater extent.

The modeled flux in the EDI energy range, $jemod$, and measured flux, $jeEDI$, are well correlated [Fig. 6(e), only antiparallel flux shown]. However, the amplitude of the modeled variations $\delta jemod/jemod\u223c0.3$ is smaller than $\delta jeEDI/jeEDI\u223c1$. For better comparison, the modeled flux has been resampled to the EDI timeline. This has removed the small scale variations related to the trapped/passing boundary seen in Fig. 5(f).

To compare the modeled electron trapping to electron distributions observed by FPI, we time-average $femod$ displayed in Fig. 6(c) and plot it together with the time-averaged $feFPI$ and *f*_{0} in Fig. 7. The purple-shaded area shows the range of observed *v _{ph}*. The green-shaded areas show the range of velocities corresponding to the border between trapped and passing electrons at the maximum potentials, $vtrap=vph\xb12e\varphi max/me$. The widths of the shaded areas represent the variations between different ESWs/EHs. We can see how

*f*

_{0}(yellow) has been significantly modified to $\u27e8femod\u27e9$ (red) in the speed range surrounding

*v*and below, $v\u2225\u2272\u22125000$ km/s. In the upper trapping range, $v\u2225\u2273\u22125000$ km/s, $\u27e8femod\u27e9\u223cf0$. The modeled distribution $\u27e8femod\u27e9$ agrees fairly well with the reduced distribution observed by FPI, $\u27e8feFPI\u27e9$ (blue). We also note that there is a hint of a dip in $\u27e8femod\u27e9$ right at

_{ph}*v*.

_{ph}## IV. DISCUSSION

In this section, we discuss the approach of employing 1D modeling of the phase space density, what 3D effects might be important, the impact of modeling parameters on $jemod$, and quantitative comparison between $jeEDI$ and $jemod$.

In a simplified 1D model, the condition for an electron to be trapped in an electrostatic potential traveling at speed *v _{ph}* along the magnetic field direction is $(me/2)(v\u2225\u2212vph)2)\u2212e\varphi <0$. In a 3D structure, finite gyroradius effects and perpendicular electric fields complicate electron trajectories and may redistribute the kinetic energy between the parallel and perpendicular directions. This means that in the general 3D case, electrons may wander between trapped and passing trajectories, leading to velocity space diffusion.

^{52,58}However, the trajectories that only temporarily satisfy the parallel trapping condition also contribute to the charge deficit required to sustain the EHs.

^{58}As a result, the 3D-effects may affect the trapped particle population, but we expect these effects to be relatively small since the perpendicular size of the EHs is significantly larger than the gyroradii of the trapped electrons.

Finite gyroradius effects such as described above may still lead to pitch angle scattering. For this event, a large part of the electron population is close to gyroresonant, i.e., the electron crossing time is comparable to the electron gyroperiod: the speed relative to the speed of the ESWs for a gyroresonant electron is $|v\u2212vph|=fcel||=1960$ km/s, where the electron cyclotron frequency is *f _{ce}* = 560 Hz. This implies that pitch angle scattering

^{21}might be significant, which was also implied when comparing the perpendicular temperatures of the warm population to those of the lobe populations in Fig. 3, $Te,\u22a5warm/Te,\u22a5lobe\u223c3$. However, we note that the heating parallel to the ambient magnetic field is dominating: $Te,||warm/Te,\u22a5warm\u223c2$.

Most existing 3D models of EHs either assume some kind of cylindrical or spherical symmetry or require detailed knowledge of the 3D EH structure.^{58–60} Given the large range of uncertainties in EH structure for this event, attempting to include them would not rend the model more reliable, but potentially the opposite. We, therefore, deem it beyond the scope of the current study.

Grabbe^{61} employed a gyrokinetic modeling approach by assuming that the perpendicular speed was well approximated by guiding center drifts. This led to a correction in the integration limits and the denominator in Eq. (2), $E+e\varphi \u2192E+e\varphi \u2212\epsilon D$, where *ε _{D}* is the energy $mevD2/2$ corresponding to the guiding center drift $vD$. Assuming $vD=vE\xd7B=E\xd7B/|B|=1500$ km/s, for

*B*=

*20 nT and $E\u22a5=30$ mV/m, we obtain $\epsilon D=6$ eV $\u226ae\varphi $. Therefore, the expected correction is small and should not alter the results much from the simpler 1D model.*

The smaller amplitude variations of $femod$ compared to $feEDI$ can be due to a number of effects. For example, modeling an EH separately resulted in larger variations than modeling them together. This can be seen by comparing Fig. 5(f) with the same EH in Fig. 6, which is the one that has the largest *δn*. This suggests that the averaging effects related to $dnetflat/d\varphi $ play a role. Other effects may be related to the 3D structure of the EHs that are not present in the 1D model we employ, uncertainties and averaging in the estimated EH parameters, the relatively simple choice of *f*_{0}, or the nonlinear and time-dependent process of EH formation in nature. For example, when filtering $E\u2225$ and detrending $\varphi $, we removed any effects related to finite potential drops,^{46} in addition to slightly lowering the amplitudes of $\varphi $. An underestimation of $\varphi $ would lead to an overestimation of the flux inside of the EHs. We also note that shifting the peak of *f*_{0} to higher velocities, closer to the EDI measuring range, increases $jemod$ between the structures but results in a lower $\delta nmod$. Also, the closer the *v _{ph}* and

*v*are, the larger the flux variations will be inside of the EHs. Nonlinear trapping effects and scattering during the wave growth are also expected to significantly affect the distribution of passing electrons both at the locations of the EHs and in between them.

_{EDI}^{14}As such, the modeled fluxes variations should be regarded as a conservative estimate.

In the upper trapping range, closer to $v||=0$ km/s, the unperturbed distribution was not significantly modified by the waves, $f0\u223cfemod$. In this speed range, *f*_{0} was chosen as to approximately correspond to the measured distribution $feFPI$. It is possible that this lower speed range of the distribution was affected by a different set of wave modes, for example, Buneman waves. Such multistage processes have been demonstrated for magnetic reconnection separatrices both at the dayside^{19} and the nightside.^{20}

Since the model is 1D, but $jeEDI$ is a field-of-view measure with limited pitch angle range, we cannot make direct quantitative comparison between $jeEDI$ and $jemod$. However, by means of the lower-cadence FPI measurements, we can relate the two measurements quantitatively in the following manner. Figures 2(f)–2(h) show that the field-of-view fluxes of EDI and FPI showed some discrepancies but were of comparable amplitude to each other. Considering the moving window average of $jeEDI$ [Fig. 2(g)], $\u27e8je,MMS1EDI/je,MMS1FPI\u27e9\u223c1$. At the same time, Fig. 7 shows that the time-averaged modeled flux also agreed fairly well with the reduced time-averaged 1D distribution of FPI. Therefore, since the field-of-view measures of $jeEDI$ agree with those of $jeFPI$, and 1D reduced measures of $jeFPI$ agree with those of $jemod$, there is an implication that there is not only qualitative, but also quantitative agreement between $jeEDI$ and $jemod$.

## V. SUMMARY AND CONCLUSION

In this paper, we have presented observational evidence of the continuous spatiotemporal modulation of electron particle flux associated with a series of EHs. Using four closely separated MMS spacecraft, we measured the speeds, length scales, and electrostatic potentials of a group of ESWs identified as EHs in the PSBL. The EHs were associated with large amplitude fluctuations in the electron particle flux as measured by EDI antiparallel to the magnetic field in the energy range 500 ± 25 eV, $\delta jeEDI/jeEDI\u223c1$. At the centers of the EHs, $jeEDI\u21920$. The EHs were observed together with an electron beam that showed a sign of thermalization.

Through 1D modeling of the phase space density of the EHs, we made quantitative predictions of the electron flux based on $\varphi $ and *v _{ph}*. While the 1D model excludes 3D and temporal effects, we could unambiguously show that the observed flux $jeEDI$ was, indeed, associated with the observed electric field. The modeled flux was well correlated with, but showed slightly smaller relative variations than, the observed EDI flux. The time-averaged modeled distribution $\u27e8femod\u27e9$ also agreed fairly well with the reduced distribution observed by FPI, $\u27e8feFPI\u27e9$.

The observed large fluctuations in electron particle flux, associated with large variations in phase space density, show the strong interaction between the waves and electrons. This strong interaction is an efficient way to thermalize beams and contributes to transforming directed drift energy to thermal energy in physical plasmas.

## ACKNOWLEDGMENTS

The authors recognize the tremendous effort in developing and operating the MMS spacecraft and instruments and sincerely thank all involved. MMS data analysis was performed using the IRFU-MATLAB package available at https://github.com/irfu/irfu-matlab. M.R.A. was supported by NASA Contract No. NNG04EB99C. C.N. and P.T. thankfully acknowledge support from the Research Council of Norway under Contract No. 300865.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts of interest to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are publicly available at https://lasp.colorado.edu/mms/sdc/public/ through the MMS Science Data Center, Ref. 62.