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

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

(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 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 ¹ 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$.Footnote ² 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 Mitchell2010a; 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 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 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 Hellinger2017a). 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 Oughton2014a; Papini et al. Reference Papini, Franci, Landi, Verdini, Matteini and Hellinger2019a). 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.

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

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

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. 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 Hellinger2019a). 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 Hellinger2019a) 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).

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

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 (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

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

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.