## 1. Introduction

The solar wind is a low-collisionality plasma produced in the solar corona (Marsch Reference Marsch2006). 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 Reference Gazis and Lazarus1982) which suggests the presence of local heating and particle-acceleration mechanisms (Goldstein *et al.* Reference Goldstein, Wicks, Perri and Sahraoui2015). 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 Reference Coleman1968; Marsch & Tu Reference Marsch and Tu1990) 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 Reference Kiyani, Osman and Chapman2015); an inertial range in which the power index varies from $-3/2$ to $-5/3$ (Iroshnikov Reference Iroshnikov1963; Marsch & Tu Reference Marsch and Tu1990; Podesta Reference Podesta2009; Boldyrev *et al.* Reference Boldyrev, Perez, Borovsky and Podesta2011); and a dissipation range in which the power index is less clearly defined (Goldstein, Roberts & Fitch Reference Goldstein, Roberts and Fitch1994; Li, Gary & Stawicki Reference Li, Gary and Stawicki2001; Howes *et al.* Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008*b*) with spectral breaks at electron scales (Alexandrova *et al.* Reference Alexandrova, Saur, Lacombe, Mangeney, Mitchell, Schwartz and Robert2009; Sahraoui *et al.* Reference Sahraoui, Goldstein, Robert and Khotyaintsev2009). 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.* Reference Tatsuno, Dorland, Schekochihin, Plunk, Barnes, Cowley and Howes2009; Eyink Reference Eyink2018; Verscharen, Klein & Maruca Reference Verscharen, Klein and Maruca2019). 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 Reference Marsch, Vocks and Tu2003; Kasper, Lazarus & Gary Reference Kasper, Lazarus and Gary2008) and non-resonant heating such as stochastic heating (Chandran *et al.* Reference Chandran, Li, Rogers, Quataert and Germaschewski2010, Reference Chandran, Verscharen, Quataert, Kasper, Isenberg and Bourouaine2013). 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 Reference Smith, Vasquez and Hollweg2011; Podesta & TenBarge Reference Podesta and TenBarge2012; Salem *et al.* Reference Salem, Howes, Sundkvist, Bale, Chaston, Chen and Mozer2012; Podesta Reference Podesta2013; Goldstein *et al.* Reference Goldstein, Wicks, Perri and Sahraoui2015). Another mechanism proposed to carry the turbulent cascade to subproton scales is magnetic reconnection (Sundkvist *et al.* Reference Sundkvist, Retinò, Vaivads and Bale2007; Franci *et al.* Reference Franci, Cerri, Califano, Landi, Papini, Verdini, Matteini, Jenko and Hellinger2017; Loureiro & Boldyrev Reference Loureiro and Boldyrev2020).

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 Reference Hesse and Schindler1988; Pontin Reference Pontin2011). 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 Reference Priest and Forbes2007; Zweibel & Yamada Reference Zweibel and Yamada2009). In the latter, reconnection events are characterised by streams of particles associated with Alfvénic disturbances and magnetic-field rotations (Gosling *et al.* Reference Gosling, Skoug, McComas and Smith2005; Davis *et al.* Reference Davis, Phan, Gosling and Skoug2006; Gosling, Eriksson & Schwenn Reference Gosling, Eriksson and Schwenn2006; Phan *et al.* Reference Phan, Gosling, Davis, Skoug, Øieroset, Lin, Lepping, McComas, Smith and Reme2006; Phan, Gosling & Davis Reference Phan, Gosling and Davis2009; Gosling Reference Gosling2012; Phan *et al.* Reference Phan, Bale, Eastwood, Lavraud, Drake, Oieroset, Shay, Pulupa, Stevens and MacDowall2020). 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 Reference Matthaeus and Lamkin1986; Strauss Reference Strauss1988; Lazarian & Vishniac Reference Lazarian and Vishniac1999; Kim & Diamond Reference Kim and Diamond2001; Servidio *et al.* Reference Servidio, Dmitruk, Greco, Wan, Donato, Cassak, Shay, Carbone and Matthaeus2011; Boldyrev & Loureiro Reference Boldyrev and Loureiro2017; Adhikari *et al.* Reference Adhikari, Shay, Parashar, Pyakurel, Matthaeus, Godzieba, Stawarz, Eastwood and Dahlin2020; Loureiro & Boldyrev Reference Loureiro and Boldyrev2020). 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.* Reference Franci, Cerri, Califano, Landi, Papini, Verdini, Matteini, Jenko and Hellinger2017; Boldyrev & Loureiro Reference Boldyrev and Loureiro2017; Cerri & Califano Reference Cerri and Califano2017; Papini, Landi & Del Zanna Reference Papini, Landi and Del Zanna2019*b*), 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.* Reference Boldyrev, Horaites, Xia and Perez2013; Mallet, Schekochihin & Chandran Reference Mallet, Schekochihin and Chandran2017). 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 Reference Boldyrev and Loureiro2017). 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.* Reference Birn, Drake, Shay, Rogers, Denton, Hesse, Kuznetsova, Ma, Bhattacharjee and Otto2001; Shay *et al.* Reference Shay, Drake, Rogers and Denton2001; Loureiro *et al.* Reference Loureiro, Uzdensky, Schekochihin, Cowley and Yousef2009; Servidio *et al.* Reference Servidio, Matthaeus, Shay, Cassak and Dmitruk2009, Reference Servidio, Matthaeus, Shay, Dmitruk, Cassak and Wan2010; Bessho *et al.* Reference Bessho, Chen, Hesse and Wang2017) and in three dimensions (Hesse, Kuznetsova & Birn Reference Hesse, Kuznetsova and Birn2001; Pritchett & Coroniti Reference Pritchett and Coroniti2001; Lapenta Reference Lapenta2003; Lapenta *et al.* Reference Lapenta, Krauss-Varban, Karimabadi, Huba, Rudakov and Ricci2006; Kowal *et al.* Reference Kowal, Lazarian, Vishniac and Otmianowska-Mazur2009; Daughton *et al.* Reference Daughton, Roytershteyn, Karimabadi, Yin, Albright, Bergen and Bowers2011; Baumann, Galsgaard & Nordlund Reference Baumann, Galsgaard and Nordlund2013; Liu *et al.* Reference Liu, Daughton, Karimabadi, Li and Roytershteyn2013; Pritchett Reference Pritchett2013; Muñoz & Büchner Reference Muñoz and Büchner2018). 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 Reference Lapenta2012; Germaschewski *et al.* Reference Germaschewski, Fox, Abbott, Ahmadi, Maynard, Wang, Ruhl and Bhattacharjee2016). 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.* Reference Told, Cookmeyer, Muller, Astfalk and Jenko2016) 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.* Reference Shay, Drake, Rogers and Denton2001; Scholer *et al.* Reference Scholer, Sidorenko, Jaroschek, Treumann and Zeiler2003; Ricci *et al.* Reference Ricci, Brackbill, Daughton and Lapenta2004; Shay *et al.* Reference Shay, Drake, Swisdak and Rogers2004; Daughton, Scudder & Karimabadi Reference Daughton, Scudder and Karimabadi2006; Daughton *et al.* Reference Daughton, Roytershteyn, Karimabadi, Yin, Albright, Bergen and Bowers2011; Leonardis *et al.* Reference Leonardis, Chapman, Daughton, Roytershteyn and Karimabadi2013; Liu *et al.* Reference Liu, Daughton, Karimabadi, Li and Roytershteyn2013; Beresnyak Reference Beresnyak2016; Goldman, Newman & Lapenta Reference Goldman, Newman and Lapenta2016). 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 Reference Howes and Nielson2013; Howes Reference Howes2015*a*), current-sheet structures (Howes Reference Howes2016), 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.* Reference Germaschewski, Fox, Abbott, Ahmadi, Maynard, Wang, Ruhl and Bhattacharjee2016) 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 (Reference Sridhar and Goldreich1994) and Goldreich & Sridhar (Reference Goldreich and Sridhar1995), 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.* Reference Verscharen, Parashar, Gary and Klein2020) 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.* Reference Franci, Cerri, Califano, Landi, Papini, Verdini, Matteini, Jenko and Hellinger2017). The r.m.s. of a quantity $\psi$ is defined as

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 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 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.* Reference Franci, Landi, Verdini, Matteini and Hellinger2018; Arzamasskiy *et al.* Reference Arzamasskiy, Kunz, Chandran and Quataert2019; Landi *et al.* Reference Landi, Franci, Papini, Verdini, Matteini and Hellinger2019).

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 2*b*). 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.

### 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.* (Reference Franci, Landi, Verdini, Matteini and Hellinger2018), we calculate the energy associated with the 3-D Fourier modes $\psi _{3D}(\boldsymbol {k})$ of a quantity $\psi$ as

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

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 Reference Cho and Vishniac2000). 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 Reference Podesta2009; Chen *et al.* Reference Chen, Mallet, Yousef, Schekochihin and Horbury2011; TenBarge *et al.* Reference TenBarge, Podesta, Klein and Howes2012; Oughton *et al.* Reference Oughton, Matthaeus, Wan and Osman2015; Gerick, Saur & von Papen Reference Gerick, Saur and von Papen2017), and the development of an anisotropic energy cascade is sufficient for the determination of reconnection events in the present study.Footnote ^{1} 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

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

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,

and the parallel 1-D reduced power spectral density,

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$.Footnote ^{2} 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.* Reference Chen, Horbury, Schekochihin, Wicks, Alexandrova and Mitchell2010*a*; Bruno, Trenchi & Telloni Reference Bruno, Trenchi and Telloni2014) as well as in hybrid simulations (Franci *et al.* Reference Franci, Landi, Verdini, Matteini and Hellinger2018; González *et al.* Reference González, Parashar, Gomez, Matthaeus and Dmitruk2019).

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 Reference Bavassano and Bruno1989; Grappin, Velli & Mangeney Reference Grappin, Velli and Mangeney1991; Wicks *et al.* Reference Wicks, Horbury, Chen and Schekochihin2010, Reference Wicks, Horbury, Chen and Schekochihin2011; Chen *et al.* Reference Chen, Mallet, Yousef, Schekochihin and Horbury2011). 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.* Reference Franci, Landi, Verdini, Matteini and Hellinger2018; Arzamasskiy *et al.* Reference Arzamasskiy, Kunz, Chandran and Quataert2019; Cerri, Groselj & Franci Reference Cerri, Groselj and Franci2019; Landi *et al.* Reference Landi, Franci, Papini, Verdini, Matteini and Hellinger2019). 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.* Reference Cerri, Franci, Califano, Landi and Hellinger2017*a*). The level of compressive fluctuations in our simulation is greater than observed in the solar wind (Chen Reference Chen2016), 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.* Reference Wan, Rappazzo, Matthaeus, Servidio and Oughton2014*a*; Papini *et al.* Reference Papini, Franci, Landi, Verdini, Matteini and Hellinger2019*a*). 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.* Reference Daughton, Roytershteyn, Karimabadi, Yin, Albright, Bergen and Bowers2011; Liu *et al.* Reference Liu, Daughton, Karimabadi, Li and Roytershteyn2013; Daughton *et al.* Reference Daughton, Nakamura, Karimabadi, Roytershteyn and Loring2014). 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 Reference Priest and Démoulin1995). 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.* Reference Eastwood, Shay, Phan and Øieroset2010, Reference Eastwood, Phan, Øieroset, Shay, Malakit, Swisdak, Drake and Masters2013). Moreover, proton temperature anisotropies in reconnection events can trigger kinetic instabilities, which then have a stabilising effect on the current sheet (Matteini *et al.* Reference Matteini, Hellinger, Goldstein, Landi, Velli and Neugebauer2013).

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 Reference Hesse, Kuznetsova and Birn2004). 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 Reference Pritchett and Coroniti2004; Pontin Reference Pontin2011). 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 Reference Somov and Titov1985; Ni *et al.* Reference Ni, Lin, Roussev and Schmieder2016). 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 Reference Hesse and Schindler1988; Priest, Hornig & Pontin Reference Priest, Hornig and Pontin2003; Pontin Reference Pontin2011), 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.* Reference Uritsky, Pouquet, Rosenberg, Mininni and Donovan2010) 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.

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

^{3}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}}$.

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

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}}$.

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.

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.* Reference Phan, Eastwood, Shay, Drake, Sonnerup, Fujimoto, Cassak, Øieroset, Burch and Torbert2018). 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.

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.* Reference Smith, Ghosh, Dmitruk and Matthaeus2004; Servidio *et al.* Reference Servidio, Dmitruk, Greco, Wan, Donato, Cassak, Shay, Carbone and Matthaeus2011). 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.* Reference Franci, Cerri, Califano, Landi, Papini, Verdini, Matteini, Jenko and Hellinger2017; Papini *et al.* Reference Papini, Franci, Landi, Verdini, Matteini and Hellinger2019*a*). In three dimensions, the reconnection rate can be computed integrating $E_{\parallel }$ along the magnetic field lines crossing the diffusion region (Schindler, Hesse & Birn Reference Schindler, Hesse and Birn1988; Pontin Reference Pontin2011). However, the complex structure of the field lines makes the application of this method to our type of simulations unclear (Liu *et al.* Reference Liu, Daughton, Karimabadi, Li and Roytershteyn2013; Daughton *et al.* Reference Daughton, Nakamura, Karimabadi, Roytershteyn and Loring2014). An extension of 2-D methods that avoid the use of the electric field (Franci *et al.* Reference Franci, Cerri, Califano, Landi, Papini, Verdini, Matteini, Jenko and Hellinger2017; Papini *et al.* Reference Papini, Franci, Landi, Verdini, Matteini and Hellinger2019*a*) 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 Reference Drake, Shay and Swisdak2008; Yamada *et al.* Reference Yamada, Yoo, Jara-Almonte, Ji, Kulsrud and Myers2014).

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 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.* Reference Phan, Eastwood, Shay, Drake, Sonnerup, Fujimoto, Cassak, Øieroset, Burch and Torbert2018; Stawarz *et al.* Reference Stawarz, Eastwood, Phan, Gingell, Shay, Burch, Ergun, Giles, Gershman and Le Contel2019; Sharma Pyakurel *et al.* Reference Sharma Pyakurel, Shay, Phan, Matthaeus, Drake, TenBarge, Haggerty, Klein, Cassak and Parashar2019; Mallet Reference Mallet2020). However, our indicators show that both ions and electrons interact with this reconnection region.

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 (Reference Petschek1964) 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.* Reference Phan, Gosling, Davis, Skoug, Øieroset, Lin, Lepping, McComas, Smith and Reme2006, Reference Phan, Gosling and Davis2009; Gosling Reference Gosling2012), 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 Reference Gosling2012; Phan *et al.* Reference Phan, Bale, Eastwood, Lavraud, Drake, Oieroset, Shay, Pulupa, Stevens and MacDowall2020), 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.* Reference Kasper, Bale, Belcher, Berthomier, Case, Chandran, Curtis, Gallagher, Gary and Golub2019; McManus *et al.* Reference McManus, Bowen, Mallet, Chen, Chandran, Bale, Larson, de Wit, Kasper and Stevens2020). 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

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 Reference Gosling2012).

Trajectory *T2* (the white line on the left in figure 9*b*) 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.* Reference Grošelj, Mallet, Loureiro and Jenko2018). 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 Reference Howes and Nielson2013; Howes Reference Howes2015*a*, Reference Howes2016). 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 Reference Boldyrev and Perez2012; Boldyrev & Loureiro Reference Boldyrev and Loureiro2019), kinetic simulations (Cerri, Servidio & Califano Reference Cerri, Servidio and Califano2017*b*; Cerri *et al.* Reference Cerri, Groselj and Franci2019), and observations in the solar wind (Wang *et al.* Reference Wang, He, Alexandrova, Dunlop and Perrone2020) 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.* Reference Franci, Landi, Verdini, Matteini and Hellinger2018; Arzamasskiy *et al.* Reference Arzamasskiy, Kunz, Chandran and Quataert2019; Landi *et al.* Reference Landi, Franci, Papini, Verdini, Matteini and Hellinger2019). 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.* Reference Howes, Dorland, Cowley, Hammett, Quataert, Schekochihin and Tatsuno2008*a*; TenBarge & Howes Reference TenBarge and Howes2012) and 3-D PIC simulations (Grošelj *et al.* Reference Grošelj, Mallet, Loureiro and Jenko2018). 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.* Reference Hesse, Kuznetsova and Birn2001; Pritchett & Coroniti Reference Pritchett and Coroniti2001; Wiegelmann & Büchner Reference Wiegelmann and Büchner2001; Lapenta *et al.* Reference Lapenta, Krauss-Varban, Karimabadi, Huba, Rudakov and Ricci2006; Liu *et al.* Reference Liu, Daughton, Karimabadi, Li and Roytershteyn2013; Vapirev *et al.* Reference Vapirev, Lapenta, Divin, Markidis, Henri, Goldman and Newman2013; Muñoz & Büchner Reference Muñoz and Büchner2018; Lapenta *et al.* Reference Lapenta, Pucci, Goldman and Newman2020), as has been the effect of turbulence on the development of reconnection events (Daughton *et al.* Reference Daughton, Nakamura, Karimabadi, Roytershteyn and Loring2014; Lapenta *et al.* Reference Lapenta, Markidis, Goldman and Newman2015; Pucci *et al.* Reference Pucci, Servidio, Sorriso-Valvo, Olshevsky, Matthaeus, Malara, Goldman, Newman and Lapenta2017; Papini *et al.* Reference Papini, Landi and Del Zanna2019*b*). 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.* Reference Daughton, Roytershteyn, Karimabadi, Yin, Albright, Bergen and Bowers2011). 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.* Reference Phan, Bale, Eastwood, Lavraud, Drake, Oieroset, Shay, Pulupa, Stevens and MacDowall2020) and different from the events observed in space. These latter events are mostly very large interface regions between plasmas (Phan *et al.* Reference Phan, Gosling, Davis, Skoug, Øieroset, Lin, Lepping, McComas, Smith and Reme2006; Gosling Reference Gosling2007). The occurrence of electron-only reconnection (Phan *et al.* Reference Phan, Eastwood, Shay, Drake, Sonnerup, Fujimoto, Cassak, Øieroset, Burch and Torbert2018; Stawarz *et al.* Reference Stawarz, Eastwood, Phan, Gingell, Shay, Burch, Ergun, Giles, Gershman and Le Contel2019) 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.* Reference Mistry, Eastwood, Haggerty, Shay, Phan, Hietala and Cassak2016), 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.* Reference Crooker, Burton, Phillips, Smith and Balogh1996; Moldwin *et al.* Reference Moldwin, Ford, Lepping, Slavin and Szabo2000).

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 Reference Gosling2012). Our artificial-spacecraft trajectory *T2* (figure 10*d*) 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.* Reference Franci, Stawarz, Papini, Hellinger, Nakamura, Burgess, Landi, Verdini, Matteini and Ergun2020), and plasma instabilities (Hellinger & Štverák Reference Hellinger and Štverák2018) 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 Reference Howes2015*b*; Lazarian *et al.* Reference Lazarian, Eyink, Jafari, Kowal, Li, Xu and Vishniac2020). 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 Reference Jacobs and Hesthaven2009) 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.

## Supplementary material and movies

Supplementary material and movies are available at https://doi.org/10.1017/S0022377821000404.

## Acknowledgements

The authors would like to thank the three anonymous reviewers whose thoughtful comments have led to significant improvements in our manuscript.

*Editor Alex Schekochihin thanks the referees for their advice in evaluating this article.*

## Funding

This work was supported by the European Space Agency's Networking/Partnering Initiative (NPI) programme and the Colombian programme Pasaporte a la Ciencia, Foco Sociedad – Reto 3 (Educación de calidad desde la ciencia, la tecnología y la innovación (CTel)), ICETEX (J.A.A.R., grant numbers 4000127929/19/NL/MH/mg, 3933061); STFC Ernest Rutherford Fellowship (D.V., grant number ST/P003826/1); STFC Consolidated Grant (D.V., grant number ST/S000240/1), (R.T.W., grant number ST/S000240/1), (C.J.O., grant number ST/S000240/1) and (G.N., grant number ST/S000240/1), and NSF grant (K.G., grant number AGS-1460190). This work was performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC Capital Grants ST/K000373/1 and ST/R002363/1, and STFC DiRAC Operations Grant ST/R001014/1. DiRAC is part of the National e-Infrastructure.

## Declaration of interests

The authors report no conflict of interest.

## Data availability statement

The data that support the findings of this study are openly available in Zenodo at reference number 10.5281/zenodo.4313310

## Appendix A. Initial conditions of the simulation

We initialise our simulation with eight anisotropic low-frequency counter-propagating Alfvén waves within a box of volume $L_{x} \times L_{y} \times L_{z}$. We take the background magnetic field to be along the $z$-axis, $\boldsymbol {B}_{0}=B_{0} \hat {z}$, and set up our fluctuations with wavevectors following the theory of critical balance by GS95. According to GS95, turbulence is isotropic at the large-scale end of the inertial range and develops an anisotropic cascade of energy with respect to the local magnetic field. The anisotropic cascade of energy is associated with a wavevector anisotropy $k_{\parallel } \propto ( |\boldsymbol {k}_{\perp }| )^{\gamma }$, where $k_{\parallel }$ and $k_{\perp }$ are the wavevector components in the directions parallel and perpendicular with respect to the local background magnetic field. The index $\gamma$ is a power index that is approximately constant in each wavevector range of the turbulent power spectrum. For the inertial range, $\gamma =2/3$. Since the fluctuations are isotropic at the large-scale end of the inertial range, we express the relation between $k_{\perp }$ and $k_{\parallel }$ as

where $C$ is a constant which is chosen so that $k_{\parallel }=k_{\perp }$ at the large-scale end of the inertial range, which we set up as $k_{\perp }d_{i}=10^{-4}$ consistent with observations (Wicks *et al.* Reference Wicks, Horbury, Chen and Schekochihin2010; Chen *et al.* Reference Chen, Mallet, Schekochihin, Horbury, Wicks and Bale2012). We define $k_{m,\perp } = \sqrt {k_{m,x}^{2} +k_{m,y}^{2}}$, where the index $m$ refers to the mode of the wave. Since we use periodic boundary conditions in our simulation, we adjust the wavelengths of our initial modes $\lambda _{m,i}$ so that $L_i$ is an integer multiple of $\lambda _{m,i}$. Then, the wavevector components are

*a–c*)\begin{equation} k_{m,x}=m\frac{2{\rm \pi}}{L_{x}} \quad k_{m,y}=m\frac{2{\rm \pi}}{L_{y}} \quad \text{and} \quad k_{m,z}=m\frac{2{\rm \pi}}{L_{z}}. \end{equation}

Since we use only $m = 1$, we drop the index $m$ for simplicity. Each wave satisfies the Alfvénic polarisation relation

where $V_{A0}=B_{0}/\sqrt {\mu _{0}n_{i}m_{i}}$ is the Alfvén speed, $n_{i}$ is the ion density and $m_{i}$ is the ion mass. Here $\boldsymbol {u}_{s,\alpha }$ is the bulk velocity of the species $s$. The index $\alpha =1,\ldots ,8$ refers to each wave. The four waves with odd $\alpha$ travel along the $z$-direction and the other four in the opposite direction. The amplitude $\boldsymbol {A}_{\alpha }$ of the perturbation $\delta \boldsymbol {B}_{\alpha }$ of each wave is perpendicular to both the background magnetic field $\boldsymbol {B}_{0}$ and to the wave's wavevector $\boldsymbol {k}_{\alpha }$. Thus, we write the components of the wavevector as

and

where $\phi _\alpha$ is the azimuthal angle of $\boldsymbol {k}_{\alpha ,\perp }$. The waves propagating in the $+z-$direction have $\phi _{\alpha }=0, {\rm \pi}, {\rm \pi}/4$ and $5{\rm \pi} /4$, whereas the waves propagating in the $-z$-direction have $\phi _{\alpha }={\rm \pi} /2, 3{\rm \pi} /2, 3{\rm \pi} /4$ and $7{\rm \pi} /4$. This distribution of azimuthal angles produces a quasi-gyrotropic distribution of fluctuations in the plane perpendicular to the background magnetic field while keeping the initial magnetic field divergence-free. The components of the fluctuating fields for each wave are given by

and

where $\psi _{\alpha }$ represents a random phase for each $\alpha$. The amplitude $|\boldsymbol {A}_{\alpha }|$, according to Chandran *et al.* (Reference Chandran, Li, Rogers, Quataert and Germaschewski2010), follows

Thus, the components of the total initial magnetic variations are

*a,b*)\begin{equation} \delta B_{T,x} = D\sum_{\alpha=1}^{8}\delta B_{\alpha,x} \quad \text{and} \quad \delta B_{T,y} = D\sum_{\alpha=1}^{8}\delta B_{\alpha,y}, \end{equation}

where $D$ is a normalisation constant defined as

which ensures that the total amplitude of all modes $|\delta \boldsymbol {B}_{T} / B_{0}| \sim 1$ at the beginning of the simulation. We assume that the nonlinear time is comparable to the linear time at the initial time, thus we initialise the simulation with strong turbulence. The nonlinearity parameter $\chi =(\delta B_{T} /B_{0})/(k_{\parallel }/k_{\perp }) \sim L_{Z}/L_{x} \sim 5.2$ at the initial time which quantitatively states that the initialised turbulence is strong. The components of the velocity fluctuations $\delta \boldsymbol {u}_{T}$ are calculated self-consistently according to (A 3).

The wavelengths of the initial waves at $k_{\perp }d_{i}=1$, are $\lambda _{\perp } = 2{\rm \pi} d_{i}$ and $\lambda _{\parallel } = 2{\rm \pi} /10^{-4/3} d_{i}$. Therefore, the size of the box required to simulate our initial ($m=1$) anisotropic Alfvén waves is $L_{z}=\lambda _{\parallel }$ and $L_{x}=L_{y}=\sqrt {2}\lambda _{\perp }$. However, we use $L_{z}=125d_{i}$, $L_{x}=L_{y}=24d_{i}$, $\lambda _{\parallel } = L_z$, and $\lambda _{\perp } = \sqrt {2}L_{x}/4$. This choice keeps the ratio $\lambda _{\perp }/\lambda _{\parallel }\approx 10^{-4/3}$ while allowing a wider spatial evolution in the perpendicular direction.

The theoretical critical-balance scaling $k_{\parallel } \sim k_{\perp }^{2/3}$ applies to Alfvén waves in the inertial range. The initial fluctuations in our simulation have $k_\perp d_{i} \sim 1$ which is at the transition scale from the inertial to the dissipation range. Natural fluctuations at this scale have an anisotropy consistent with the critical-balance scaling based on the size of the inertial range (Wicks *et al.* Reference Wicks, Horbury, Chen and Schekochihin2010). The scale dependence of the anisotropy in the inertial range also varies when considering dynamic alignment and intermittency (Cho & Lazarian Reference Cho and Lazarian2004; Boldyrev *et al.* Reference Boldyrev, Perez, Borovsky and Podesta2011; Chandran, Schekochihin & Mallet Reference Chandran, Schekochihin and Mallet2015; Chen Reference Chen2016). We assume a critical-balance scaling over an inertial range of four decades to capture the relative amplitude of the anisotropy without simulating the true evolution of the inertial-range turbulence. Therefore, we initialise with fluctuations at $k_{\perp }d_{i} \sim 1$ that have such an anisotropy. The wavevector anisotropy in the dissipation range is less well understood and, at kinetic scales, it is not clear whether the turbulence is mostly carried by KAWs, whistler waves, or a combination of compressive and non-compressive modes (Schekochihin *et al.* Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009; Chen *et al.* Reference Chen, Wicks, Horbury and Schekochihin2010*b*; Boldyrev & Perez Reference Boldyrev and Perez2012). Moreover, pressured-balanced structures also contribute to the turbulent cascade (Verscharen *et al.* Reference Verscharen, Marsch, Motschmann and Müller2012; Narita & Marsch Reference Narita and Marsch2015; Verscharen, Chen & Wicks Reference Verscharen, Chen and Wicks2017). Nevertheless, our anisotropic initialisation is supported by solar-wind measurements (Horbury, Forman & Oughton Reference Horbury, Forman and Oughton2008; Alexandrova *et al.* Reference Alexandrova, Saur, Lacombe, Mangeney, Mitchell, Schwartz and Robert2009; Wicks *et al.* Reference Wicks, Horbury, Chen and Schekochihin2010, Reference Wicks, Horbury, Chen and Schekochihin2011) and allows a kinetic cascade to develop self-consistently as the simulation evolves.

## Appendix B. Second-order structure function

Following Cho & Vishniac (Reference Cho and Vishniac2000), we define the local magnetic field between two points $\boldsymbol {r}_{1}$ and $\boldsymbol {r}_{2}$ as

We define the coordinate parallel to the local magnetic field $\boldsymbol {B}_{l}$ as $r_{\parallel }=\hat {z}\boldsymbol {\cdot } (\boldsymbol {r}_{2} - \boldsymbol {r}_{1})$ and the coordinate perpendicular as $r_{\perp }=|\hat {z} \times (\boldsymbol {r}_{2} - \boldsymbol {r}_{1})|$, where $\hat {z}=\boldsymbol {B}_{l}/|\boldsymbol {B}_{l}|$. With these definitions, we calculate the second-order structure function of the magnetic fluctuations $\boldsymbol {b}(\boldsymbol {r_{1}}) = \boldsymbol {B}_{l} - \boldsymbol {B}(\boldsymbol {r}_{1})$ as

where $\langle \ \rangle$ represents the average over the spatial domain. In order to discretise the $r_{\perp }r_{\parallel }$-plane, we calculate the values of $r_{\perp }$, $r_{\parallel }$, and $Fb2$ for each pair of points $\boldsymbol {r}_{1},\boldsymbol {r}_{2}$. Then, for each pixel, we calculate the mean value as the sum of all $Fb2$ divided by the number of combinations ($\boldsymbol {r}_1$,$\boldsymbol {r}_2$) in each pixel. We apply a filter to remove the pixels with less than $\sqrt N$ combinations, where $N$ is the total number of combinations in the $r_{\perp }$, $r_{\parallel }$ space.

Figure 11 shows $\log (Fb2)$ in the $r_{\perp },r_{\parallel }$-plane for the time steps $t=0, t=12/\omega _{\textrm {pi}}, t=t_{\mathrm R}=120/\omega _{\textrm {pi}}$ and $t=240/\omega _{\textrm {pi}}$. At $t=12/\omega _{\textrm {pi}}$, the structure function indicates a perpendicular cascade of the magnetic energy. On the other hand, the structure function does not give evidence of a strong parallel cascade and is, instead, still consistent with our initial conditions in terms of the parallel extent of the magnetic-field fluctuations. At $t=t_{\mathrm R}=120/\omega _{\textrm {pi}}$, the green horizontal structure suggests that the magnetic energy has been redistributed and cascaded to smaller parallel scales. The analysis of the structure functions is consistent with our analysis of the Fourier spectra in figure 5.