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.

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 established12 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 rapid13 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 agyrotropy18,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.

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 Rem3/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.

The simulations were performed using the fully kinetic 2.5D (X, Y in real space, VX, VY, VZ 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 Δx=0.02di0.73λD and Δt=0.0015Ωci1, 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 β=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Ωci1. This is slightly less than 2 τnl 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.

The second simulation is (102.4di)2 in size with a grid spacing of Δx=0.0125di and Δt=0.0025Ωci1 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 β=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Ωci1 (where Ωci1 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λD where NλD=πλD2×ppc/(Δx)2.

TABLE I.

Parameters for different turbulent simulations. The simulation length in di (lx, ly), the grid spacing in di (Δx), the time step (Ωci1 based on a uniform magnetic field of B0=1), the speed of light ccA, the ion to electron mass ratio (mime), the number of particles per cell (ppc), the number of particles per Debye circle (NλD=πλD2×ppc/(Δx)2), and the ion and electron plasma beta based on the mean magnetic field and initial mean temperatures.

Namelx = lyΔxΔtccAmimeppcNλDβi=β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 
Namelx = lyΔxΔtccAmimeppcNλDβi=β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=v×B. 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.

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 Rm10000) 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 ∼104, 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ωpe1, 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 space37,38 or sharp Fourier space filtering.22,43 In this work, we will employ a Gaussian filter for our real space filtering, defined as
    g(Bi,j)=k,lBk,l2πλd/Δxe(ik)2+(jl)22(λd/Δx)2,
    (1)

    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 ks. For both Gaussian and Fourier filtering, we applied the filters at the scale λ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.

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=2axy. 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.

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 kde2. 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.

FIG. 1.

Omnidirectional spectrum of the magnetic field (a + b) and out-of-plane electric field (c) for OTV simulations taken at 15Ωci1. Cases are shown with different particles per cell (a) and employing time averaging, Gaussian filtering, and Fourier filtering for the 3200 ppc simulation (b + c). The subplot inside (a) shows the value of the magnetic spectrum at the Debye length (EB(kλd=1)) as a function of particles per cell. The slope of the best fit line is approximately −1 which implies that the noise floor is inversely proportional to the number of particles per cell). The vertical dashed lines correspond to the wave number of 3 major lengths scales: the ion inertial length, electron internal length, and the Debye length from left to right, respectively.

FIG. 1.

Omnidirectional spectrum of the magnetic field (a + b) and out-of-plane electric field (c) for OTV simulations taken at 15Ωci1. Cases are shown with different particles per cell (a) and employing time averaging, Gaussian filtering, and Fourier filtering for the 3200 ppc simulation (b + c). The subplot inside (a) shows the value of the magnetic spectrum at the Debye length (EB(kλd=1)) as a function of particles per cell. The slope of the best fit line is approximately −1 which implies that the noise floor is inversely proportional to the number of particles per cell). The vertical dashed lines correspond to the wave number of 3 major lengths scales: the ion inertial length, electron internal length, and the Debye length from left to right, respectively.

Close modal

The lower two panels of Fig. 1 show the magnetic and out-of-plane electric field (Ez) 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 Ez. Time averaging on a few plasma frequency time scales reduces spectral power even at scales ∼di. 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.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.

FIG. 2.

Omnidirectional spectrum of the magnetic field (a) and out-of-plane electric field (b) respectively for time averaged, Gaussian filtered, Fourier filtered, and unfiltered fields for the large turbulent simulation taken at 206.25Ωci1. The three vertical dashed lines correspond to the wave numbers associated with the Debye length, the electron inertial length, and the ion inertial length.

FIG. 2.

Omnidirectional spectrum of the magnetic field (a) and out-of-plane electric field (b) respectively for time averaged, Gaussian filtered, Fourier filtered, and unfiltered fields for the large turbulent simulation taken at 206.25Ωci1. The three vertical dashed lines correspond to the wave numbers associated with the Debye length, the electron inertial length, and the ion inertial length.

Close modal

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 (Jz) 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)].

FIG. 3.

PDFs of the out-of-plane current density (Jz) for different ppc OTV simulations taken at 15Ωci1. In (a), we plot the PDFs of Jz with different ppc values ranging between 50 and 3200. Panel (b) shows the PDFs of the time averaged, Gaussian filtered, Fourier filtered, and unfiltered Jz for the 3200 simulation. Panel (c) shows the effect of Fourier filter on the PDF of the 50 ppc simulation and compares it with the unfiltered 50 and 3200 ppc simulation PDFs.

FIG. 3.

PDFs of the out-of-plane current density (Jz) for different ppc OTV simulations taken at 15Ωci1. In (a), we plot the PDFs of Jz with different ppc values ranging between 50 and 3200. Panel (b) shows the PDFs of the time averaged, Gaussian filtered, Fourier filtered, and unfiltered Jz for the 3200 simulation. Panel (c) shows the effect of Fourier filter on the PDF of the 50 ppc simulation and compares it with the unfiltered 50 and 3200 ppc simulation PDFs.

Close modal

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

SDK(r)=δBx(r)4δBx(r)22,
(2)

where the angled brackets denote a spatial average and δBx(r) is the increment of the magnetic field in the x-direction, defined as δBx(r)=Bx(x+r)Bx(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 Bx. 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 δBx 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 δBx 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.43 

FIG. 4.

Scale dependent Kurtosis of the magnetic field (Bx) for the OTV simulations with different numbers of ppc (a) and with the field time averaged, Gaussian filtered, Fourier filtered, and unfiltered for the ppc = 3200 simulation, all taken at 15Ωci1. The three vertical dashed lines correspond to the Debye length, the electron inertial length, and the ion inertial length.

FIG. 4.

Scale dependent Kurtosis of the magnetic field (Bx) for the OTV simulations with different numbers of ppc (a) and with the field time averaged, Gaussian filtered, Fourier filtered, and unfiltered for the ppc = 3200 simulation, all taken at 15Ωci1. The three vertical dashed lines correspond to the Debye length, the electron inertial length, and the ion inertial length.

Close modal

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 δBx for different averaging/filtering techniques. In all cases, the scale dependent kurtosis matches very well down to ∼0.5di, 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.5de. This implies that a Gaussian filter applied to field data cannot capture the smallest scale intermittent structures as well as a Fourier filter.

FIG. 5.

The scale dependent kurtosis of Bx for the 102.4 × 102.4 di simulation with different smoothing techniques applied at 206.25Ωci1.

FIG. 5.

The scale dependent kurtosis of Bx for the 102.4 × 102.4 di simulation with different smoothing techniques applied at 206.25Ωci1.

Close modal

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

PVI=|ΔB|2|ΔB|2,

and a scalar measurement of the electron pressure tensors agyrotropy normalized over the unit interval, defined in Swisdak19 as

Qe=Pe122+Pe232+Pe312Pe2+2PePe,

where Pe 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., Pe11 = Pe∥ and Pe22=Pe33=Pe).

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 di, well separated from the Debye length.

FIG. 6.

Horizontal cut of the partial variance increment (PVI) of the magnetic field with a lag of 1di in the OTV simulations taken at 15Ωci1. The cut is taken through the middle of the box (Y=5.12di) and shows the effect of (a) ppc and (b) different filtering techniques.

FIG. 6.

Horizontal cut of the partial variance increment (PVI) of the magnetic field with a lag of 1di in the OTV simulations taken at 15Ωci1. The cut is taken through the middle of the box (Y=5.12di) and shows the effect of (a) ppc and (b) different filtering techniques.

Close modal

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.

FIG. 7.

Horizontal cut of the electron pressure agyrotropy (Qe) in the OTV simulations taken at 15Ωci1. The cut is taken through the middle of the box (Y=5.12di) and shows the effect of (a) ppc and (b) different filtering techniques. In (c), the effects of filtering are compared between the 50 and 3200 ppc simulations.

FIG. 7.

Horizontal cut of the electron pressure agyrotropy (Qe) in the OTV simulations taken at 15Ωci1. The cut is taken through the middle of the box (Y=5.12di) and shows the effect of (a) ppc and (b) different filtering techniques. In (c), the effects of filtering are compared between the 50 and 3200 ppc simulations.

Close modal

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 ∼104. 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 (>105). 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.

FIG. 8.

Overview of the 3200 particles per cell Orszag-Tang simulation at 15Ωci1 where we plot the Fourier filtered magnetic vector potential a and with “x” denoting the location of identified critical saddle points and o's and triangles denoting minimum and maximum, respectively. From this overview, it is clear that there should only be 4 X-points, and 4 max/minimums.

FIG. 8.

Overview of the 3200 particles per cell Orszag-Tang simulation at 15Ωci1 where we plot the Fourier filtered magnetic vector potential a and with “x” denoting the location of identified critical saddle points and o's and triangles denoting minimum and maximum, respectively. From this overview, it is clear that there should only be 4 X-points, and 4 max/minimums.

Close modal
FIG. 9.

Number of critical points vs particles per cell (ppc) in Orszag-Tang simulations with different field filtering included. Red, green, and black lines correspond to the number of maximum, minimum, and saddle critical points identified in each simulation, respectively. The solid line is for the critical points identified from the unfiltered magnetic fields, the dashed lines for the time averaged fields, the dot dashed line for the Gaussian filtered fields, and the dotted line for Fourier filtered fields.

FIG. 9.

Number of critical points vs particles per cell (ppc) in Orszag-Tang simulations with different field filtering included. Red, green, and black lines correspond to the number of maximum, minimum, and saddle critical points identified in each simulation, respectively. The solid line is for the critical points identified from the unfiltered magnetic fields, the dashed lines for the time averaged fields, the dot dashed line for the Gaussian filtered fields, and the dotted line for Fourier filtered fields.

Close modal

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.

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Ωci1). For each time sample, we apply a Fourier filter to the magnetic field data for kλd>1 (λ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.

TABLE II.

Number of critical points at each different time.

Time Ωci1MinMaxX-pointsTotal
206.25 110 116 226 452 
250.00 144 165 309 618 
292.50 159 159 318 636 
Time Ωci1MinMaxX-pointsTotal
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=77di,y=22di). 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.

FIG. 10.

(a) Overview of the large turbulent simulation at 206.25Ωci1 where we plot colored contours of the Fourier filtered magnetic vector potential a and with “x” denoting the location of identified critical saddle points. (b) is an enlarged subsection from (a) as identified by the black box. It is clear from this figure that these topological structures exist at different scales throughout the entire simulations and clearly correspond to apparent coherent structures.

FIG. 10.

(a) Overview of the large turbulent simulation at 206.25Ωci1 where we plot colored contours of the Fourier filtered magnetic vector potential a and with “x” denoting the location of identified critical saddle points. (b) is an enlarged subsection from (a) as identified by the black box. It is clear from this figure that these topological structures exist at different scales throughout the entire simulations and clearly correspond to apparent coherent structures.

Close modal

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 Ez and |Ez| shown as the black triangles in the two panels of Fig. 11. We normalize Ez to the root mean square of the in plane magnetic field (note, in this simulation Bx,y,rms1) 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 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.

FIG. 11.

PDF of the reconnection rates. (a) Ez and (b) |Ez| at the X-points identified in the three different times of the large PIC simulation (black triangles) and the values found in a 20482 MHD turbulence simulation (Run 6 in Wan et al.,25 cyan squares).

FIG. 11.

PDF of the reconnection rates. (a) Ez and (b) |Ez| at the X-points identified in the three different times of the large PIC simulation (black triangles) and the values found in a 20482 MHD turbulence simulation (Run 6 in Wan et al.,25 cyan squares).

Close modal

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 simulations12,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 configurations30 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 Howes47) 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 instabilities48,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.

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.

1.
P. J.
Coleman
, Jr.
,
Astrophys. J.
153
,
371
(
1968
).
2.
J.
Saur
,
H.
Politano
,
A.
Pouquet
, and
W. H.
Matthaeus
,
Astron. Astrophys.
386
,
699
(
2002
).
3.
E.
Marsch
,
Living Rev. Sol. Phys.
3
,
1
(
2006
).
4.
A.
Retinò
,
D.
Sundkvist
,
A.
Vaivads
,
F.
Mozer
,
M.
André
, and
C. J.
Owen
,
Nat. Phys.
3
,
236
(
2007
).
5.
R.
Bruno
and
V.
Carbone
,
Living Rev. Sol. Phys.
2
(1),
4
(
2005
).
6.
M.-M.
Mac Low
,
Astrophys. J.
524
,
169
(
1999
).
7.
N.
Banerjee
and
P.
Sharma
,
Mon. Not. Roy. Astron. Soc.
443
,
687
(
2014
).
8.
M.
Yamada
,
R.
Kulsrud
, and
H.
Ji
,
Rev. Mod. Phys.
82
,
603
(
2010
).
9.
W.
Matthaeus
and
S. L.
Lamkin
,
Phys. Fluids
29
,
2513
(
1986
).
10.
11.
W.
Matthaeus
and
M.
Velli
,
Space Sci. Rev.
160
,
145
(
2011
).
12.
S.
Servidio
,
W.
Matthaeus
,
M.
Shay
,
P.
Cassak
, and
P.
Dmitruk
,
Phys. Rev. Lett.
102
,
115003
(
2009
).
13.
A.
Lazarian
and
E. T.
Vishniac
,
Astrophys. J.
517
,
700
(
1999
).
14.
S.
Servidio
,
F.
Valentini
,
D.
Perrone
,
A.
Greco
,
F.
Califano
,
W. H.
Matthaeus
, and
P.
Veltri
,
J. Plasma Phys.
81
,
325810107
(
2015
).
15.
C. K.
Birdsall
and
A. B.
Langdon
,
Plasma Physics vis Computer Simulation
(
McGraw-Hill Book Company
,
New York
,
1985
), pp.
20
22
.
16.
L.
Franci
,
A.
Verdini
,
L.
Matteini
,
S.
Landi
, and
P.
Hellinger
,
Astrophys. J. Lett.
804
,
L39
(
2015
).
17.
E.
Camporeale
,
G.
Delzanno
,
B.
Bergen
, and
J.
Moulton
,
Comput. Phys. Commun.
198
,
47
(
2016
).
18.
A.
Greco
,
P.
Chuychai
,
W. H.
Matthaeus
,
S.
Servidio
, and
P.
Dmitruk
,
Geophys. Res. Lett.
35
,
L19111
, doi: (
2008
).
19.
M.
Swisdak
,
Geophys. Res. Lett.
43
,
43
, doi: (
2016
).
20.
M.
Wan
,
W. H.
Matthaeus
,
H.
Karimabadi
,
V.
Roytershteyn
,
M.
Shay
,
P.
Wu
,
W.
Daughton
,
B.
Loring
, and
S. C.
Chapman
,
Phys. Rev. Lett.
109
,
195001
(
2012
).
21.
M.
Wan
,
W. H.
Matthaeus
,
S.
Servidio
, and
S.
Oughton
,
Phys. Plasmas
20
,
042307
(
2013
).
22.
M.
Wan
,
S.
Oughton
,
S.
Servidio
, and
W.
Matthaeus
,
Phys. Plasmas
17
,
082308
(
2010
).
23.
M.
Wan
,
S.
Oughton
,
S.
Servidio
, and
W. H.
Matthaeus
,
Phys. Plasmas
16
,
080703
(
2009
).
24.
A.
Greco
,
S.
Perri
,
S.
Servidio
,
E.
Yordanova
, and
P.
Veltri
,
Astrophys. J. Lett.
823
,
L39
(
2016
).
25.
D.
Montgomery
and
C. W.
Nielson
,
Phys. Fluids
13
,
1405
(
1970
).
26.
H.
Karimabadi
,
V.
Roytershteyn
,
M.
Wan
,
W. H.
Matthaeus
,
W.
Daughton
,
P.
Wu
,
M.
Shay
,
B.
Loring
,
J.
Borovsky
,
E.
Leonardis
,
S. C.
Chapman
, and
T. K. M.
Nakamura
,
Phys. Plasmas
20
,
012303
(
2013
).
27.
A.
Zeiler
,
D.
Biskamp
,
J. F.
Drake
,
B. N.
Rogers
,
M. A.
Shay
, and
M.
Scholer
,
J. Geophys. Res.
107
,
1230
, doi: (
2002
).
28.
S.
Orszag
and
C.
Tang
,
J. Fluid Mech.
90
,
129
(
1979
).
29.
R.
Dahlburg
and
J.
Picone
,
Phys. Fluids B: Plasma Phys.
1
,
2153
(
1989
).
30.
T. N.
Parashar
,
M. A.
Shay
,
P. A.
Cassak
, and
W. H.
Matthaeus
,
Phys. Plasmas
16
,
032310
(
2009
).
31.
B.
Vasquez
and
S.
Markovskii
,
Astrophys. J.
747
,
19
(
2012
).
32.
P.
Wu
,
S.
Perri
,
K.
Osman
,
M.
Wan
,
W. H.
Matthaeus
,
M. A.
Shay
,
M. L.
Goldstein
,
H.
Karimabadi
, and
S.
Chapman
,
Astrophys. J. Lett.
763
,
L30
(
2013
).
33.
S.
Servidio
,
W. H.
Matthaeus
,
M. A.
Shay
,
P.
Dmitruk
,
P. A.
Cassak
, and
M.
Wan
,
Phys. Plasmas
17
,
032315
(
2010
).
34.
D.
Biskamp
,
Phys. Fluids
29
,
1520
(
1986
).
35.
W.
Daughton
,
V.
Roytershteyn
,
B. J.
Albright
,
H.
Karimabadi
,
L.
Yin
, and
K. J.
Bowers
,
Phys. Rev. Lett.
103
,
065004
(
2009
).
36.
W.
Matthaeus
,
S.
Dasso
,
J.
Weygand
,
L.
Milano
,
C.
Smith
, and
M.
Kivelson
,
Phys. Rev. Lett.
95
,
231101
(
2005
).
37.
M.
Germano
,
J. Fluid Mech.
238
,
325
336
(
1992
).
38.
G. L.
Eyink
,
Phys. D: Nonlinear Phenom.
207
,
91
(
2005
).
39.
H.
Aluie
and
G. L.
Eyink
,
Phys. Fluids
21
,
115108
(
2009
).
40.
Y.
Yang
,
W. H.
Matthaeus
,
T. N.
Parashar
,
C. C.
Haggerty
,
V.
Roytershteyn
,
W.
Daughton
,
M.
Wan
,
Y.
Shi
, and
S.
Chen
,
Phys. Plasmas
24
(
7
),
072306
(
2017
).
41.
U.
Frisch
,
Turbulence
(
Cambridge University Press
,
1996
).
42.
T. N.
Parashar
,
W. H.
Matthaeus
,
M. A.
Shay
, and
M.
Wan
,
Astrophys. J.
811
,
112
(
2015
).
43.
M.
Wan
,
W.
Matthaeus
,
V.
Roytershteyn
,
H.
Karimabadi
,
T.
Parashar
,
P.
Wu
, and
M.
Shay
,
Phys. Rev. Lett.
114
,
175002
(
2015
).
44.
J.
Scudder
and
W.
Daughton
,
J. Geophys. Res.: Space Phys.
113
,
A06222
(
2008
).
45.
M. A.
Shay
,
J. F.
Drake
,
B. N.
Rogers
, and
R. E.
Denton
,
Geophys. Res. Lett.
26
,
2163
, doi: (
1999
).
46.
C.
Chen
,
A.
Mallet
,
A.
Schekochihin
,
T.
Horbury
,
R.
Wicks
, and
S.
Bale
,
Astrophys. J.
758
(
2
),
120
(
2012
).
47.
48.
N.
Loureiro
,
A.
Schekochihin
, and
S.
Cowley
,
Phys. Plasmas
14
,
100703
(
2007
).
49.
F.
Pucci
and
M.
Velli
,
Astrophys. J. Lett.
780
,
L19
(
2014
).