Catalog of fine-structured electron velocity distribution functions $-$ Part 1: Antiparallel magnetic-field reconnection (Geospace Environmental Modeling case)

To understand the essential physics needed to reproduce magnetic reconnection events in 2.5-D particle-in-cell (PIC) simulations, we revisit the Geospace Environmental Modeling (GEM) setup. We set up a 2-D Harris current sheet (that also specifies the initial conditions) to evolve the reconnection of antiparallel magnetic fields. In contrast to the GEM setup, we use a much smaller initial perturbation to trigger the reconnection and evolve it more self-consistently. From PIC simulation data with high-quality particle statistics we study a symmetric reconnection site, including separatrix layers, as well as the inflow and the outflow regions. The velocity distribution functions (VDFs) of electrons have a fine structure and vary strongly depending on their location within the reconnection setup. The goal is to start cataloging multidimensional fine-structured electron velocity distributions showing different reconnection processes in the Earth's magnetotail under various conditions. This will enable a direct comparison with observations from, e.g., the NASA Magnetospheric MultiScale (MMS) mission, to identify reconnection-related events. We find regions with strong non-gyrotropy also near the separatrix layer and provide a refined criterion to identify an electron diffusion region in the magnetotail. The good statistical significance of this work for relatively small analysis areas reveals the gradual changes within the fine structure of electron VDFs depending on their sampling site.


Introduction
The notion of magnetic reconnection was originally introduced by Giovanelli (1946) to the space and astrophysical plasma physics community in order to explain violent energy releases, such as solar flares and coronal mass ejections at the Sun. Nowadays, magnetic reconnection is known to occur also in the Earth magnetosphere, in particular at the dayside magnetopause, the cusp region, and in the magnetotail. Magnetic reconnection is the most likely mechanism to drive auroral substorms (Russell and McPherron, 1973). Various theoretical models have been proposed to explain the mechanism operating in magnetic reconnection. Major examples are the viscous-type reconnection by Sweet (1958) and Parker (1957), the slow-shock acceleration model by Petschek (1964), and the discontinuity-compound model by Sonnerup (1970). A number of numerical simulations as well as remote and in situ observations have also been performed to understand the spatial and temporal development of the reconnection process (Paschmann et al., 2013;Karimabadi et al., 2013;Baumjohann, 2013, 2015). Magnetic reconnection is also observed in laboratory plasmas, yet it is difficult to reach the non-collisional regime in reconnection experiments and transfer the obtained results to a spaceplasma environment; see reviews by Zweibel and Yamada (2009) and Yamada et al. (2010).
Magnetic reconnection requires the breakdown of the frozen-in magnetic field from the magnetohydrodynamic point of view. The motion of the magnetic field lines is described by the induction equation, and solving this equation requires detailed knowledge of the electric fields in the plasma, particularly on the kinetic scales where individual particle motions (gyration, drift, wave-particle reso- Ph.-A. Bourdin: Catalog of electron velocity distributions -Part 1: GEM case nance) will be effective. Therefore, the reconnection region is divided into distinct scales: macroscopic magnetohydrodynamic scales and microscopic kinetic scales around the reconnection site where individual particle species need to be treated. In an approximation based on the two-fluid model of plasma, the electric field on the kinetic scales is evaluated by the generalized Ohm's law, although it is a rather simplified picture, neglecting the wave-particle interactions such as the cyclotron or Landau resonances (Pritchett, 2001). We evaluate the induction equation with the help of a reduced protonto-electron mass ratio (which is about 1836 in reality) and the generalized Ohm's law.
In particular, when the current sheet thickness (measured by the gradient scale or inhomogeneity of the magnetic field) becomes comparable or even smaller than the particle gyroradius, the velocity distributions are no longer Maxwellian nor gyrotropic (Hoshino et al., 2001). Recent kinetic simulations show that the velocity distributions are indeed unique with various realizations: two-sided together with a triangular distribution (Ng et al., 2012) and also including a swirl distribution (Bessho et al., 2014;Shuster et al., 2015). These non-gyrotropic velocity distributions are obtained from a particle-in-cell (PIC) simulation in a two-dimensional reconnection setup. Qualitatively, the non-gyrotropic distributions can be understood as effects of electron meandering motions (Horiuchi and Sato, 1994), acceleration through electric fields, and deflection by the magnetic field -similar to the situation observed for ions by Nagai et al. (2015).
Earlier and recent kinetic studies indicate that the electron stress must be the most relevant effect at the center of the reconnection region in a steady state (Horiuchi and Sato, 1994;Hesse and Winske, 1998;Pritchett, 2001;Hesse et al., 2011). We aim to associate the electron velocity distribution functions with various regions of magnetic reconnection. To this end, we run a numerical experiment to generate magnetic reconnection following Geospace Environmental Modeling (GEM) (Birn et al., 2001;Pritchett, 2001) and systematically characterize the 3-D electron velocity distribution functions, in particular how the distribution functions are non-gyrotropic, and where they occur.
Here we present a comprehensive catalog of nongyrotropic electron velocity distribution functions that are relevant to magnetic reconnection. Such simulation results are to be compared to Magnetospheric MultiScale (MMS) spacecraft observations from the magnetotail, similar to what was done on the dayside by . The catalog can be used to sum up electron velocity distributions along the trajectory of a spacecraft in order to characterize the magnetic-field configuration that was crossed. In particular, we describe a characteristic feature of the electron diffusion region within antiparallel field reconnection. Additional such catalogs for different configurations are needed for a better understanding of observed electron velocity distribution functions (VDFs).

PIC simulation and analysis
To compare the GEM setup with magnetotail observations, we basically require an antiparallel magnetic-field configuration, which is implemented by a Harris current sheet (Harris, 1962). We add an initial perturbation to the in-plane magnetic-field components B x and B z to trigger the reconnection in the center of the simulation domain.
We use the open-source code "iPic3D" (Markidis et al., 2010), which implements the GEM setup together with our changes 1 described in Sect. 2.1. With a large number of super-particles in our simulation we have access to good statistics about particle parameters even for small analysis regions. We use a total of 13.1 million super-particles for each species (electrons and ions) per d 2 i , where d i is the inertial length of background ions (see Sect. 2.1).
The simulation results comprise the particle and bulk velocities, the magnetic and electric fields in three components, as well as the full pressure tensor of electrons and protons.

Initial parameters and boundary conditions
We use 25 as a proton-to-electron mass ratio. We perform another simulation with a mass ratio of 100 and an unchanged background magnetic field (B 0 ). The main differences we find are as follows: (1) electron velocities are scaled up proportionally to the square root of the mass ratio; (2) the electron diffusion region shrinks in its z extent with the same proportionality; and (3) the electron velocity shear layer is located closer to the separatrix. As a result, we find similar shapes of the electron velocity distributions when we consider the slight spatial displacement of the electron velocity shear layer (see Sect. 2.3) towards the almost unchanged location of the separatrix. Using the same computational resources, a higher mass ratio implies stronger PIC noise, and the fine structure in the electron VDFs becomes less significant. Still, one should continue this study and provide catalogs for comparison with more realistic mass ratios.
The computational domain covers 25.6 × 12.8 d 2 i with 512 × 256 grid points, where d i is the ion inertial length for the background number density n 0 = 0.2 and the initial background magnetic field B 0 = 0.05477. The resulting grid spacing is ∆r = 0.05 d i . The initial condition follows B x (z) = B 0 tanh(z/D 0 ) for the magnetic field and n(z) = n 0 (1 + 5 cosh −2 (z/D 0 )) for the number density with D 0 = 0.5447 d i as the half-thickness of the current sheet.
The initial thermal velocity is v th,e ≈ 0.072 c for electrons and v th,i ≈ 0.032 c for ions. Together with B 0 = 0.05 and the mass ratio of 25, this implies a plasma beta of β e = 1/6 and β i = 5/6 for electrons and ions, respectively. The Debye length λ D relates to the grid spacing as λ D ≈ 0.29∆r. The time step remains fixed at ∆t = 0.5/Ω i , where Ω i = e B 0 /m is the ion gyrofrequency determined by the charge e, the ion mass m, and the initial background magnetic field amplitude B 0 . In the x direction we use a periodic boundary condition, as well as conducting walls at the z boundaries. y is degenerated in our case and is hence quasiperiodic or invariant.
We tested to what extent guide fields of B y = 0.1, 1, and 10%B 0 influence the results obtained from the GEM case that has no guide field. We find that guide fields of up to 1%B 0 do not significantly change the obtained electron VDFs presented in this work. For a larger guide field of 10%B 0 , we do see significant variations and hence propose to continue this cataloging study with guide fields of 1%B 0 and more, departing from the exactly antiparallel magnetic-field configuration in a stepwise way.
We use another simulation run with a quadruple box size of 51.2 × 25.6 d 2 i to verify that we have no significant influence due to the domain boundaries in the results of this work. Reflected particles from the z boundaries mainly propagate along the magnetic field, which is horizontal (mainly oriented along the x direction near the z boundaries) and does not yet connect to any of our regions of interest. Particles that cross the x boundaries penetrate only the outer simulation domain and do not reach closer than about x = ±9 d i to the reconnection center in significant quantities. An example of such particles traveling inwards is visible as a minor third peak in v x > 0.1 c at the position of x = 9 and z = 0.6 d i ; see the video in the Supplement. 2 Because we need to provide the statistical significance in order to obtain low-noise fields and fine-structured electron velocity distributions, we require an unprecedentedly large number of particles (in total 8.6 billion). When scaling the mass ratio to higher values, typical spatial scales become smaller for the electrons, which requires having smaller analysis areas. Hence, this either requires a currently unfeasible increase in computational demands while trying to maintain the data quality from this work for much higher mass ratios, or one has to change the original GEM parameters (like B 0 ), which would on the other hand prevents us from comparing this catalog directly with earlier GEM simulation results. Also, the simulation run with a larger box size needs to be conducted with less particles per d 2 i . Therefore, for now we retain the original GEM settings, including the mass ratio and simulation box size to obtain the lowest possible PIC noise level.

Comparison of reconnection rates
While the original GEM setup proposes a large perturbation (covering the whole simulation domain), we trigger the reconnection with a perturbation that is 10 times smaller in its spatial extent. The effect of this modification is a smaller gasto-magnetic pressure disequilibrium in the initial condition and a later onset of the reconnection. The reconnected field topology becomes more symmetric because both peaks in the perturbation are closer to the middle of the simulation do-  (2001), Birn et al. (2001), and Schmitz and Grauer (2006). The vertical gray dashed line indicates the time of the data snapshot we analyze in this work.
main. Therefore, we evolve the reconnection in a more selfconsistent way.
The reconnection rate is the slope of the reconnected flux. We find that our setup generates a similar reconnection rate as Pritchett (2001), while the reconnected flux is different due to the initial condition; see Fig. 1. If we subtract the difference of both initial conditions and consider that the onset of the reconnection is about ∆t = 2 Ω −1 i later in our simulation, both curves become very similar. A later onset was also observed by Birn et al. (2001), while the reconnection rate evolves differently in their work. Schmitz and Grauer (2006) used a Vlasov code to model the GEM setup, and they also find a later onset. Their reconnection rate lies in between the ones observed by Birn et al. (2001) and Pritchett (2001).
To check if the simulation domain size has a significant impact on the reconnected flux, we compare it with a simulation run that has a quadruple-sized box. All other parameters and the initial conditions are kept identical. We find an almost identical evolution of the reconnected flux until t = 22 Ω −1 i . After that time, the reconnected fluxes start to deviate from each other (see the dashed red and orange lines in Fig. 1). We find that the plateau seen in Pritchett (2001), Birn et al. (2001), and Schmitz and Grauer (2006) after t = 27 Ω −1 i is caused by the box size and is hence due to influence from the boundary conditions. Until t = 20.54 Ω −1 i , we see no significant influence from the box boundaries. Therefore, we use this snapshot for our analysis and refrain from using data at later times.
From our lager-box simulation run, we see that the magnetic influx and hence the reconnection processes are slowly being suppressed and come to a halt after t = 33 Ω −1 i ; see the differences between the red solid and the orange dashed line in Fig. 1.
We note that the comparability between different simulations is limited because we basically watch at different evolution stages of the reconnection when we check for identical times; see the vertical gray dashed line in Fig. 1. As compared to our results, the amount of total reconnected flux at this time is about 52 % higher in Birn et al. (2001); Pritchett (2001) and 21 % higher in Schmitz and Grauer (2006). Therefore, not the time of the data snapshot but instead the amount of reconnected flux (without the initial perturbation) is the better quantity to compare between simulation works.

Evolved reconnection
We display a snapshot of the reconnection in Fig. 2, where the amount of flux that has reconnected after the initial condition is 1.39 B 0 . The simulation time here is t = 20.54 Ω −1 i . We obtain the reconnection center exactly in the middle of the simulation box, where the out-of-plane component of the reconnection electric field To estimate the non-gyrotropy of electrons we use the index proposed by Swisdak (2016) (Eqs. A5-A8) that we compute as where I 1 is the trace of the pressure tensor P, its field-parallel component is P = e B T Pe B with a unit vector along the magnetic field e B , and I 2 is The region with a high non-gyrotropy index Q e roughly follows the reconnection current sheet, which is strongest at x = 0 and z = ±0.15 d i and is elongated along the mean background magnetic field direction x.
Also, we find an increased value of Q e that follows the electron velocity shear layer close to the separatrix; see green and yellow encoded u e,y located around z = ±1.4 d i in Fig. 2a. This shear layer is enclosed by an electron bulk that is accelerated away from (towards) the reconnection site along the x direction; see the black/red (white/blue) contour lines in Fig. 2 that correspond to downstreaming (upstreaming) electrons, respectively. Exactly in between these x shear flows we see a strongly enhanced out-of-plane bulk velocity along the positive y direction. There, the non-gyrotropy index Q e is strongly enhanced even outside the electron diffusion region and in the absence of a significant out-of-plane electric field; see Fig. 2b and c.
In Fig. 3 we show cuts of Q e and E rec together with the electron and ion bulk velocities u e and u i along the x and z axes. In the cut along x we find a peak in the non-gyrotropy index Q e in the reconnection center, where E rec is strongest. However, also in regions away from the center √ Q e is significantly enhanced and reaches values above 0.3 where E rec is practically not present. Therefore, we need another indicator besides the non-gyrotropy index in order to identify an electron diffusion region.
Along the z direction, the bulk velocity u e and √ Q e both have a double-peak shape, while the reconnection electric    Figure 3. Cuts for the non-gyrotropy index Q e , the reconnection electric field E rec , the scalar product the of total electric field and current density E rec · j, and the bulk velocities u, along the z = 0 line (a) and along x = 0 (b).
field E rec clearly has a single peak. The separation distances of the double peaks in √ Q e and u e,y differ by 0.1 d i . We compare √ Q e (Swisdak, 2016) with the scalar product of the reconnection electric field and the current density E rec · j . This latter quantity was also used to identify the electron diffusion region in Burch et al. (2016). We find that these two quantities anticorrelate well along the x axis (z = 0) within |x| < 2d i near the electron diffusion region; see the orange dashed and the black lines in Fig. 3a. This means that Q e is large along z = 0 and near the reconnection center, where the current density j is antiparallel to E so that their scalar product becomes negative. This correlation breaks down for |x| > 5d i , where E rec · j vanishes but √ Q e is significantly high. Also, the strong peaks in E rec · j at x = ±4d i are not seen in √ Q e . On the other hand, along the z axis we simply do not see a similar anticorrelation as along the x axis; see Fig. 3b.
For completeness, we checked that the non-gyrotropy index AØ e (as defined in Scudder and Daughton, 2008) gives similar results to Q e , except that AØ e shows a stronger relative enhancement near the separatrices as compared to the current sheet in the central reconnection region. The electron diffusion region is also indicated as being slightly larger in AØ e along the z direction than compared to Q e . Altogether, we do not find a clear correlation of Q e with any other quantity plotted in Fig. 3. This suggests that the off-diagonal elements of the pressure tensor increase also through processes that are subsequent to the electron acceleration from the reconnection electric field E rec -or that are, in other words, not directly induced in the central electron diffusion region.

Electron velocity distribution functions
We select small regions of interest with a size of 0.2 × 0.2 d 2 i , where we bin the electron velocities of particles contained in the region. This size also reflects the electron gyro-radius where the magnetic field reaches B 0 /4, like near the reconnection center. Zenitani and Nagai (2016) average the electron velocity distributions over regions of 0.5 × 0.5 d 2 i , which results in less fine-structured electron velocity distributions; see Fig. 4. Hence, we need to average over areas of 0.2 × 0.2 d 2 i in order to capture distinct electron populations -or to stay within about one electron gyration radius in a weak-field regime, like near the reconnection site. Even smaller analysis regions of 0.1 × 0.1 d 2 i would not reveal significant additional fine structures above the noise level; see row (c) in Fig. 4.
The differences between the row (a) in our Fig. 4 and the original Fig. 4, panels (e1)-(e3) in Zenitani and Nagai (2016) are due to the different time in the evolution of the reconnection process. We see that the fine structure formed during the free evolution of the reconnection (t = 20.54 Ω −1 i ), and it slowly decays or washes out when reaching the plateau phase (due to numerical constraints; see We also average the non-gyrotropy index Q e proposed by Swisdak (2016) within the analysis regions. Each region contains about half a million super-particles per species.
In Figs. 5 to 11 we provide a comprehensive catalog of 3-D electron velocity distribution functions as 2-D cuts along the simulation coordinate directions and with integrated distributions along the direction orthogonal to each 2-D cut. The text label color indicates the location of the analyzed regions group, cf. the colored dots in Fig. 2d. The 2-D velocity distributions are along the mean initial magnetic field (x), the outof-plane direction (y), and along the initial magnetic-field gradient (z, perpendicular to the central current sheet). Fig. 5 contains the reference area (top row) that shows a clearly Maxwellian shape, represented by Gaussian distributions in all three directions. The other panels in Fig. 5 are samples within the inflow region, where we find rectangular shapes in the v x -v z and the v x -v y plane, while the v y -v z remains Maxwellian, as also reported by Schmitz and Grauer (2006). The region labeled "late inflow" is already close to the reconnection center and features a gradually appearing triangular shape in the v x -v y cut.
When leaving the inflow region and approaching the reconnection site, we sample a single-peak shape along v z in   Fig. 6. At the reconnection center, we find that the inflowing electrons are accelerated by the reconnection electric field E rec along the y direction, which forms the elongated tip of the triangular-shaped distribution in the v x -v y cut; see "reconnection center" row in Fig. 6. Also, we can confirm the increasing tilt of that tip when we consecutively sample regions on the x axis but with increasing distance from the reconnection center, as shown in Shuster et al. (2015). This increasing tilt is due to the electric field vector that points perpendicular to the x − y plane only in the exact reconnection center, but has a growing in-plane component when going away from the reconnection center. The fine structure in the v x -v y panels (Fig. 6b, e, h, k, and n) is caused by the number of electron meandering motions through oppositely oriented magnetic fields above and below the reconnection center. This is very similar to the behavior of electrons for an antiparallel field plus a guide field, as described by Ng et al. (2012). We also find gradual changes from x = z = 0 (Fig. 6, middle row) to x = 1, 1.5, 2; see the three upper rows in Fig. 7 (visible even better in the online movie). It becomes clear that individual populations, like the red tip of the downwards-pointing triangle and the red vshaped population above this tip in Fig. 6h (two lowermost  peaks), are the same populations as the green and orange distributions in the v x -v y panel shown in Fig. 7h (two leftmost peaks) at x = 2, z = 0.
The strong red peaks in the center of the three lower rows of Fig. 7 are therefore from electrons that have not completed any full meandering motion, but that have basically evolved from the inflowing velocity distribution population directly; see row (c) in Fig. 5.
After crossing the reconnection center, we find a doublepeak shape along v z that is maintained until the plasma reaches the "acceleration center" (Fig. 7) which has a similar spiral shape in the v x -v y plane, as also found by Bessho et al. (2014). Figure 8 contains samples of the plasma exhaust downstream of the reconnection and acceleration regions. The velocity distributions gradually become more gyrotropic again, as seen from the "late downstream" to the "early outflow" regions. In particular, the differences in the work of Shuster et al. (2015) between their panels rows F (t = 20) and G (t = 29 Ω −1 i ) in their Fig. 3 are explicable by the suppression of the reconnection process due to the box size (or the boundary conditions) at later times (cf. our Sect. 2.2). While at t = 20 their and our distributions are of course similar, we see that the previously created non-gyrotropic distributions downstream of the reconnection site decay with time. In particular, these non-gyrotropic distributions are again approaching a more Maxwellian-like shape after t = 20, which means that the actual process that initially created those non-           Figure 11. Same as Fig. 5 but following the separatrix layer in the upstream direction towards the inflow. gyrotropic distributions has either stopped or has at least become significantly weaker. In Fig. 9 we highlight some additional regions of enhanced and unexpectedly high non-gyrotropy index Q e values in the reconnection outflow and around the secondary peaks in the reconnection electric field E rec .
We find electrons that are strongly accelerated along the local magnetic-field direction, visible as multiple distinct stripes of enhanced probability at negative v x ; see the rows "second reconnection" and "second acceleration" in Fig. 9c and d. This indicates that such electrons were undergoing multiple acceleration processes. One possible cause is that these electrons were performing multiple meandering motions within the electron diffusion region, which may indeed explain mostly equidistant stripes that are roughly orthogonal to the background magnetic field. Because we also find a region with an enhanced outwards acceleration at x = ±3.2 and z = ±0.9 d i , together with a significant reconnection electric field E rec , we identify a small secondary reconnection site at the location x = ±2.6 and z = ±0.6 d i . This finding is underpinned by the fact that the secondary acceleration region is clearly distinct from the main acceleration region that surrounds the reconnection center; see the contour lines in Fig. 2. In earlier works this feature may not have been observed as clearly because of a higher PIC noise level.
Of particular interest regarding the non-gyrotropy are the samples across the separatrix layer that we show in Fig. 10. We also find double-peak shapes (in the v x -v z distribution) at and below the separatrix layer, which basically comes from an electron velocity shear caused by nearby upstream and downstream flows. The orientation angle of these double peaks is well aligned with the local magnetic-field vector; see the green line in the "downstream shear" row ( Fig. 10e).
In the region "current sheet" we see again stripes mostly parallel to v x , which is explicable here by multiple meandering motions within the electron diffusion region. We do not see a significant bulk motion with a negative v x here.
When we follow the separatrix layer along the upstream direction, we see strongly non-gyrotropic distributions in Fig. 11 that might misleadingly be interpreted as being close to (or within) the electron diffusion region; see the "separatrix upstream" row ( Fig. 11a) with a non-gyrotropy index above √ Q e ≥ 0.3. It is important to note that the Q e parameter is indicative of crossing the electron diffusion region only with some additional criterion. For example, we find that one should also see a distribution with a double peak oriented along the z direction (or the background magnetic-field gradient) in the v y -v z components of the electron distribution function (reflecting a meandering motion) in order to identify the electron diffusion region; see the "reconnection" regions in Fig. 6c, f, i, l, and o.

Discussion regarding MMS observations
A complete set of all electron velocity distribution functions within the antiparallel reconnection site (for the original GEM case) discussed in this work is available as a movie online. 3 One should still note that MMS observations are typically time integrations that represent trajectories through the simulation domain. Hence, one probably needs to sum up multiple electron VDFs in order to match observations of antiparallel field reconnection. In follow-up publications we plan to expand this catalog with guide fields, different plasma densities and temperatures, and more realistic mass ratios.
With respect to the recently observed and discussed "crescent"-shaped electron VDFs (see Hesse et al., 2014;Burch et al., 2016) it is worth noticing that we find no such distribution. This is expected, because the antiparallel field configuration of this catalog does not fit the dayside magnetosphere.
In this work we find that non-gyrotropic velocity distribution functions for the electrons can occur not only in the electron diffusion region, but also in extended regions, in particular at the electron velocity shear layer close to the separatrix. Recent MMS observations within the dayside of Earth's magnetosphere revealed non-gyrotropic distributions that are associated with asymmetric reconnection .
The MMS mission is going to detect non-gyrotropic distributions also at the nightside of the magnetosphere and one may be misled to a wrong interpretation of the reconnection process because the association between the non-gyrotropic distributions and the spatial regions at or around the reconnection center is difficult. For an unambiguous identification of the electron diffusion region, we suggest looking for a double-peak electron velocity distribution in the yz plane along the magnetic-field gradient. This distribution should also be symmetric with respect to v z = 0 (where z is along the background magnetic-field gradient), together with a nongyrotropy index Q e of 0.3 or higher ( √ Q e ≥ 0.55).

Outlook for solar physics
Non-Maxwellian electron distributions have been predicted theoretically (Roussel-Dupré, 1980) and recently observed (Lee et al., 2017) in the solar atmosphere. Unstable solar magnetic-field configurations (e.g., triggering flares or coronal mass ejections) imply that magnetic reconnection takes place and hence currents exist that may be dissipated to heat the corona (Bourdin et al., 2013(Bourdin et al., , 2014(Bourdin et al., , 2015. Magneticfield parallel electric fields explain the localized acceleration of individual particles (Threlfall et al., 2016). From this work we see that a "reconnection-induced current" is often not a Maxwellian distribution of electrons that is shifted to-wards the direction of their center-of-mass motion. Instead, such currents have non-gyrotropic electron velocity distributions caused by the magnetic reconnection processes. These distributions may feature additional instabilities, within current sheets and when propagating into background plasma (Maneva et al., 2016), which would allow for a better understanding of (or new) onset mechanisms of solar eruptive events.

Outlook for future simulations
This particular work was intentionally performed with a rather simplified PIC setup. It is obvious that future simulations could be performed with a more realistic mass ratio (at least 10 times larger) and a larger box size allowing for a reconnection that may evolve for longer without influence from any boundary conditions. Both approaches will result in a significant increase of computational demands. The GEM parameter settings can be improved with respect to better applicability to the Earth's magnetotail by changing the density, and hence the plasma beta, to more realistic values. Also the influence of weaker and stronger guide fields should be investigated further.
For a better understanding of the plasma-kinetic processes involved, it is a good idea to repeat this experiment while tracking specific particles that resemble certain populations of interest and to inspect their individual trajectories in order to gain insights into the physical processes involved.
Recent 2-D and 3-D kinetic simulations demonstrate that nonsteady turbulent features arise when considering a more realistic large system size and/or 3-D space (Daughton et al., 2006(Daughton et al., , 2011Fujimoto, 2011;Lapenta et al., 2015). In this study, we did not treat such nonsteady features. A necessary future research topic would be to provide a catalog including nonsteady regions.
We also suggest adding a much smaller perturbation in the initial condition for similar simulations because this helps to trigger the reconnection more precisely in the box center and allows us to evolve the reconnection more self-consistently.

Outlook for future theoretical work
While we show in our catalog that distribution functions gradually change while we follow the bulk plasma through reconnection, we still find quite characteristic distribution functions for specific locations, like the inflow region, the reconnection center, the acceleration region, and the outflow. A fundamental physics questions is as follows: can we decompose any distribution found in our model as a superposition of individual transformations that are specific to distinct physical process involved in magnetic reconnection? In a sense, one could answer this by finding a fundamental and complete set of distributions (or transformations) that would allow us to construct any observed velocity distribution, where one could give "coefficients" that represent the influence of each distinct physical process that was involved in forming an observed distribution function. In return, one would gain insights into which kinetic processes the plasma was undergoing in its history before an in situ measurement.
Data availability. The simulation code, including the changes and input parameters used for this work, can be obtained from https://github.com/IWF-Graz/iPic3D/; check out the release tag "GEM-2D_2016-v1" and use the "GEM-smallpert.inp" input file. The data from this work may be provided on request.
Competing interests. The author declares that he has no conflict of interest.