Electron acceleration and thermalization at magnetotail separatrices

In this study we use the Magnetospheric Multiscale (MMS) mission to investigate the electron acceleration and thermalization occurring along the magnetic reconnection separatrices in the magnetotail. We find that initially cold electron lobe populations are accelerated towards the X line forming beams with energies up to a few keV's, corresponding to a substantial fraction of the electron thermal energy inside the exhaust. The accelerated electron populations are unstable to the formation of electrostatic waves which develop into nonlinear electrostatic solitary waves. The waves' amplitudes are large enough to interact efficiently with a large part of the electron population, including the electron beam. The wave-particle interaction gradually thermalizes the beam, transforming directed drift energy to thermal energy.


Introduction
Magnetic reconnection is a universal process where magnetic energy is often explosively released, leading to particle acceleration and heating. Observations suggest that magnetic reconnection and subsequent processes can accelerate electrons to energies of tens or even hundreds of keVâĂŹs [e.g. Hoshino et al., 2001;Øieroset et al., 2002;Vaivads et al., 2011;Fu et al., 2019]. The particle energization associated with magnetic reconnection is known to take place in several regions: in the inflow region and along the separatrices [e.g. Nagai et al., 2001;Egedal et al., 2008;Hesse et al., 2018a;Eriksson et al., 2018], inside the ion and electron diffusion regions [e.g. Torbert et al., 2018;Hesse et al., 2018b;Wang et al., 2018], in the magnetic reconnection exhaust [e.g. Bessho et al., 2015], in the vicinity of magnetic islands [Chen et al., 2008;Huang et al., 2012], both during island coalescence [Pritchett, 2008] and contraction [Drake et al., 2006], and at dipolarization fronts [e.g. Fu et al., 2011;Vaivads et al., 2011]. Where, how, and to what extent the particles are accelerated depend not only on fundamental properties such as the particle species and the relative composition of species, but also on changing properties, such as the particle's velocity. Two examples of the former is that the presence of heavier ions or cold ionospheric ions can act as energy sinks in addition to reducing the rate at which magnetic flux is being reconnected [e.g. Toledo-Redondo et al., 2017;Tenfjord et al., 2019]. One example of the latter is Fermi type B acceleration in the reconnection exhaust or in magnetic islands where the energization is more efficient if the initial velocity is higher [Northrop, 1963;Drake et al., 2006]. A clear indication of the non-uniform energization of a particle species is the fact that the energy partition is generally not uniform, with some particles being accelerated to superthermal energies, while some remain thermal [e.g. Hoshino et al., 2001]. How this energy-dependent energiza-tion affects the bulk energization of a species is unclear. For example, a study of the change in electron temperature between the magnetosheath and the reconnection exhaust during reconnection at the dayside magnetopause did not show any strong dependence on the initial electron temperature in the magnetosheath [Phan et al., 2013]. The acceleration mechanisms can vary between direct acceleration by electric fields, for example the reconnection electric field inside the diffusion regions [Bessho et al., 2015], the already mentioned Fermi acceleration [Fermi, 1949], and betatron acceleration [e.g. Northrop, 1963]. Ultimately, due to the conservation of energy, the final plasma energies must depend on the amount of available magnetic energy compared to the amount of plasma to be reconnected, which varies during the reconnection process .
In addition, energy transfer does not always occur directly between the magnetic field and the particles, but often in steps, between different plasma populations, mediated by electromagnetic fields. One such example is the Hall magnetic fields, which are due to the different motions of ion and electrons. Observations from the terrestrial magnetotail show that at the separatrices of magnetic reconnection, lower-energy field aligned electrons flow into the reconnection region while higher-energy electrons flow out of the reconnection region [Øieroset et al., 2001;Nagai et al., 2001;Asano et al., 2008] carrying the outward and inward Hall currents, respectively [Nagai et al., 2003]. The acceleration leading to the formation of these electron flows and associated currents has by some authors been suggested to be a necessity to maintain quasi-neutrality inside the ion diffusion region [Uzdensky and Kulsrud, 2006;Egedal et al., 2008]. It has been explained as following: inside the ion diffusion region, the demagnetized ions are free to move across the magnetic field while the magnetized electrons are tightly bound to the magnetic field lines. As a magnetic flux tube expands while convecting inwards, the ion density can thus remain close to constant while the electron density must decrease. The resulting charge separation produces an electric field that accelerates electrons inward, which can lead to the formation of beams and temperature anisotropies [Egedal et al., 2005]. In some cases, the electric field can become localized leading to the formation of double layers [Ergun et al., 2009;Wang et al., 2014;Egedal et al., 2015]. The ultimate effect of electron acceleration along the separatrices remains disputed. For example, Bessho et al. [2015] found that the separatrix acceleration occurring at the inbound leg of an electron trajectory was mostly negated by the decelerating effect when the same electron arrived close to the separatrix on the opposite side of the neutral sheet. Meanwhile, Egedal et al. [2015] argued that the confining nature of the potential could lead to more efficient energization within the exhaust by the reconnection electric field. Furthermore, as mentioned above, an initially higher electron velocity would also lead to more efficient Fermi acceleration.
On a more local scale, electromagnetic waves can mediate energy transfer between different plasma populations. For example, the counter-streaming hot and cold electron populations occurring at reconnection separatrices has been studied extensively with numerical simulations. They have been shown to be unstable to the generation of electrostatic waves, leading to the thermalization of the cold electron beam [Divin et al., 2012;Fujimoto, 2014;Huang et al., 2014;Egedal et al., 2015;Chen et al., 2015]. Depending on the velocity at which the waves are generated, they can interact with different parts of the electron distributions [Omura et al., 1996;Graham et al., 2015], and transfer energy between them. In the nonlinear stages of instabilities, it is common that electron trapping by the strong wave potential leads to electron phase space holes and electrostatic solitary waves (ESWs) [e.g. Mozer et al., 2018]. At reconnection separatrices, the interface between the inflowing and outflowing electrons represents a velocity shear. In such an environment, the instabilities developing may lead to transfer of energy not only between different energy ranges, but also across the boundary [Hesse et al., 2018a]. Although electrostatic waves and ESWs are commonly observed at reconnection separatrices in conjunction with electron beams [Cattell et al., 2005;Viberg et al., 2013] or plateau distributions associated with significant drift speed [Ergun et al., 1998], their effect on plasma populations has not been firmly established.
In order to determine the importance of the separatrix acceleration and subsequent waveparticle interaction for the overall electron energization during magnetic reconnection, it is necessary to observe these phenomena in space. In this study we will do so, using highcadence plasma and field measurements by the four closely separated Magnetospheric Multiscale Mission (MMS) spacecraft. We are able to make detailed measurements of both the electron acceleration and subsequent wave-particle interaction at separatrix regions in the magnetotail.

Observations
In this section, we report MMS observations from the plasma sheet boundary layer. The electric field is from the Electric field Double Probes (EDP) [Lindqvist et al., 2014;Ergun et al., 2014], the magnetic field is from the FluxGate Magnetometer (FGM) [Russell et al., 2014], the plasma distributions and moments are from the Fast Plasma Investigation (FPI) [Pollock et al., 2016]. All times are given in Coordinated Universal Time (UTC). Unless otherwise stated, positions and vectors are given in Geocentric Solar Ecliptic (GSE) coordinates.
We will first present observations of relatively thin channels of electron jets directed opposite to the broader ion and electron flows that is the exhaust flow of magnetic reconnection. For a few events, we will quantify the level of acceleration and compare it to the thermal energy of the lobe population and the plasma sheet population. We will then investigate the wave activity within these regions of accelerated electrons to infer the wave-particle interaction. We focus on electrostatic waves that propagate predominantly along the ambient magnetic field.

Electron acceleration channel
In this section we present an event from the magnetotail observed on June 7, 2018 at [−18, 4, 2] Earth radii, investigated previously by Huang et al. [2019]. Figure 1a shows the electron differential energy flux (DEF), in which we can identify the lobe at lower energies (E e 1 keV), and the plasma sheet at higher energies (E e 1 keV). The spacecraft are initially located in the southern lobe before they enter the plasma sheet boundary layer and the outer edges of the plasma sheet. The spacecraft make partial exits into the lobe two more times before residing in the plasma sheet until the end of the displayed time interval. During this time, the magnetic field is predominantly tailward (B x < 0). However, at 00:54:20 when the spacecraft encounters the plasma sheet for the first time, a significant northward component (B z > 0) appears, closely associated with changes in B x and B y . Huang et al. [2019] interpreted this as a passing flux rope. The ion flow is Earthward (v ix > 0), indicating the spacecraft are located in the Earthward exhaust of a magnetic reconnection X line. We note that although v ix maximizes at ∼ 800 km/s during the shown interval, the ion distribution consists of two populations: one cold population with bulk speed close to v ix = 0 km/s, and another hotter population streaming Earthward at speeds > 1000 km/s (not shown). At later times (not shown), the ion flow reverses, indicating that the X line is moving tailward. The Earthward ion flow is matched by an Earthward electron flow (v ex > 0). However, at the edges of the Earthward flow, three shorter intervals of larger amplitude tailward (v ex < 0) flows are observed. The electron flows are consistent with the current derived from the magnetic field (not shown). These regions are associated with density cavities where n e ∼ 0.01 cm −3 (Figure 1e). In comparison, the density in the lobes (before 00:54:10) is n lb e ∼ 0.05 cm −3 , and the largest density during the interval, which we associate with a plasma sheet population, is n sh e ∼ 0.15 cm −3 . We shall henceforth refer to the regions of low densities and large amplitude electron flows as acceleration channels.  Figure 1f shows the phase space density (PSD) of the reduced electron distribution projected onto the magnetic field: Inside the three acceleration channels, we can clearly see the accelerated population for v e > 0. These acceleration channels with electron flows directed towards the magnetic reconnection X line, opposite to the exhaust flow, are prominent features seen in the separatrix regions in numerical simulations of magnetic reconnection [e.g. Divin et al., 2012;Fujimoto, 2014;Egedal et al., 2015;Hesse et al., 2018a]. The reduced distributions presented in the simulations bear strong resemblance to the reduced electron distributions observed by MMS presented here.
All the acceleration channels occur at the edges of the reconnection outflow, where we expect the separatrices to be. Therefore, consistent with numerical simulations, we identify the acceleration channels as being located at the separatrices of a magnetic reconnection site. However, from numerical simulations, we know that density cavities associated with the accelerated populations do not always extend over the whole length of the separatrices [e.g. Egedal et al., 2015]. In addition, in the presence of a guide field the acceleration regions are to some degree suppressed at two opposing of the four separatrices [e.g. Pritchett and Coroniti, 2004]. Therefore, while accelerations channels and reconnection separatrices are closely related, they are not always coincident. Also, importantly, electron acceleration channels are not exclusively related to magnetic reconnection but can occur in a multitude of plasma environments.
Since we do not observe any reconnection outflow reversals in the immediate proximity of the acceleration channels, we have no straightforward means of determining how far away from the X line they are located.
The acceleration channels are associated with large amplitude parallel (Figure 1g-1h) and perpendicular (shown for a shorter time interval in Figure 5) electric field fluctuations. Since instabilities driven by parallel beams often result in large amplitude parallel electric fields, we will in sections 5 and 6 focus on investigating the relation between the field-aligned electric fields and the parallel streaming populations. First, however, we will quantify in more detail the electron acceleration.

Acceleration potential
To obtain a quantitative estimate of the acceleration potential that the electrons have passed through, we investigate the reduced electron distribution in more detail. Figure 2 again shows the reduced electron distribution, now for a slightly shorter time interval. The thinner black line shows the parallel electron bulk speed v bulk e, | | , and the thin dashed line shows v bulk e | | ± v te, | | , where v te, | | = 2k B T e, | | /m e is the electron thermal speed based on the parallel temperature.
To estimate the acceleration potential we use the Liouville approach and assume that a lobe population initially at rest has been accelerated to the energy corresponding to the energy where the maximum phase space density is observed, i.e. eψ = m e v 2 acc /2, where v acc = v e, | | ( f e = f local e,max ). To avoid picking up small variations in f e we have first applied a running average over three full 3-D distributions (the averaged distribution in Figure 2 can be compared to the original distribution in Figures 1f or 6a), and thereafter again averaged v acc and the corresponding ψ over three time steps. The obtained speeds v acc are shown as thick solid black lines in Figure 2. We note that the phase space density between the lobe and the acceleration channels is not completely conserved, f lb e > f acc e . This indicates that non-adiabatic processes are at work, for example wave-particle interaction. Although beam thermalization through wave-particle interaction will decrease the average drift velocity of the beam population, it can initially tend to shift the peak phase space density to higher energies [see e.g. Figure 2 in Che et al., 2009]. When we find v acc , we therefore exclude instances where the beam is thermalized beyond a certain threshold, even if it is possible to find a f local e,max . This is for example the case during the last part of the first acceleration channel. We will discuss the implications of this later.
The accelerated populations have larger speeds than the moments calculated from the entire . This is due to the presence of an additional electron population close to v e, = 0 inside the acceleration channels. Inside the acceleration channels, the temperature (as indicated by the distance between the two dashed lines marking v e, ± v te, ) is increased relative to the lobes. This initial jump in temperature is also largely due to the presence of the additional population at v e, = 0. The population with low parallel speed could be due to wave-particle interactions, or leakage from the plasma sheet or even the lobes.
The maximum potentials associated with v acc for the three acceleration channels are ψ = [1800, 2400, 1400] V, respectively (also written on top of Figure 2). In comparison, characteristic temperatures in the lobe and plasma sheet are T lb e = 220 eV and T sh e = 3700 eV (obtained from the blue and red intervals shown in Figure 1). In terms of these characteristic energies, eψ ≈ [8, 11, 6] T lb e ≈ [0.5 0.6, 0.4] T sh e . Note that the thermal energy of the plasma sheet is usually larger closer to the neutral sheet [Baumjohann et al., 1989]. Because |B| > 15 nT during the entire shown interval, the spacecraft stay relatively close to the plasma sheet boundary layer, and therefore T sh e should be considered as a lower bound.
We have performed the same analysis as described above for a few other acceleration channels from a total of five burst intervals during two days in July, 2017, listed in Table 1. The results are shown in Figure 3 (blue circles). As before, all events we have included from MMS show clear features of a cold population being accelerated through a potential drop. We have not included flat-top distributions, or events in which the entire electron beam has likely been thermalized already. However, many events show evidence of some degree of thermalization. We show the results as a function of electron beta in the lobe β lb e . Note that the acceleration channels that are from the same burst interval can correspond to the same lobe and/or plasma sheet intervals. It is also possible that acceleration channels that are observed in close proximity to each other can be channels that are crossed multiple times. For these MMS events we find that ψ = 1 − 8 keV, eψ/T lb e = 1 − 15 (with one value at 35), and eψ/T sh e = 0.1 − 1.7. However for the last seven acceleration channels in Table 1, T sh e is likely underestimated, as MMS only skirted the plasma sheet boundary layer. Regardless of this, electrons passing these acceleration channels have already reached a substantial fraction Figure 3. Summary of acceleration potential for a few electron acceleration channels, compared to previous events observed by Cluster [Borg et al., 2012;Egedal et al., 2015]. Cluster results are adapted from Table 1 in Egedal et al. [2015]. (a) ψ is inversely proportional to β lb e , and (b) many times the electron thermal energies per charge in the lobe T lb e /e. (c) eψ is comparable to plasma sheet thermal energies T sh e . Note that in events where the spacecraft only stay at the edge of the boundary layer, T sh e is likely underestimated, see Table 1. of their final energy before entering the magnetic reconnection exhaust proper. In agreement with previous results obtained by Cluster [Borg et al., 2012;Egedal et al., 2015] (red circles), the acceleration potentials show an inverse dependence on β lb e . We note that these two studies seem to cover different ranges of β lb e . This could be due to selection bias, or instrumental differences related to the accuracy to which the densities and temperatures can be determined. Another statistical study of electron distributions in magnetic reconnection regions by Cluster found electrons beams directed towards the X line with an occurrence frequency of about 20% in the off-equatorial region [Asano et al., 2008]. The beams had energies of 4-10 keV.
of their events the spacecraft observed electron beams propagating inward towards the X line. They found that the beams had a higher occurrence frequency Similar to previous observational studies of both dayside [e.g. Vaivads et al., 2004] and nightside [e.g. Lu et al., 2010;Wang et al., 2012] magnetic reconnection, the acceleration channels we study in this paper are associated with density cavities. Figure 4a shows the relation between the lobe densities n lb e and the minimum densities inside the acceleration channels n sep e . For all events, n sep e < n lb e . We do not show it here, but for all the events, the densities on the plasma sheet side of the acceleration channels are larger than both n lb e and n sep e . Figure 4b shows that the ratio of densities between the lobes and the acceleration channels, n sep e /n lb e , becomes smaller with increasing acceleration potential ψ. The decrease in density between the lobe and the acceleration channels is expected from the conservation of phase space density of an accelerated plasma population [e.g. Schamel, 1982]. The existence of density cavities at the separatrices is also in agreement with numerical simulations of symmetric antiparallel [e.g. Shay et al., 2001;Lu et al., 2010;Egedal et al., 2015], and guidefield [e.g. Pritchett and Coroniti, 2004] magnetic reconnection. In the case of guide-field reconnection, the reconnection electric field has a component parallel to the magnetic field. This results in enhanced parallel acceleration by the reconnection electric field, and larger density cavities at two opposing of the four separatrices. In this study, we have not differentiated between strictly antiparallel and guide-field reconnection. However, since all the events are from the tail, it is likely that any guide field, if present, is low or moderate.

Width of acceleration region
Before we continue, we shall make a rough estimate of the width of the first acceleration channel in Figure 2, which had ψ = 1800 V. For this event, it is not possible to reliably determine the spacecraft trajectory relative to the boundary layer from timing analysis. We  Table 1. Acceleration potential ψ obtained from a total of five burst intervals in July 2017. T sh e is chosen as the largest temperature observed in the proximity of the acceleration channel. In events where the spacecraft only stay at the edge of the boundary layer, T sh e is underestimated, this is likely the case for the seven last events. therefore take a different approach using the perpendicular electric field. We make the following assumptions: (1) The acceleration potential is electrostatic in nature, and the parallel potential drop is therefore accompanied by a perpendicular potential drop. This is consistent with the divergent electric field E ⊥,z centered around the electron flow shown in Figures 5a and 5b. We show the original field and the field downsampled to 3 Hz, to highlight the DC variations. We can determine that the field is divergent because MMS cross the southern separatrix from the lobe to the plasma sheet and observes a negative-positive polarity of E ⊥,z . A divergent electric field is associated with a positive electrostatic potential, consistent with the acceleration of electrons in towards the X line. This can also be seen in numerical simulations of magnetic reconnection [e.g. Figures 4b-4c in Divin et al., 2012]. However, some simulations also show that the electron acceleration along separatrices are in part due to inductive electric fields [Egedal et al., 2015;Bessho et al., 2015]. Figure 5c shows the obtained potential ψ at the original cadence and downsampled to 3 Hz, like the electric field. (2) The perpendicular profile of the acceleration channel is Gaussian: ψ ⊥ (z) = ψ 0 exp(−z 2 /2l 2 z ), where z ∼ z GSE is the coordinate perpendicular to both B and the main electron flow v e , and ψ 0 = ψ(z = 0) is the potential in the center of the acceleration channel. While the Gaussian shape is somewhat arbitrary, we have no way to better determine the exact shape. The perpendicular electric field associated with this potential structure has peaks values |E max z | = l −1 z ψ 0 exp(−1/2) at z = ±l z , where the potential is ψ(z = ±l z ) = ψ 0 exp(−1/2). The half width is thus given by (1) We now choose two time steps from where the electric field is the strongest, marked by yellow squares in Figures 5b and 5c. At these times, because the observed electric field maximizes here, the spacecraft are presumably located at an intermediate distance from the center of the acceleration channel, close to l z . For the two points |E max z | = [25, 19] mV/m, and ψ(z = ±l z ) = [700, 1300] V, giving estimated half widths l z = [25, 65] km, respectively. The estimated thickness of the acceleration channel is thus L = 2l z = 50 − 130 km. In comparison, the ion and electron thermal gyroradii ranges between 100 − 400 km, and 2 − 8 km, respectively, where the smaller (larger) values are taken at the lobe (sheet) side of the acceleration channel. The Debye length is ∼ 0.5 km in the lobes and reaches 4 km inside the acceleration channel.

Wave activity
Inside, and in the vicinity of the acceleration channels, large amplitude parallel electric fields are typically observed. An example is shown in Figure 1h, where the largest amplitude fields form bipolar pulses, often termed electrostatic solitary waves (ESW) (see also Figure 5 in Huang et al. [2019]).
To quantify to what extent the electric field can affect the electrons and modify their velocity, it is helpful to look at the sum of the kinetic and potential electron energy in the frame of the wave traveling at speed v ph : which is a constant of motion. If U < 0, the electron is following a trapped trajectory, and if U > 0, the electron is following a passing trajectory [Bernstein et al., 1957]. The electrons can transition from passing to trapped trajectories (or vice versa) if the wave field is growing (or decaying) -the electrons become trapped (or released). The limiting velocity that separates trapped and passing electron trajectories at the point where the potential is the largest, called the trapping velocity. The trapping velocity is a good indication of what part of the electron distribution is likely to interact efficiently with the wave. To find the trapping velocities, we need to find the propagation velocities v ph , and electrostatic potentials φ of the waves. When the same wave structure is observed by two or more of the spacecraft, we can perform interferometry measurements to obtain the propagating velocity. That is, we measure the delay between the times the structure is observed by the different spacecraft, and compare it to the spacecraft separations. This is possible because the spacecraft separation of about 10 km is comparable to the typical length scale of the wave forms ∼ 10λ De , where λ De = ( 0 k B T e /ne 2 ) 1/2 is the Debye length [Graham et al., 2016]. Inside the acceleration channel λ De = 3.3 km (using T e ∼ 2000 eV and n e = 0.01 cm −3 ). We find that v ph is loosely proportional to the velocity of the drifting electron population. Figure 6b shows the normalized power spectrum of E as a function of frequency f and parallel wavenumber k obtained from four spacecraft interferometry for the time interval shown in Figure 6a. The method is described for two points of measurement by Graham et al. [2016], but here generalized to four points. This removes the need to assume a given propagation direction. The resolvable k 's are related to the inter spacecraft separation as k ,max = π/max(∆l i j, ), where ∆l i j, is the distance between the individual spacecraft pairs (denoted by indices i and j) parallel to the ambient magnetic field. In our case max(∆l i j, ) = ∆l 14, = 25 km, giving |k | max ≈ 0.125 km −1 . The lower power found at k 0.05 km −1 and f 0.4 kHz might be due to spatial aliasing. The black lines mark the phase velocity of the individual ESWs (as shown in Figure 6a). While the general trend is that the the phase velocities increase towards the plasma sheet, it is possible to roughly divide the ESWs into two groups based on their speeds: one earlier group with lower speeds and lower f and one later group with higher speeds and higher f .
Using v ph , we can obtain the distances between the positive and negative peaks of the bipolar electric fields, which we call the peak-to-peak length scale. We find that they vary be-  tween l pp = 20 km and l pp = 120 km with a mean value of l pp = 57 km (Figure 7). The wavenumbers corresponding to wavelengths λ = 2l pp of the individual ESWs, k • ∼ π/l pp , are marked by circles in Figure 6. For some of the ESWs, the electric field perpendicular to the ambient magnetic field was comparable or even larger than the parallel field (not shown). This suggest the perpendicular length scales can be comparable to the parallel length scales [Franz et al., 2000]. A large range of estimated l pp 's are comparable to the estimated width of the acceleration channel, which was L ≈ 50 − 130 km.
The potential of the waves along the trajectory of the spacecraft are calculated by integrating the parallel electric field, using the parallel component of the measured phase velocity, dl = −v ph dt: The maximum electrostatic potential of each ESW along the spacecraft trajectory φ max varies between φ max = 300 V and φ max = 4000 V with a mean value of φ max = 1500 V (Figure 7). The corresponding trapping speeds v tr (Eq. 3) are shown in Figure 6a as two black lines bounding v ph . At earlier times where v ph < 5 × 10 3 km/s, the trapping range encompasses significant parts of the beam, indicating favorable conditions for strong waveelectron interaction. At later times where 10 × 10 3 < v ph < 35 × 10 3 km/s, the beam is not as apparent and has likely become significantly thermalized, also indicative of strong waveparticle interaction.
Note that at the later part of the interval, where v ph are larger, the beam has become significantly thermalized. Since we did not take into account significantly thermalized beams when determining ψ, the maximum acceleration speed obtained for this acceleration channel (v acc ∼ 25 × 10 3 km/s, corresponding to ψ = 1800 V), is smaller than the observed phase velocities, v acc < v ph .
We have performed this analysis for a number of events, and present two more of them in Figure 8. We observe both similarities and differences between the different cases. For example, in Figure 8a, although there is an asymmetry in f e (v ) close to the lobe towards the end of the interval, the electron bulk velocity is close to zero, and no distinct beam is observed. However, the phase velocities are proportional to the energy at which the phase space density begins to decrease rapidly (we refer to this energy as the shoulder energy). This might indicate that the beam has already been destroyed by the wave-electron interaction, and that what we observe is the thermalized beam. This is supported by the fact that the trapping range v ph ± v tr covers a large part of the electron distribution. Because no distinct beam is observed, this time interval is not included in Table 1 or Figure 3. Figure 8c shows an acceleration channel with a distinct beam, which is included in Table 1 and Figure 3. Again we find that the phase velocities are proportional to the beam speed and that the wave interaction range covers the beam. Figures 8b and 8d show that the phase velocities obtained from timing analysis of individual ESWs correspond well to the maximum wave power in the dispersion relation obtained form four spacecraft spectral analysis. In Figure 8d, we can see clear indications of spatial aliasing.

Spatiotemporal evolution and instability analysis
To investigate whether the observed plasma distributions can account for the generation of the ESWs, we solve the unmagnetized, electrostatic dispersion equation: for the event presented in Figure 6. Here, Z is the plasma dispersion function, and 's' denotes the different plasma populations. Based on the observed ESW characteristics, the ESWs in Figure 6 can be divided roughly into two groups. One group with lower v ph observed at earlier times and one group with larger v ph observed at later times. We will therefore investigate different combinations of plasma populations.
We expect the accelerated lobe electron population to be the main driver of the instability. This beam, however, can be in different stages of evolution. Further downstream of the acceleration channel, observations and simulations show that the peak phase space density is shifted toward larger speeds, but also that the beam is weaker, i.e. it has a lower density relative to the electron population at lower speeds (compare the blue stars and purple circles in Figure 9a). To the zeroth order, the shift to larger speed is due to the acceleration. However, the wave-particle interaction can also contribute to this as follows: electron trapping removes phase space density from the lower speed edge of the beam, decreasing its density and at the same time shifting the beam peak to higher speeds [e.g. Che et al., 2009]. Recall that when acquiring the acceleration potential ψ in section 4, we did not include significantly thermalized beams. The distribution shown with purple circles in Figure 9 is an example of such a distribution that was considered too thermalized. Note that a beam drift speed of 35000 km/s would correspond to an acceleration potential of 3500 eV. The electron population at lower speeds can be plasma sheet electrons that enter the acceleration channel during part of their gyromotion. It can also be the trapped electrons that were originally part of the beam. A large part of the low-velocity electrons are within the trapping range of the ESWs. The ions can be both cold lobe and hotter plasma sheet populations. However, the ion thermal speeds of both lobe and plasma sheet ions are both low in comparison to electron and phase speeds, so in the dispersion analysis we will only consider a single (medium hot) ion population with a temperature of T i = 5 keV, corresponding to a thermal speed of v ti = 980 km/s. An ion temperature of 10 keV would correspond to a thermal speed of 1380 km/s, which is not a significant difference considering the phase speeds and electron speeds.
We consider two different electron distributions, based on observed distributions at slightly earlier (blue) and later (purple) times where the observed v ph are slower and faster, respectively (Figure 9a)  Figure 9b. For both cases we obtain good matches to the real frequencies f r , while the maximum growth for the faster ESWs are shifted towards higher k 's. The growth rate for the faster beam is about ten times larger than that for the slower beam. The phase speeds at maximum growth rate are v slow ph = 2500 km/s and v f ast ph = 27000 km/s, respectively, shown as vertical dashed-dotted lines in Figure 9a. For the slower ESWs, depending on small variations of the input parameters, either the ionelectron or electron-electron modes are dominating. The ion-electron mode is essentially a Buneman type mode with a hot electron background that does not interact with the drifting electron population [Norgren et al., 2015]. The slow electron-electron mode is an electronacoustic wave. The faster ESWs are generated by an electron beam-mode instability with close to constant phase speed regardless of wavenumber. The evolution of instabilities from Buneman to electron beam-mode is similar to what was described by Che et al. [2009] for For the faster ESWs, the peak growth rate is shifted towards larger k , or equivalently smaller wavelengths. guide field reconnection. Although the electron beams we study here are not located inside the EDR, the local dynamics can be the same.
Since we observe large amplitude highly nonlinear localized structures, we are not observing the waves in the linear stage of instability. The predicted growth is therefore not necessarily expected to coincide with the range of k 's where the linear growth rate maximizes. Although some simulations do show good correspondence between observed wave characteristics and linear instability growth rates [e.g. Fujimoto, 2006;Chen et al., 2015], ESWs are for example known to merge with each other and grow in size [e.g. Mottez et al., 1997], which would correspond to a shift to smaller k 's. If we would assume that the wave power at the initial stages of the beam-mode instability peaked at 0.15 km −1 but due to coalescence has shifted to 0.05 km −1 , this would correspond to a change from a wavelength λ = 40 km to a corresponding length scale of the ESWs of 120 km. Another effect that may play a role in modifying the wave growth is the velocity shear and the perpendicular structure of the flow channel inside which they grow [Che et al., 2011]. The waves can not be considered as plane waves with infinite extent in the perpendicular direction, which Eq. 4 assumes. However, extending the wave analysis to include these effects is beyond the scope of this paper. ESW's are also known to be limited in size by transverse instabilities [Muschietti et al., 2000;Graham et al., 2016]. If the bounce frequency of electrons trapped in the potential well of the ESWs exceeds the gyrofrequency, the ESWs tend to become unstable and dissipate. The bounce frequency for a given potential increases with decreasing length scale, or correspondingly larger k's. We have examined this relation here (not shown) and find that the limiting k's are about three times larger than the k's where the wave power maximizes. We therefore do not think this is the deciding factor in determining the range of observed length scales.

Discussion
We have investigated the electron acceleration and wave-particle interaction in acceleration channels located at magnetic reconnection separatrices in the magnetotail. Generally, we found that the lobe populations were gradually accelerated up to a significant fraction of the thermal energies in the plasma sheet. Nonlinear ESWs that were observed at the same time as the accelerated populations had large enough potentials to interact with a significant part of the electron distribution, including the beam. Here we will discuss how the wave-particle interaction and spatial effects are expected to alter the beam and thereby our estimate of the acceleration potential ψ. We will also discuss how the wave properties can be related to the thermalized beam.
The beam represents the free energy of the system. When the waves have grown to amplitudes such that they are able to trap the beam, the beam can become thermalized, and the free energy source is removed. The wave trapping range at saturated wave growth is therefore expected to encompass the beam responsible for the instability growth. This is consistent with our observation: when a distinct beam could be observed, such as seen at earlier times in Figure 6a or in Figure 8c, the wave trapping range encompassed a large part or the entirety of the beam: The waves observed at the same time as the beams were highly nonlinear solitary structures, indicating they could be close to saturation.
When we estimated the acceleration potential ψ, we studied the beam component of the reduced electron distribution. However, the wave-particle interaction can affect the beam and thereby affect our estimate of ψ. For example, electrons will locally be accelerated to higher energies in the potential φ of the wave. Since the sampling rate of FPI is lower than the time scale of an individual ESW, the beam will appear spread out in velocity space, and the peak of the spread out beam energy will appear at higher energies than what was accounted for by the accelerating potential ψ. The waves can also eventually trap electrons, and thus gradually remove electrons from the lower speed edge of the beam. This will further spread out the beam in velocity space and shift the peak phase space density of the beam to higher energies. However, these thermalization effects occur gradually, and can have affected the beam significantly even before the beam has reached its final energy as governed by the larger scale potential ψ.
For the case presented in Figure 6a, we could clearly see that the closer to the plasma sheet the accelerated population were observed, the weaker, or more thermalized, the beam became. In the region with higher velocity ESWs, it is hard to perceive the beam at all. Since we did not take into account significantly thermalized beams when determining ψ, the maximum acceleration speed obtained for this acceleration channel (v acc ∼ 25 × 10 3 km/s, corresponding to ψ = 1800 V) was smaller than many of the observed phase velocities, v acc < v ph . However, due to the nature of the generating instabilities, where the phase speeds fall in between the two populations, it is reasonable to assume that the speed of the unthermalized beam was larger than v ph as shown in the inequality in Eq. 5. For example, for the instability analysis of the faster ESWs in Section 6, the drifting component of the electron distribution had a speed of 35 × 10 3 km/s. Therefore, although the beam-wave interaction can make the beam appear at higher energies than what is accounted for by ψ, the observed phase speeds should give a minimum value. This is an indication that the maximum potential drop for this acceleration channel was in fact larger than what we estimated in Figure 2. In Table 1 and Figure 3 ψ could therefore be considered as lower bounds.
Like mentioned above, for the beam in Figure 6a, the beam thermalization became more prominent when the beam had already been accelerated to larger energies. The point in time when the beam seemed to become significantly thermalized coincided roughly with the appearance of high velocity ESWs. It is possible that the increased thermalization could be related to the generating instabilities. The electron beam-mode instability, which we found could be responsible for the wave generation in the large velocity region had a growth rate roughly ten times larger than the instabilities that were active when the beam had lower speeds.
To some extent, the increased thermalization should also be due to the integrated effect of wave-particle interactions along the acceleration channel. That is, the further down the acceleration channel the beam has progressed, the longer distance the wave-particle interaction has had the time to affect the beam.
In other cases where no distinct beam could be observed (we presented one such case in Figure 8a), the phase speeds were instead proportional to the energy at which the phase space density started to decrease rapidly. Due to the similar wave behaviour between the cases with and without distinct beams, it is likely that the presumed beam in the latter case had already been completely thermalized. In such cases where no beam can be identified, but the wave characteristics can be obtained, it could be possible to use the inequality in Eq. 5 to indirectly estimate the beam characteristics and therefore the amount of electron acceleration as well as thermalization. Inferring the characteristics of the generating electron beam after it has already been thermalized has been done in previous studies [e.g. Norgren et al., 2015].
Here, however, we have here tracked the continuous change in the beam, accompanied by the continuous change in wave characteristics, certifying the credibility of such an approach.
We also consider the case where the weakened beam could be partially due to spatial effects perpendicular to the magnetic field. For example, if we assume that the initial lobe population is isotropic, electrons with pitchangles close to 90 • would cover a larger perpendicular distance throughout their gyromotion than electrons with pitchangles closer to 0 • or 180 • . The beam would therefore be weaker at the edges of the acceleration channel than at the center. However, the gyroradii of lobe electrons, based on T lb e = 220 eV and B lb = 20 nT is ρ e = 2.5 km, which is significantly smaller than the estimated thickness of the acceleration channel L ≈ 50 − 130 km. For this gyration effect to be important, it is likely that some prior heating and pitch angle scattering would have to had taken place.

Conclusions
In this study we investigated the electron acceleration and subsequent wave generation and wave-particle interaction at magnetic reconnection separatrices in the magnetotail. We summarize our conclusions below.
• Adjacent to the reconnection exhaust, we found relatively thin regions of electron lobe populations accelerated towards the X line. The electrons were accelerated to energies of 1-7 keV, several times the thermal energies within the lobe, and a significant fraction of the thermal energies inside the outflow. • All acceleration regions were associated with density cavities. The decrease in densities from the lobes to the acceleration channels increased with increasing acceleration potential, consistent with theoretical predictions. • For two acceleration channels presented in more detail, we could observe how the lobe populations were gradually accelerated. For one of them, the resulting beam became significantly weaker closer to the plasma sheet. • Electrostatic solitary waves observed in the acceleration regions had phase speeds proportional to the beam speeds. The potentials of the waves were large enough such that the waves could interact efficiently with a large part of the electron population, including the beam. This indicates that the waves play an important role in controlling the evolution of the beam, aiding to thermalize it. • For one acceleration channel we investigated the instability of the evolving electron distribution and found that it could account for the observed wave properties. When the beam had been accelerated to moderate speeds, it was unstable to a combination of competing Buneman and electron acoustic instabilities, generating waves at low phase speeds. When the beam had been accelerated to larger speeds and had become weaker, the distribution was unstable to an electron beam-mode instability, generating waves at larger speeds.
For the events that we have examined, these wave-particle interaction ranges presented above represent the typical conditions when ESWs are observed. This shows the efficiency with which waves are able to turn directed drift energy into thermal energy of the plasma. These results are not only applicable to magnetic reconnection, but to any process in which superthermal electron beams form.