Ion reflection and acceleration near magnetotail dipolarization fronts associated with magnetic reconnection

Dipolarization fronts (DFs) are often associated with the leading edge of earthward bursty bulk flows in the magnetotail plasma sheet. Here multispacecraft Time History of Events and Macroscale Interactions during Substorms (THEMIS) observations are used to show that a spatially limited region of counterpropagating ion beams, whose existence is not evident in either the plasma moments or the electric field, is observed on the low‐density side of DFs. The THEMIS magnetic field data are used to establish appropriate comparison cuts through a particle‐in‐cell simulation of reconnection, and very good agreement is found between the observed and simulated ion distributions on both sides of the DF. Self‐consistent back tracing shows that the ion beams originate from the thermal component of the preexisting high‐density plasma into which the DF is propagating; they do not originate from the inflow region in the traditional sense. Forward tracing shows that some of these ions can subsequently overtake the DF and pass back into the high‐density preexisting plasma sheet with an order‐of‐magnitude increase in energy; this process is distinct from other ion reflection processes that occur directly at the DF. The interaction of the reconnection jet with the preexisting plasma sheet therefore occurs over a macroscopic region, rather than simply being limited to the thin DF interface. A more general consequence of this study is the conclusion that reconnection jets are not simply fed by plasma inflow across the separatrices but are also fed by plasma from the region into which the jet is propagating; the implications of this finding are discussed.

preexisting plasma observed prior to the DF arriving at the satellite [Zhou et al., 2011]. At higher latitudes near the plasma sheet boundary layer, the reflected ions can stream along field lines and are associated with transient brightening in the proton aurora Zhou et al., 2012]. If the DF is not planar, then this reflection process may cause dawn-dusk asymmetry [Zhou et al., 2014]. An important feature of this model is that the ions gyrate (and are reflected) solely in the plane of the current sheet (i.e., x GSM -y GSM ). However, Birn et al. [2013] examined the behavior of ion test particles in a 3-D MHD simulation of magnetotail reconnection and did not find the ion reflection reported by Zhou et al. [2010], pointing out that Zhou et al. assumed both 2-D structures and that the DF propagates at constant speed.
While the test particle approach can reveal important physics, it does not include self-consistent feedback, and so for a more complete treatment of the problem, a particle-in-cell (PIC) or other kinetic approach is required. Using PIC simulations, Wu and Shay [2012] showed that reflected ions are observed on the high-density side of the DF but did not examine the reflection mechanism in more detail. Hybrid simulations, which also capture ion kinetic physics, were used to explore the behavior of ions at the leading edge of reconnection jets, corresponding to the low-density side of the DF, reporting that counterpropagating beams sourced from the lobe are present, but this work focused on the physics of reconnection jet formation rather than the DF interface itself [Fujimoto et al., 1996;Nagai et al., 2002; M. S. Nakamura et al., 1998]. Overview of the THEMIS satellite data. THEMIS P1 measurements of (a) ion energy flux, (b) ion velocity, (c) ion density, (d) ion temperature, (e) electron temperature, (f) magnetic field components, (g) magnetic field strength, and (h) specific entropy are shown. (i-p) THEMIS P2 measurements in the same format are also presented. In each case, 7 min of data, centered on the DF, is shown. Note that the time interval shown for each satellite is therefore different.
To better understand the relationship between DFs and reconnection and more specifically the nature of the ion kinetic physics in the vicinity of the DF, here we examine and compare multipoint experimental observations of a DF with a PIC simulation of reconnection onset. This particular DF was observed by Time History of Events and Macroscale Interactions during Substorms (THEMIS) [Angelopoulos, 2008] on 27 February 2009 and has been extensively documented [e.g., Deng et al., 2010;Ge et al., 2012;Runov et al., 2011aRunov et al., , 2009, meaning that it can be regarded as a clean and uncontroversial example, which, in the first instance, provides a good baseline comparison with the simulation. We concentrate on the THEMIS P1 (THEMIS-B) and THEMIS P2 (THEMIS-C) data, since our goal is to examine the structure of the DF before it reaches the dipolar magnetotail. This ensures that the physics of the event is not strongly influenced by processes in the flow braking region  as the Earth's dipolar field is not modeled in the simulation.

Overview of Spacecraft Observations
In this study data from the fluxgate magnetometers [Auster et al., 2008], the electrostatic analyzers (ESAs) [McFadden et al., 2008], the electric field instruments [Bonnell et al., 2008], and the solid state telescopes (SSTs) [Angelopoulos, 2008] are used. Spacecraft data are shown in the GSM coordinate system. During the time interval under study, both probes were operating in burst mode and returning a full 3-D plasma distribution every spacecraft spin (3 s). The ion moments were calculated by combining the burst mode ESA distributions with the full mode SST data so as to more accurately measure the temperature.
An overview of the data is presented in Figure 1. The DF, corresponding to the increase in B zGSM (at the time indicated by the vertical black lines), was observed by the THEMIS P1 satellite at t P1 = 07:51:26 UT and x P1 = À20.1 R E and by the THEMIS P2 satellite at t P2 = 07:52:35 UT and x P2 = À16.7 R E . Prior to the arrival of the DF, P1 was located in the northern hemisphere of the magnetotail (B xGSM > 0) and P2 in the southern hemisphere (B xGSM < 0); both observed relatively high density plasma, indicating that they were located in the plasma sheet, and the overall conditions were relatively stable, with T i~3 keV. Detailed analysis of the magnetic field measured by both satellites prior to the arrival of the DF indicates that the current sheet is well aligned with the GSM coordinate system (to within 10°) as is often observed in this region of the magnetotail [Eastwood et al., 2010]. The current sheet normal points in the z GSM direction, and the guide field is very small in this event; if present, it is on the order of 1 nT in the Ày GSM direction.
The ion density dropped as the DF arrived at each satellite, and the earthward ion flow speed subsequently increased to its maximum value over a period of about a minute. At both satellites, this roughly corresponds to the interval over which B zGSM (and therefore |B|) was enhanced. Both the ions and electrons were hotter in the low-density region after the DF arrived, with the ions exhibiting a larger increase in the perpendicular temperature. In the case of P2, the ion temperature subsequently became isotropic, whereas at P1, the parallel ion temperature became larger. This difference is associated with the different cuts that the satellites made through the structure. After 07:52:12 UT, B xGSM increased at P1, and so P1 moved toward the edge of the jet (for example due to current sheet flapping or some other relative motion), at which point, the parallel ion temperature was observed to dominate. In contrast, at P2, B xGSM remained relatively small, and so more isotropic ion temperatures were observed. The electrons were in both cases observed to be more isotropic, with some perpendicular heating seen only at P2 immediately adjacent to the DF. In general, the similarities between the two encounters immediately indicate the extent to which DFs remain The P1 data (red) has been shifted by +69 s so that the increase in B zGSM occurs at the same time as P2.
Journal of Geophysical Research: Space Physics 10.1002/2014JA020516 relatively stable as they propagate earthward, as has previously been noted . The specific flux tube entropy pV 5/3 , calculated according to the method described by Wolf et al. [2006], is thought to characterize the degree to which BBFs may easily propagate earthward into the more dipolar regions of the magnetotail [Pontius and Wolf, 1990]. At both satellites, this parameter is reduced to less than 0.3 for an interval of about 40 s following the arrival of the DF.
We now examine in more detail the magnetic field time series in the vicinity of the DF. Figure 2 shows the magnetic field data from THEMIS P1 and P2 overlaid, where the THEMIS P1 data have been time shifted by +69 s. The B zGSM profiles ( Figure 2c) are indeed very similar, and the B xGSM profiles (Figure 2a) are roughly mirror images of each other from 07:51:00 to 07:52:50 UT. If the magnetic field can be used as a guide to determine the distance to the current sheet, the similarity in |B xGSM | implies that the two satellites were approximately equidistant from the B xGSM = 0 plane. The B yGSM time series are also largely similar, but we note that just prior to the arrival of the DF, B yGSM is negative at P1 and positive at P2.
The observations can be used to establish characteristic lengths, times, and speeds for this event. Prior to the arrival of the DF, at P1, <n> = 0.72 ± 0.02 cm À3 , and <B xGSM > = À8.4 nT ± 0.3 (averaged over 07:50:21-07:50:41 UT). At P2, <n> = 0.99 ± 0.03 cm À3 , and <B xGSM > = +8.4 ± 0.6 nT (averaged over 07:51:30-07:51:50 UT). Assuming that the observed plasma is in simple pressure balance with the lobe (where plasma pressure is negligible), the lobe field strength is estimated to be~25-30 nT, which corresponds to a characteristic ion cyclotron frequency of f ci = 0.42 Hz. Similarly, the preexisting current sheet density (specifically, we use n~0.8 cm À3 as the average of the values measured on opposite sides by the two spacecraft) corresponds to a characteristic ion inertial length of d i~2 55 km. Finally, a characteristic Alfvén speed based on these density and magnetic field strength values is~670 km s À1 . We therefore find that the spacecraft separation in x GSM was~3.4 R E~8 5 d i and the temporal separation was Δt~69 s~29 f ci À1 (ion cyclotron periods). The spacecraft separation in z GSM was 0.7 R E , corresponding to~17 d i and implying, given its apparent orientation, that the layer was relatively thick; the reconnection onset generating the fast earthward flow occurred further downtail in a region where the current sheet is presumed to have been thinner. The DF itself propagated between the satellites at v DF~3 14 km s À1 .

Observed Distribution Functions
During this event, the satellites were in "burst mode," meaning that full 3-D distributions were downlinked every spin (3 s). Figure 3 shows several distributions observed by THEMIS P2 during the DF encounter, together with time series of the ion density, magnetic field, and the y GSM component of the electric field (as measured and as calculated from Àv × B). The distribution functions are shown as cuts in the v-B plane, as this most clearly shows the overall structure in the distribution.
Immediately prior to the arrival of the DF, a single ion population was observed in the high-density stagnant current sheet (Figure 3a). It is asymmetric with a net flux in the +x GSM direction, corresponding to high-energy ions moving earthward from the DF. In a previous study, these ions were identified as so-called "reflected" ions and linked to enhanced proton precipitation and enhanced proton aurora through the use of global MHD simulations . Figure 3b shows the distribution observed just after the peak in B zGSM in the low-density plasma immediately adjacent to the DF. Two distinct and well-resolved populations are observed, separated in the ±v par direction but displaced in the +v perp direction. They are bisected by the black line representing x GSM . In fact, here the magnetic field is dominated by B zGSM (Figure 3f), so the parallel direction is approximately equivalent to the +z GSM direction, and the perpendicular direction is equivalent to +x GSM . The parallel speed of the two populations is several hundred km/s, a significant fraction of the characteristic Alfvén speed for this event.
These counterpropagating beams, moving north and south parallel to the magnetic field which is dominated by B zGSM , are a persistent feature and are observed for several tens of seconds after the passage of the DF, as illustrated by Figure 3c, which shows the distribution midway through the interval. The counterpropagating beams were observed at P2 until 07:53:17 UT (duration 42 s), which corresponds to the end of elevated B zGSM values, and when both the peak ion velocity was observed and the specific entropy was no longer reduced. Figure 3d shows the ion distribution function at this time: the counterpropagating beams are no longer visible, and the distribution is similar to those observed in the plasma sheet for smaller values of B x [Raj et al., 2002]. Figure 4 shows the same data for P1, and similar behavior was observed. Prior to the arrival of the DF, a single population, with an earthward flux of more energetic ions, was observed. Interpenetrating ion beams were observed after the passage of the DF until 07:52:15 UT (duration 49 s), at which time B zGSM was no longer elevated, the jet speed was maximized, the specific entropy increased, and the counterstreaming beams were no longer observed ( Figure 4d). We note that the distribution is different from that observed in Figure 3d and this is associated with the larger value of B x (as mentioned in the context of Figure 1). This difference is reflected in the fact that the parallel temperature dominates at P1 closer to the edge of the current sheet, whereas the plasma is more isotropic at P2 where B x was smaller.
The counterstreaming distributions in fact appear to be similar to those reported by Fujimoto et al. [1996] and Nagai et al. [2002] based on single-spacecraft Geotail observations. By contrast, here we now use the two spacecraft THEMIS measurements to establish spatial properties and examine their context in terms of the plasma moments and the electric field measurements. Using the observed speed of the DF (314 km/s), these durations correspond to a spatial scale of~2.1-2.4 R E or~52-60 d i in x GSM . Although we note that the region is slightly smaller at P2, it is unclear as to whether this represents a temporal evolution or is simply due to the spacecraft following slightly different trajectories through the structure, but in any case, it would appear that this region of counterpropagating beams changes little in size and thus survives in quasi-dynamic equilibrium for tens of ion gyroperiods as it travels past P1 and P2. Direct confirmation that this is a discrete region attached to the DF with scale size less than the P1-P2 separation in x GSM (85 d i ) comes from the fact that P1 exits this interpenetrating beam region before the DF reaches P2.
Figures 3g and 4g show the y GSM component of the electric field. The blue time series shows E yGSM measured by the Electric Field and Waves instrument at 128 samples per second. The orange time series shows the y GSM component of E vxb = Àv × B. Strong short-period fluctuations of the electric field with amplitudes 20-30 mV/m were measured at the DF itself, but in the longer interval where the counterpropagating beams were observed, the ion plasma appears to be largely frozen-in. Whilst these electric field data are consistent with previous observations of other DFs Fu et al., 2012], it does not provide any obvious hint of the underlying structure in the ion plasma. Similarly, the temperature anisotropy does not reflect the existence of the counterpropagating beams. Referring to Figure 1, T perp > T par after the passage of the DF. However, the two interpenetrating beams are separated along B, i.e., the v par direction.

Simulations
To better understand the observations, we now examine particle-in-cell (PIC) simulations of reconnection onset and associated DF formation. The simulations have been performed using the Parsek2d implicit particle-in-cell code [Markidis et al., 2009]. The simulation is two dimensional in space (∂/∂z SIM = 0), and the box size (x SIM and y SIM ) is 200 and 30 d i or 5120 and 768 cells and an average of~400 particles per cell. The system is initialized with a Harris sheet equilibrium: B xSIM (y) = B 0 tanh[(y SIM -y 0 )/(0.5d i0 )], where the current sheet width is 0.5 d i0 . A small guide field B zSIM = B g = 0.1 is also included to break unphysical symmetries and to yield more realistic fields, flow velocities, and second moments. The density profile is given by n(y SIM ) = n 0 sech 2 [(y SIM -y 0 )/(0.5 d i0 )] + n b , where the background density n b = 0.1 n 0 , but the inertial length d i is based on the unperturbed Harris current sheet density, n 0 . The ion to electron mass ratio m i /m e = 256, the electron Journal of Geophysical Research: Space Physics 10.1002/2014JA020516 thermal velocity v (e, th) /c = 0.045, and c/V A0~1 00, which for physical electrons gives T e~1 keV and T i~5 keV. The simulation coordinate system is such that the current sheet normal points in the +y SIM direction and the reconnecting magnetic field component points in the +x SIM direction above the current sheet. The conversion between the GSM and SIM coordinate systems is as follows: x GSM = x SIM , y GSM = Àz SIM , and z GSM = y SIM . Figure 5a shows the simulated density (logarithmic scale) at t = 30/Ω ci . The center of the preexisting current sheet is located along y 0 = 15 d i . At this time, reconnection has started at an X point at (x, y) SIM~( 100, 15), but the Alfvénic jets are far from the edge of the box at x SIM = 200 d i . A sharp pressure front has formed at x SIM~1 24 d i , separating the high-density preexisting remainder of the torn Harris plasma sheet (dark red region to the right) from the low-density reconnection jet (blue through light red region to the left).
Overplotted magnetic field lines show that at this interface, there is a strong y SIM component to the magnetic field, which corresponds to the z GSM signature of the DF in the satellite data. The overall stability of the DF in the experimental data is confirmed by the simulation as the DF smoothly propagates from left to right away from the X line as the simulation progresses. In comparing the satellite data with the simulation, we first use the magnetic field data to determine qualitatively the most consistent trajectory through the simulation. While it is possible to generate time series at a specific location from the simulation data (i.e., making a true virtual spacecraft observation), the stability of the DF suggests that spatial cuts can be used for the comparison. Figure 6 shows cuts through the simulation domain along x SIM at five different values of y SIM (y SIM = 12.5, 14, 15, 16, and 17.5 d i ) at t = 30/Ω ci . The locations of these cuts are shown in Figure 5a.
We first examine the behavior in B xSIM (Figure 6c), which corresponds to B xGSM . Referring to Figure 2, it can be seen that the simulation cuts at y SIM = 14 and y SIM = 16, in green and cyan, respectively, are most similar to the trajectories of THEMIS P2 and P1. The cuts at y SIM = 12.5 (blue) and y SIM = 17.5 (purple), further from the current sheet, do not show an increase and decrease in the magnitude of B x as is observed. The cut along y SIM = 15 (red) corresponds to the midplane of the current sheet. Figure 6e shows the B ySIM component which exhibits the expected increase in this field component that is fundamentally characteristic of DFs. The cuts at y SIM = 14 (green) and y SIM = 16 (cyan) are again consistent with the THEMIS data. Figure 6d shows the ÀB zSIM component (corresponding to B yGSM the out of plane component). The simulation cut along y SIM = 16 (cyan) captures the negative deviation in the out of plane component on the high-density side that was observed by P1 (Figure 2b). The deviations between y SIM = 14 (green) and y SIM = 16 (cyan) data are qualitatively similar to the observations. Finally, cuts of the electric field components are shown in   Figure 6i shows the out of plane, or reconnection, electric field whose large-scale structure is similar to that observed (Figures 3g and 4g).
Consequently, we conclude that the observed magnetic field data are most consistent with cuts made just above and just below the midplane with P1~y SIM = 16 (cyan) and P2~y SIM = 14 (green). However, we emphasize that care is required when comparing experimental and simulation data. First, the available simulation domain is smaller than the observed structures, so the simulated DF is very much closer to the X line than it is in reality. Second, the thickness of the simulated Harris current sheet is significantly less than that observed, so the simulated DF is relatively curved and not extended in the y SIM direction, unlike the observations, which indicate that the DF is extended in z GSM direction. Third, the simulation is 2.5 dimensional (i.e., ∂/∂z SIM = 0) and so is not able to capture the development of, e.g., instabilities which may grow in the out of plane direction. These last two points may explain in part the differences between the observed and measured electric fields near the magnetic field ramp associated with the DF. We note that the separation of the two satellites in y GSM is~1 R E , which is less than the expected scale length of variability in BBFs [Nakamura et al., 2004], and in fact, the qualitative agreement between the experimental data and the simulations extends to the ion distribution functions as is described in the next section.

Simulation Distribution Functions
Having established that the appropriate comparison to be made with the THEMIS P1 and P2 data is with cuts along y SIM = 16 and y SIM = 14 through the DF, Figures 5b and 5c show the simulation ion distribution functions from the cut above the current sheet at y SIM = 16, and Figures 5d and 5e show similar data from the cut at y SIM = 14. Figure 5a shows where along the cut the ion distributions are extracted. The ion distributions are shown in the x-y SIM plane, with the x SIM axis horizontal and y SIM axis vertical. It should be noted that the simulation distributions are "instantaneous" (averaged over space at one time), whereas experimentally observed distributions are accumulated over time (averaged over space and time). The simulation distributions are also integrated along the third dimension, whereas cuts through the experimental distributions were shown in Figures 3 and 4. We have extensively analysed the simulation data and the experimental data, using a variety of visualisations that employ cuts, integrations, and different coordinate systems, and these formats reveal the structure most cleanly.
In the high-density region to the right of the DF, a single distribution is observed with asymmetry in the +v xSIM direction. In the low-density region to the left of the DF in the leading edge of the reconnection jet, multiple ion populations are observed, including peaks at +v xSIM /+v ySIM and +v xSIM /Àv ySIM with speeds of the order of the Alfvén speed. The magnetic field is dominated by B ySIM , and so the beams are separated in velocity space approximately along the magnetic field.
Comparison of the observed and simulated ion distribution functions must also be treated with care. To the first order, the simulation data can be compared to the satellite data by rotating the simulated distributions 90°counterclockwise, since v xSIM~vxGSM and v ySIM~vzSIM . These populations are therefore consistent with the experimental data shown in Figures 3 and 4, with features appearing in the appropriate locations in both velocity and real space, and with similar spread in velocity. It should be noted that because of the domain size, the simulation does not reproduce the observed separation of the macroscopic region of counterstreaming ions from the X line itself; no bulk flow distributions similar to those presented in Figures 3d and 4d are present. This also means that the features in the background are somewhat different. However, the goal of the present study is to characterize the ion distributions associated with the DF, and careful examination of the simulation therefore confirms the existence of the different main types of population observed in the data, on either side of the DF.

Time History of Ion Distributions
To better understand the dynamics associated with the observed ion distributions, we use the simulation data to trace backward and forward in time the counterpropagating beams found in the simulation in Figure 5 at t = 30. To do this, we concentrate on the distribution observed at x = 120 and y = 14 shown in Figure 5d. Figures 7a-7c show the three orthogonal 2-D views of this velocity distribution. As before, each view is integrated over the third orthogonal velocity component. Two small regions (cubes in velocity space) were chosen to represent the two dominant beam-like components (black: beam A and red: beam B), with the edges of the cube taken to contain the maxima in the v xSIM -v ySIM and v ySIM -v zSIM plots (Figures 7a and 7c).
The resulting regions enclosed in v xSIM -v zSIM (Figure 7b) do not appear to correspond to obvious maxima, but these boundaries are dictated by the boundaries already chosen. The ions composing the two beams are clearly distinguishable, with beam A moving primarily downward and to the right, while beam B moves primarily upward and to the right (as labeled in Figure 7d). Movies S1 and S2 in the supporting information show the animations of these particle trajectories superposed on E zSIM and E ySIM , respectively.
The ions do not gain energy by crossing the separatrices from the lobe into the exhaust but instead are initially part of the unperturbed thermal population on the flanks of the Harris sheet and get swept up by the front, which is moving past them. At the starting point of the trajectories (t = 23.28), the ions are to the right of the front (which is then at x~115 as indicated by the maximum in |E zSIM | in Figure 7d), although the "Hall" electric field (E ySIM ) already extends well past the front (Figure 7e) (this corresponds to the rapid propagation pointing flux in the vicinity of the separatrices that has previously been reported [Eastwood et al., 2013;Shay et al., 2011]). At this time, the ions constituting beam A form part of the thermal plasma on the high-density side of the DF, where E zSIM (the reconnection electric field, in the Àz SIM direction) is weak, and the magnetic field is dominated by B xSIM .
As the DF approaches, at t ≈ 28, the reconnection field (E zSIM ) strengthens, and ions in beam A are accelerated in Àz SIM . The beam is driven toward the midplane of the current sheet by the Lorentz force associated with B xSIM , causing the velocity to rotate from the Àz SIM to the Ày SIM direction. Closer to the midplane, the magnetic field is dominated by B ySIM , and this causes the velocity of the ions in beam A to be rotated from the Àz SIM to the +x SIM direction by the v zSIM × B ySIM force, which is strong because B ySIM is enhanced. It should be noted that throughout this process, the ions are not frozen in. The acceleration of the ions constituting beam A in the x direction from rest means that their average speed is less than that of the DF, and so by t ≈ 30, the ions in beam A are in the low-density region where they meet the ions constituting beam B (Figures 7f and 7g). Ions associated with beam B move the opposite way since they start from a region of negative B xSIM . The ions in beams A and B are accelerated by the DF and become entrained in the reconnection jet. By t ≈ 34, the ions in beam A have caught up with the DF. Figure 8 shows that by t ≈ 40, some ions from beam A have passed back into the upstream region moving ahead of the DF (cf., Figure 5c), whereas others remain part of the reconnection jet to the left of the DF.
To better understand the relationship between these ions and the ions which are reflected and accelerated by gyrating in the strong normal magnetic field at the DF near the B xSIM = 0 plane, Figure 8 also shows beam C (green).  Figure 7). At the time the snapshot of E zSIM is taken, the ions of beam A correspond to the particles now to the right of the DF at x > 140, whereas the ions of beam B are located over a range of y values immediately left of the DF. Beam C, previously reflected at the DF, has escaped beyond the right-hand edge of the figure. See also Movies S1 and S2 in the supporting information.

10.1002/2014JA020516
The ions constituting beam C are located near B xSIM = 0 in the high-density plasma adjacent to the DF at x and y = 125 and 14.5 and at t ≈ 30 (the time when beams A and B cross in physical space). These ions form part of the initial Harris current sheet. When the DF reaches these ions, they experience a rapid impulsive acceleration in the Àz SIM direction and gyrate so as to subsequently move in the +x SIM direction with a velocity greater than the DF. Thus, the acceleration and entrainment of ions in beams A and B coexist with ion reflection at the DF, and both processes contribute to the more general population of so-called reflected ions in the high-density region. Movies S1 and S2 in the supporting information show the full trajectories of the ions contained all three beams.
The time history of the ion kinetic energy (KE) in the various beams is shown in Figure 8. The kinetic energy is scaled such that an ion moving with velocity V A0 would have KE = 1. An ion at the initial thermal velocity of the simulation (5 keV) would have KE = 0.42, consistent with the range of ion energies at the start of the trajectories. When beams A and B cross at t ≈ 30, both have gained significant energy, with beam A having the higher energy. After crossing, both beams continue to gain energy, and by t ≈ 40, the distribution of beam A has been rapidly thermalized, with individual particle energies ranging from 1X to 10X their original values. Beam C is more energetic, but beam A is more dispersed in energy. The peak energies of beam A, an order of magnitude greater than the initial thermal energy, i.e., several tens of keV, are slightly less than the energy of beam C but comparable to the reported energies of earthward moving ions reflected by DFs . Finally, we note that at t ≈ 40, beam B appears to remain on the low-density side of the DF, where it continues to form part of the reconnection jet, and its energy is reduced.

Discussion
Close examination of the ion distribution functions observed in the vicinity of a DF in the Earth's magnetotail reveals a spatially limited region of counterpropagating ion beams moving parallel/antiparallel to the magnetic field, which is dominated by B zGSM . These counterpropagating beams are a persistent feature, seen first by P1 and then by P2, separated by 69 s (~29 ion cyclotron periods) and 3.4 R E (~85 ion inertial lengths). The simultaneous P1 and P2 observations show that this region is confined to the leading edge of the jet and does not extend all the way back to the vicinity of the X line. The region is macroscopic, with the observed durations corresponding to spatial scales of~2.1-2.4 R E or~52-60 d i in x GSM . For completeness, we have also examined the ion distribution plots for this event at P3 and P4, which were located closer to the Earth; these spacecraft also show similar counterpropagating ion beams when B zGSM is large. Furthermore, this behavior does not appear to be restricted to this event. Examining six other previously published DF events observed in the midmagnetotail, five are found to display signatures which are similar to this event (the other appears to go out of the ESA energy range) as described in Table 1.
The existence of the counterpropagating beams is masked by the moments. In this event, it is observed that T i, perp > T i, par immediately after the arrival of the DF, even though the parallel and antiparallel ion beams may introduce an increased effective T i, par . More generally, it has been found in flow bursts that T i, perp /T i, par > 1 in 63% of the flow burst events and that T i, perp /T i, par is enhanced when compared to the unperturbed plasma sheet observed prior to the arrival of the flow burst [Kim et al., 2010]. Since the DF exists in approximate steady state, this requires the electric field to remain approximately uniform, and so the decrease in the plasma jet bulk speed v xGSM closer to the DF simply corresponds to the increase in B zGSM near the DF, and the increase in perpendicular temperature might naturally also be associated with the increase in B zGSM . The existence of the counterpropagating ion beams and the associated rich kinetic physics is thus not evident from consideration of the fluid properties alone. One implication of this concerns calculations of depleted specific entropy in magnetotail flow bursts (i.e., reduced PV γ ), which is considered crucial to explaining how and why they can penetrate deep into the dipolar tail region. In fact, this parameter is reduced precisely in the region where the counterpropagating ion beams are observed. Since the ion distribution remains non-Maxwellian as it propagates a significant distance in the magnetotail, and given the complex ion kinetic physics at work, entropy calculations must be treated extremely carefully, particularly in global MHD simulations of reconnection onset and substorms.
We now discuss the source of the beams and the implications for reconnection as a whole. Counterpropagating beams are a recognized signature of reconnection exhausts seen both in the solar wind [Gosling et al., 2005] and in the magnetotail [Raj et al., 2002]. For example, in magnetotail bursty bulk flows, it is not uncommon to see two beams of different temperatures, especially at higher latitudes away from B xGSM = 0 [Raj et al., 2002]. These beams are due to ions entering the exhaust from opposite sides of the current sheet, and differences in the temperature and speed of the beams are thought to be due to whether the beams have crossed the jet and the current sheet to the observation point or not. We also note that ion beams have also been observed and simulated in the ion diffusion region [Aunai et al., 2011;Wygant et al., 2005].
However, the ion beams we observe here adjacent to the DF, and which were also observed in Geotail data [Fujimoto et al., 1996;Nagai et al., 2002], are different from these previous observations. Even though |B| is large near the DF, comparison of the observed magnetic field with the simulated field shows that the spacecraft are at "low latitude" because of how the magnetic field geometry is distorted by the DF. The PIC simulation of reconnection leading to DF formation reproduces many of the key observed features of the counterpropagating beams, with agreement in the speed and typical temperature, and enables the history of the beams to be studied. In particular, the simulation shows that the ions constituting the beams come from the preexisting population at the edges of the preexisting Harris current sheet and not from the lobe as previously thought [Nakamura et al., 1998], with important conceptual implications.
In classic pictures of well-developed, steady state reconnection, plasma from the inflow regions on the two sides of the current sheet is reconnected and that plasma alone forms the jet. However, the counterpropagating ion beams observed in this event do not cross the separatrices from the inflow region in the traditional sense but come from the preexisting plasma sheet. As such, the leading edge of the jet is not fed simply by plasma crossing the separatrices from the inflow region, it is also continuously fed by the preexisting plasma of the current sheet that is being reconnected. While the DF may locally appear to be a tangential discontinuity [e.g., Fu et al., 2012], and thus a well-defined boundary separating the preexisting current sheet from the reconnection jet, its limited extent perpendicular to the reconnection jet means that in fact the two plasmas can mix. The present work therefore shows that when studying reconnecting current sheets, the nature of the plasma in the preexisting current sheet must be considered as it can be entrained and become part of the leading edge of the reconnection jet itself.
Furthermore, ions from the preexisting plasma sheet can be accelerated and escape back into the high-density side by overtaking the DF. This mechanism appears to exist in parallel with the more well-known ion reflection that occurs in the vicinity of the B x = 0 plane due to the strong normal magnetic field at the DF. The consequences for this particular event have been studied by Ge et al. [2012], who showed that earthward streaming ions were present in the high-density preexisting plasma adjacent to the DF and used global MHD simulations to show this was linked to enhanced proton precipitation and enhanced proton aurora. Interestingly, Figure 8b shows that the ions from beam A, which have escaped the DF, are on field lines away from the midplane, and so in a geophysical context would presumably more easily propagate earthward toward auroral latitudes.

Conclusions
In this paper we have presented new THEMIS observations of ion distribution functions associated with DFs in the Earth's magnetotail. A region of counterpropagating ion beams is observed in a macroscopic volume spanning~2.1-2.4 R E or 52-60 d i in x GSM on the low-density side of the DF. This corresponds to the leading edge of the BBF/reconnection jet. Comparing the observations of P1 and P2 (separated by 69 s (~29 ion gyroperiods) and 3.4 R E (~85 d i ) in x GSM ), we find that this region and the associated DF appear to exist in a Journal of Geophysical Research: Space Physics 10.1002/2014JA020516 state of approximate dynamic equilibrium (d/dt = 0 in the DF frame) and is indeed spatially limited. The electric field is uniform, and calculations suggest that the ion plasma remains largely frozen-in except in the DF itself. Thus, such calculations do not reflect the very significant ion kinetic behavior that is observed. We have examined several other events observed by THEMIS, which are consistent with this picture.
Analyses of the experimental and simulated magnetic fields were used to determine appropriate cuts through the simulation for comparison. The simulated ion distribution functions reproduce the observations, particularly the counterpropagating beams. By studying the time history of the particles, we find that the beams are sourced from the thermal population in the preexisting plasma sheet and do not come from the inflow region or the lobe; the beams do not cross the separatrices. In fact, they are accelerated by the reconnection electric field and gyrate around the edge of the DF into the jet. Once entrained into the jet, they may then subsequently overtake the DF and leak back into the high-density region ahead of the DF. This leakage occurs in parallel with the more well-known ion reflection process that occurs in the strong normal magnetic field at DF near the B x = 0 symmetry plane. Their energies are consistent with previous observations of earthward moving ions reflected by DFs.
A very important consequence of this analysis is that simple fluid-based cartoons of a reconnection jet interacting with the preexisting plasma sheet, causing it to be pushed out of the way, or piled up, do not capture important aspects of the physics. In fact, some of the thermal ions initially in the preexisting plasma sheet are swept up in the low-density, fast-moving reconnection jet, and indeed become part of the leading edge of the jet itself. Thus, it would appear that the leading edge of the reconnection jet does not solely push the preexisting current sheet out of the way, but also entrains some of this preexisting plasma, meaning that the region over which the jet interacts with the preexisting plasma is not simply defined by the thin DF, but is macroscopic.
This also means it is not correct to think of the BBF as an isolated bubble of discrete plasma consisting of separate low-entropy plasma propagating earthward. The PIC simulations presented here show that nonfrozen-in ion behavior allows particles to enter and leave the low-density "bubble" region continuously, as it propagates a significant distance. This may have important implications for understanding the energy balance at the leading edge of reconnection jets. Finally, a limitation of the present work is that the simulation domain is insufficiently large to reproduce the observed separation of this macroscopic region of counterstreaming ions from the X line itself, and further work is required using simulations with larger domains to fully explore these processes.