Exploring the Statistics of Magnetic Reconnection X-points in Kinetic Particle-in-Cell (PIC) Turbulence

TL;DR

This paper analyzes magnetic reconnection in turbulent plasmas using kinetic PIC simulations, revealing broader reconnection rate distributions and identifying optimal noise reduction techniques for data analysis.

Contribution

It introduces a generalized analysis method for PIC data and demonstrates the effectiveness of Fourier filtering at the Debye scale for reconnection studies.

Findings

Reconnection rates can reach up to 0.5 in PIC turbulence.

Fourier filtering at the Debye scale optimizes noise reduction.

Reconnection rate distribution is broader than in MHD simulations.

Abstract

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 generalized to fully kinetic 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. Finlay, 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.

Tables

Table 1. Table 1: Parameters for different turbulent simulations. The simulation length in d i subscript 𝑑 𝑖 d_{i} ( l x subscript 𝑙 𝑥 l_{x} , l y subscript 𝑙 𝑦 l_{y} ), the grid spacing in d i subscript 𝑑 𝑖 d_{i} ( Δ ​ x Δ 𝑥 \Delta x ), the time step ( Ω c ​ i − 1 superscript subscript Ω 𝑐 𝑖 1 \Omega_{ci}^{-1} based on a uniform magnetic field of B 0 = 1 subscript 𝐵 0 1 B_{0}=1 ), the speed of light c c A 𝑐 subscript 𝑐 𝐴 \frac{c}{c_{A}} , the ion to electron mass ratio ( m i m e subscript 𝑚 𝑖 subscript 𝑚 𝑒 \frac{m_{i}}{m_{e}} ), the number of particles per cell (ppc) and the number of particles per Debye circle ( N λ D = π ​ λ D 2 × ppc / ( Δ ​ x ) 2 subscript 𝑁 subscript 𝜆 𝐷 𝜋 superscript subscript 𝜆 𝐷 2 ppc superscript Δ x 2 N_{\lambda_{D}}=\pi\lambda_{D}^{2}\times\rm{ppc}/(\Delta x)^{2} ).

Name	$l_x=l_y$	$\Delta x$	$\Delta t$	$\frac{c}{c_A}$	$\frac{m_i}{m_e}$	ppc	$N_{{\lambda }_D}$
OTV1	10.24	.02	.0015	20	25	12	71
OTV2	10.24	.02	.0015	20	25	50	295
OTV3	10.24	.02	.0015	20	25	200	1178
OTV4	10.24	.02	.0015	20	25	800	4712
OTV5	10.24	.02	.0015	20	25	3200	18850
PWu1	102.4	.0125	.0025	30	25	400	11170

Table 2. Table 2: Number of critical points at each different time.

Time ${\Omega }_{ci}^{-1}$	Min	Max	X-points	total
206.25	110	116	226	452
250.00	144	165	309	618
292.50	159	159	318	636

Equations

$$g(B_{i,j})=\sum _{k,l}\frac{B_{k,l}}{\sqrt{2\pi }{\lambda }_d/\Delta x}{e}^{-\frac{(i-k{)}^2+(j-l{)}^2}{2({\lambda }_d/\Delta x{)}^2}}$$

$$\mathrm{SDK}(\mathrm{r})=\frac{⟨\delta {\mathrm{B}}_{\mathrm{x}}(\mathrm{r}{)}^4⟩}{{⟨\delta {\mathrm{B}}_{\mathrm{x}}(\mathrm{r}{)}^2⟩}^2}$$

Full text

††preprint: APS/123-QED

I Introduction

Most naturally occurring plasmas are observed (e.g., 1; 2; 3; 4; 5) or believed to be (e.g., 6; 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 Yamada, Kulsrud, and Ji (2010). 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 (Matthaeus and Lamkin, 1986a; Lapenta, 2008; Matthaeus and Velli, 2011). It has been established (Servidio et al., 2009) 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 (Lazarian and Vishniac, 1999) large scale reconnection. Retinò et al. (2007) 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. (2009) 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 Servidio et al. (2015). A less computationally demanding approach to simulating the fully kinetic model is particle-in-cell (PIC) Birdsall and Langdon (1985) 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) (Birdsall and Langdon, 1985). 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 (Franci et al., 2015; Camporeale et al., 2016) 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 and the number of X-points). 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, while 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 (Matthaeus and Lamkin, 1986b; Wan et al., 2012; Matthaeus and Velli, 2011). To quantify the role of reconnection in kinetic turbulence a 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, were 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 $\sim\ Re_{m}^{3/2}$) (Wan et al., 2013), 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 converge to a stable value when the resolution is approximately 3 times smaller than the Kolmogorov dissipation scale (Wan et al., 2010a). It is important to note that while 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 (Wan et al., 2010a). A sequence of studies Wan et al. (2009, 2010a, 2013) 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 (Greco et al., 2016).

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 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 (Birdsall and Langdon, 1985). 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 (Montgomery and Nielson, 1970). 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 (Wan et al., 2012; Karimabadi et al., 2013). 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. (2010a) 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 Zeiler et al. (2002). Two sets of simulations are studied. The first set is the Orszag-Tang vortex (OTV) setupOrszag and Tang (1979); Dahlburg and Picone (1989); Parashar et al. (2009); Vasquez and Markovskii (2012), performed in a doubly periodic domain of $(10.24d_{i})^{2}$ with a grid spacing and time step of $\Delta x=.02d_{i}$ and $\Delta t=.0015\Omega_{ci}^{-1}$ respectively. The simulations uses an artificial mass ratio of $m_{i}/m_{e}=25$ and speed of light $c/c_{A0}=20$ where $c_{A0}$ 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 $T_{i}=T_{e}=.3m_{i}c_{A0}^{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 were performed to times just after the peak of dissipation (peak in $|J_{z}^{2}|$) at $15\Omega_{ci}^{-1}$. This is slightly less than $2\tau_{nl}$ where $\tau_{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.

The second simulation is $(102.4d_{i})^{2}$ in size with a grid spacing of $\Delta x=0.0125d_{i}$ and $\Delta t=0.0025\Omega_{ci}^{-1}$ and $c/c_{A0}=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. This simulation was preformed to study von-Kárman energy decay in kinetic plasmas and more details about the simulation can be found in Wu et al. (2013). Three different time snapshots of the simulation at times $t=206.25,250.0,292.5\ \Omega_{ci}^{-1}$ (where $\Omega_{ci}^{-1}$ 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 1 along with the number of macro-particles in a Debye circle $N_{\lambda_{D}}$ where $N_{\lambda_{D}}=\pi\lambda_{D}^{2}\times\rm{ppc}/(\Delta x)^{2}$

IV Analysis Methods

IV.1 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 Servidio et al. (2010); Wan et al. (2010b). 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 Wan et al. (2013), 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, $R_{m}<1000$), one would not expect physical secondary islands to form dynamically, using the well-known empirical criterion (threshold of $R_{m}\sim 10,000$) suggested by Biskamp (1986). As far as we are aware, for kinetic plasma a similar threshold is obtained (see, e.g. Daughton et al. (2009)) with the system size determining the effective Reynolds number Matthaeus et al. (2005).

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 $\sim 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 $\sim 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 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/\sqrt{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}^{-1}$, which is roughly 200 times smaller than the typical nonlinear time (for the $(102.4d_{i})^{2}$ simulation.

• •

Gaussian/ Fourier Spatial Filtering: A standard technique used in studying the scale to scale transfer of energy is real space Germano (1992); Eyink (2005) or sharp Fourier space filtering Aluie and Eyink (2009); Wan et al. (2012) In this work we will employ a Gaussian filter as for our real space filtering, defined as:

[TABLE]

where each i, j grid point is smoothed by neighbors with $k$ and $l$ truncated at $\pm 4$. For Fourier filtering, we apply a sharp cutoff in Fourier space at a scale of interest $k_{s}$. For both Gaussian and Fourier filtering we applied the filters at the scale $\lambda_{d}$ for each simulation.

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.

IV.2 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. (2009, 2010). 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, $M_{i,j}=\frac{\partial^{2}a}{\partial x\partial y}$. 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. Fig  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 particles, and this is reflected in the spectra. The discrepancies in magnetic spectrum for different ppc simulations can be seen at scales as large as $kd_{e}\sim 2$. The insert in the top panel shows the power in magnetic field at the Debye scale as a function of ppc. As the noise in a variable is expected to go down as $\sqrt{N}$ with 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 spectrum very close to the Debye scale, reducing the noise by a factor of a few. The Gaussian filtered data also matches the unprocessed data up to the Debye scale. The Fourier filter by definition matches identically until its sharp cutoff at 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 $\sim d_{i}$. This implies that the time averaging has 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 Yang et al. (tted). 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. Fig 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. 3c).

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

[TABLE]

where the angled brackets denote a spatial average and $\delta B_{x}(r)$ is the increment of the magnetic field in the x-direction, defined as $\delta B_{x}(r)=B_{x}(x+r)-B_{x}(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 can not have large kurtosis Parashar et al. (2015). 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 $\delta 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 $\delta B_{x}$ for the ppc 3200 OTV run for different smoothing techniques. The problem of small scale gaussianization are alleviated by almost all of the processing techniques and the kurtosis saturates to a constant value at the smallest scales Wan et al. (2015).

The scale dependent kurtosis for the larger PIC run (Pwu1 in Table 1), 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 (Parashar et al., 2015). Fig 5 shows the scale dependent kurtosis for $\delta B_{x}$ for different averaging/filtering techniques. In all cases, scale dependent kurtosis matches very well down to $\sim 0.5d_{i}$, at which point the Gaussian filtered data starts to diverge from the other curves. Unfiltered, Fourier filtered and time averaged data sets match with each other down to scales $\sim 0.5d_{e}$. This implies that a Gaussian filter applied to field data cannot capture the smallest scale intermittent structures as well as a Fourier filter.

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 mean square electric current density. At this time the number of reconnection sites can be visually counted to be 4 as shown in Fig. 6. If the critical point finding algorithm is applied directly to the unprocessed simulation data, the noise introduces artificially large number of minima and maxima. Fig  7 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 $\sim 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 brings down the number of X-points significantly (an order of magnitude less) but the number still is significantly large and these cases also follow a power-law decrease 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 three 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 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. (2012). 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 1. We collect statistics from three different times in the same simulation ($t=206.25$, $250.0$ and $292.5\,\Omega_{ci}^{-1}$) For each time sample we apply a Fourier filter to the magnetic field data for $k\lambda_{d}>1$ ($\lambda_{d}=.0375d_{i}$). The results from the topological analysis of the filtered data for these three sets of output can be found in Table 2.

Across the three times we find 853 critical saddle points. The first snapshot of the three is shown in Fig. 8. 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\,d_{i}$, $y=22\,d_{i}$). There are, however, numerous regions with many X-points clustered together. These frequently corresponds 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 $|E_{z}|$ shown as the black triangles in the two panels of Fig 9. We normalize $E_{z}$ to the root mean square of the in plane magnetic field (note, in this simulation $B_{x,y,rms}\approx 1$) and so Fig 9 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 .5, with an average magnitude of about 0.1. Fig 9 also includes cyan squares that represent the reconnection rates found using the same procedure applied to MHD. This data is from Fig 10 in Wan et al. (2010a). While 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 Shay et al. (1999).

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 Servidio et al. (2009, 2010) 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), while 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 simulation for both Orszag Tang vortex configurations Parashar et al. (2009) and larger broad band turbulence simulations Wu et al. (2013).

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 Wan et al. (2010a). 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 processesShay et al. (1999).

While 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 Wan et al. (2010a). 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 Wan et al. (2013), where for large systems, at high magnetic Reynold 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 Loureiro, Schekochihin, and Cowley (2007); Pucci and Velli (2014) 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.

Acknowledgements.

Research is supported by NSF AGS-1460130 (SHINE), and NASA grant NNX15AW58G, NNX14AI63G (Heliophysics Grand challenge Theory), NNX08A083G–MMS, 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.