Three-dimensional magnetic reconnection in particle-in-cell simulations of anisotropic plasma turbulence

Abstract

We use three-dimensional (3-D) fully kinetic particle-in-cell simulations to study the occurrence of magnetic reconnection in a simulation of decaying turbulence created by anisotropic counter-propagating low-frequency Alfvén waves consistent with critical-balance theory. We observe the formation of small-scale current-density structures such as current filaments and current sheets as well as the formation of magnetic flux ropes as part of the turbulent cascade. The large magnetic structures present in the simulation domain retain the initial anisotropy while the small-scale structures produced by the turbulent cascade are less anisotropic. To quantify the occurrence of reconnection in our simulation domain, we develop a new set of indicators based on intensity thresholds to identify reconnection events in which both ions and electrons are heated and accelerated in 3-D particle-in-cell simulations. According to the application of these indicators, we identify the occurrence of reconnection events in the simulation domain and analyse one of these events in detail. The event is related to the reconnection of two flux ropes, and the associated ion and electron exhausts exhibit a complex 3-D structure. We study the profiles of plasma and magnetic-field fluctuations recorded along artificial-spacecraft trajectories passing near and through the reconnection region. Our results suggest the presence of particle heating and acceleration related to small-scale reconnection events within magnetic flux ropes produced by the anisotropic Alfvénic turbulent cascade in the solar wind. These events are related to current structures of the order of a few ion inertial lengths in size.

1. Introduction

The solar wind is a low-collisionality plasma produced in the solar corona (Marsch 2006). It expands across the solar system exhibiting spatial and temporal variations in composition, density, velocity and temperature as well as in the electric and magnetic fields. The solar wind shows a non-adiabatic temperature profile with distance from the Sun (Gazis & Lazarus 1982) which suggests the presence of local heating and particle-acceleration mechanisms (Goldstein et al. 2015). Unlike in collisional plasmas, in the solar wind, the energy dissipation cannot be attributed to the viscous interaction due to binary particle collisions nor to any process that depends directly on collisions, such as the collisional electric resistivity for instance. In the solar wind, the magnetic-field fluctuations exhibit a power-law distribution of the magnetic energy across a large range of spatial scales from 0.1 au to subproton scales (Coleman 1968; Marsch & Tu 1990) which indicates the presence of turbulence in the solar wind. The energy cascade has three regimes: the so-called injection range in which the power index of the magnetic-field fluctuations is $-1$ (Kiyani, Osman & Chapman 2015); an inertial range in which the power index varies from $-3/2$ to $-5/3$ (Iroshnikov 1963; Marsch & Tu 1990; Podesta 2009; Boldyrev et al. 2011); and a dissipation range in which the power index is less clearly defined (Goldstein, Roberts & Fitch 1994; Li, Gary & Stawicki 2001; Howes et al. 2008b) with spectral breaks at electron scales (Alexandrova et al. 2009; Sahraoui et al. 2009). The transport of energy between scales is known as the energy cascade. At subproton scales, kinetic dissipation mechanisms become important, particles are energised, and the entropy of the system irreversibly increases (Tatsuno et al. 2009; Eyink 2018; Verscharen, Klein & Maruca 2019). The nature of the fluctuations at subproton scales and the properties of the plasma determine whether the turbulent energy is mainly dissipated by ions or whether it cascades to electron scales at which it is ultimately dissipated by electrons. In the framework of wave turbulence, the energy-dissipation mechanisms are classified into two main categories: resonant heating such as Landau damping and ion-cyclotron damping (Marsch, Vocks & Tu 2003; Kasper, Lazarus & Gary 2008) and non-resonant heating such as stochastic heating (Chandran et al. 2010, 2013). In this framework, turbulent fluctuations with polarisation properties consistent with kinetic Alfvén waves (KAWs) and whistler waves are often evoked as the mechanisms that carry the turbulent cascade to electron scales. In general, observations more often find evidence for KAW-like fluctuations than for whistler-wave-like fluctuations (Smith, Vasquez & Hollweg 2011; Podesta & TenBarge 2012; Salem et al. 2012; Podesta 2013; Goldstein et al. 2015). Another mechanism proposed to carry the turbulent cascade to subproton scales is magnetic reconnection (Sundkvist et al. 2007; Franci et al. 2017; Loureiro & Boldyrev 2020).

Magnetic reconnection is a process in which particles are heated and accelerated while the magnetic-field topology changes. It takes place when magnetic structures form a region in which the frozen-in condition is locally broken allowing the exchange of particles between the magnetic structures and leading to a change in the magnetic connectivity (Hesse & Schindler 1988; Pontin 2011). Magnetic reconnection is a multiscale phenomenon that appears in both space and laboratory plasmas under conditions reaching from fully collisional to effectively collisionless. It has been predicted to occur in coronal mass ejections, solar flares, explosive events in planetary magnetospheres, accretion discs, star-formation regions, fusion plasmas and in the solar wind (see Priest & Forbes 2007; Zweibel & Yamada 2009). In the latter, reconnection events are characterised by streams of particles associated with Alfvénic disturbances and magnetic-field rotations (Gosling et al. 2005; Davis et al. 2006; Gosling, Eriksson & Schwenn 2006; Phan et al. 2006; Phan, Gosling & Davis 2009; Gosling 2012; Phan et al. 2020). These structures are interpreted as the so-called ‘exhaust regions’ of reconnection events. Although magnetic reconnection has been studied for over 60 years, there is still no consensus in terms of a complete theory to describe magnetic reconnection at all scales involved. The problem is rooted in the fact that the range of spatial ($L$) and temporal ($\tau$) scales involves fluid-like behaviour at $L \gg \rho _{i}, d_{i}$, where $\rho _i$ is the ion gyroradius and $d_{i}$ is the ion inertial length, as well as kinetic behaviour and energy dissipation at subproton scales, $L \ll \rho _{i}, d_{i}$. In addition, since plasmas are often in a turbulent state, the presence of a turbulent field alters the onset and evolution of reconnection events (Matthaeus & Lamkin 1986; Strauss 1988; Lazarian & Vishniac 1999; Kim & Diamond 2001; Servidio et al. 2011; Boldyrev & Loureiro 2017; Adhikari et al. 2020; Loureiro & Boldyrev 2020). It is unclear how turbulence and reconnection affect each other and how the energy is partitioned between particles and fields through both processes. For instance, although the role of reconnection in the small-scale turbulent cascade has been studied previously (Franci et al. 2017; Boldyrev & Loureiro 2017; Cerri & Califano 2017; Papini, Landi & Del Zanna 2019b), it is still unclear how 3-D reconnection proceeds in the turbulent solar wind. It is not well understood whether 3-D reconnection disrupts current sheets and coherent magnetic-field structures associated with intermittency at small scales in the same way as it disrupts these structures at large scales (Boldyrev et al. 2013; Mallet, Schekochihin & Chandran 2017). Moreover, it is unclear how reconnection changes the turbulent cascade as the wavevector anisotropy increases with decreasing scale and how turbulence affects the reconnection process itself (Boldyrev & Loureiro 2017). Therefore, it is necessary to study the energy partition as well as the links between turbulence and reconnection at small scales in order to fully understand the mechanisms of energy dissipation and plasma heating in the solar wind.

The use of numerical simulations has been proved to be an invaluable tool to test existing theories over a wide range of parameters. Moreover, using simulations, we self-consistently explore nonlinear problems which lie beyond analytical theory. Simulations expand our knowledge regarding magnetic reconnection processes in two dimensions (Birn et al. 2001; Shay et al. 2001; Loureiro et al. 2009; Servidio et al. 2009, 2010; Bessho et al. 2017) and in three dimensions (Hesse, Kuznetsova & Birn 2001; Pritchett & Coroniti 2001; Lapenta 2003; Lapenta et al. 2006; Kowal et al. 2009; Daughton et al. 2011; Baumann, Galsgaard & Nordlund 2013; Liu et al. 2013; Pritchett 2013; Muñoz & Büchner 2018). The use of high-performance computing facilities and increasing computational capabilities facilitate the study of plasmas from first principles using particle-in-cell (PIC) simulations (Lapenta 2012; Germaschewski et al. 2016). These simulations are able to resolve proton and electron scales and to account for phenomena that only reveal themselves using kinetic theory. For instance, electron-kinetic effects can affect ion-scale processes (Told et al. 2016) even in linear theory. These effects may be even enhanced in nonlinear processes. Currently, full PIC simulations are unable to cover the whole range of scales involved in natural plasma turbulence and reconnection since they are expensive in terms of computing memory and require small time steps to satisfy stability criteria. However, their ability to model the physics behind the energy partition at small scales makes PIC the most appropriate method to address subproton and electron-scale phenomena as well as collisionless energy dissipation.

Kinetic simulations of magnetic reconnection are often based on idealised conditions, such as the Harris current-sheet configuration (Shay et al. 2001; Scholer et al. 2003; Ricci et al. 2004; Shay et al. 2004; Daughton, Scudder & Karimabadi 2006; Daughton et al. 2011; Leonardis et al. 2013; Liu et al. 2013; Beresnyak 2016; Goldman, Newman & Lapenta 2016). In this work, we study the formation of current structures and the occurrence of 3-D magnetic reconnection as a result of turbulent dynamics in PIC simulations of collisionless anisotropic Alfvénic turbulence. We initialise our simulation with counter-propagating Alfvén waves that then self-consistently interact and generate turbulence (Howes & Nielson 2013; Howes 2015a), current-sheet structures (Howes 2016), and regions of magnetic reconnection. The overall objective of this work is to discover the properties of reconnection events that terminate the inertial-range cascade of solar-wind turbulence and define criteria that identify such features in future 3-D simulations and in spacecraft data. These results will allow future work to advance the study of linked reconnection and turbulence based on a solid and consistent framework of observable features. In § 2, we describe our initial conditions for the simulation as well as our numerical set-up. We present our results in § 3 and our conclusions in § 4.

2. Simulation

We use the explicit Plasma Simulation Code (known as PSC) (Germaschewski et al. 2016) to simulate eight anisotropic counter-propagating Alfvén waves in an ion–electron plasma. Since the theories of turbulence dissipation through reconnection in the solar wind are intrinsically connected to anisotropy through the generation of thin structures that form the precursors of current sheets, our initial waves are anisotropic. The anisotropy of the initial fluctuations is set up according to the theory of critical balance by Sridhar & Goldreich (1994) and Goldreich & Sridhar (1995), henceforth GS95. A detailed explanation of the initial conditions is presented in Appendix A. The normalisation parameters are the speed of light $c$, the vacuum permittivity $\epsilon _{0}$, the magnetic permeability $\mu _{0}$, the Boltzmann constant $k_{b}$, the elementary charge $q$, the ion mass $m_{i}$, the initial density of ions and electrons $n_{i}=n_{e}$, and the ion inertial length $d_{i}=c/\omega _{\textrm {pi}}$, where $\omega _{\textrm {pi}}=\sqrt {n_{i}q^{2}/m_{i}\epsilon _{0}}$ is the ion plasma frequency. We set $\beta _{s,\parallel }=1$ and $T_{s,\parallel }/T_{s,\perp }=1$, where $\beta _{s,\parallel }=2 n_s \mu _{0} k_{B}T_{s,\parallel }/B_{0}^{2}$, $\boldsymbol {B}_{0}$ is the background magnetic field, and the index $s$ indicates the plasma species. Here $T_{s,\parallel }$ and $T_{s,\perp }$ are the parallel and perpendicular temperatures, respectively. The magnetic field is normalised to the value of the constant background field $B_0$ and the Alfvén speed ratio is $V_{A}/c=0.1$, where $V_{A}=B_{0} / \sqrt {\mu _{0}n_{i}m_{i}}$ is the ion Alfvén speed. We use $100$ particles per cell ($100$ ions and $100$ electrons) and a mass ratio of $m_{i}/m_{e} = 100$ so that $d_e = 0.1 d_{i}$, where $m_{e}$ is the electron mass and $d_{e}$ is the electron inertial length. The simulation box size is $L_{x} \times L_{y} \times L_{z} = 24d_{i}\times 24d_{i}\times 125d_{i}$, and the spatial resolution is $\Delta x =\Delta y = \Delta z = 0.06d_{i}$. We use a time step of $\Delta t =0.06/ \omega _{\textrm {pi}}$. In our normalisation, the Debye length $\lambda _{D}=d_{i}\sqrt {\beta _{i}/2}V_{A}/c$ defines the minimum spatial distance that must be resolved in the simulation and $\lambda _D=0.07d_i$. Although our numerical parameters $V_A/c$ and $m_i/m_e$ are not identical to the corresponding parameters in the solar wind, they allow us to perform simulations within the computational limitations. With these parameters, the simulated electrons are mildly relativistic, which they are not in the real solar wind. However, the effect of mildly relativistic electrons on the propagation and damping of kinetic-scale normal modes, including KAWs, Alfvén/ion-cyclotron (known as A/IC) waves and fast-magnetosonic/whistler (known as FM/W) waves, is negligible (Verscharen et al. 2020) and not important for the evolution of the turbulent cascade, regardless of the processes that carry the cascade to subproton scales.

3. Results

In this section, we discuss the time evolution (§ 3.1) and the spectral properties (§ 3.2) of the turbulence in our simulation. We then define a new set of indicators of reconnection based on two-dimensional (2-D) and 3-D reconnection models and study a self-consistently formed reconnection region in detail (§ 3.3). We record and discuss the plasma properties that an artificial spacecraft observes in the spacecraft frame as it passes through our simulation box (§ 3.4).

3.1. Time evolution and formation of current structures

We first identify a representative time $t_R$ for our subsequent analysis of the turbulence properties. The root mean square (r.m.s.) of the current density $J^{\textrm {rms}}$ is an indicator commonly used to identify the time at which the system reaches a quasi-stationary state. At this time, the generation of current sheets by waves is balanced by their decay so that the growth of $J^{\textrm {rms}}$ saturates, which marks the time of maximum turbulent activity in the simulation (Franci et al. 2017). The r.m.s. of a quantity $\psi$ is defined as

(3.1)\begin{equation} \psi^{\textrm{rms}} = \sqrt{\langle \psi^{2} \rangle - \langle \psi \rangle^{2}}, \end{equation}

where $\langle \cdots \rangle$ represents the spatial average over the whole simulation domain. Figure 1 shows the time evolution of the r.m.s. of the current density $\boldsymbol {J}$ (blue), the magnetic field $\boldsymbol {B}$ (black) and the ion velocity $\boldsymbol {v}_{i}$ (red) in our simulation. Since we start our simulation under the assumption that the linear time $\tau _{l}$ is approximately equal to the nonlinear time $\tau _{\textrm {nl}}$, we estimate $\tau _{\textrm {nl}} \sim \tau _{l} \sim 1/k_{\parallel }V_{A} \sim L_{z}/2{\rm \pi} V_{A} \approx 200/\omega _{\textrm {pi}}$. This estimate for the nonlinear time $\tau _{\textrm {nl}}$ is, therefore, related to the scale of the initial fluctuations and represents an upper limit. We observe a peak in $J^{\textrm {rms}}$ at $t=12/\omega _{\textrm {pi}}$ which is due to the self-consistent formation of current structures as a response to the initial magnetic-field fluctuations. The variation in $B^{\textrm {rms}}$ and $J^{\textrm {rms}}$ during the initial phase, between $t=12/\omega _{\textrm {pi}}$ and $t=96/\omega _{\textrm {pi}}$, suggests that the system is still in a phase of self-adjustment. The formation of the plateau in $J^{\textrm {rms}}$ at $t \approx \tau _{\textrm {nl}}/2 \approx 100/\omega _{\textrm {pi}}$ indicates that the system has reached a quasi-stationary state. Therefore, we expect the formation of current structures such as current sheets and current filaments by this time. The vertical dashed line marks the time $t = 120/\omega _{\textrm {pi}}$ at which $J^{\textrm {rms}}$ begins to decrease monotonically until the simulation ends. In this sense, the time $t=120/\omega _{\textrm {pi}}$ represents the beginning of the decaying phase in our system. As the system evolves in time, current and magnetic structures dissipate, and we expect an exchange of the energy stored in the magnetic field with the kinetic energy of the particles. Based on these considerations, we use the time $t_{R}=120/\omega _{\textrm {pi}}$ to study the spectral properties of the turbulence in our system.

Figure 1. Time evolution of the r.m.s. of the current density $\boldsymbol {J}$ (blue), magnetic field $\boldsymbol {B}$ (black) and ion velocity $\boldsymbol {v}_{i}$ (red). The vertical dashed line marks the time $t_{R} = 120/\omega _{\textrm {pi}}$ at which $J^{\textrm {rms}}$ begins to decrease.

Figure 2 shows a 3-D rendering of the magnitude of the transverse magnetic field $|\boldsymbol {B}_{xy}| = \sqrt {B_{x}^{2}+B_{y}^{2}}$ at two different time steps: $t=0$ (panel (a)) and $t=t_{R}$ (panel (b)). Figure 2(a) shows the anisotropic interference pattern of the linear superposition of Alfvén waves at $t=0$. Initially, there are no coherent eddies present because no nonlinear interaction has taken place yet. However, the initial magnetic-field fluctuations are already anisotropic. Figure 2(b) shows that at time $t=t_{R}$, there is a clear presence of magnetic eddies with varying cross-section diameters $L_{D}$ and elongations $L_{\parallel }$, where $L_{\parallel }$ represents the length of these eddies along the local magnetic field. Even though we start with a superposition of only eight waves, nonlinear interactions generate magnetic eddies of different shapes and anisotropies. At this time, the magnetic-field structures consist of a combination of linear fluctuations and magnetic eddies. To estimate the shape of the magnetic structures at $t=0$ and $t=t_{R}$, we calculate $\Delta B=\sqrt {B_{x}^{2}+B_{y}^{2}+(B_{z}-B_{0})^{2}}$ and use an intensity threshold defined as $\Delta B > \langle \Delta B \rangle +2 \Delta B^{\textrm {rms}}$. We define a magnetic structure as the combination of those cells in our simulation that are connected as next neighbours and fulfil this threshold condition. The exact value of the threshold is chosen to improve the performance of the algorithm in the identification of these structures. After the identification of the structures, we calculate their principal axes. We define $L_D=\sqrt {L_{\perp 1}^2+L_{\perp 2}^2}$, where $L_{\perp 1}$ and $L_{\perp 2}$ are the two orthogonal diameters in the plane perpendicular to the local magnetic field and $L_{\parallel }$ is the elongation of the structure along the local magnetic field.

Figure 2. Three-dimensional rendering of the transverse magnetic field magnitude $|\boldsymbol {B}_{xy}| = \sqrt {B_{x}^{2}+B_{y}^{2}}$ at $t=0$ (a) and $t=t_{R}$ (b). The colour bar ranges from the minimum magnitude (black) to the maximum magnitude (yellow) throughout the simulation domain at $t=t_{R}$. We use the same colour bar in both panels for a direct comparison. The initial background magnetic field is directed along the $z$-direction. At the initial time, the fluctuations are anisotropic and elongated along the $z$-direction. At $t=t_{R}$, small-scale magnetic eddies have formed and interact nonlinearly with each other. The eddies present varying cross-section diameters $L_{D}$ and lengths $L_{\parallel }$.

Figure 3(a) shows the probability distribution function (PDF) of $L_{\parallel }$, $L_{D}$, and the aspect ratio $L_{\parallel }/L_{D}$ at $t=0$ and figure 3(b) at $t=t_{R}$. The mean value and standard deviation of the distributions of $L_{\parallel }$, $L_D$, and $L_{\parallel }/L_D$ at $t=0$ are ${L}_{\parallel } = (16.33 \pm 8.32)d_{i}$, ${L}_{D} = (1.55 \pm 0.95)d_{i}$, and ${(L_{\parallel }/L_{D})} = (11.01 \pm 7.06)d_{i}$. At $t=t_{R}$, we find ${L}_{\parallel } = (2.16 \pm 5.08)d_{i}$, ${L}_{D} = (0.62 \pm 0.72)d_{i}$, and ${L_{\parallel }/L_{D}} = (2.55 \pm 1.94)d_{i}$. According to this analysis, the nonlinear interaction has formed magnetic structures with smaller elongations and cross-section diameters continuously distributed between $L_{D}=1d_{i}$ and $8 d_{i}$. The distribution of aspect ratios is less uniform at $t=t_{R}$ than at $t=0$. The number of magnetic structures with nearly isotropic aspect ratios is greater at $t=t_{R}$. To study the distribution of the large-scale structures at $t=t_{R}$, we further apply a filter to remove all regions with an equivalent volume $V \leqslant 1d_{i}^{3}$, where $V$ is defined as the space filled by the sum of all contiguous cells associated with a given magnetic structure. For all structures with $V>d_{i}^{3}$, we find ${L}_{\parallel } = (14.97 \pm 9.01)d_{i}$, ${L}_{D} = (3.14 \pm 2.25)d_{i}$, and ${L_{\parallel }/L_{D}} = (5.46 \pm 2.48)d_{i}$. The distribution of the large-scale magnetic structures maintains an anisotropy consistent with our initial conditions. Figure 3(c) shows the scaling between $L_{\parallel }$ and $L_{D}$ for the magnetic structures at $t=0$. The linear fit to these structures, dashed line, reveals the scaling $L_{\parallel } \sim L_{D}^{0.66}$, which is consistent with our initial anisotropy, i.e. $L_{\parallel } \sim L_{D}^{2/3}$. Figure 3(d) shows the scaling between $L_{\parallel }$ and $L_{D}$ for the magnetic structures at $t=t_{R}$. The orange dots represent structures satisfying $V > d_{i}^{3}$ while the blue dots show structures satisfying $V \leqslant d_{i}^{3}$. The linear fit to the former population, top black dashed line, reveals the scaling $L_{\parallel } \sim L_{D}^{0.7}$. In contrast, the linear fit to the latter population, bottom red dashed line, shows an isotropic scaling, $L_{\parallel } \sim L_{D}$. Around $L_{D} \sim d_{i}$, we find a transition and mixing between structures with both scalings. This suggests that the large-scale structures tend to maintain the initial anisotropy while the small-scale structures become more isotropic. This isotropic scaling at subproton scales has also been observed in hybrid simulations (Franci et al. 2018; Arzamasskiy et al. 2019; Landi et al. 2019).

Figure 3. (a,b) Probability distribution functions (PDF) of elongations $L_{\parallel }$ (top), cross-section diameters $L_D$ (middle), and aspect ratios $L_{\parallel }/L_D$ (bottom) of the magnetic structures at $t=0$ (a) and $t=t_R$ (b). (c) Scaling between $L_{\parallel }$ and $L_D$ at $t=0$. The black dashed line shows the linear fit. (d) Scaling between $L_{\parallel }$ and $L_D$ of the large-scale population (orange) and small-scale population (blue) at $t=t_{R}$. The top black dashed line shows the linear fit to the former while the bottom red dashed line shows linear fit to the latter.

Figure 4 shows 3-D renderings of $B_{z}$ and $|\boldsymbol {J}|$ at $t=t_{R}$. Figure 4(a) shows $B_{z}$, from the same vantage point as figure 2(b). Although the initial $\boldsymbol {B}_{0}$ is uniform and points in the $+z$-direction, nonlinear interactions generate regions in which $B_{z}$ is negative. These regions are mostly localised in the centres of the small eddies in figure 2(b). Figure 4(b) shows that the locations of the most intense current filaments coincide with the centres of the magnetic eddies. Current filaments are intense quasi-cylindrical current structures. Similar to the case of the magnetic structures, we apply the threshold $|J| \geqslant \langle |J|\rangle + 4(|J|)^{\textrm {rms}}$ to determine the shape of the current filaments. The mean cross-section diameter of these current filaments is $\hat {L}_{D} = (1.94 \pm 0.84)d_{i}$. Their mean elongation is $\hat {L}_{\parallel } = (12.32 \pm 6.70)d_{i}$, and their mean aspect ratio is $\hat {L}_{\parallel }/\hat {L}_{D} = (6.84 \pm 3.48)$. The filaments are mostly elongated along the ${z}$-direction. Some filaments have undergone bending and twisting due to the nonlinear interactions. The elongations of the current filaments are distributed similarly to the elongations of the magnetic eddies (not shown here) and vary in the range of scales from ${\sim } 4 d_{i}$ to ${\sim } 30 d_{i}$. Panel (b) shows in addition the formation of thin current-sheet-like structures at the edges of the eddies where the perpendicular component of the magnetic field is nearly zero (see figure 2b). We define current sheets as current structures in which $L_{cs} \gg \delta _{cs}$ and $\varDelta _{cs} \gg \delta _{cs}$, where $L_{cs}$ is the current-sheet length along the local magnetic field, $\varDelta _{cs}$ is the current-sheet width tangential to the magnetic eddies, and $\delta _{cs}$ is the current-sheet thickness normal to the edge of the eddies. The formation of these current sheets is due to the turbulent motions that squeeze the eddies together. In the supplementary material available at https://doi.org/10.1017/S0022377821000404 we provide a movie that shows the time evolution of the volume rendering of $J_{z}$ in the $zx$-plane. We observe the tearing and breaking up of current sheets as well as the onset of instabilities arising from the nonlinear interactions and of jets oblique to the major axes of the current sheets as a result of the turbulent evolution. However, a detailed study of these phenomena is beyond the scope of this work.

Figure 4. Visualisations of the simulation domain at $t=t_{R}$. (a) Three-dimensional rendering of the magnetic-field component $B_{z}$. Blue represents the negative, red the positive and white the zero values of $B_{z}$. The eddies’ centres present different values of $B_z$ with either positive or negative polarity. (b) Three-dimensional rendering of the magnitude of the current density $| \boldsymbol {J} |$ from the same vantage point as panel (a). The colour represents in blue (red) the smallest (largest) values of $|\boldsymbol {J}|$. Filaments of intense current density are aligned with the eddies’ centres. Current filaments and extended current-sheet-like structures are mainly elongated along the $z$-direction.

3.2. Evidence of turbulence

A broad power-law spectrum of the fluctuations indicates the presence of turbulence as the energy cascades from large to small scales. To analyse the spectral properties of the system, following Franci et al. (2018), we calculate the energy associated with the 3-D Fourier modes $\psi _{3D}(\boldsymbol {k})$ of a quantity $\psi$ as

(3.2)\begin{equation} \psi_{3D}(\boldsymbol{k}) = \tilde{\psi}(\boldsymbol{k})\tilde{\psi}^{*}(\boldsymbol{k}), \end{equation}

where $\boldsymbol {k}$ is the wavevector, $\tilde {\psi }(\boldsymbol {k})$ is the 3-D spatial Fourier transform of $\psi$ and $\tilde {\psi }^{*}(\boldsymbol {k})$ represents its complex conjugate. If $\psi$ is a vector quantity, the 3-D Fourier transform is taken over each component, and the product is defined as

(3.3)\begin{equation} \tilde{\psi}(\boldsymbol{k})\tilde{\psi}^{*}(\boldsymbol{k}) = \sum_{i}\tilde{\psi}_{i}(\boldsymbol{k})\tilde{\psi}_{i}^{*}(\boldsymbol{k}), \end{equation}

where the index $i$ represents the components $x,y$, and $z$. Since our system does not include any anisotropy within the plane perpendicular to the background magnetic field on average, we assume that the energy distribution in the turbulent fluctuations remains axially symmetric on average. Thus, the wavevector can be expressed, without loss of generality, as its perpendicular and parallel components $(k_{\perp }, k_{\parallel })$. We note that the local (rather than the global) average magnetic field defines the cylindrical symmetry axis for the turbulent fluctuations (Cho & Vishniac 2000). However, we use the global background magnetic field as a proxy. This simplification is motivated by the strong alignment of the eddies with the background magnetic field at this time in our simulation (see figures 2 and 4). Moreover, the definition of the local magnetic field is a matter of ongoing research and debate (Podesta 2009; Chen et al. 2011; TenBarge et al. 2012; Oughton et al. 2015; Gerick, Saur & von Papen 2017), and the development of an anisotropic energy cascade is sufficient for the determination of reconnection events in the present study.¹ Thus, we calculate the perpendicular and parallel components of the wavevector as $k_{\perp }=\sqrt {k_{x}^{2}+k_{y}^{2}}$ and $k_{\parallel }=k_{z}$, respectively, and assume that the fluctuations are statistically independent of the azimuthal angle. We integrate $\psi _{3D}$ over concentric rings in $k_\perp$-space. The energy associated with the $j$th-ring is

(3.4)\begin{equation} \psi^{j}_{2D}(k_{{\perp}},k_{{\parallel}}) = \int_{k_{{\perp}}^{j}}^{k_{{\perp}}^{j}+ \textrm{d}k_{{\perp}}} \psi_{3D}(k_{{\perp}}, k_{{\parallel}}) 2{\rm \pi} k_{{\perp}}^{j\prime}\,\textrm{d}k_{{\perp}}^{\prime}, \end{equation}

where the thickness $\textrm {d}k_{\perp }$ of these rings is taken as the magnitude of the smallest perpendicular wavevector in our system $\textrm {d}k_{\perp } = 2{\rm \pi} / \sqrt {2}L_{x}$. To visualise the energy cascade in $k$-space as well as the level of anisotropy in the system, we compute the reduced 2-D power spectral density $P^{\psi }_{2D}(k_{\perp },k_{\parallel })$ as

(3.5)\begin{equation} P^{\psi}_{2D}(k_{{\perp}},k_{{\parallel}}) = \sum_{j} \frac{1}{k_{{\perp}}}\psi^{j}_{2D}(k_{{\perp}},k_{{\parallel}}). \end{equation}

Figure 5 shows the logarithm of the 2-D reduced power spectral density of the magnetic-field fluctuations $P^{\boldsymbol {B}}_{2D}$ normalised to $\max {P^{\boldsymbol {B}}_{2D}}$ in the $k_{\parallel }$–$k_{\perp }$ plane at $t=0$ (a), $t=12/\omega _{\textrm {pi}}$ (b), $t=t_{R}$ (c) and $t=240/ \omega _{\textrm {pi}}$ (d). The horizontal dashed line marks $k_{\perp }d_{e}=1$ which corresponds to $k_{\perp }d_{i} = 10$ owing to our mass ratio of $m_i/m_e=100$. The vertical dashed line marks $k_{\parallel }d_{i}=1$. At $t=0$, the energy is entirely stored in the initial modes. At $t=12/\omega _{\textrm {pi}}$, the isocontours show that the energy has already cascaded to $k_{\perp } d_e >1$ whereas the parallel cascade has not yet reached the kinetic range. At $t=t_{R}$, the perpendicular cascade has not proceeded any farther but the parallel energy transport reached $k_{\parallel } d_i >1$. At $t=240/\omega _{\textrm {pi}}$, the energy distribution has not considerably changed compared with the distribution at $t=120/\omega _{\textrm {pi}}$. For comparison with analytical predictions, we overplot the expected critical-balance scaling of $k_{\perp }\sim k_{\parallel }^{3/2}$ as a dashed line at small $k_{\perp }$. We note, however, that $P_{2D}^{\boldsymbol {B}}$ exhibits a broad distribution in $\boldsymbol {k}$-space around this prediction. In order to explore the anisotropy of the cascade in more detail, we compute the perpendicular one-dimensional (1-D) reduced power spectral density,

(3.6)\begin{equation} P^{\psi}_{1D_{{\perp}}}(k_{{\perp}}) = \int_{0}^{\infty} P^{\psi}_{2D}(k_{{\perp}},k_{{\parallel}})\,\textrm{d}k_{{\parallel}}, \end{equation}

and the parallel 1-D reduced power spectral density,

(3.7)\begin{equation} P^{\psi}_{1D_{{\parallel}}}(k_{{\parallel}}) = \int_{0}^{\infty} P^{\psi}_{2D}(k_{{\perp}},k_{{\parallel}})\,\textrm{d}k_{{\perp}} ,\end{equation}

of multiple fluctuating quantities $\psi$. Figure 6(a) shows the perpendicular 1-D reduced power spectral density of the magnetic-field fluctuations $P^{{B}}_{1D_{\perp }}$ (black line), of the ion velocity fluctuations $P^{{v}_{i}}_{1D_{\perp }}$ (red line), and of the ion density fluctuations $P^{{n}_{i}}_{1D_{\perp }}$ (blue line) at $t=t_{R}$. The vertical dashed lines mark $k_{\perp }d_{i}=1$, $k_{\perp }d_{e}=1$, and $k_{\perp }\lambda _{D}=1$. The enhancement in $P^{{v}_i}_{1D_{\perp }}$ at $k_{\perp }d_{i} = 17$ is an artefact created by Debye-length effects and the finite spatial resolution of the system. The scale of the initial waves in the perpendicular direction coincides with the transition point of the energy cascade from inertial to kinetic scales, i.e. $k_{\perp }d_{i}=1$. Therefore, our simulations do not describe the cascade at $k_{\perp }d_{i}\leqslant 1$. During the first nonlinear time, the system develops a broadband spectrum of perpendicular density fluctuations in the kinetic range. Furthermore, $P^{{B}}_{1D_{\perp }}$ and $P^{{v}_{i}}_{1D_{\perp }}$ exhibit similar spectral indices in part of the kinetic range between $k_{\perp }d_{i} \sim 3$ and ${\sim } 6$. Within the same interval, $P^{{n}_{i}}_{1D_{\perp }}$ follows a steeper spectrum. These features suggest the presence of both Alfvénic and compressive fluctuations, consistent with the presence of KAWs. Also $P^{{B}}_{1D_{\perp }}$ in the interval $k_{\perp }d_{i} \sim 1.8$ to ${\sim } 7$ follows a power-law scaling with a spectral slope of $-3$. In the range between $k_{\perp }d_{i} \sim 7$ and ${\sim } 20$, the slope is slightly steeper with a power index of approximately $-4$.² Although we calculate the energy spectrum of the magnetic-field fluctuations using the global background magnetic field, these values are within the range of slope variability measured in the solar wind (Chen et al. 2010a; Bruno, Trenchi & Telloni 2014) as well as in hybrid simulations (Franci et al. 2018; González et al. 2019).

Figure 5. Isocontours of $\log _{10}P^{\boldsymbol {B}}_{2D}$ of the fluctuating magnetic field as a function of $k_{\parallel }$ and $k_{\perp }$ at different time steps. The dashed lines provide a reference for the scaling of $k_{\perp }$ and $k_{\parallel }$. The horizontal (vertical) dashed line marks $k_{\perp }d_{e}=1$ ($k_{\parallel }d_{i}=1$). At $t=0$, the spectrum shows the modes of our initialisation and their Fourier harmonics. At $t=12/\omega _{\textrm {pi}}$, the cascade in the perpendicular direction (vertical axis) has proceeded beyond electron scales ($k_{\perp }d_{i} \geqslant 10$). At $t=t_{R}$, although the perpendicular cascade has not proceeded significantly farther, the cascade in the parallel direction (horizontal axis) has reached the kinetic range ($k_\parallel d_i\approx 1$) up to ion scales but not to electron scales. At $t=240/\omega _{\textrm {pi}}$, the distribution has not considerably changed compared with $t=t_{R}$.

Figure 6. (a) Perpendicular and (b) parallel reduced 1-D power spectral densities $P^{B}_{1D_{\parallel , \perp }}$ (black), $P^{v_{i}}_{1D_{\parallel , \perp }}$ (red), and $P^{n_{i}}_{1D_{\parallel , \perp }}$ (blue) at $t=t_{R}$. The vertical dashed lines indicate $k_{\parallel , \perp }d_{i} = 1$, $k_{\parallel , \perp }d_{e} = 1$, and $k_{\parallel , \perp }\lambda _{D} = 1$.

Figure 6(b) shows the parallel 1-D reduced power spectral density of the magnetic-field fluctuations $P^{{B}}_{1D_{\parallel }}$ (black line), ion velocity fluctuations $P^{{v}_{i}}_{1D_{\parallel }}$ (red line), and ion density fluctuations $P^{{n}_{i}}_{1D_{\parallel }}$ (blue line) at $t=t_{R}$. The vertical dashed lines mark $k_{\parallel }d_{i}=1$, $k_{\parallel }d_{e}=1$, and $k_{\parallel }\lambda _{D}=1$. At $k_{\parallel }d_{i} \leqslant 1$, $P^{{B}}_{1D_{\parallel }}$ and $P^{{v}_{i}}_{1D_{\parallel }}$ follow a similar trend as expected for Alfvénic turbulence. The spectral slope for $P^{B}_{1D_{\parallel }}$ is close to $-2$ between $k_{\parallel }d_{i} \sim 0.1$ and ${\sim }0.3$ which is in agreement with the magnetic-field power spectrum $k_{\parallel }^{-2}$ observed in the solar wind (Bavassano & Bruno 1989; Grappin, Velli & Mangeney 1991; Wicks et al. 2010, 2011; Chen et al. 2011). At smaller parallel scales, the spectrum steepens to $-2.5$ between $k_{\parallel }d_{i} \sim 0.4$ and ${\sim } 2$ and farther towards $-4$ between $k_{\parallel }d_{i} \sim 2$ and ${\sim } 4$. Both the perpendicular and parallel spectral indices have values of -4. The equality of these exponents has been observed in 3-D hybrid PIC simulations and has been suggested to be a consequence of the anisotropy being frozen at subproton scales (Franci et al. 2018; Arzamasskiy et al. 2019; Cerri, Groselj & Franci 2019; Landi et al. 2019). Although we initialise the system with non-compressive waves, the simulation swiftly develops a cascade of density fluctuations which suggests that compressive modes form self-consistently in the energy cascade. The development of compressive fluctuations has been suggested to depend on the plasma parameters rather than the initial conditions (Cerri et al. 2017a). The level of compressive fluctuations in our simulation is greater than observed in the solar wind (Chen 2016), but the reasons for the creation of such strong compressive fluctuations is unknown. At $k_{\parallel }d_{i} \approx 1.4$, the slope of $P^{{v}_{i}}_{1D_{\parallel }}$ separates from the slope of $P^{B}_{1D_{\parallel }}$ and approaches the slope of $P^{{n}_{i}}_{1D_{\parallel }}$. The flattening of $P^{{n}_{i}}_{1D_{\parallel }}$ at $k_{\parallel }d_{i} \approx 4$ is due to finite particle noise.

3.3. Reconnection sites

In this section, we confirm that magnetic reconnection occurs in our simulation domain. Methods to find reconnection sites in 2-D simulations are based on the identification of magnetic islands and their closest $x$-point within a current sheet (Wan et al. 2014a; Papini et al. 2019a). However, the interaction of magnetic structures such as flux tubes, which are the 3-D equivalent of 2-D magnetic islands, is more complex than in the 2-D case, and magnetic reconnection does not happen at a single point but in an extended region (Daughton et al. 2011; Liu et al. 2013; Daughton et al. 2014). In 2-D and 3-D theories of reconnection, strong current sheets are often associated with reconnection events as the key locations of energy dissipation. However, there are events in which the $x$-points are not placed exactly within the current sheet (Priest & Démoulin 1995). The presence of a strong guide magnetic field and asymmetries of the reconnection event can shift the position of the $x$-point and even preclude the reconnection event (Eastwood et al. 2010, 2013). Moreover, proton temperature anisotropies in reconnection events can trigger kinetic instabilities, which then have a stabilising effect on the current sheet (Matteini et al. 2013).

In our turbulent simulation set-up, we expect that once the reconnection events have occurred, most of them exhibit local asymmetries due to the turbulent nature of the domain. Moreover, the background magnetic field acts as a guide field in our reconnecting flux ropes. Therefore, in order to capture all reconnection events in such a complex and asymmetric field geometry, we require a new method to determine reconnection sites in our 3-D simulations. Strong gradients in at least one component of the magnetic field as well as magnetic null points are common features of both 2-D and 3-D reconnection events. Strong gradients directly relate to the presence of current sheets according to Ampère's law. The presence of magnetic null points is not a requirement for reconnection though. In 2-D reconnection, for instance, the presence of a guide field removes this requirement (Hesse, Kuznetsova & Birn 2004). 3-D reconnection, on the other hand, can take place in collapsing structures that form current sheets related to quasi-separator lines, which do not require magnetic null points (Pritchett & Coroniti 2004; Pontin 2011). Exhaust regions in which particles are accelerated to velocities near the Alfvén speed are another common feature. Magnetic reconnection not only accelerates particles but also increases their thermal energy. Hence, an enhancement in the population of heated particles is a further indicator of reconnection as long as it occurs near a region in which accelerated particles and magnetic field gradients are present.

During magnetic reconnection, the electric field is responsible for the energy exchange between particles and fields in the current sheet. The associated energy exchange is quantified by $\boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E}$ (Somov & Titov 1985; Ni et al. 2016). We expect to find coherent regions in the simulation domain in which $\boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E}$ is non-zero. According to 3-D steady-state theories of magnetic reconnection (Hesse & Schindler 1988; Priest, Hornig & Pontin 2003; Pontin 2011), when a magnetic field line enters a diffusion region, the integral of the parallel electric field ($E_{\parallel }=\boldsymbol {E}\boldsymbol {\cdot } \boldsymbol {B}/| \boldsymbol {B}|$) along the magnetic field line within the diffusion region must be different from zero. Since a non-zero $E_{\parallel }$ can indicate the presence of non-vanishing diffusive terms in Ohm's law, we use the presence of non-zero $E_{\parallel }$ as a possible indicator for a diffusion region located within a finite volume. Although $E_{\parallel }$ is not a good indicator in the absence of a guide magnetic field, we expect to find coherent regions in the simulation domain with non-zero $E_{\parallel }$.

In summary, we identify the following indicators that we consider essential for the presence of reconnection in a region of our simulation domain. We adopt a clustering detection method (Uritsky et al. 2010) based on the mean value of each quantity $\psi$, its r.m.s. value $\psi ^{\textrm {rms}}$, and a threshold value $N_{\textrm {th}}$. Thus, we search the simulation domain for regions in which $\psi \geqslant \langle \psi \rangle +N_{\textrm {th}}(\psi )^{\textrm {rms}}$. Our indicators for magnetic reconnection are the following.

1. C1 Current-density structures, $|\boldsymbol {J}| \geqslant \langle |\boldsymbol {J}| \rangle +N_{\textrm {th}}(|\boldsymbol {J}|)^{\textrm {rms}}$.³

2. C2 Fast ions and electrons, $|\boldsymbol {v}_{i,e}| \geqslant \langle |\boldsymbol {v}_{i,e}| \rangle +N_{\textrm {th}}(|\boldsymbol {v}_{i,e}|)^{\textrm {rms}}$.

3. C3 Heated particles, $T_{i,e} \geqslant \langle T_{i,e} \rangle +N_{\textrm {th}}(T_{i,e})^{\textrm {rms}}$.

4. C4 Energy transfer between fields and particles, $| \boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E} - \langle |\boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E}| \rangle | \geqslant N_{\textrm {th}}(|\boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E}|)^{\textrm {rms}}$.

5. C5 Non-zero parallel electric fields, $| E_{\parallel } - \langle |E_{\parallel }| \rangle | \geqslant N_{\textrm {th}}(|E_{\parallel }|)^{\textrm {rms}}$.

To find the number of events satisfying these conditions, we use the first-neighbour volumetric method described in § 3.1. We apply the algorithm to identify clusters of contiguous cells fulfilling each condition separately as well as combinations of them. Afterwards we apply a filter to remove all regions with an equivalent volume $V \leqslant 1d_{i}^{3}$, where $V$ is defined as the sum of the volumes of all contiguous cells associated with the cluster. This is motivated by the fact that we are mostly interested in events in which both ions and electrons experience reconnection. Therefore, we expect to find coherent regions with a size of at least $d_i$. We analyse two values for the threshold: $N_{\textrm {th}}=3$ and $N_{\textrm {th}}=4$. We present our results in table 1, where C2$_{i}$ and C2$_{e}$ refer to the separate application of criterion C2 to ions and to electrons, respectively. The same definitions apply to C3. Next, C4$_{+}$ and C4$_{-}$ refer to the application of condition C4 separated by cases in which $\boldsymbol {J}\boldsymbol {\cdot } \boldsymbol {E} >0$ ($+$) and $\boldsymbol {J}\boldsymbol {\cdot } \boldsymbol {E} <0$ ($-$). The same definitions apply to C5. As expected, a larger number of locations fulfil these conditions if the threshold is lower. Moreover, all events detected with $N_{\textrm {th}}=4$ are also detected when using $N_{\textrm {th}}=3$. There are no events that fulfil our condition C5. The reason for this result is that, although local regions fulfil C5, the volume of contiguous cells fulfilling C5 is never greater than $1d_i^3$. We attribute this effect to particle noise, which has a strong effect on parallel electric fields in PIC simulations. If we reduce the threshold to $N_{\textrm {th}}=2$, the algorithm is also unable to define clusters of cells, because our method is based on intensity thresholds which perform well for quantities with heavy tail distributions. The distribution of $E_{\parallel }$ in our simulation is spread with $\langle |E_{\parallel }| \rangle = 2.0 \times 10^{-3} B_{0}c$ and standard deviation $(|E_{\parallel }|)_{\textrm {std}} = 1.5 \times 10^{-3} B_{0}c$. The same argument applies to $\boldsymbol {J} \cdot \boldsymbol {E}$. Despite detecting at least 17 regions fulfilling C4 with $N_{\textrm {th}}=4$, there are no regions that satisfy all conditions C1 to C4 within a volume greater than $1d_{i}^3$. However, if we reduce the equivalent volume threshold to $0.3d_{i}^3$, we find six regions that fulfil conditions C1 to C4. We mark the corresponding numbers with an asterisk in table 1.

Table 1. Number of events in our simulation domain at time $t=t_R$ fulfilling each condition.

In figure 7, we visualise our indicators for magnetic reconnection. We use a 2-D projection on the $zx$-plane of a part of our simulation domain, $50 d_{i} < L_{z} <100 d_{i}$. Figure 7(a) shows the isosurfaces of $|\boldsymbol {J}| = \langle |\boldsymbol {J}| \rangle +3(|\boldsymbol {J}|)^{\textrm {rms}}$ (indicator C1) colour-coded in light blue. The selected structures mainly correspond to current filaments. Figure 7(b) shows regions in which $|\boldsymbol {v}_{i}| = \langle \boldsymbol {v}_{i} \rangle +3(\boldsymbol {v}_{i})^{\textrm {rms}}$ (green) and $|\boldsymbol {v}_{e}| = \langle \boldsymbol {v}_{e} \rangle +3(\boldsymbol {v}_{e})^{\textrm {rms}}$ (purple), our indicator C2. The locations of fast electrons according to C2 coincide with the locations of large currents according to C1, since the electrons are the main carriers of the electric current. This electron behaviour is consistent with observations in space plasma and reproduced in simulations (Phan et al. 2018). We identify five structures in which accelerated ions coincide with our condition C1. Figure 7(c) shows isosurfaces of $T_{i} = \langle T_{i} \rangle +3(T_{i})^{\textrm {rms}}$ (gold) and $T_{e} = \langle T_{e} \rangle +3(T_{e})^{\textrm {rms}}$ (pink), according to our indicator C3. Although the electric current is mostly carried by electrons, we find current structures that are not associated with high-temperature electrons and vice versa. The structures associated with heated electrons have mostly filamentary shapes. Figure 7(d) shows the application of our indicator C4. The regions in which $\boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E} = \langle \boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E} \rangle \pm 3(\boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E})^{\textrm {rms}}$ is positive (negative) are colour-coded in red (blue). There are large and diffuse clusters of positive and negative $\boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E}$ between $z=55d_{i}$ and $z=85d_{i}$. We also locate filamentary structures of positive $\boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E}$ which partially coincide with the regions fulfilling C3. Figure 7(e) shows our indicator C5. The regions in which $E_{\parallel } = \langle |E_{\parallel } |\rangle \pm 2(|E_{\parallel }|)^{\textrm {rms}}$ is positive (negative) are colour-coded in orange (blue). The effect of particle noise on the electric field leads to difficulties in the determination of the associated clusters. Figure 7(f) shows the combination of our indicators C1 to C4. We define two regions, $R_1$ and $R_2$, as the regions in which our indicators C1 to C4 are fulfilled. According to our assumptions, magnetic reconnection is taking place in the vicinity of these regions.

Figure 7. Reconnection indicators projected onto a 2-D cut in the $zx$-plane at $y=21d_{i}$. (a) Indicator C1: isosurfaces of $|\boldsymbol {J}| = \langle |\boldsymbol {J}| \rangle +3(|\boldsymbol {J}|)^{\textrm {rms}}$ (light blue). (b) Indicator C2: isosurfaces of $|\boldsymbol {v}_{i,e}| = \langle \boldsymbol {v}_{i,e} \rangle +3(\boldsymbol {v}_{i,e})^{\textrm {rms}}$ for ions (green) and for electrons (purple). (c) Indicator C3: isosurfaces of $T_{i,e} = \langle T_{i,e} \rangle +3(T_{i,e})^{\textrm {rms}}$ for ions (gold) and for electrons (pink). (d) Indicator C4: isosurfaces of $\boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E} = \langle \boldsymbol {J} \cdot \boldsymbol {E} \rangle \pm 3(\boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E})^{\textrm {rms}}$ for positive $\boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E}$ (red) and negative $\boldsymbol {J} \boldsymbol {\cdot } \boldsymbol {E}$ (blue). (e) Indicator C5: isosurfaces of $E_{\parallel } = \langle |E_{\parallel } |\rangle \pm 2(|E_{\parallel }|)^{\textrm {rms}}$ for positive $E_{\parallel }$ (orange) and negative $E_{\parallel }$ (blue). Panel (f) shows, on top of the isosurfaces related to indicators C1 to C4, magnetic field lines colour-coded with $|\boldsymbol {B}|$. The magnetic field lines suggest the reconnection of a twisted flux rope with an adjacent flux rope. The white sphere of radius $1 d_i$ at $(z,x)=(77,13.5)d_{i}$ in panel (f) is a reference point that marks the position of a reconnection site. In panel (f), we also indicate the regions $R_1$ and $R_2$ defined in the text. We provide a movie of the evolution of the magnetic field lines in the supplementary material.

To visualise the change of magnetic connectivity, we trace magnetic field lines in our simulation domain. The region of most intense $|\boldsymbol {B}|$ is collocated with $R_{1}$. The magnetic field lines suggest the reconnection of a twisted flux rope with an adjacent flux rope. The white sphere of radius $1 d_i$ at $(z,x)=(77,13.5)d_{i}$ is a reference region that marks the position at which the magnetic field lines associated with the flux ropes exchange connectivity. We provide a movie to support this claim in the supplementary material. The change of connectivity between the flux ropes lasts for $\sim 96/\omega _{\textrm {pi}} \sim 0.46 \tau _{\textrm {nl}}$, which is a long time compared with the time the turbulent cascade requires to develop. The long existence of connectivity exchange and of the current structure can be associated with the suppression of nonlinearities in the current sheet. In 2-D geometries, the rate of magnetic-flux change between two magnetic islands, the so-called reconnection rate, is determined by the electric field at the $x$-point (Smith et al. 2004; Servidio et al. 2011). It can also be computed as the difference in the out-of-plane component of the magnetic vector potential between the $x$-point and the $o$-point (Franci et al. 2017; Papini et al. 2019a). In three dimensions, the reconnection rate can be computed integrating $E_{\parallel }$ along the magnetic field lines crossing the diffusion region (Schindler, Hesse & Birn 1988; Pontin 2011). However, the complex structure of the field lines makes the application of this method to our type of simulations unclear  (Liu et al. 2013; Daughton et al. 2014). An extension of 2-D methods that avoid the use of the electric field (Franci et al. 2017; Papini et al. 2019a) to the 3-D case requires the calculation of the vector potential which (i) is elaborate in 3-D PIC simulations of the type used in this study and (ii) impractical in the comparison with spacecraft data.

As the flux rope twists, it bends towards the region of changing magnetic connectivity, henceforth we refer to this region as the ‘$x$-region’. During the flux-rope bending, plasma ions are accelerated towards the $x$-region. To illustrate this behaviour, we visualise the streamlines of the ion and electron bulk velocities that leave the reconnection region. Figure 8(a) shows a view over an $xy$-plane cut of $J_{z}$. Grey colour represents negative values, red colour represents positive values, and white indicates a value of zero for $J_z$. The displayed streamlines of the ion bulk velocity emerge from the centre of the $x$-region. The streamlines are colour-coded with $v_{ix}$. The dark-blue segment near the dark-grey region indicates that the ions primarily move towards the reconnection site in the negative $x$-direction. As the ions approach the $x$-region, their speed decreases and their trajectories are deflected into the $y$-direction. The displayed streamlines maintain a coherent shape of width ${\sim } 2 d_{i}$ along the $z$-direction. Figure 8(c) shows the same ion velocity streamlines but over an $zx$-plane cut of $J_{z}$. The region where ions have large $|v_{ix}|$ coincides with the core of the twisted flux rope in Figure 8(f) (black region) which suggests that they are accelerated by the bending of the flux rope. Considering that the ion velocity streamlines are indicative of the shape of the exhaust region associated with the $x$-region, the branch of the stream lines on the right-hand side in figure 8(a) represents the reconnection exhaust of the event. It is 3-D and asymmetric. Likewise, the electron motion associated with the $x$-region is asymmetric. However, it differs considerably from the ion motion. Figure 8(b) shows the electron velocity streamlines colour-coded with $v_{ez}$ in the same view as in figure 8(a). These streamlines remain contained within a smaller region compared with the ion streamlines. They are mainly aligned with the $z$-direction. On the left-hand side of the reconnection site in figure 8(d), the electron streamlines are directed along the $J_{z}$ structure as expected since the current is mostly (but not entirely) carried by electrons. In contrast, on the right-hand side of the reconnection site, the electrons move in directions towards and away from the reconnection site as is shown by the arrows. Considering the electron velocity streamlines, the electron exhaust is also asymmetric and 3-D but smaller than the ion exhaust. The diffusion region associated with the $x$-region of the reconnection event is likely to be the large structure of positive $J_{z}$ crossing the $x$-region in the $z$-direction in figure 8(c). The shape of the electron streamlines suggests a diffusion region that resembles the distorted diffusion region observed in 3-D Hall magnetic reconnection (Drake, Shay & Swisdak 2008; Yamada et al. 2014).

Figure 8. Streamlines of the ion and electron bulk velocities over 2-D cuts of the simulation plane showing $J_{z}$. (a,b) The view over the $xy$-plane in which the $x$-direction points downward and the $y$-direction points towards the right-hand side. (c,d) The view over the $zx$-plane in which the $x$-direction points downward and the $z$-direction points towards the left-hand side. Panels (a,c) show ion bulk velocity streamlines colour-coded with $v_{ix}$. Panels (b,d) show electron velocity streamlines colour-coded with $v_{ez}$. The arrows indicate the direction of the ion bulk motion and of the electron bulk motion.

In summary, our set of indicators suggests the presence of multiple reconnection sites in our simulation domain. Our automated identification based on our indicators allows for a detailed inspection of the magnetic-field connectivity of each event. Our method searches for clusters of cells fulfilling all conditions. This approach misses events in which ions and electrons are accelerated and heated in different locations near the reconnection site. If the event is large enough to affect both ions and electrons, we expect streams of accelerated particles in both species related to the reconnection event. Given the variability in the shape and size of these particle outflows, the volume threshold must be adjusted depending on the problem at hand in different simulation set-ups.

3.4. One-dimensional trajectories across the reconnection region

In-situ measurements of spacecraft typically record the plasma and magnetic-field fluctuations along the spacecraft trajectory. In order to compare such measurements with our 3-D simulations, we ‘fly’ an artificial spacecraft through our simulation box along three trajectories, T1, T2, and T3, and record the plasma and magnetic-field fluctuations along these trajectories. According to Taylor's hypothesis, we assume that the plasma structures are static as they are convected over the spacecraft with the average solar-wind bulk speed. The trajectories are taken within the $xy$-plane and are shown as the white lines in figure 9. The trajectory T1, shown in figure 9(a), passes close to the reconnection site when it crosses the white square although it does not carry the spacecraft right through the centre of the $x$-region.

Figure 9. Trajectories of an artificial spacecraft crossing our simulation domain. (a) Trajectory T1. The spacecraft moves from the top-left corner to the bottom-right corner. This trajectory crosses a region that we identify as a reconnection exhaust. (b) Trajectories T2 and T3 are parallel toeach other. The former crosses through the reconnection site while the latter passes out of the reconnection site.

Figure 10(a) shows the plasma and magnetic-field fluctuations for our trajectory T1. We normalise these quantities to their initial values at the beginning of the simulation. Thus, the ion and electron temperatures are normalised to $T_{0}$. The magnetic field and its components are normalised to the initial background magnetic field $B_{0}$. The ion density is normalised to the initial density $n_{0}$. The ion and electron velocities are normalised to the initial Alfvén speed $V_{A0}$. The shaded area in figure 10(a) represents the region delimited by the white square in figure 9(a). The ion and electron temperatures are positively correlated with each other as well as with the density across this trajectory. The magnetic-field and ion-density fluctuations exhibit mainly anticorrelation with each other across the trajectory. This anticorrelation directly reflects the presence of slow-mode-like compressive fluctuations. The electron speed shows local peaks at $r \sim 11 d_{i}$ and $r \sim 15 d_{i}$ with no associated peaks in ion speed. This behaviour suggests the presence of local mechanisms that accelerate electrons only. This behaviour resembles electron-only reconnection events (Phan et al. 2018; Stawarz et al. 2019; Sharma Pyakurel et al. 2019; Mallet 2020). However, our indicators show that both ions and electrons interact with this reconnection region.

Figure 10. (a–c) Plasma and magnetic-field fluctuations associated with our trajectory $\textit {T1}$. (d–f) Plasma and magnetic-field fluctuations associated with our trajectory $\textit {T2}$. Panels (a,d) show the particle temperature $T_{i,e}$, magnetic field $B$, ion density $n_{i}$, and particle speed $v_{i,e}$ normalised as described in the text. The shaded areas mark the data recorded within the white squares in panel (a,b) of figure 9, respectively. Panels (b,e) show the components of the magnetic field (black) and ion velocity (red) for $\textit {T1}$ and $\textit {T2}$, respectively. Panels (c,f) show the derivative correlations $\rho _{v_{i}B}$ and $\rho _{|v||B|}$ for trajectory $\textit {T1}$ and for trajectories $\textit {T2}$ and $\textit {T3}$, respectively.

When the artificial spacecraft trajectory T1 enters the region marked with the white square in figure 9(a), it encounters a coherent structure which exhibits an enhancement in the ion and electron temperatures by a factor of approximately 1.5 to 2 compared with the background level at $r \sim 20 d_{i}$. At this position, the spacecraft observes a decrease in the magnetic field associated with an increase in the particle speed as well as an increase in the particle density. These are characteristic features associated with slow-mode-like fluctuations and shocks. Since in the trajectories shown in this section, the particle bulk speed is always less than the local magnetosonic speed, these events are not slow-mode shocks but rather fluctuations with a slow-mode-like polarisation. At $r \sim 22 d_{i}$, there is a slight enhancement in the electron speed which corresponds to the spike within the two large eddies seen in the white square in figure 10(a). At $r \sim 23 d_{i}$, the spacecraft observes another slow-mode-polarised region which corresponds to the large structure in the middle of the square. According to the Petschek (1964) model of magnetic reconnection, the exhaust of particles is limited by a pair of slow-mode shocks. However, in recent studies of reconnection in the solar wind (Phan et al. 2006, 2009; Gosling 2012), the boundaries of reconnection exhausts often lack these features. Instead, exhausts are typically characterised through a rotation in the magnetic field along with a change in the sign of the correlation between the particle speed and the magnetic field (Gosling 2012; Phan et al. 2020), consistent with our simulation results. Figure 10(b) shows from top to bottom $B_{x}$, $B_{y}$, $B_{z}$, and $|B|$ in black as well as $v_{ix}$, $v_{iy}$, $v_{iz}$, and $|v_{i}|$ in red for trajectory T1. In the shaded area (i.e., near the reconnection site), the velocity component $v_{ix}$ changes its sign between $r \sim 19 d_{i}$ and $r \sim 25 d_{i}$ while $\boldsymbol {B}$ undergoes a partial rotation. During the same interval, $v_{iy}$ shows little variation and $B_{iy}$ reverses its sign. Since the background magnetic field dominates $B_z$, the variations in the magnetic components $B_{x}$ and $B_{y}$ are more pronounced than the variations in $B_{z}$. As seen in the profile of $v_{iz}$, although ions are mostly stationary in the direction parallel to the background magnetic field, they are accelerated in the parallel direction near the slow-mode-like fluctuations. We note that the velocity spikes and magnetic-field drop-offs as seen in the $z$-component of $\boldsymbol {B}$ to a certain degree resemble the properties of the magnetic-field switchbacks observed in the solar wind (Kasper et al. 2019; McManus et al. 2020). Moreover, the blue regions in figure 4(a) suggest the presence of magnetic reversals within the simulation domain. A comparison and further study is required to establish a potential correspondence between our simulation and observational data.

To visualise the correlation between the magnetic-field and velocity components we define the derivative correlation between the two variables $v_{j}$ and $B_{j}$ as

(3.8)\begin{equation} \rho_{vBj} = \frac{\Delta v_{j}}{\Delta r }\frac{\Delta B_{j}}{\Delta r }, \end{equation}

where $\Delta r$ is a distance increment, $\Delta v_{j} = v_{j}(r + \Delta r) - v_{j}(r)$, and $\Delta B_{j} = B_{j}(r + \Delta r) - B_{j}(r)$. We use $\Delta r=0.6d_{i}$ to reduce the effect of noise when calculating the derivative while keeping the spatial step small to cover small-scale fluctuations. Figure 10(c) shows from top to bottom $\rho _{v B x}$, $\rho _{v B y}$, $\rho _{v B x}$, and $\rho _{|v||B|}$ for trajectory T1, where $\rho _{|v||B|}$ is defined accordingly with the magnitudes of $\boldsymbol {v}$ and $\boldsymbol {B}$. The $v_{ix}$ and $B_{x}$ components exhibit mostly positive correlation along the trajectory. However, there are two strong peaks of anticorrelation within the shaded area. Likewise, the $v_{iy}$ and $B_{y}$ components show more variability in the correlations from positive and negative derivative correlations within the shaded area than outside the area. This is due to the transit of the artificial spacecraft through the slow-mode-like fluctuations. In particular, around $r = 23 d_{i}$, all three components present a change from anticorrelation to positive correlation. The presence of a pair of slow-mode-like fluctuations along with a magnetic-field rotation suggests that this region is indeed an exhaust region similar to those reported in previous observational studies in the solar wind (Gosling 2012).

Trajectory T2 (the white line on the left in figure 9b) carries the spacecraft right through the centre of the $x$-region. In figure 10(d), at $r\sim 5 d_i$ and $r\sim 9 d_i$, the artificial spacecraft records particle temperature minima associated with density cavities as well as local peaks in the magnetic field. As the spacecraft moves towards the $x$-region, within the shaded region, the particle temperature remains approximately constant. There is a local minimum in the magnetic field which corresponds to the centre of the $x$-region at $r = 14 d_{i}$. On either side of the $x$-region, we find small enhancements in the electron speed. These peaks, in addition to the electron streams in figure 8(d), suggest the presence of electron-only streams in the vicinity of the $x$-region. The ion speed decreases as the spacecraft enters the $x$-region and increases as the spacecraft leaves the $x$-region. After leaving this region, the spacecraft encounters the highly twisted flux rope at $r = 16 d_{i}$ where it records an enhancement in all bulk quantities and in the magnetic field. The pair formed by the $x$-region and the closest twisted flux rope resembles the known pairs of $x$-points and magnetic islands known from 2-D models of reconnection. At the end of the trajectory, at $r\sim 23 d_i$, the spacecraft encounters a slow-mode-polarised structure which corresponds to the bright structure in the right-bottom-corner of figure 9(b). Figure 10(e) shows the components of the magnetic field and ion bulk velocity for trajectory T2. From $r =10 d_{i}$ to $r =16 d_{i}$, $B_{x}$ changes polarity and, from $r =8 d_{i}$ to $r =15 d_{i}$, $B_{y}$ undergoes a partial rotation. The change in the sign of $v_{iy}$ at the point where the spacecraft enters the shaded area and its value of approximately zero at the point where it leaves the shaded area in T2 shows a local stream of particles leaving the region along the $y$-direction. This corresponds to the right-hand side branch of the ion velocity streamline in figure 8(a). At $r =13 d_{i}$, $v_{iz}$ presents a mild peak corresponding to a weak current sheet. Entering the shaded area and up to $r \sim 19 d_{i}$, $v_{ix}$ is negative along T2 consistent with the stream of ions described in figure 8.

Trajectory T3 is parallel to trajectory T2, and the separation of these trajectories is $1.5d_{i}$. Along trajectory T3, $B_{z}$ and $B_{y}$ as well as $v_{iz}$ and $v_{iy}$ follow approximately similar behaviours (not shown here). However, the local variations are more pronounced along T2 as this trajectory crosses through the centres of multiple structures. Figure 8(f) shows the derivative correlation of the magnetic field and velocity components for trajectories T2 (black) and T3 (cyan). Trajectory T2 shows stronger positive and negative correlations in all components due to the transit through the structures. For the $x$-component, the peak of positive correlation corresponds to the transit through the flux rope which is associated with particle acceleration.

4. Discussion and conclusions

We simulate plasma turbulence created by the collision of counter-propagating Alfvén waves with a wavevector anisotropy consistent with the GS95 theory of critical balance at the small-scale end of the inertial range. Our initial waves have wavenumbers near the spectral breakpoint from the inertial to the kinetic range of turbulence. This choice allows us to set up the system with Alfvén waves and let the system develop kinetic and compressive fluctuations in the kinetic range self-consistently and with an anisotropy reminiscent of the solar wind, with the aim of developing reconnection features consistent with solar-wind turbulence. The use of a strong anisotropy in the initial waves allows the system to undergo nonlinear interactions and to create flux ropes during the first nonlinear time, which is in agreement with earlier simulation work (Grošelj et al. 2018). Our initial anisotropic set-up reduces the simulation time that a fully 3-D PIC simulation of turbulence without this imposed anisotropy would require in order to develop reconnection as a product of anisotropic turbulence.

The nonlinear interaction of the anisotropic waves self-consistently creates Alfvénic turbulence and generates magnetic-field and current-density structures such as current filaments and current sheets as part of the turbulent cascade (Howes & Nielson 2013; Howes 2015a, 2016). The initial scaling between $L_{\parallel }$ and $L_{D}$, for the magnetic structures, is $L_{\parallel } \sim L_{D}^{2/3}$. At $t=t_{R}$, the magnetic structures satisfying $V > d_{i}^{3}$ maintain an anisotropy consistent with the initial conditions and follow a $L_{\parallel } \sim L_{D}^{0.7}$ scaling. Although theoretical predictions including those based on intermittency (Boldyrev & Perez 2012; Boldyrev & Loureiro 2019), kinetic simulations (Cerri, Servidio & Califano 2017b; Cerri et al. 2019), and observations in the solar wind (Wang et al. 2020) suggest the scaling $L_{\parallel } \sim L_{D}^{2/3}$ at subproton scales, our analysis of structures with $V \leqslant d_{i}^{3}$ is more consistent with an isotropic scaling $L_{\parallel } \sim L_{D}$ which has also been observed in hybrid simulations (Franci et al. 2018; Arzamasskiy et al. 2019; Landi et al. 2019). The change of anisotropy over time (figure 3) is also observed in the evolution of the 2-D reduced power spectral density (figure 5). The anisotropy initially decreases due to the change in the mean value of the distribution of cross-section diameters and of the elongations of the magnetic structures.

The spectral index of the corresponding perpendicular 1-D power spectrum of the magnetic-field fluctuations in the kinetic range varies between $-3$ and $-4$. Meanwhile, the spectral index of the parallel power spectrum of the magnetic-field fluctuations varies from $-2$ in the interval $0.1\lesssim k_{\parallel } d_i \lesssim 0.3$ to $-4$ at subproton scales. These results show that our simulation develops an anisotropic turbulent cascade and the associated 3-D structures predicted to contribute to reconnection as a dissipation mechanism for turbulence.

The critical-balance theory of Alfvénic turbulence has been tested using gyrokinetic simulations (Howes et al. 2008a; TenBarge & Howes 2012) and 3-D PIC simulations (Grošelj et al. 2018). The evolution and morphology of 3-D reconnection events, starting from a Harris current-sheet configuration, have been studied at kinetic scales (Hesse et al. 2001; Pritchett & Coroniti 2001; Wiegelmann & Büchner 2001; Lapenta et al. 2006; Liu et al. 2013; Vapirev et al. 2013; Muñoz & Büchner 2018; Lapenta et al. 2020), as has been the effect of turbulence on the development of reconnection events (Daughton et al. 2014; Lapenta et al. 2015; Pucci et al. 2017; Papini et al. 2019b). However, little attention has been given to the occurrence of small-scale reconnection as a product of the turbulent cascade in a fully 3-D geometry. Our study contributes to the understanding, identification, and geometry of these reconnection events.

We establish a set of indicators to find regions in which magnetic reconnection takes place in 3-D PIC simulations consistent with magnetic reconnection theories. These indicators are based on the presence of current-sheet structures (C1), fast particles (C2), heated particles (C3), diffusion regions marked by energy transfer between fields and particles (C4), and non-zero parallel electric fields (C5). Since our method is based on thresholds for the bulk quantities, the selected regions correspond to high-intensity structures. Our method uses fast ions as an indicator (C2). Thus, it  does not identify all reconnection events, especially not those related to electron-only reconnection. In a follow-up study, it is worthwhile to investigate the role of the threshold level for the identification of reconnection sites and the relaxation of ion-based conditions to enable the identification of electron-only reconnection events. Our method is applicable as a first approach in the exploration of reconnection events in large 3-D PIC simulations in which the handling of the kinetic particle information is computationally expensive due to the large number of particles.

We identify three regions that fulfil our set of indicators C1 to C4 for $N_{\textrm {th}}=3$ and have an equivalent volume larger than $1d_{i}^{3}$. We also illustrate the working of our method in a subset of our simulation domain. We inspect the time evolution of the magnetic field lines and observe the change of connectivity between a highly twisted flux rope and a less twisted flux rope. We find a good agreement between the geometry of the flux ropes formed by turbulence in our simulation with the flux ropes formed by the turbulent disruption of a Harris current-sheet (Daughton et al. 2011). We observe the occurrence of a complex reconnection event in which the region of changing connectivity ($x$-region) has a volume of $\sim 12.5 d_{i}^{3}$. This event dissipates turbulent fluctuations in current structures of the order of a few $d_{i}$, which are smaller than the smallest events recently observed in the solar wind (Phan et al. 2020) and different from the events observed in space. These latter events are mostly very large interface regions between plasmas (Phan et al. 2006; Gosling 2007). The occurrence of electron-only reconnection (Phan et al. 2018; Stawarz et al. 2019) and electron-scale turbulent fluctuations suggests that events as the one we describe take place in the solar wind.

Although there is a good agreement between studies using the Harris configuration and solar-wind observations (Mistry et al. 2016), our event is considerably more complex than the idealised steady and non-turbulent Harris current-sheet configuration often invoked to study magnetic reconnection. The shape of our reconnection region is asymmetric, and the regions in which particle heating and acceleration occur are mostly associated with current filaments rather than current sheets. This suggests that the twist of the flux ropes plays a crucial role for the particle heating in our simulation. In addition, this finding supports the notion that reconnection events occur in the solar wind through small-scale flux ropes (Crooker et al. 1996; Moldwin et al. 2000).

We trace 1-D artificial-spacecraft trajectories across the simulation domain to study the fluctuations in the quantities $n_{i}$, $\boldsymbol {v}_{i,e}$, $T_{i,e}$, and $\boldsymbol {B}$. These samplings may facilitate direct comparisons between our simulations and spacecraft observations in the solar wind. Our trajectories T1 and T3 pass near the identified reconnection region, and our trajectory T2 crosses through the centre of the $x$-region. We observe the presence of slow-mode-polarised fluctuations as anticorrelated fluctuations in $n_{i}$ and $|B|$, rotations in the magnetic field, and changes in the sign of the correlation between the magnetic field and the ion velocity consistent with reconnection exhausts observed in the solar wind (Gosling 2012). Our artificial-spacecraft trajectory T2 (figure 10d) shows an enhancement in all bulk quantities, which may be associated with a reconnecting flux rope. Moreover, this trajectory suggests that the encounter of a magnetic minimum followed by an enhancement in all bulk quantities may be associated with the encounter of an $x$-region and a flux rope. Such a pair $x$-region–flux-rope corresponds to the traditional pair $x$-point–$o$-point in 2-D models of reconnection. It would be worthwhile to compare our simulated spacecraft trajectories with spacecraft observations of small-scale reconnection events and reconnection exhausts in the solar wind. The instrumentation on-board Solar Orbiter and Parker Solar Probe has the appropriate time resolution for such a comparison.

In our reconnection event, ions and electrons behave differently as shown in figure 8. Both ions and electrons move towards and away from the $x$-region but in different directions. Our trajectories in the vicinity of the reconnection event suggest that the slow-mode-like features associated with the partial rotation in the magnetic field and the change in the $\boldsymbol {v}_i$–$\boldsymbol {B}$ correlation are also present in these spontaneously created small-scale events.

The finite number of particles per cell has an important effect on the determination of coherent regions of strong $E_{\parallel }$, our indicator C5. Therefore, C5 is not a good indicator when the number of particles per cell is ${\lesssim }100$. Although 2-D studies of turbulence, magnetic reconnection (Franci et al. 2020), and plasma instabilities (Hellinger & Štverák 2018) are able to use considerably larger numbers of particles per cell (${\sim }1000$), our work requires the third dimension in order to model the turbulence and the complex reconnection geometry more appropriately (Howes 2015b; Lazarian et al. 2020). Nonetheless, the increasing computational power of high-performance-computing facilities will allow us to perform increasingly more accurate 3-D PIC simulations and to test all of our indicators over a wider range of parameters. Before these methods become computationally viable, divergence-cleaning of the electric field (Jacobs & Hesthaven 2009) is a possible route to reduce the effect of particle noise.

Our data set possibly includes further reconnection sites that can be studied in more detail in the future. In this project, we use bulk quantities to study the reconnection events. In our future work, it will be interesting to study the changes in the particle distribution functions as a result of the identified small-scale reconnection events. Such a more detailed study of the associated particle kinetics will allow us to understand the energy exchange between fields and particles and the details of the energy dissipation through small-scale reconnection events in the solar wind.