Magnetic reconnection is a ubiquitous phenomenon in turbulent plasmas. It is an important part of the turbulent dynamics and heating of space and astrophysical plasmas. We examine the statistics of magnetic reconnection using a quantitative local analysis of the magnetic vector potential, previously used in magnetohydrodynamics simulations, and now employed to fully kinetic particle-in-cell (PIC) simulations. Different ways of reducing the particle noise for analysis purposes, including multiple smoothing techniques, are explored. We find that a Fourier filter applied at the Debye scale is an optimal choice for analyzing PIC data. Finally, we find a broader distribution of normalized reconnection rates compared to the MHD limit with rates as large as 0.5 but with an average of approximately 0.1.

## I. INTRODUCTION

Most naturally occurring plasmas are observed (e.g., Refs. 1–5) or believed to be (e.g., Refs. 6 and 7) in a turbulent state. Turbulent plasma dynamics create strong gradients in the magnetic field, leading to conditions in which the magnetic topology may change at kinetic scales. This can produce fast, bursty outflows associated with magnetic reconnection.^{8} Magnetic reconnection not only mediates the development of turbulence but also is very efficient at converting magnetic field energy into kinetic energy, both in flows and in thermal degrees of freedom. It is, therefore, of great importance to quantify the role of reconnection in turbulence, an issue addressed here for collisionless kinetic plasma, which is particularly relevant for space and astrophysical applications.

The interplay of reconnection and turbulence can be seen as a two-way interaction: on the one hand, turbulence can establish and sometimes control the conditions for reconnection, and on the other hand, phenomena associated with reconnection, such as exhaust jets, can drive turbulence.^{9–11} It has been established^{12} that MHD turbulence causes magnetic flux tubes to interact and reconnect, leading to a broad statistical distribution of reconnection rates. Similarly, external driving of turbulence can induce rapid^{13} large scale reconnection. Retinò *et al.*^{4} showed how the dynamics of reconnection and turbulence are closely intertwined in Earth's magnetosheath. In this study, we are interested in the former problem, reconnection in the midst of broadband plasma turbulence, extending the Servidio *et al.*^{12} study to the case of collisionless plasma.

Large, noise-free, fully kinetic Eulerian Vlasov simulations of turbulence are at present extremely computationally intensive. These are essentially out of reach of present day computers although the hybrid Vlasov variation, with fluid electrons, is tractable.^{14} A less computationally demanding approach to simulating the fully kinetic model is the particle-in-cell (PIC)^{15} model, which we employ here. As we will see below, to assess statistical properties of reconnection requires identification of physically correct critical points (here, X-points). For the PIC method, this involves not only the possibility of numerical issues associated with use of finite differences in space, but also additional subtleties connected with finite numbers of macro-particles per cell (or, per Debye sphere).^{15} Understanding how this numerical issue can affect the sub-ion scale dynamics of PIC simulations is an important part of accurately simulating collisionless plasma turbulence. Several works have shown that a large number of marco-particles per cell are required to capture these dynamics using an explicit scheme.^{16,17} In this paper, to achieve reliable and physically correct results, we study the effects of changing the counting statistics as well as the effects of filtering the electromagnetic fields. The influence of these variations on different-order statistical quantities will be considered [i.e., the spectrum, the scale dependent kurtosis (SDK), and the number of X-points]. The approach described in Servidio *et al.*^{12} is only applicable in 2D simulations, so these filters are also applied to cuts through the simulations of the partial variance increment (PVI) and the electron agyrotropy^{18,19} to show that the filtering techniques employed in this work could be used in identifying X-points in 3D geometries. Using the appropriate filtering, we analyze the number of X-points generated in a large fully developed decaying turbulence simulation and thereby arrive at a physically reliable assessment of the associated reconnection rates in the many flux tube interactions that occur.

## II. ACCURACY OF TURBULENCE IN SIMULATIONS

Turbulence involves dynamical activity that spans broad ranges of temporal and length scales, a property that makes accurate simulation inherently difficult. Computing dynamics of large structures requires long simulation times, whereas appropriately simulating the much faster activity at small dissipation scales requires fine spatial resolution and well resolved short time scales. In the midst of these complex dynamics, understanding the role of magnetic reconnection at inertial range scales is important for understanding both the topological and the energetic features of the turbulent cascade.^{9,11,22} To quantify the role of reconnection in kinetic turbulence, the first step must be to locate putative reconnection sites. In two dimensional turbulence, this means finding the magnetic X-type critical points, namely saddle points of the magnetic potential, where the in-plane magnetic field is null. One needs also to understand whether the identified X-points are physical, or if they are numerical artifacts. This is not a trivial issue. In 2D MHD, it has been found that the number of X-points generated during turbulence depended on the magnetic Reynold's number (Number of X-points $\u223c\u2009Rem3/2$),^{21} even when care is taken to control for numerical errors. It was also previously shown that inadequate spatial resolution results in the generation of spurious X-points, but that the number of critical points converges to a stable value when the resolution is approximately three times smaller than the Kolmogorov dissipation scale.^{22} It is important to note that although the energy spectra in an under-resolved simulation might display the correct characteristic turbulent features, the value of different higher order statistics (e.g., the scale dependent kurtosis or the number of X-points) can be incorrectly calculated due to the grid scale Gaussian fluctuations.^{22} A sequence of studies Wan *et al.*^{21–23} has led to a reasonable level of understanding of the generation of X-points due to both physical and numerical effects for the case of MHD turbulence. However, open questions remain about proliferation of X-points in kinetic plasmas and in observations.^{24}

The above-mentioned background makes it clear that even though a simulation might qualitatively appear to be properly resolved, the smallest scale dynamics, and the topology of the magnetic field, might not be accurately accounted for. An understanding of this issue is potentially critical for studying kinetic physics in turbulent PIC simulations. PIC is a Monte-Carlo technique for a numerical solution of the Vlasov equation that breaks up the distribution function into a large number of “macro particles” and then traces their trajectory in both real and velocity space. Fields in PIC are susceptible to noise issues as well as entropy conservation violation associated with poor counting statistics.^{15} In a PIC simulation with fixed spatial resolution, these statistics are ultimately tied to the number of macro particles per grid cell (ppc). Simulations with an inadequate number of ppc can lead to the gaussianization of real and velocity space gradients on time scales comparable to those of dynamical interest.^{25} It has previously been shown in turbulent PIC simulations that coherent structures ranging in lengths as small as the electron scale are generated and that the plasma becomes intermittent.^{20,26} The present work is geared towards beginning the process of understanding what controls the proliferation of X-points in kinetic simulations. However, before addressing that question, or even the question of how the simulation resolution affects this process, we need to understand if and how the number of macro particles affects the higher order statistics of smaller scale processes. Here, we study this by using some of the tests outlined in Wan *et al.*,^{22} namely the number of X-points and the behavior of the scale dependent kurtosis.

## III. SIMULATION DETAILS

The simulations were performed using the fully kinetic 2.5D (*X*, *Y* in real space, *V _{X}*,

*V*

_{Y},

*V*

_{Z}in velocity space), electromagnetic PIC code P3D.

^{27}Two sets of simulations are studied. The first set is the Orszag-Tang vortex (OTV) setup,

^{28–31}performed in a doubly periodic domain of $(10.24di)2$ with a grid spacing and time step of $\Delta x=0.02di\u22480.73\lambda D$ and $\Delta t=0.0015\Omega ci\u22121$, respectively. The simulations use an artificial mass ratio of $mi/me=25$ and speed of light $c/cA0=20$ where $cA0$ is the Aflvén speed based on a reference density and magnetic field of 1 in simulation units. The simulation begins with a uniform density, a mean out-of-plane magnetic field of 1, and in-plane magnetic and velocity fluctuation of r.m.s. strength 0.2. The ion and electron temperatures are initially uniform with $Ti=Te=0.3micA02$ and a total plasma beta of $\beta =1.2$. Five simulations were performed, varying the number of particles per cell (ppc =12, 50, 200, 800, 3200), but keeping all other parameters fixed. The simulations run until just after the peak of dissipation (peak in $|Jz2|$) at $15\Omega ci\u22121$. This is slightly less than 2

*τ*where

_{nl}*τ*is the large scale nonlinear time (eddy turnover time) of the system. This system at this time is ideal to test the reconnection site finding algorithm, as the number of physical reconnection sites can be visually identified and counted.

_{nl}The second simulation is $(102.4di)2$ in size with a grid spacing of $\Delta x=0.0125di$ and $\Delta t=0.0025\Omega ci\u22121$ and $c/cA0=30$. The simulation was initialized with an MHD like initial condition with “Alfvénic perturbations” in the in-plane *B* and *V* with an initially specified spectrum. The initial r.m.s of the *B* and *V* fluctuations have a strength of 1. There is an out-of-plane mean field of 5 and density of 1. Initially the density and temperatures are uniform and have $Ti=Te=1.25micA02$ with a total plasma beta of $\beta =0.2$. It is interesting to note that in this simulation, the grid scale relative to the Debye length is about a factor of 3 smaller than the OTV simulations. This simulation was originally performed to study von-Kárman energy decay in kinetic plasmas and this grid spacing and time step were required for robust energy conservation. However, by reducing the grid spacing in this simulation with only a few hundreds of particles per cell, the number of effective macro particles per Debye circle is the same order as that of OTV5. More details about the simulation can be found in Wu *et al.*^{32} Three different time snapshots of the simulation at times $t=206.25,250.0,292.5\u2009\Omega ci\u22121$ (where $\Omega ci\u22121$ is based on the mean out-of-plane magnetic field of 5) are used in this study for better statistics. Details for all of the simulations are presented in Table I along with the number of macro-particles in a Debye circle $N\lambda D$ where $N\lambda D=\pi \lambda D2\xd7ppc/(\Delta x)2$.

Name . | l_{x} = l_{y}
. | $\Delta x$ . | $\Delta t$ . | $ccA$ . | $mime$ . | ppc . | $N\lambda D$ . | $\beta i=\beta e$ . |
---|---|---|---|---|---|---|---|---|

OTV1 | 10.24 | 0.02 | 0.0015 | 20 | 25 | 12 | 71 | 0.6 |

OTV2 | 10.24 | 0.02 | 0.0015 | 20 | 25 | 50 | 295 | 0.6 |

OTV3 | 10.24 | 0.02 | 0.0015 | 20 | 25 | 200 | 1178 | 0.6 |

OTV4 | 10.24 | 0.02 | 0.0015 | 20 | 25 | 800 | 4712 | 0.6 |

OTV5 | 10.24 | 0.02 | 0.0015 | 20 | 25 | 3200 | 18850 | 0.6 |

PWu1 | 102.4 | 0.0125 | 0.0025 | 30 | 25 | 400 | 11170 | 0.1 |

Name . | l_{x} = l_{y}
. | $\Delta x$ . | $\Delta t$ . | $ccA$ . | $mime$ . | ppc . | $N\lambda D$ . | $\beta i=\beta e$ . |
---|---|---|---|---|---|---|---|---|

OTV1 | 10.24 | 0.02 | 0.0015 | 20 | 25 | 12 | 71 | 0.6 |

OTV2 | 10.24 | 0.02 | 0.0015 | 20 | 25 | 50 | 295 | 0.6 |

OTV3 | 10.24 | 0.02 | 0.0015 | 20 | 25 | 200 | 1178 | 0.6 |

OTV4 | 10.24 | 0.02 | 0.0015 | 20 | 25 | 800 | 4712 | 0.6 |

OTV5 | 10.24 | 0.02 | 0.0015 | 20 | 25 | 3200 | 18850 | 0.6 |

PWu1 | 102.4 | 0.0125 | 0.0025 | 30 | 25 | 400 | 11170 | 0.1 |

In both the Orszag-Tang simulations and the Wu *et al.*^{32} simulation, the ion particle velocities are set to the initial bulk flow and the electron velocities are set to supply the current. The ion and electron densities and temperatures are initially constant and uniform and so pressure balance is not initially enforced. Similarly, the initial electric field was set to 0 and not $E=\u2212v\xd7B$. In the first few time steps of each simulations, these imbalances are corrected through the generation of transient waves throughout the system which decay on a time scale much shorter than a non-linear time.

## IV. ANALYSIS METHODS

### A. Noise reduction

As alluded to above, any attempt to identify and tabulate X-points, active reconnection sites, and associated reconnection rates in two dimensional turbulence models must confront both physical and numerical issues, and importantly, must distinguish between them. In this regard, a useful observation is that the onset of numerical resolution issues in MHD are signified by the appearance of phase errors and “gaussianization” of noisy fluctuations at the smallest scales.^{25,36} Apart from proper resolution to determine the physical reconnection rate, an added complication is that at high magnetic Reynolds numbers, the number of X-points increases due to *physical* effects,^{21} as the most intense current sheets themselves become turbulent and unstable at high local Reynolds numbers.

For the OTV configuration in MHD and at moderate large scale Reynolds numbers (say, $Rm<1000$), one would not expect physical secondary islands to form dynamically, using the well-known empirical criterion (threshold of $Rm\u223c10\u2009000$) suggested by Biskamp.^{34} As far as we are aware, for kinetic plasma, a similar threshold is obtained (see, e.g., Daughton *et al.*^{35}) with the system size determining the effective Reynolds number.^{36}

Carrying the above ideas over to the present case of kinetic simulation, we expect that the number of reconnection sites in a moderate size (Reynolds number) OTV simulation at the peak of current to be ∼4. However, if a snapshot is taken from a PIC simulation at this time and analyzed directly, the number of X-points may be found to be as large as ∼10^{4}, depending on parameters. It transpires, as we will show below, that particle noise has to be removed before the physical X-points can be unmasked from the “noise X-points.” For the results discussed here, we will consider three different approaches to reducing the particle noise:

*Increasing the number of particles per cell:*Ideally PIC simulations should be run with as many particles per cell (ppc) as possible. However, computational limits often restrict the choice to a few hundred for typical simulation. The particle noise reduces as $1/N$ where*N*is the number of ppc. In this study, we varied ppc from 12 to 3200.*Time averaging:*A technique for reducing particle noise in PIC simulation output data is to time average the data over time scales much shorter than those of dynamical interest. This averaging can reduce plasma oscillations associated with finite particle number charge imbalances. In this study, we employed time averaging of the data over a period of $3\omega pe\u22121$, which is roughly 200 times smaller than the typical nonlinear time for the $(102.4di)2$ simulation.*Gaussian/Fourier Spatial Filtering:*A standard technique used in studying the scale to scale transfer of energy is real space^{37,38}or sharp Fourier space filtering.^{22,43}In this work, we will employ a Gaussian filter for our real space filtering, defined as(1)$g(Bi,j)=\u2211k,lBk,l2\pi \lambda d/\Delta xe\u2212(i\u2212k)2+(j\u2212l)22(\lambda d/\Delta x)2,$where each i, j grid point is smoothed by neighbors with

*k*and*l*truncated at ±4. For Fourier filtering, we apply a sharp cutoff in Fourier space at a scale of interest*k*. For both Gaussian and Fourier filtering, we applied the filters at the scale_{s}*λ*for each simulation._{d}

We apply these techniques in concert with the development of the specific diagnostics that will be used to analyze reconnection properties. Subsequently, we turn towards the physical properties revealed by the statistics of reconnection. It should be noted that out of the three noise reduction procedures mentioned above, only the increase in number of particles per cell will reduce the effects of noise in the time evolution of the simulation. Time averaging and Gaussian/Fourier filter are post-processing procedures that do not remove the problems associated with particle noise during the time evolution of the simulation.

### B. Identifying critical points in 2D PIC turbulence

To identify the critical points in 2.5D simulations, we follow the procedure devised by Servidio *et al.*^{12,33} We examine the 1st and 2nd derivatives of the magnetic vector potential. Magnetic null points are identified as zeros of the 1st derivatives. At the null points, we consider the sign of the product of the eigenvalues of the second derivative (Hessian) Matrix, $Mi,j=\u22022a\u2202x\u2202y$. Two negative eigenvalues indicate a maximum of the magnetic vector potential, and two positive eigenvalues indicate a minimum. Null points where the product of the eigenvalues is negative are identified as saddle points, and possible physical X-points.

## V. EFFECTS OF TIME AVERAGING AND SPATIAL FILTERING

We start by comparing the results of smoothing on turbulent statistics. First, we examine power spectra of electromagnetic fields. Figure 1 shows the omnidirectional spectra of the magnetic field and out-of-plane electric field for the OTV simulations. The top panel shows the magnetic field spectra for varying number of particles per cell. The result, as expected, is that the particle noise goes down with increasing number of ppc, and this is reflected in the spectra. The discrepancies in the magnetic spectrum for different ppc simulations can be seen at scales as large as $kde\u223c2$. The insert in the top panel shows the power in the magnetic field at the Debye scale as a function of ppc. As the noise in a variable is expected to go down as $N$ with a number of particles, the energy is expected to go down linearly. This is indeed the case because the slope of the line in the figure is −1.

The lower two panels of Fig. 1 show the magnetic and out-of-plane electric field (*E _{z}*) spectra for the largest ppc run for time averaging as well Gaussian and Fourier filtering. Time averaging the magnetic field only changes the spectrum very close to the Debye scale, reducing the noise by a factor of a few. The Gaussian filtered data also match the unprocessed data up to the Debye scale. The Fourier filter by definition matches identically until its sharp cutoff at the Debye scale.

The most drastic effect of time averaging is visible on *E _{z}*. Time averaging on a few plasma frequency time scales reduces spectral power even at scales ∼

*d*. This implies that the time averaging has a drastic effect on electric fields. The adverse effect of time averaging on the electric field is also apparent in the spectra of the large PIC run as shown in Fig. 2. Broadband reduction of the electric field by time averaging can suppress the computed estimates of scale to scale transfer of energy.

_{i}^{40}An understanding of whether this is a consequence of the numerical method, or if it relates to high frequency physics remains to be investigated.

To look at the effects of smoothing on statistics of turbulent fluctuations, we plot probability distribution functions (PDFs) for the out-of-plane electric current density (*J _{z}*) in the OTV simulations. Figure 3 shows the PDFs for various ppc and smoothing techniques. The smallest ppc simulations are dominated by random fluctuations associated with poor counting statistics and the corresponding PDF of current density resembles a Gaussian PDF. However, when the ppc = 50 simulation is smoothed with the Fourier filter, the shape of the current density PDF converges towards the ppc = 3200 case [see Fig. 3(c)].

Next, we turn our attention to the scale dependent kurtosis (SDK), one of the common measures of intermittency in a turbulent system.^{41} The scale dependent kurtosis is defined as

where the angled brackets denote a spatial average and $\delta Bx(r)$ is the increment of the magnetic field in the x-direction, defined as $\delta Bx(r)=Bx(x+r)\u2212Bx(x)$. In principle, the scale dependent kurtosis can be calculated for the increment of any field and any vector component; however, in this work, we elect to only present the scale dependent kurtosis of *B _{x}*. It should be kept in mind that the OTV simulation, being relatively small in size, has a very small effective Reynolds number and hence cannot have large kurtosis.

^{42}Moreover, the large scale inhomogeneity of OTV makes the kurtosis drop below 3 at larger scales. However, the scale dependent kurtosis can still be computed and its convergence to a stable value for different noise reduction techniques can be studied.

The top panel of Fig. 4 shows the scale dependent kurtosis of *δB _{x}* for different ppc cases of the OTV simulation. At large scales, the scale dependent kurtosis matches for all simulations. However, at smaller scales, higher ppc simulations have significantly larger kurtosis. Particle noise in lower ppc simulations evidently randomizes (gaussianizes) the smaller scale structures, decreasing the kurtosis at smaller scales. The bottom panel of Fig. 4 shows the scale dependent kurtosis of

*δB*for the ppc 3200 OTV run for different smoothing techniques. The problem of small scale gaussianization is alleviated by almost all the processing techniques and the kurtosis saturates to a constant value at the smallest scales.

_{x}^{43}

The scale dependent kurtosis for the larger PIC run (Pwu1 in Table I), however, tells a slightly different story. This simulation's larger size allows a greater separation between the energy containing scales and the “dissipative” scales (i.e., kinetic ion and electron scales) and thus has a larger effective Reynolds number. The larger Reynolds number allows the generation of stronger small scale coherent structures, and thus the energy cascade of this simulation more closely resembles the energy transfer in the turbulent MHD limit.^{42} Figure 5 shows the scale dependent kurtosis for *δB _{x}* for different averaging/filtering techniques. In all cases, the scale dependent kurtosis matches very well down to ∼0.5

*d*, at which point the Gaussian filtered data start to diverge from the other curves. Unfiltered, Fourier filtered, and time averaged data match with each other down to scales ∼0.5

_{i}*d*. This implies that a Gaussian filter applied to field data cannot capture the smallest scale intermittent structures as well as a Fourier filter.

_{e}Although in 2D simulations, it is straightforward to identify x-points by examining the magnetic scalar potential and looking for critical points, identifying x-lines in a 3D system is a significantly more challenging problem since a single magnetic potential can no longer be defined. However, magnetic reconnection is a non-linear, boundary layer process and generates distinct features that can help to characterize and identify it in three dimensions. The small scale structure of reconnection implies that it must occur in regions of intermittency. One simple measurement of intermittent and thus one potential identifying measurement of reconnection is the partial variance increment (PVI) of the magnetic field, defined in Greco *et al.*^{18} as

and a scalar measurement of the electron pressure tensors agyrotropy normalized over the unit interval, defined in Swisdak^{19} as

where *P _{e}* has been rotated into a frame where the 1st axis is parallel to the local magnetic field and the 2nd and 3rd axes are rotated such that the 2,2 and 3,3 diagonal components of the tensor are equal (i.e.,

*P*

_{e11}=

*P*

_{e∥}and $Pe22=Pe33=Pe\u22a5$).

The PVI is calculated for a single horizontal cut in the center of the OTV simulations for a separation of $1di$ and is shown in Fig. 6. This location has the benefit of crossing through two current sheets at the edge and the center of the simulation domain. In Fig. 6(a), the effects of ppc on the PVI is shown. Although ppc reduces the noisy fluctuations of this quantity, the general structure of the signal is unchanged. In Fig. 6(b), we apply the three filtering techniques to the PVI and find that they are mainly in accord with one another. The difference between the 3200 ppc simulation and the Gaussian or Fourier filtered data, however, is almost undetectable. This suggests that PVI calculated across an ion inertial length is insensitive to random fluctuations associated with ppc. In this simulation, this makes sense because most of the noise associated with PIC counting statistics should occur at Debye length scales and PVI is mainly sensitive to signals at the chosen lag, which here is *d _{i}*, well separated from the Debye length.

The nonlinear processes that occur in magnetic reconnection generate non-zero off diagonal terms in the electron pressure tensor which ultimately result in electron agyrotropy. Agyrotropy been proposed as an identifier of reconnection,^{44} this quantity requires an accurate measurement of the off diagonal elements of the second moment of the distribution function, which are even more susceptible to particle noise than the magnetic fields or the current densities. This is evident in Fig. 7 where a cut of the electron pressure agyrotropy is taken along X at $Y=5.12di$. In Fig. 7(a), the agyrotropy of the 50 ppc OTV simulation is entirely dominated by noise, and no clear signal is evident. It is only by increasing the ppc or filtering [Fig. 7(a)] can the signal be seen. The peak values of the agyrotropy are an order of magnitude smaller than the noise amplitude in the 50 ppc simulation. Again, all three filtering techniques provide similar refinement of the signal. Finally, in Fig. 7(c), the Fourier filtered 50 ppc simulations are compared with the Fourier filtered and unfiltered 3200 ppc simulation. Some of the strongest agyrotroic features survive the filtering process (e.g., the features near $X=0.5di$ and $X=10di$ in Fig. 7). However, there are still regions where the noise associated with poor counting statistics is stronger than the signal of interest. This analysis reaffirms that higher order moments of the distribution function are sensitive to counting statistics and it shows that even the agyrotropy from simulations with many particles per cell can benefit from filtering.

Finally, we study the number of identified X-points for simulations with different numbers of ppc as well as the effects of time averaging, Gaussian filtering and Fourier filtering. For this purpose, we choose to work with the OTV near the peak of root mean square electric current density. At this time, the number of reconnection sites can be visually counted to be 4 as shown in Fig. 8. If the critical point finding algorithm is applied directly to the unprocessed simulation data, the noise artificially introduces a large number of minima maxima and saddle points. Figure 9 shows the number of critical points as a function of number of particles per cell. For very small number of particles, the number of critical points is ∼10^{4}. The number of critical points decreases as a power law with increasing ppc. However, the number of particles required to achieve convergence to the physical number of critical points is very large (>10^{5}). Time averaging reduces the number of X-points significantly (an order of magnitude less). However, each time averaged simulation has many more than 4 x-points and follows a decreasing power law with increasing ppc. On the other hand, Fourier filter and Gaussian filter (both at the Debye scale) remove the particle counting noise to reveal the physical critical points even for rather small ppc. Even for ppc = 12, the number of identified critical points is less than a factor of 3 too large.

Combining the results of analysis of turbulence quantities, and the results of critical point finding, we can conclude the following:

Ideally one would run a simulation with a large number of particles per cell. However, the number of particles required to reduce noise significantly and hence capture physical reconnection sites is prohibitively large. This would require more computational time and so restricts the size of the simulation.

Time averaging, although a common and simple technique fails to capture the physical reconnection sites and also adversely affects the electric fields for analysis purposes.

Gaussian filtering appears to capture the physical reconnection sites and spectra very well but has a slight negative effect on the scale dependent kurtosis.

Fourier filtering appears to reduce the noise effects while minimally interfering with the physical effects discussed above.

Based on the above considerations, we conclude that the optimum method to analyze PIC simulations for reconnection studies is via the application of a Fourier filter as was done by Wan *et al.*^{20} We now identify the reconnection sites in the large PIC simulation to study their statistics.

## VI. RECONNECTION IN TURBULENCE

To examine the statistics of reconnection in turbulence, we analyze the larger 2.5D turbulent PIC simulation again carried out with P3D (PWu1 in Table I). We collect statistics from three different times in the same simulation (*t* = 206.25, 250.0 and $292.5\u2009\Omega ci\u22121$). For each time sample, we apply a Fourier filter to the magnetic field data for $k\lambda d>1$ ($\lambda d=0.0375di$). The results from the topological analysis of the filtered data for these three sets of output can be found in Table II.

Time $\Omega ci\u22121$ . | Min . | Max . | X-points . | Total . |
---|---|---|---|---|

206.25 | 110 | 116 | 226 | 452 |

250.00 | 144 | 165 | 309 | 618 |

292.50 | 159 | 159 | 318 | 636 |

Time $\Omega ci\u22121$ . | Min . | Max . | X-points . | Total . |
---|---|---|---|---|

206.25 | 110 | 116 | 226 | 452 |

250.00 | 144 | 165 | 309 | 618 |

292.50 | 159 | 159 | 318 | 636 |

Across the three times, we find 853 critical saddle points. The first snapshot of the three is shown in Fig. 10. The first panel shows the whole domain and the X-points identified within it. At first inspection, there are only a handful of locations that resemble the classical picture of reconnection synonymous with PIC reconnection simulations (long straight oppositely directed field lines with large coherent current sheets e.g., $x=77\u2009di,\u2009y=22\u2009di$). There are, however, numerous regions with many X-points clustered together. These frequently correspond to “secondary islands.” The second panel shows the zoomed-in region denoted by the black square in the first panel, and it becomes clear that the X-points marked in the simulation do in fact correspond to critical points of the magnetic vector potential.

Because of the 2D nature of this simulation, we know that reconnection must occur in the X-Y plane, and so the reconnecting electric field must point out-of-plane (Z direction). So for each identified saddle point, we interpolate the Fourier filtered, out-of-plane electric field. We generate the PDFs for *E _{z}* and $|Ez|$ shown as the black triangles in the two panels of Fig. 11. We normalize

*E*to the root mean square of the in plane magnetic field (note, in this simulation $Bx,y,rms\u22481$) and so Fig. 11 can be interpreted as the PDF of reconnection rates in our simulation. From this, we find the reconnection rates in our kinetic simulation can be as large as 0.5, with an average magnitude of about 0.1. Figure 11 also includes cyan squares that represent the reconnection rates found using the same procedure applied to MHD. These data are from Fig. 10 in Wan

_{z}*et al.*

^{22}Although the MHD and the PIC results have some similarities, it is clear that the PIC distribution is broader than the MHD distribution, in both the range of reconnection rates and the shape of the PDF. It is also worth noting that the average magnitude in the MHD case (0.044) is a little more than a factor of 2 less than the PIC case (0.10). This result reaffirms the idea that turbulence can potentially lead to enhanced reconnection rates and is in apparent agreement with the idea that reconnection rates in kinetic plasmas tend to be larger than in MHD and of order 0.1.

^{45}It is also worth while noting that the present (filtered) PIC results for the distribution of reconnection rates is quite similar to the distribution of rates seen in 2.5D Hybrid Vlasov (HVM) simulation, as seen in Fig. 5 of Ref. 15. The distributions agree mainly at low values of reconnection rate, but the present PIC results, which include kinetic electron physics, extend higher reconnection rate values.

## VII. CONCLUSIONS

In this work, we have begun to examine the statistics of x-type critical points (X-points) in fully kinetic 2.5D PIC simulations. This work extends the procedures applied to MHD simulations^{12,33} to PIC. We find that noise fluctuations in the magnetic field associated with the counting statistics corresponding to the number of particles per cell (ppc) result in a noisy magnetic vector potential, and ultimately spurious numbers of X-points. Increasing the number of ppc helps to lessen this effect. Other noise reduction techniques examined include post processing the simulation magnetic field data by using either time averaging over a plasma frequency, or spacial filtering using a Fourier or Gaussian method. We showed how each of these affected different statistics of the turbulence, including the omni-directional energy spectrum, the probability density function, the scale dependent kurtosis (SDK) and the number of X-points. We find that the number of ppc would need to be increased several orders of magnitude to have accurate enough field data to identify the correct number of X-points (a prospect that is currently computationally intractable), whereas time averaging significantly alters the spectrum of the electric field. However, imposing a spatial filter at the Debye scale stabilizes the number of X-points, regardless of the number of ppc, while not dramatically changing the spectrum or the scale dependent kurtosis. The above tests were carried out using kinetic PIC simulations for both Orszag Tang vortex configurations^{30} and larger broad band turbulence simulations.^{32}

A clear limitation of this x-point identifying procedure given the importance of including the third dimensions in turbulence (e.g., Chen *et al.*,^{46} Howes^{47}) is its applicability to only 2D simulations. This approach to finding x-type critical points from the magnetic scalar potential cannot be used in 3D data. Although this work does not directly address identifying x-lines in 3D simulations, it does show that the filtering techniques applied in the 2D case can be used to better resolve indicators of 3D reconnection. In particular, quantities that depend on the second moment of the Boltzmann distribution are very susceptible to noise from poor counting statistics, yet this effect can be mediated by the Debye scale Fourier filtering.

With these filtering techniques, we identify the X-points in a large decaying turbulent PIC simulation, and we calculate the reconnection rate at each of these points. We find that the magnitude of the reconnection rates range between 0 and 0.5 in standard normalized Alfvén units and have an average value 0.1 which is approximately a factor of 2 larger than the MHD result.^{22} Note that reconnection rates of this magnitude are ordinarily associated with a “fast reconnection” process although in the present case, the normalization is with respect to the global r.m.s fields rather than the local upstream quantities. The PDF of the PIC reconnection rates has a broader shape than in MHD which implies that there are a larger fraction of X-points that have large reconnection rates. This is consistent with idea that reconnection rates can be boosted by kinetic plasma processes.^{45}

Although this work demonstrates a clear procedure to identify X-points in a kinetic PIC simulation, it does not address questions about how varying the number of particle per cell affects the proliferation of X-points during a simulation. It is clear in this work that the number of ppc is an important quantity for the accuracy of the fields, and it has been shown in the MHD case that spatial resolution has a dramatic effect on the number of X-points generated during a simulation.^{22} It is clearly plausible that PIC simulations could be susceptible to a similar issue related to poor counting statistics and spatial resolution. This is an important question for the PIC modeling community that should be addressed in greater detail in the future. In addition, there remain important questions regarding the *physical*, rather than numerical proliferation of X-points in turbulent plasma. This phenomenon has been previously demonstrated for MHD,^{21} where for large systems, at high magnetic Reynold's number, the expected number of islands can increase dramatically, even when numerical inaccuracies are carefully controlled. The turbulent proliferation of reconnection sites is clearly a nonlinear dynamical effect, but is likely related to the family of linear instabilities known as plasmoid instabilities^{48,49} that are initiated from equilibrium field configurations. Accurately tracking the putative physical proliferation of X-points in kinetic turbulence would require simultaneously computing the dynamics of a large, high effective Reynolds number plasma, while respecting the numerical issues we have explored in the present analysis. In this regard, the present study is only an initial step in trying to answer these larger questions regarding reconnection in kinetic plasma turbulence.

## ACKNOWLEDGMENTS

Research is supported by NSF AGS-1460130 (SHINE), and NASA Grant Nos. NNX15AW58G, NNX14AI63G (Heliophysics Grand challenge Theory), NNX08A083GMMS, the MMS mission through grant NNX14AC39G, and the Solar Probe Plus science team (IOSIS/Princeton subcontract No. D99031L). Simulations and analysis were performed at NCAR-CISL and at NERSC. The data used are listed in the text, references, and are available by request.