Nonthermal electron velocity distribution functions due to 3D kinetic magnetic reconnection for solar coronal plasma conditions

Magnetic reconnection can convert magnetic energy into kinetic energy of non-thermal electron beams. Those accelerated electrons can, in turn, cause radio emission in astrophysical plasma environments such as solar flares via micro-instabilities. The properties of the electron velocity distribution functions (EVDFs) of those non-thermal beams generated by reconnection are, however, still not well understood. In particular properties that are necessary conditions for some relevant micro-instabilities. We aim at characterizing the EVDFs generated in 3D magnetic reconnection by means of fully kinetic particle-in-cell (PIC) code simulations. In particular, our goal is to identify the possible sources of free energy offered by the generated EVDFs and their dependence on the strength of the guide field. By applying a machine learning algorithm on the EVDFs, we find that: (1) electron beams with positive gradients in their 1D parallel (to the local magnetic field direction) velocity distribution functions are generated in both diffusion region and separatrices. (2) Electron beams with positive gradients in their perpendicular (to the local magnetic field direction) velocity distribution functions are observed in the diffusion region and outflow region near the reconnection midplane. In particular, perpendicular crescent-shaped EVDFs (in the perpendicular velocity space) are mainly observed in the diffusion region. (3) As the guide field strength increases, the number of locations with EVDFs featuring a perpendicular source of free energy significantly decreases. The formation of non-thermal electron beams in the field-aligned direction is mainly due to magnetized and adiabatic electrons, while in the direction perpendicular to the local magnetic field it is attributed to unmagnetized electrons.


I. INTRODUCTION
Magnetic reconnection is a fundamental mechanism of energy conversion from magnetic energy into plasma heating, bulk flow kinetic energy and particle acceleration in astrophysical, space and laboratory plasmas [1][2][3] . In particular, magnetic reconnection causes electron acceleration and non-thermal velocity electron distribution functions (EVDFs) in collisionless plasmas, for example, stellar coronae. The non-thermal features of EVDFs are not only caused by local plasma processes but they can also carry signatures of processes such as particle's acceleration and their interaction with turbulence in reconnection. Measurement of non-thermal EVDFs can be used to investigate local processes in magnetic reconnection. During most of the space age, ion velocity distribution functions are measured with the limited time resolution by instruments on-board spacecrafts, for example, in the solar wind [4][5][6] and in the Earth's magnetosphere [7][8][9] . Recently the availability of high-resolution instruments allows in-situ observation of EVDFs, for example, in the solar wind by the 3DP instrument onboard the WIND spacecraft 10,11 and in the Earth's magnetosphere by the multi-spacecraft Magnetospheric Multiscale Mission (MMS) allows to obtain highly-resolved EVDFs [12][13][14] . However, different from in-situ measured magnetospheric EVDFs, it is still impossible to measure in-situ non-thermal EVDFs at the reconnection sites in the solar corona.
Features and consequences of EVDFs formed by magnetic reconnection have been theoretically and numerically investigated under a wide variety of conditions and approaches since decades ago [15][16][17][18][19][20][21][22][23] . The most commonly found EVDFs of beams moving along the reconnection separatrices are due to acceleration by the reconnection electric field 21,24,25 . There are other anisotropic and non-gyrotropic EVDFs including but not limited to triangular-shaped with striations near the X-point 18 , arcs 20 , flat-top 26,27 , crescent and U-shaped distribution functions 28 . They are attributed to the meandering motion of electrons near the X-point or other processes like pitch angle scattering or magnetic gradients in the outflow region. Here we focus only on those features that are a necessary condition for the microscopic plasma instabilities that can cause radio emission, e.g., Type III solar radio bursts (SRBs), and thus allow remote diagnostics of magnetic reconnection.
The Type III SRBs are thought to be caused by electron beams accelerated in flare magnetic reconnection. The energetic electron beams then propagate outwards along the open magnetic field, e.g., into the solar wind, and emit radio emission along their way [29][30][31][32] . Non-thermal EVDFs in the solar corona are not directly measurable but their consequences in the form of electromagnetic 2 waves are, for example, by means of ground-based radio-telescopes 33,34 .
The relevant feature of EVDFs that can generate the micro-instabilities leading to formation of radio emission is related to their positive velocity gradients. In more precise terms, the so-called Penrose criterion 35 establishes a necessary and sufficient condition for electrostatic instabilities. It is based on an integral of the derivative of the velocity distribution function in the (parallel) velocity space. It represents a weighted average of the positive gradients of the distribution functions, which contribute to the instability criterion, compared to the rest of the distribution function. For an electron beam enough separated from the main electron distribution function, it is equivalent to say that a positive gradient makes the distribution unstable to electrostatic waves, in particular Langmuir waves. The physics behind this criterion is just attributed to the inverse Landau damping. There are similar but much less known sufficient and necessary conditions for transverse waves and instabilities propagating along the magnetic field due to positive velocity gradients in the parallel direction as well 36 . For purely electromagnetic instabilities due to positive gradients in the perpendicular velocity direction, and in particular for non-gyrotropic distribution functions, there seems to be a lack of a general theorem equivalent to the Penrose criterion to our knowledge.
But for the specific case of the electron cyclotron maser instability (explained in detail below), the positive gradient(s) in the perpendicular velocity direction is a necessary condition for the instability, representing the population inversion of the maser mechanism. Note that positive velocity gradients in the electron distribution functions will always generate unstable waves because of the inverse Landau damping. But their growth rate may be too small compared to the frequency of the unstable waves if those velocity gradients are relatively weak. Such unstable waves, and so the instability itself, would be negligible in comparison with the normal plasma modes or the surrounding thermal noise. So in the remainder of this paper we search for positive velocity gradients in the electron distribution functions as a necessary condition for the instabilities discussed in the following.
In the field-aligned direction, EVDFs with parallel velocity space gradients, i.e., v · ∂ f /∂ v > 0, can cause streaming-like instabilities, which in turn usually generate unstable electrostatic Langmuir waves at the local plasma frequency and its harmonics, i.e., ω = nω pe (n = 1, 2, ...) (see Refs. 30,[37][38][39][40] and references therein). Therefore the typical frequencies of this emission mechanism depend on the local plasma density. As proposed by the widely accepted plasma emission mechanism, these beam-generated Langmuir waves can further interact with ion-sound waves to produce electromagnetic emission via a multistage nonlinear wave-wave process. These electromagnetic waves can eventually escape from the ambient plasma and be remotely observed, provided their frequencies are above the (cutoff) local plasma frequency ω pe and their phase speed is equal to or greater than the speed of light ω/k ≥ c (see Refs. 30,[41][42][43][44][45][46][47] and references therein). The plasma emission mechanism due to EVDFs with parallel velocity gradients has been investigated by kinetic simulations (see Refs. 48,49 and references therein). This kind of EVDFs leading to streaming instabilities are often found in simulations of magnetic reconnection, their formation is mainly due to the parallel reconnection electric field 16,22,50 . However, the properties of the emitted waves due to those distribution functions with positive sources of free energy have been comparatively much less studied, with a strong focus on the streaming instabilities in the separatrix region of reconnection 19,51-54 . In contrast, EVDFs with a positive velocity gradient in the direction perpendicular to the local magnetic field, i.e., ∂ f /∂ v ⊥ > 0, can cause the so-called electron cyclotron maser instabilities (ECMIs, see Refs. [55][56][57][58] and references therein). The typical frequency of electromagnetic waves caused by ECMIs lies at the electron cyclotron frequency and its harmonics, i.e., ω = nΩ ce (n = 1, 2, ...). The conventional electron cyclotron maser emission mechanism operates when the electron cyclotron frequency is larger than the plasma frequency, i.e., Ω ce ≥ ω pe . This condition occurs in region with relatively low density and strong magnetic field, for example, density cavities in magnetic reconnection (see Refs. 59,60 and references therein). Those generated electromagnetic waves can then directly emit out as X polarized waves from the ambient plasma.
Treumann and Baumjohann 58 pointed out that electromagnetic waves can escape from plasmas with the opposite frequency condition Ω ce < ω pe , if the harmonics of Ω ce are excited strongly enough. In such a case, the remotely observed radio emission is possibly not at the fundamental frequency Ω ce but at its harmonics.
Analytical studies and kinetic PIC simulations of the ECMI are often relied on the so-called ring EVDF, which resembles a ring in the 2D velocity space perpendicular to the local magnetic field and satisfies ∂ f /∂ v ⊥ > 0 40,61-65 . Although a very idealized model of EVDF, ring distributions of electrons have been obtained in kinetic magnetic reconnection simulations 66,67 . However, what it is practically more often observed is partial-ring EVDFs in the 2D velocity space perpendicular to the local magnetic field, they are named as perpendicular crescent-shaped EVDFs or electron crescents. Another type of crescent along the parallel velocity direction is often discussed as well 12 .
Perpendicular crescent-shaped EVDFs have been observed by the MMS mission in the Earth's magnetopause where magnetic reconnection takes place in asymmetric configurations 12,68-71 and 4 in the Earth's magnetotail where magnetic reconnection is symmetric 72 . Kinetic PIC simulations of asymmetric reconnection with conditions similar to those in magnetopause reconnection have generated the perpendicular crescent-shaped EVDFs 20,28,67,73-77 . There are different proposed formation mechanisms behind those perpendicular crescentshaped EVDF in magnetic reconnection. In general, those EVDFs are associated with finite gyroradius effect, which preferentially occurs in the neighborhood of the X-point of reconnection where the magnetic field strength is weak and particles perform meandering motion, or near the separatrices where there are steep gradients 71 . Analytical and numerical analysis show meandering electrons have indeed led to crescent-shaped EVDFs 73,77 and even earlier to similar crescent-shaped velocity distribution of ions [78][79][80] . In particular, some mechanisms attribute the formation of those EVDF to the E × B drift caused by the electric field associated to asymmetric reconnection 20,75 , while others have claimed that crescent-shaped EVDFs do not depend on those asymmetric electric fields but rather only on magnetic field gradients 81 .
Numerical studies of electron bunches crossing tangential discontinuities have indeed revealed formation of crescent-shaped EVDFs due to gradient-B drift 82 . It is therefore expected that under strong guide magnetic fields, where magnetic field strength gradients are weaker, crescent-shaped EVDFs would tend to be suppressed. Bessho et al. 74 investigated the guide field effects on the EVDFs formed in kinetic PIC simulations of asymmetric reconnection. They identified two effects associated with a weak guide field on the crescent-shaped EVDFs: a widening of the opening angle and a cutoff in a reduced EVDF along a direction that is oblique to the magnetic field.
Egedal et al. 67 found that the observed crescent-shaped EVDFs in asymmetric reconnection can be accounted for by an extension of their trapping model 83 , basically considering an additional population of electrons coming from the magnetosheath side of the current sheet. The physical mechanism is due to an energy cutoff in the distribution functions that can be traced back to gradients and diamagnetic drifts on scales of the order of an electron gyroradius.
Note that most of the aforementioned studies are based on either a test particle approach or 2.5D PIC simulations of kinetic magnetic reconnection. But the features of 3D magnetic reconnection can be very different to its 2D counterpart. The additional degree of freedom (let us say, z-direction) could allow for a more efficient release of the free energy. Indeed, it is expected that EVDFs with positive gradients in the out-of-plane direction (i.e., along v z direction) will release their energy via streaming instabilities via unstable waves with a k z vector. This process is not possible in a 2D configuration. This is particularly relevant under the influence of a (relatively strong) guide field, because the reconnection electric field will accelerate electrons along the zdirection, which direction will be mostly aligned to the total magnetic field pointing also very near the z−direction (at least near the reconnection midplane). So the evolution of those electron distribution functions will be clearly different whether the z-direction is allowed (3D) or not (2D). 3D magnetic reconnection kinetic simulations can become more turbulent as a result of field-aligned streaming instabilities 84,85 . But a systematic characterization of the EVDFs in 3D kinetic magnetic reconnection is still lacking. And more so the identification of the EVDFs features relevant to the instabilities that can cause radio emission: i.e., velocity space gradients in EVDFs along the parallel or perpendicular directions to the local magnetic field. Only a very limited number of 3D kinetic magnetic reconnection simulations have shown evidence of, e.g, crescent-shaped EVDFs 86,87 . Moreover guide field effects on all those processes remain also unexplored, a parameter that plays a critical role in solar coronal conditions where those instabilities presumably develop.
In addition to all the previous physical considerations, the systematic analysis of the availability of sources of free energy offered by EVDFs generatd by magnetic reconnection is a challenging problem. Only very recently the first attempt for such a diagnostic has been carried out 88 . They implemented an unsupervised machine learning technique to the EVDFs formed by 2D magnetic reconnection simulations in different guide field strengths, non-Maxwellian features of EVDFs can be automatically detected. By using this method they identified the reconnection region providing an additional signature beyond more traditional criteria like those based on the non-gyrotropy of the electron pressure tensor 89 or the dissipation measure in the electron frame of reference 90 .
What is more interesting of this machine learning technique is its capability to automatically identify non-Maxwellian features such as formation of non-thermal electron beams and their possibility to offer sources of free energy causing micro-instabilities, their temperatures and associated anisotropies.
In order to fill the gaps from previous studies, we thus aim at a systematic characterization of the possible sources of free energy from the EVDFs due to 3D kinetic magnetic reconnection.
We restrict ourselves in particular to positive velocity gradients (along the parallel or parallel or perpendicular direction) in those EVDFs since they are necessary conditions for the microinstabilities eventually leading to radio emission. Note that this does not necessarily imply that the reconnection sites are unstable and emit radio waves. Rather that radio emission is due to the electron beams generated by reconnection, as per the standard model of solar flares. The electron 6 beams will eventually emit radiation further away from the reconnection region as they propagate in the solar corona and solar wind. A proper description of this radio emission would require a transport model for the electron propagation at AU scales. The objective of our paper is to characterize just the first step in this chain of events, the formation of electron beams.
For this purpose, we extend and apply a similar unsupervised machine learning algorithm already developed by 88 to determine the those features of the EVDFs in reconnection. We in particular assess the guide field effects on those sources of free energy.
The organization of this paper is as follows: in the Sec. II, simulations setup, parameters and diagnostics are described. In the Sec. III, the results of formation of non-thermal electron beams, possible sources of free energy and their dependence on the strength of guide field are discussed.
In the Sec. IV, we summarized the formation of EVDFs and the expected sources of free energy in reconnection plane in dependence on the guide field strength.

A. Simulation setup
In order to analyze the EVDFs formed in 3D kinetic magnetic reconnection, we carried out numerical simulations using the 3D fully kinetic particle-in-cell (PIC) code ACRONYM 91 (see also http://plasma.nerd2nerd.org/papers.html). We initialized our simulations with a Harris current sheet equilibrium 92 with different guide field strengths. Note that some parameters and initial conditions of the simulations are similar to those in Refs. 21,85 .
The simulations were initialized with two current sheets sufficiently separated to avoid their interaction at short timescales. The main reason to apply two current sheets is due to the implement of periodic boundary conditions in all three directions. Note that we only concentrate on one of the CSs for all investigations of EVDFs.
The initial magnetic field of the two current sheets is expressed as follows: here L = 0.25d i is the current sheet halfwidth with d i the ion skin depth. B ∞ is the asymptotic reconnection magnetic field, B g is the strength of the guide field in the out-of-the-reconnectionplane direction (or simply the out-of-plane direction, here it is the z direction) and b g = B g /B ∞ is defined as the relative guide field. L x is the simulation box size along x direction (while L y and L z 7 are the sizes along y and z directions respectively). We chose homogeneously distributed plasma and equal ion and electron temperatures, i.e., T i /T e = 1.0. The ion-to-electron mass ratio is set to be µ = m i /m e = 100. Provided the Harris sheet equilibrium and the associated equilibrium between magnetic and thermal pressures, the asymptotic magnetic field is constrained by the electron thermal speed v the = k B T e /m e (here v the = 0.1c) in the following way: here ω pe = 4πn 0 e 2 /m e is the electron plasma frequency calculated based on n 0 : the central peak density of each species defined via the following density profile derived from the Harris sheet equilibrium for the two current sheets: here n bg is the background population density added to avoid vacuum out of the current sheets.
Note that the expression of particle number density Eq. The expression is valid for the left half of the domain in the x direction (i.e., for the left current sheet), while the right half-domain (with the right current sheet) has opposite drift velocities for both species. In our simulation U z = 0.04c. This implies that the total out-of-plane current density is oppositely directed for either current sheet in the following way: the electron and ion current densities equally contribute to the total current density in each current sheet.
The systems are initialized with a perturbation in the magnetic field components B x and B y to accelerate the reconnection onset by seeding an O and X point at either current sheet. It can be derived from the following out-of-reconnection-plane vector potential: on the reconnection plane z = 8d i for each simulation.
with an amplitude of δ P = 0.07 and σ z = 0.25d i . In particular, the localization in the z direction is to trigger reconnection at the middle of the simulation box (in contrast to trigger reconnection at all points along the z direction).
All the simulations are calculated on a 3D mesh with We carried out three simulations of magnetic reconnection with varying external magnetic fields (see Table I). The initial relative guide field b g is determined by the critical magnetic field b c 93 . This quantity, which plays an important role on the linear tearing mode instability, represents the relative guide field that is necessary to fully magnetize the electrons with thermal speed v the in the central layer of the current sheet and is equal where r i = v thi /Ω ci is the thermal ion gyroradius derived from B ∞ . In this study, the critical magnetic field value is b c = 0.376. The regime b g < b c is denoted as the weak guide field limit in which electrons trajectories are dominantly affected by the Harris magnetic field B ∞ (x), while the opposite limit b g b c refers to the strong guide field regime, in which electrons are fully magnetized by the guide field. The transitional regime between the unmagnetized and fully magnetized electron trajectories is denoted in the following as the intermediate range of guide fields 93 .

B. Diagnostics
In this study, resulting EVDFs are evaluated in the following local velocity reference frame at each location in the reconnection plane 94 : here b = B/|B| is the unit vector of the local magnetic field B and u is the local electron bulk flow velocity. These three basis unit vectors determine a local right-handed orthogonal velocity coordinate system. e is the unit vector along (or parallel) to the local magnetic field, while e ⊥1 and e ⊥2 are other two unit vectors perpendicular to the local magnetic filed.
The EVDFs to be investigated are in essence the probability density functions, namely, they are estimated as the frequency of particle number per bin in the 2D or 1D velocity space and normalized by the total number of particles. As a result, the integral of the EVDF over the whole space where it is estimated will yield 1. For example, for the 2D EVDF calculated in Here v and v ⊥ = v 2 ⊥1 + v 2 ⊥2 are the local parallel and perpendicular electron velocities in the local reference frame determined by Eq. (6).
As discussed above, a possible source of free energy and necessary condition for microinstabilities eventually leading to radio emission is the presence of positives gradients in the 1D EVDFs f (v ) and f (v ⊥ ) separately. However, different from investigating positive gradients in the field-aligned EVDF f (v ), we look for positive gradients in the perpendicular EVDF In this study, an unsupervised machine learning algorithm, namely, the Bayesian Gaussian mixture model (BGMM) 95 , is extended to fit the 1D EVDFs and to identify the possible sources of free energy generated in simulations of magnetic reconnection.
We now briefly describe why we used this approach. The appropriate fitting distributions to EVDFs should take both physical and mathematical considerations into account. In low-energy regions, like those away from the current sheet, the plasma is Maxwellian (or bi-Maxwellian) distributed. Some common deviations from this distribution, prone to happen in collisionless magnetic reconnection, are electron beams with eventually anisotropic temperatures, which are usually associated to more high-energy regions. Their total EVDFs can then be represented as a superposition of Maxwellian EVDFs, at least to a first order approximation, although more complex features can also exist. This simplified assumption could be justified based on the fact that large deviations from Maxwellian EVDFs in the form of distribution functions with positive gradients are not sustainable for a long time in collisionless plasmas. This is because those gradients tend to be reduced due to the inverse Landau damping effect and consequently the deviations from a Maxwellian are reduced as well.
A distribution that combines both low-energy Maxwellian and high-energy non-Maxwellian features is the Kappa distribution, which has been extensively applied to collisionless space and astrophysical plasmas 96,97 . Mathematically speaking, a single Kappa distribution population can be approximately fitted by the sum of a central Gaussian distribution function to fit the core and a sum of Gaussian distributions function with large widths to fit its wide tail. In this way, many EVDFs in space and astrophysical plasma can be fitted by a sum of Gaussian distributions 88 .
Based on this assumption, we firstly fit the velocity distribution function of electrons by a sum of Gaussians by the so-called Gaussian mixture model (GMM) 95 as follows: here A k corresponds to the weight of the k-th component of the Gaussian distribution, N indicates the multivariate Gaussian distribution parameterized by the mean vector µ k and the covariance the parameters of the GMM. The quantity |A i | 2 represents the intensity of the k-th component of the Gaussian mixture, and µ k and Σ k indicate the bulk flow velocity and thermal speed, respectively.
The GMM model is ideal for our purposes because it allows to fit the most common distribution functions expected in reconnection and easily determine their properties such as mean bulk flow speed or temperature.
In order to compensate for the overfitting of GMM models, the second step is to assess the appropriate number of sub-populations based on the model selection technique named Bayesian information criterion (BIC). By introducing a penalty term associated with the number of parameters in the model, this model selection method can automatically determine the appropriate number of sub-populations 95 . See a detailed discussion about the BGMM algorithm applied to the EVDFs in Appendix A.
The formation of EVDFs for different guide field regimes critically depend on the electron dynamics. During the non-linear process of magnetic reconnection, the electrons crossing the magnetic field reversal can be characterized by a curvature parameter in the following form 98,99 : here

A. Global evolution and selection of time snapshots for analysis
As mentioned above, we focused on one of the current sheets in each simulation in this study. Fig. 1 (a,c,e) show the 3D spatial configuration of the out-of-plane current density j z within the 25 Ω −1 ci for Run1 (the antiparallel magnetic reconnection case with b g = 0). They roughly have a similar structure on each reconnection plane, namely, the X-point varies around x = d i , y = 4d i to some degree and the Fig. 1(b,d,f) show the corresponding 3D spatial structure of magnetic field lines near the reconnection plane at that on each reconnection plane, the topology of magnetic field lines in most regions especially in the inflow and separatrix regions behaves like the typical 2D reconnection structure. There are variations in the topology of magnetic field lines along the z direction, but they can still be roughly approximated by a 2D reconnection structure. Those variations can be seen in Fig. 10 showing the structure of the current density component j z in the out-of-the-reconnection plane direction, which modulates the variations of the magnetic field at the calculated plane (near the X-point).   Fig. 2(a)) between the X-line, O-line and boundaries of the simulation box, namely, dΨ/dt (see solid curves in Fig. 2(b1)). Note that for this calculation we assumed a standard 2D reconnection geometry for the location of the X and O points. This ignores the variation in the out-of-plane direction and the more complicated structure of the magnetic field due to the wavy structure of the current sheet in the out-of-plane direction (see Fig. 10). But in average (along the z-direction) the X and O points are approximately located at the points required by our calculations of the reconnection rate as explained above.
Due to the initial perturbation and reconnected magnetic flux available in the simulation box, reconnection saturates at t = 6.5Ω −1 ci for Run1, t = 7.5Ω −1 ci for Run2 and after t = 8.2Ω −1 ci for Run3 respectively. Those times are relatively short in units of the inverse cyclotron times in comparison with other simulations. This mainly has to do with the relatively small simulation box size along the x direction and the consequent small available magnetic flux to be converted by magnetic reconnection. The larger the simulation box the larger the available magnetic flux which implies that magnetic reconnection can be sustained for longer times. A second reason is the relatively thin initial current sheet. That enhances the growth rates of the tearing mode instability triggering magnetic reconnection. As a result simulations of magnetic reoconection will take less time to evolve and to reach saturation in comparison with equivalent simulations using an initially thicker current sheet.
The maxima of reconnection rates are larger than 0. Our larger values can be attributed to numerical effects due to the reduced size of the simulations box and the interaction of the second current sheet (not shown in those plots).
The corresponding 1D parallel EVDF f (v ), which are obtained by integrating the 2D EVDFs along the perpendicular direction (see Fig. 3(c1)), also shows locations of two separate beams at about v = ±0.2c in the parallel direction. The 2D EVDF of background electrons (see Fig. 3(b2)) are mainly Maxwellian distributed around v = 0, v ⊥ = 0 in the v − v ⊥ plane, and its 1D EVDF Fig. 3(c2)) shows a nearly bell-shaped distribution function around v = 0. In this way, the 1D EVDF of total electrons (as a sum of both Harris and background electrons) shows three separated maxima with comparable values along the parallel direction, i.e., two from the Harris electrons at v = ±0.2c and one from the background electrons at v = 0 (see Fig. 3(c3)). Fig. 4 shows the velocity distribution of electrons collected from a spherical region of radius ≤ 0.1d i centered at point B (see Fig. 2 (a2)) a bit off the X point into the inflow region on the reconnection plane at z = 8d i and t = 5.25Ω −1 ci for Run1. Fig. 4(a) shows the distribution of total electrons in the 3D velocity space: the thermal core (see dots with cold colors) and non-thermal electron beams off the local magnetic field B axis and to the positive v ⊥1 direction (see dots with warm colors). The beam is better seen in the 2D EVDF f (v ⊥1 , v ⊥2 ) of Harris electrons shown in the v ⊥1 − v ⊥2 plane (see Fig. 4(b1)), which are obtained by integrating along the parallel direction.
The Harris electrons population shows a crescent-shaped feature in the v ⊥1 − v ⊥2 plane, while the azimuthally integrated 1D EVDF f (v ⊥ )/2πv ⊥ shows that a local maximum occurs at v ⊥ ≈ 0.16c (see Fig. 4(c1)). The latter implies a positive gradient in the perpendicular EVDF. On the other hand the bulk flow speed of the background electrons is near zero and these electrons exhibit a Maxwellian distribution (see 2D EVDF in Fig. 4(b2) and its 1D EVDF in Fig. 4(c2)). As a result, the resulting total 2D EVDF shows both thermal core and the crescent-shape feature (see Fig. 4(b3)), while its 1D EVDF, on the other hand, shows no positive gradient in the f (v ⊥ )/2πv ⊥ (see Fig. 4(c3)).
Note that at the points A and B, even though positive gradients exist in the 1D EVDFs calculated from Harris electrons (see Fig. 3(c1) and Fig. 4(c1)), the gradients are not present in the 1D EVDFs from total electrons (see Fig. 3(c3) and Fig. 4(c3)). This is because the background electrons are mainly thermally distributed and contribute more significantly to the total electron density than the Harris electrons, which has a significant influence on the formation of positive gradients in 1D EVDFs of total electrons. As a result, the locations where sources of free energy are available is a subset of those where non-thermal electron beams exist.
The results show that the Harris electron population tend to be the population responsible for the non-thermal features in the total EVDFs (in both parallel and perpendicular direction to the local magnetic field), while the background population deviates slightly from a Maxwellian distribution.
In this study, we concentrated on the non-thermal electron beams which can offer sources of free energy in the form of positive gradients in their EVDFs as a necessary conditions for microinstabilities. For this purpose, we identify and focus on two types of EVDFs: (1) those with a positive velocity gradient in 1D EVDFs from the Harris electron population, which is useful to identify the formation mechanism of non-thermal electron beams. (2) those with a positive gradient in the 1D EVDF from the total electrons (as a sum of Harris and background electrons), which determines the existence of parallel/perpendicular sources of free energy. which are computed from the out-of-plane vector potential to a first order approximation.

C. Harris electron population EVDFs
( Fig. 5(b1)) to plasma frequency Ω ce /ω pe and the curvature parameter κ (Fig. 5(b2)) calculated by Eq. (8) on the plane z = 8d i at t = 5.25Ω −1 ci for Run1. In this case, EVDFs with a positive parallel velocity gradient in the Harris electron population mainly form in the diffusion region below the X point and two bottom-branches of separatrices, while the corresponding perpendicular velocity gradients are generally distributed in the diffusion region and outflow region near the midplane (roughly located at x = d i ) of reconnection plane. The cyclotron to plasma frequency ratio Ω ce /ω pe is generally greater than 0.5 in the inflow regions, while it is less than 0.5 elsewhere.
In the diffusion region and outflow region near the midplane, the curvature parameter κ is roughly about 1, namely, κ ≈ 1. This implies electrons here are mainly unmagnetized (or weakly magnetized) and non-adiabatic, they mainly contribute to the formation of EVDFs with a positive perpendicular velocity gradient. While EVDFs with a positive parallel velocity gradient are mainly due to magnetized electrons with κ ≥ 3 in the separatrices. Fig. 6 shows 2D EVDFs in the v − v ⊥ plane as well as the associated integrated 1D EVDFs f (v ) at four different points (i.e., points P i (i = 1, 2, 3, 4) in Fig. 5(a1)). These EVDFs are derived from Harris electrons within a spherical region of radius ≤ 0.1d i at each location on the reconnection plane z = 8d i and time t = 5.25Ω −1 ci for Run1. Let us discuss firstly the EVDFs at points P 2 and P 3 in the diffusion region. Fig. 6 (a2) shows that at point P 2 (see Fig. 5(a1)) two counterpropagating beams with parallel bulk flow speed ≈ v ≈ ±0.18c in the v − v ⊥ plane, while Fig. 6 (b2) clearly shows corresponding 1D EVDF f (v ) having two peaks located at v ≈ ±0.18c. A rough estimation of the uncertainty of the fitting is given by the standard deviation σ between the data f i and the fitted curve f (v ), i.e., quantities, obtained from the fitting data, should be compared with the standard deviation Eq. 9 in order to assess how significant are they compared to the input data. We can clearly see that they are significantly larger by at least a factor of 5, so this confirms the reliability of our fittings and in particular the calculated gradients. This is of course only a simple estimation with significant drawbacks that will be discussed later.
All this implies this EVDF offer sources of free energy which could possibly cause counterstreaming instabilities, provided that the background population does not significantly contribute to the total electron density. Similar results are observed in the EVDFs at point P 3 (see Fig. 6 (a3,b3)).
Let us focus now on the EVDFs in the separatrix region farther away from the diffusion region.  Fig. 6(b1)). The value differences (along vertical axis) for theses two gradients is about ∆ f = 1.85c −1 and 0.76c −1 separately, they are more than the deviation σ = 0.13c −1 . Their velocity gradient slopes are k = −11.81 and k = 5.86 (in units of c −1 ) respectively. At point P 4 (see Fig. 5(a1)), Fig. 6 (a4) exhibits an electron beam located at v = 0.1c and v ⊥ = 0.15c, while its 1D EVDF features a peak located at about v = 0.15c (see Fig. 6(b4)). The value difference (along vertical axis) for this gradient is about ∆ f = 0.74c −1 being more than the deviation σ = 0.14c −1 . Since this EVDF has a positive gradient at v = 0.15c, it offers a source of free energy which could possibly cause bump-on-tail instabilities, as long as the background population does not significantly contribute to the total electron density. Fig. 7 shows 2D EVDFs in the v ⊥1 − v ⊥2 plane, these electrons are calculated from the Harris 23 electron population at six different points (i.e., points Q i (i = 1, 2 . . . , 6) in Fig. 5(a2)) on the reconnection plane z = 8d i and at time t = 5.25Ω −1 ci for Run1. The points located near and below the diffusion region, i.e., Q 3 , Q 4 , Q 5 and Q 6 (see Fig. 5(a2)) feature perpendicular crescent-shaped EVDFs, which are partial rings centered at v ⊥1 = 0, v ⊥2 = 0 with a radius of about 0.2c. Note that the EVDFs at the points Q 4 , Q 5 and Q 6 (see Fig. 7(b1,b2,b3)) feature an additional population in the quadrant v ⊥1 < 0, v ⊥2 < 0. This implies their origin is different from the main partial ring population. Because those points are located in regions with κ ≥ 3, it is plausible to think their origin is due to the magnetized and adiabatic electrons taking gyrotropic electron motion. These EVDFs can offer perpendicular sources of free energy to cause ECMIs since positive velocity gradients exist in their corresponding integrated 1D perpendicular EVDFs f (v ⊥ ), assuming that the contribution of the background electron population to the total electron density is not significant.
The EVDFs in the outflow region, i.e., at points Q 1 and Q 2 (see Fig. 5(a2)) are dominated by a thermal core (see Fig. 7(a1,a2)). This is possibly attributed to the thermalization of non-thermal electron beams generated near the X point as they move away from the diffusion region into the outflow region. The perpendicular EVDF at point Q 2 has a positive gradient in its 1D EVDF f (v ⊥ ) and thus it could offer a source of free energy to cause ECMIs, provided again that the background electron population does not significantly contribute to the total electron density. However, the perpendicular EVDFs at point Q 1 are unable to offer sources of free energy since the electrons are mainly Maxwellian distributed around v ⊥1 = 0, v ⊥2 = 0 with a thermal spread width of about 0.1c. Fig. 8(b) shows examples of typical parallel and perpendicular 1D EVDFs at four different points (i.e., T i (i = 1, 2, 3, 4) in Fig. 8(a)) in the diffusion region and separatrices of the current sheet on the reconnection plane at z = 8d i and t = 5.25Ω −1 ci for Run1. The 1D EVDF at each point T i (i = 1, 2, 3, 4) is calculated from the total electrons within a spherical region of radius ≤ 0.1d i centered at the point on the reconnection plane. Fig. 8(b1) shows the 1D parallel EVDF f (v ) at the point T 1 . The EVDF is practically a Maxwellian distribution with a mean velocity v ≈ −0.01c. The fitting by the BGMM method allows to determine the width of this thermal EVDF as v the ≈ 0.14c. There is no positive velocity gradient in the 1D EVDF and thus it cannot offer any source of free energy. Fig. 8(b2) shows the 1D parallel EVDF f (v ) at point T 2 . Its standard deviation (an estimation of the uncertainty of the fitting) from the fitted curve, using the definition of Eq. 9 is σ = 0.11c −1 . Besides a peak centered at v = 0, there are other two peaks located at v ≈ ±0.2c. These two peaks indicate two possible non-thermal field-aligned electron beams. The left-most beam, centered at v ∼ −0.2c constitutes 36% of the total electron population. It features a positive velocity gradient with slope k = −4.3 (in units of c −1 , see the red fitting slopes). Its associated difference (along the vertical axis) is ∆ f = 0.42c −1 is greater than the standard deviation σ = 0.11c −1 , but not by the same (larger) amount as in Fig. 6. Although this could indicate the reliability of the fitting and the calculated gradient or slope, it is important to notice that the standard deviation of the input data is not homogeneous, but small at the tails of the VDFs and larger near the maxima, where the slopes are actually calculated. So it is necessary to be cautious with the interpretation of the goodness of fit of our calculations, keeping in mind that the determined uncertainties provide just a very rough estimation. We have also carried out similar calculations for different bin widths of the input VDF data (not shown here). From this we can conclude that the deviations between the input data and the GMM fittings decrease as the bin width of the input data increases, and it is good enough in such a way that the slopes of the positive velocity gradients can be sufficiently distinguished, in the sense that their vertical differences (∆ f ) are larger that the (global) standard deviation between the data and the fitting. More work would be needed to make more accurate statements about those uncertainties. This all implies that source of free energy, i.e., v · ∂ f /∂ v > 0, could be available to eventually cause streaming-like instabilities. Fig. 8(b3) shows the 1D perpendicular EVDF f (v ⊥ )/2πv ⊥ at point T 3 . The EVDF has a maximum near v ⊥ = 0 and monotonically decreases to zero. There is no source of free energy. Fig. 8(b4) shows the 1D perpendicular EVDF f (v ⊥ )/2πv ⊥ at the point T 4 . The standard deviation between data and fitting from Eq. 9 is about σ = 1.06c −2 . There is a positive velocity gradient with a slope of k = 17.4 (in unit of c −2 ) in f (v ⊥ )/2πv ⊥ between v ⊥ ≈ 0.05c and v ⊥ ≈ 0.16c. Its associated vertical difference is ∆ f = 1.58c −2 , which is just barely about σ = 1.06c −2 . This indicates that the uncertainty of this gradient with respect to the input data is relatively large, so that the positive velocity gradient is not that robust and prone to noise, but it exists. The associated effective slope in f (v ⊥ ) in f (v ⊥ ) can be estimated as k e f f = 5.6 (in unit of c −1 ). We apply the BGMM method to fit the 1D EVDF f (v ⊥ )/2πv ⊥ , the non-thermal beam constitutes about 61% of the total electron population. As a result, perpendicular sources of free energy could be available to cause ECMIs. from total electrons within a spherical region of radius ≤ 0.1d i centered at this point. In this way, not only the spatial distribution of the resulting EVDFs with positive velocity gradient(s) and thus the possible sources of free energy for micro-instabilities can be visualized, but it also allows to partly understand their formation mechanism by comparing them with the spatial distributions of electromagnetic fields and other quantities. Fig. 9 shows the distribution of possible sources of free energy, i.e., EVDFs with positive gra-dients in the parallel and perpendicular direction to the local magnetic field separately along the out-of-plane direction at z = 12d i , 8d i , 4d i respectively, at time t = 5.25 Ω −1 ci for Run1. Along the field-aligned direction, sources of free energy at z = 12d i are distributed in the left-top and rightbottom separatrix branches, as well as in the diffusion region below the X point (see Fig. 9(a1)).

D. Identification of sources of free energy
Most of the sources of free energy at z = 8d i are observed in the diffusion region below the X point and bottom separatrix regions along the magnetic field lines (see Fig. 9(b1)). Fig. 9(c1) shows sources of free energy at z = 4d i mainly form in the diffusion region above the X point and top separatrix regions along the magnetic field lines. In the direction perpendicular to local magnetic field, sources of free energy at z = 4d i are mainly distributed in the outflow regions near the midplane of reconnection at different heights (see Fig. 9(c2)), while those at z = 12d i and z = 8d i are observed in the diffusion region and outflow regions near the midplane (see Fig. 9(a2,b2)).   (Fig. 10(a1)), the current density structure is wavy along the z direction, probably due to the density gradients at the edge of the current sheet driving the so-called lower-hybrid drift instability. We observed that this current sheet "kinking" weakens with increasing guide field (see Fig. 10(a2,a3)). Fig. 10(a3) shows reconnection with a stronger guide field actually makes the current sheet structure to behave more 2D-like. This kinking of the current sheet modulates the distribution of sources of free energy (distribution functions with positive gradients) along the z-direction. Such a modulation would be absent in 2D reconnection, because there is not such a kinking instability in a 2D geometry. This current sheet oscillation causes a complicated flow pattern with alternate and different shear and counter-streaming bulk plasma flows on different z−planes. And this, in turn, explains the asymmetry of the distribution of sources of free energy, in particular, in antiparallel magnetic reconnection Run1 (see Fig. 9(left)).  (see Fig. 11(a1)) and the strong guide field reconnection Run3 (see Fig. 11(c1) Eq. (8) (right column) in the reconnection plane z = 8d i for Run1 (b g = 0, t = 5.25Ω −1 ci ), Run2 (b g = b c , t = 5.5Ω −1 ci ), Run3 (b g = 2b c , t = 6.75Ω −1 ci ), respectively. Other quantities are same to those in Fig. 9. energy mainly form in the diffusion region off the X point and separatrices, while sources of free energy are mainly located in the top two separatrices for weak guide field reconnection Run2 (see Fig. 11(b1)). Along the perpendicular direction, sources of free energy are generally distributed in the diffusion region and the outflow regions near the midplane for the antiparallel reconnection Run1 (see Fig. 11(a2)). However, few sources of free energy form on the reconnection plane for Run2 (see Fig. 11(b2)) and Run3 (see Fig. 11(c2)). This implies the guide field strength has a significantly negative influence on the formation of the sources of free energy in the direction perpendicular to the local magnetic field: as the guide field strength increases, the spatial extent of perpendicular sources of free energy decreases significantly. The curvature parameter κ is also estimated based on Eq. (8) at reconnection plane z = 8d i for simulations of magnetic reconnection in different guide fields strengths (see Fig. 11(a3,b3,c3)). Along the field-aligned direction, sources of free energy are mainly due to magnetized and adiabatic electrons with κ ≥ 3 in the separatrices. While in the perpendicular direction to the local magnetic field, sources of free energy are generally observed in the regions with κ ≤ 3, this implies formation of perpendicular sources of free energy is mainly due to the unmagnetized and non-adiabatic electrons. For the antiparallel reconnection Run1, the curvature parameter κ is less than 3 near the midplane (i.e., x ≈ d i ) on reconnection plane, and the perpendicular sources of free energy are mainly generated in the region with κ ≤ 3. However, for the reconnection cases with finite guide-field Run2 and Run3, the value of the curvature parameter is practically always larger than 3 everywhere. In these cases, few perpendicular sources of free energy are observed, which proves a significantly negative correlation between the electron magnetization and perpendicular sources of free energy.

IV. CONCLUSIONS AND DISCUSSIONS
We carried out 3D PIC code simulations of magnetic reconnection in order to investigate the formation of non-thermal electron beams and EVDFs for plasma conditions suitable to the solar corona. We focused on the EVDF features that can be considered as necessary conditions for microscopic plasma instabilities that could play a role in radio emission processes. This requires the identification of possible source of free energy, namely, positive velocity space gradient(s) in the EVDFs in the direction parallel and perpendicular to the local magnetic field, respectively. For this sake we implemented a machine learning algorithm to fit the reduced 1D EVDFs and searched for their velocity gradient(s). The effects of a guide field on those feature of non-thermal electron beams were also investigated.
Our results are summarized as follows: • Possible parallel (or field-aligned) sources of free energy for streaming-like instabilities, are mainly generated in the separatrices and diffusion region. This kind of source of free energy is determined by positive gradient(s) in the 1D parallel EVDF, i.e., v ·∂ f (v )/v > 0.
Those EVDFs are a necessary condition for the instabilities leading to the plasma emission mechanism, which cause waves at the plasma frequency and possibly their harmonics 40 .
• Possible perpendicular sources of free energy for cyclotron maser instabilities are mainly formed in the diffusion and outflow region near the midplane of reconnection. This kind of source of free energy is determined by the positive gradient in the perpendicular 1D EVDF, i.e., v ⊥ · f (v ⊥ ) > 0. These non-thermal electrons are non-adiabatic and their EVDF is mostly characterized by a crescent-shaped feature. Those EVDFs represent a necessary condition that can cause ECMIs and generate waves at the harmonics of the electron cyclotron frequency 40 .
• As the strength of guide field increases, 1D EVDFs with positive velocity gradient(s) in perpendicular direction to the local magnetic field are less likely to be generated by reconnection.
Dupuis et al. 88 introduced a machine learning method to detect the formation of non-thermal electron beams by comparing the 2D EVDFs, i.e., f (v , v ⊥ ) and f (v ⊥1 , v ⊥2 ), generated by 2D magnetic reconnection with Maxwellian distributed EVDFs. In this study, we extended this method to 1D EVDFs, i.e., f (v ) and f (v ⊥ ) generated by 3D magnetic reconnection. We investigated the positive velocity gradient(s) in the 1D parallel and perpendicular EVDFs in order to determine whether sources of free energy are available at each locations in the reconnection region.
Our results depend on the reliability of the fitting, for which an estimation of its accuracy was not performed. Only a very rough estimate of the uncertainties of the EVDFs fitting by the BGMM method were carried out, in particular to assess the reliability of the positive velocity gradients compared to those fittings. This does not guarantee the accuracy our results, but it provides some kind of positive evidence. A more rigorous goodness-of-fit statistics still needs to be performed to guarantee the reliability of our method but that is deferred to a future work.
An investigation of microscopic plasma instabilities, such as the streaming-like instabilities and ECMIs, due to those nonthermal EVDFs are beyond the scope of this study. We can only speculate here that it is unlikely that this kind of instabilities can be detected in these simulations, because positive gradients in EVDFs are relatively weak, the possible micro-instabilities caused by them are thus relatively weak in comparison with other macroscopic instabilities and plasma flows. We plan to investigate the micro-instabilities and properties of resulting wave emission due to non-thermal EVDFs featuring possible sources of free energy in a forthcoming publication, in which instabilities associated to those EVDFs are separately investigated and they could become 32 more relevant.
A 3D magnetic reconnection configuration is essential for our purposes since EVDFs evolve differently comparing to their counterparts in 2D magnetic reconnection. For example, it is expected that EVDFs with positive gradients in the out-of-plane direction (i.e., along v z direction) will release their energy via streaming-instabilities via unstable waves with a k z vector. This process is not possible in a 2D configuration. This is particularly relevant for reconnection under the influence of a (relatively strong) guide field, because the total magnetic field nearly points to the z−direction (at least near the reconnection midplane), the reconnection electric field will accelerate electrons along this direction. So the evolution of those electron distribution functions will be clearly different in 3D and 2D reconnetion.
The formation mechanism of those non-thermal electron beams and resulting velocity distribution functions is also valuable to investigate in a future work. In particular, for the field-aligned non-thermal electron beams, the reconnection electric field and parallel electric field have significant influence on the electron motion and acceleration 16,50,67 . For the non-thermal electron beams that forms in the direction perpendicular to the local magnetic field, there are several proposed mechanisms (as those discussed in the introduction), including the E × B or gradient-B drifts. It is still not clear which process dominantly affects the electron acceleration in our 3D kinetic magnetic reconnection simulations. For this purpose, a study of dynamic orbits of electrons and ions, their interaction with reconnection electromagnetic fields and their dependence on the strength of guide field is necessary, which is beyond the scope of this work.

V. ACKNOWLEDGEMENTS
We gratefully acknowledge the developers of the ACRONYM code, the Verein zur Förderung kinetischer Plasmasimulationen e.V. and the financial support by the German Science Foundation (DFG), projects MU-4255/1-1 and BU 777/15-1. We also gratefully acknowledge the possibility of using the computing resources of the Max Planck Computing and Data Facility (MPCDF, formerly known as RZG) at Garching and of the Max-Planck-Institute for Solar System Research at Göttingen as well as of the Technical University Berlin, Germany. We also thank the referees for their comments and suggestions that allowed us to improve the presentation of our results. the condition distribution in terms of variable v with a given z thus is The joint distribution function of variables x and z is obtained by the product rule as follows: Then the marginal distribution function with respect to variable v Eq. (A1) is obtained by summing (or integrating) the joint distribution function Eq. (A7) over all possible states of the latent variable z, namely It is helpful to estimate the parameters of the Gaussian mixture involved the variable z by the above marginal distribution function Eq. (A8).
Given a set of data of the velocity, namely, v 1 , v 2 , · · · , v n . Thus each data v i is determined by a specific value of the latent variable z k . It is convenient to introduce the condition probability with respect to variable z using the Bayes's theorem 95 , namely The mixing probability A k is thus the prior probability of z k = 1, and the quantity γ ki can be treated as the corresponding posterior probability for a particular observed v.
Assuming the observed data points v 1 , v 2 , · · · , v n are independent from the distribution, the log likelihood of the Gaussian mixture Eq. (A1) can be constructed by the data v as follows: The parameters of the mixture Gaussian model are thus corresponding to a solution of the maximum likelihood, i.e., max L (v|Φ).
Letting the gradient of the likelihood function to be zero with respect to the mean µ k , the covariance Σ k and the weight (mixture property) A k , yields, for arbitrary k ∈ [1, 2, · · · , K], namely, Given the observed data of v, the posterior probability γ ki appears in the expressions above, thus it is viewed as the responsibility of the kth component relate to the observation v 95 .
Provided a Gaussian mixture model, an expectation-maximization algorithm 104,105 is developed to maximize the likelihood expression L as follows: • Step 1. Initializing the mean µ k , the covariance Σ k and the mixing probability A k . Evaluating the initial value of the log likelihood L .
• Step 3. Maximum step (M-step): re-estimating the value of the mean µ k , the covariance Σ k and the weight (mixture property) A k using Eq. (A11)-Eq. (A13) based on the γ ki obtained in the E step.
• Step 4. Estimating the log likelihood L using Eq. (A10). Then checking the convergence criterion, if the convergence criterion is not satisfied, then go to the E-step iteratively.
In order to compensate for the overfitting of GMM models, it is necessary to assess the appropriate number of sub-populations from the K-component Gaussian mixture models based on the model selection technique named Bayesian information criterion (BIC) 106 , namely where N is the sample size of the observation v, K denotes K components Gaussian mixture model obtained from training of the data by the EM algorithm mentioned above. By introducing this penalty term associated with the number of parameters of the GMM model, this model selection method can automatically determine the appropriate number of sub-populations 95 .