3 Results

The evolution of the simulations in time can be summarized as follows (Burgess et al. Reference Burgess, Gingell and Matteini2016): (i) linear growth of the tearing instability from the symmetric, highly ordered initial state; (ii) nonlinear growth of the tearing instability, characterized by island merging and reconnection across several current sheets; (iii) saturation of the tearing instability, and relaxation towards a chaotic, ‘turbulent’ state. The magnetic structure at the end of the simulation for the AP, GF and SH cases is shown in figure 1. These figures visualize the field line topology (but not field strength or direction) using a line integral convolution, for which a randomly generated field of noise is smeared out over magnetic field lines. This provides a dense representation of the most important features we discuss in this paper, i.e. the development of magnetic islands of varying scales, without the clutter of following multiple field lines in three dimensions. The time evolution of the tearing instability and associated magnetic structure for these multiple current sheet systems is described in detail by Burgess et al. (Reference Burgess, Gingell and Matteini2016). Discussion of the effects of other instabilities on the evolution of single three-dimensional ion-scale current sheets, including the drift-kink, firehose and ion cyclotron instabilities, can be found in Gingell et al. (Reference Gingell, Burgess and Matteini2015).

We find that growth and merging of magnetic islands by reconnection in the AP case generates an apparently three-dimensional and isotropic cascade of magnetic vortices. In the GF case, the magnetic structure remains largely two-dimensional, and the $x$ – $z$ and $y$ – $z$ planes in figure 1(b) are therefore dominated by the $z$ -directed mean field.

The introduction of super-Alfvénic shear in the SH case stabilizes the tearing instability (e.g. Chen & Morrison Reference Chen and Morrison1990; Cassak & Otto Reference Cassak and Otto2011; Landi & Bettarini Reference Landi and Bettarini2012; Doss et al. Reference Doss, Komar, Cassak, Wilder, Eriksson and Drake2015), leading to persistence of some of the initial current sheets longer than in the zero shear, anti-parallel case. Thus, the inverse cascade by merging of magnetic islands is less advanced in the super-Alfvénic shear case, and the scales of the magnetic islands are consequently smaller in figure 1(c) than in figure 1(a). From other simulations (not shown) we note that the reconnection rate is significantly reduced only when super-Alfvénic shear velocities are imposed. This dependence of the reconnection rate on shear velocity therefore enables local modulation of the magnetic structure by introducing significant differences in the shear at each current sheet. For example, including a single super-Alfvénic shear layer with a set of zero or sub-Alfvénic shear layers leads to a magnetic structure in which a single, persistent current sheet is embedded within a background of magnetic islands.

Figure 2. (a) Magnetic power spectra for the AP run, at different times, with the power-law fit. (b) The time evolution of the spectral index $\unicode[STIX]{x1D6FC}$ obtained from the fits for the three runs. The horizontal line marks the scaling exponent value $8/3$ .

In order to characterize the apparent cascade seen in the magnetic structure in figure 1, we examine the time evolution of the spectrum of magnetic fluctuations. These spectra are shown for the anti-parallel case in figure 2(a). Each spectrum shown is the trace power spectrum of the full three-dimensional box, integrated over shells of constant $k$ . At early times, strong peaks in the power spectrum are associated with spatial harmonics of the current sheet separation scale. The evolution from $t=0$ to $t\unicode[STIX]{x1D6FA}_{i}\approx 150$ consists of in-filling of the spectral peaks, generating a power-law spectrum approximating $E_{B}(k)\sim k^{-8/3}$ . For inertial magnetohydrodynamic (MHD) turbulence, one might expect a Kolmogorov power-law exponent $-5/3$ . Here, the steeper slope may be a consequence of the inverse cascade associated with merging of magnetic islands.

In order to quantitatively follow the time evolution of the spectral properties, it is possible to use the value of the exponent $\unicode[STIX]{x1D6FC}$ of the isotropically integrated magnetic power spectra, obtained through a fit with the power law $E_{B}(k)\sim k^{-\unicode[STIX]{x1D6FC}}$ . Examples of such a fit over an appropriate range of wave vectors are shown in figure 2(a), at some different times in the simulation, highlighting the variation of the scaling exponent with time. The range of wave vectors is chosen to approximate the inertial range. Given the ion kinetic scales and the constraints of the size of the simulation domain, the inertial range is therefore limited to approximately half a decade in $k$ -space. We note that at high- $k$ , close to the grid scale, the apparent rise in power is an artefact of the isotropic integration of the full three-dimensional magnetic spectra, and is not associated with particle noise. The evolution of the spectral index for different simulation cases is shown in figure 2(b). The initial values are estimated by fitting the envelope of the spectral peaks observed at early times, and do not carry any significance in terms of nonlinear energy cascade. At $t\unicode[STIX]{x1D6FA}_{i}\simeq 25$ , the spectral peaks have already merged and the spectrum is flat. At following times, it starts to broaden and steepen, with the spectral index increasing from zero to approximately $\unicode[STIX]{x1D6FC}\simeq 3$ at $t\unicode[STIX]{x1D6FA}_{i}\simeq 150$ for the AP run, or to $\unicode[STIX]{x1D6FC}\simeq 2.7$ for the SH and GF runs. At later time, $t\unicode[STIX]{x1D6FA}_{i}\gtrsim 250$ , an approximate steady state is reached at $\unicode[STIX]{x1D6FC}\simeq 2.6$ for the AP and SH runs, while $\unicode[STIX]{x1D6FC}\simeq 2.1$ for the GF run.

Although the magnetic spectra have been isotropically integrated over the three-dimensional $k$ -space, the initial conditions may introduce sources of spectral anisotropy. The multiple current sheet system is itself anisotropic, however in the AP case fast growth of the tearing instability and associated fluctuations evolve the system towards an approximately isotropic state by $t\unicode[STIX]{x1D6FA}_{\text{i}}\gtrsim 25$ . However, the guide field is known to introduce anisotropy in the magnetic spectrum of an MCS system (Burgess et al. Reference Burgess, Gingell and Matteini2016). In the GF case, a relative reduction in the power in the parallel direction compared to the perpendicular direction, visible also in the magnetic structure in figure 1(b), leads to an effective reduction of the spectral index for the isotropically integrated spectrum. This difference is reflected in the reduced $\unicode[STIX]{x1D6FC}$ of the GF run compared with the AP and SH runs at late times, visible in figure 2(b). Velocity shear is also known to introduce spectral anisotropy (Wan et al. Reference Wan, Servidio, Oughton and Matthaeus2009). However, in the SH case presented in this study, for which we have a set of oppositely directed velocity shear layers with a range of magnitudes and net zero bulk velocity, the system nevertheless evolves towards an isotropic state over a time scale determined by the saturation of the tearing instability, approximately $t\unicode[STIX]{x1D6FA}_{i}\gtrsim 150$ . Hence, towards the end of the simulation for the AP and SH cases, spectral anisotropy is not expected to affect any statistics for which isotropy is assumed.

The description of turbulence requires a more detailed description than simply the spectral analysis. In fact, fully developed turbulence is characterized by intermittency (Frisch Reference Frisch1995), or lack of scale invariance of the field fluctuations, which in turn leads to non-uniform distribution of energy dissipation at small scales. For this reason, it is customary to estimate the scale dependence of the statistical properties of the two-point field increments $\unicode[STIX]{x0394}B(x,l)=B(x+l)-B(x)$ , where $l$ is the scale parameter. A complete description is provided by the probability distribution functions (PDFs) of the scale-dependent, standardized field increments. In turbulent flows these are usually Gaussian at large scale, and have increasing tails as the scale decreases, because of the formation of large amplitude structures that concentrate on smaller and smaller scales. This is the signature of intermittency, a process naturally and universally arising in fully developed Navier–Stokes or MHD turbulence (Frisch Reference Frisch1995). PDFs of longitudinal magnetic field increments of the $x$ component are shown in figure 3 at two times ( $t\unicode[STIX]{x1D6FA}_{i}\simeq 150$ and $t\unicode[STIX]{x1D6FA}_{i}\gtrsim 300$ ), and for three different scales. It is evident that PDFs at different scales collapse on a roughly similar (Gaussian) shape, shown with a dashed line, with weak deviation of the tails. At large scale, the PDF is slightly leptokurtic, i.e. tails are lower than Gaussian, while at smaller scale they progressively increase toward weakly hyperkurtic values. This is an indication of weak scale dependence, or intermittency, so that the field is nearly self-similar. Similar behaviour is observed for all times, for all field components, for all simulation cases. However, we note that for small scales at $t\unicode[STIX]{x1D6FA}_{i}\simeq 150$ , the tail regions of the PDF are more hyperkurtic than at $t\unicode[STIX]{x1D6FA}_{i}\simeq 300$ , as a remnant of the fluctuations generated by the tearing instability during the period of fast, nonlinear reconnection, which saturates at $t\sim 100$ .

Figure 3. Examples of the PDFs of the longitudinal field increments of the $x$ component, $\unicode[STIX]{x0394}B_{x}$ , at three different scales (different colours), estimated at two times in the simulation, for the anti-parallel case. The weak scale dependence is highlighted by the comparison with a reference Gaussian (dashed lines). A similar behaviour is observed at all times, for all the field components, and for all simulation cases (not shown).

A more quantitative estimate is given by means of the PDF structure functions, i.e. the high-order moments of the PDFs, $S_{q}(\unicode[STIX]{x0394}B)=\langle |\unicode[STIX]{x0394}b|^{q}\rangle$ , with the brackets indicating spatial averages (Frisch Reference Frisch1995). In turbulent flows, it is expected that the structure functions scale as power laws of the separation $l$ , $S_{q}(l)\sim l^{\unicode[STIX]{x1D701}_{q}}$ , in the inertial range. The behaviour of the scaling exponents with the moment order $q$ is indicative of the presence of intermittency. In particular, the linear dependence $\unicode[STIX]{x1D701}_{q}\propto q$ indicates self-similarity, while deviation from such linearity indicates intermittency. Due to the limited extension of the inertial range in the numerical simulations, extended self-similarity (ESS) (Benzi et al. Reference Benzi, Ciliberto, Tripiccione, Baudet, Massaioli and Succi1993) has been used, as customary. In this case, the structure functions $S_{q}$ are plotted as a function of the second-order moment $S_{2}$ . The corresponding extended-range fit with the power law $S_{q}(S_{2})\sim S_{2}^{\unicode[STIX]{x1D709}_{q}}$ gives the scaling exponents $\unicode[STIX]{x1D709}_{p}\propto \unicode[STIX]{x1D701}_{p}$ , which can therefore be used to measure deviation from self-similarity. Intermittency models can be used to fit the function $\unicode[STIX]{x1D709}_{q}$ , providing a quantitative estimate. Here we use the $p$ -model (Meneveau & Sreenivasan Reference Meneveau and Sreenivasan1987), whose prediction for the structure functions anomalous scaling is $\unicode[STIX]{x1D709}_{q}=1-\log _{2}\left[p^{hq}+(1-p)^{hq}\right]$ . Here $p\in [0.5,1]$ is the parameter that quantifies deviation from self-similarity, with $p=0.5$ indicating absence of intermittency and $p>0.5$ indicating increasing intermittency. Figure 4(a) shows one example of ESS structure functions for AP case at the end of the simulation, for the $B_{x}$ component. Power-law fits are also indicated. Figure 4(b) shows the order dependency of the scaling exponents $\unicode[STIX]{x1D709}_{q}$ , along with the $p$ -model fit, for the same case as in the figure 4(a), again taken the end of the simulation ( $t\unicode[STIX]{x1D6FA}_{i}=300$ ) and at earlier times. In figure 4(c), the time evolution of the parameter $p$ obtained from the fit is shown, for the $x$ component, and for all three cases. For the AP and GF cases the intermittency parameter reaches a stationary value $p\simeq 0.6$ at early times in the simulation, around $t\unicode[STIX]{x1D6FA}_{i}\simeq 50$ . Such value indicates a very weak degree of intermittency, in agreement with the observation of the scale-dependent PDFs. On the contrary, no intermittency ( $p\simeq 0.5$ ) is observed for the SH case.

Figure 4. (a) One example of structure functions (ESS is used), for the $B_{x}$ increments in the AP run, at the final time of the simulation. Power-law fits are indicated. (b) The anomalous scaling of the structure functions exponents $\unicode[STIX]{x1D709}_{q}$ at three different times, along with the $p$ -model fit, for the $B_{x}$ increments in the AP run. (c) The time evolution of the parameter $p$ for the $B_{x}$ components, for all three runs.

An alternative use of the structure functions to determine intermittency is to evaluate the kurtosis, or the normalized fourth-order moment $K(l)=S_{4}(l)/S_{2}(l)^{2}$ , which describes the scale-dependent flatness of the distribution functions. For Gaussian field increments PDFs (e.g. at large scale) the kurtosis has a definite value $K_{Gauss}=3$ . Again, for intermittent fields the kurtosis increases as the scale decreases, indicating the presence of enhanced tails of the distributions. Conversely, constant kurtosis is found for self-similar, non-intermittent fields. An example of scaling dependence of the kurtosis is shown in figure 5(a), where the weak increase towards small scales is visible. At large scales, a value $K<3$ is observed, confirming the presence of low-tailed PDFs. This is consistent with the observation of low-tailed PDFs for large scales at mid and late times in figure 3. These sub-Gaussian statistics are generally associated with large-scale, coherent equilibria. In our case, this indicates that a remnant of the initial current sheet configuration is still present. At smaller scales, after an approximately steady increase in the inertial range, the kurtosis reaches values slightly larger than the Gaussian value $K=3$ , indicating the presence of weakly higher tails, although the bulk of the PDF remains Gaussian. The differences among the three runs are recovered in the kurtosis values, again showing a less developed turbulence in the SH run. In the figure 5(b), the time evolution of the maximum value of the kurtosis over the whole range of scales is shown. This indicates the time of maximum deviation from Gaussian, which was confirmed as occurring at small scales. The deviation from $K=3$ is weakly variable after $t\unicode[STIX]{x1D6FA}_{i}=100$ , and reaches a steady value close to Gaussian at $t\unicode[STIX]{x1D6FA}_{i}=250$ .

Figure 5. (a) The scale dependence of the kurtosis for the $B_{x}$ increments for all three 3-D simulations, at the final time of the simulation. (b) The time evolution of the maximum of the scale-dependent kurtosis across all scales, again for the $B_{x}$ increments.

Figure 6. Time dependence of the scale-dependent kurtosis for the AP case (a) and SH case (b). Note that at early times in the AP case high kurtosis is associated with multiples of the sheet separation scale.

The full dependence of the kurtosis $K$ on both scale and time is shown for the AP and SH cases in figure 6. For the AP case, the regions of high kurtosis $K>6$ during the earliest period of significant reconnection, $t\unicode[STIX]{x1D6FA}_{i}\sim 20$ , are associated with multiples of the current sheet separation scale $7.5d_{i}$ . This is consistent with the growth of a tearing instability on each current sheet. In contrast, for the SH case we find that the breaking of the symmetry due to variable velocity shears between current sheets removes the peaks associated with the separation scale, and instead high kurtosis is found at the smallest length scales. This case is consistent with a turbulent state with high intermittency. However, as also shown in figures 4 and 5, this state is not maintained beyond $t\unicode[STIX]{x1D6FA}_{i}\sim 50$ .

Even in those cases where the turbulence is not fully developed, the properties of the flow can be described in terms of the chaotic nature of its fluctuations using the so-called cancellation analysis. This was first introduced in the framework of fluid turbulence and for the magnetohydrodynamic dynamo (Ott et al. Reference Ott, Du, Sreenivasan, Juneja and Suri1992). Subsequently, it has been used to describe the scaling properties of MHD structures in two-dimensional simulations (Sorriso-Valvo et al. Reference Sorriso-Valvo, Carbone, Noullez, Politano, Pouquet and Veltri2002; Graham, Mininni & Pouquet Reference Graham, Mininni and Pouquet2005; Martin et al. Reference Martin, De Vita, Sorriso-Valvo, Dmitruk, Nigro, Primavera and Carbone2013) and for solar photospheric active regions (Abramenko, Yurchishin & Carbone Reference Abramenko, Yurchishin and Carbone1998a ,Reference Abramenko, Yurchishin and Carbone b ; Sorriso-Valvo et al. Reference Sorriso-Valvo, Carbone, Veltri, Abramenko, Noullez, Politano, Pouquet and Yurchyshyn2004; Sorriso-Valvo et al. Reference Sorriso-Valvo, De Vita, Kazachenko, Krucker, Primavera, Servidio, Vecchio, Welsch, Fisher and Lepreti2015; De Vita et al. Reference De Vita, Vecchio, Sorriso-Valvo, Briand, Primavera, Servidio, Lepreti and Carbone2015). It was pointed out that the cancellation analysis could represent a phenomenological measure of the fractal dimension of the field structures (Sorriso-Valvo et al. Reference Sorriso-Valvo, Carbone, Noullez, Politano, Pouquet and Veltri2002). In analogy to the fractal singularity analysis of the measure of a positive defined signal, cancellation analysis is based on the study of the singularity of a signed field. In particular, if a suitably defined measure of a field changes sign on arbitrarily fine scale, then the measure is said to be sign singular (Ott et al. Reference Ott, Du, Sreenivasan, Juneja and Suri1992). The quantitative description of such a singularity gives information on sign changes, which is relevant in the presence of sign-defined, smooth coherent structures, such as the ones that arise in a turbulent cascade.

For a scalar field $f(\boldsymbol{r})$ with zero mean, defined on a $d$ -dimensional domain $Q(L)$ of size  $L$ , its signed measure is the normalized field across scale-dependent subsets $Q(l)\subset Q(L)$ of size  $l$ ,

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D707}(l)=\frac{\displaystyle \int _{Q(l)}\,\text{d}\boldsymbol{r}\;f(\boldsymbol{r})}{\displaystyle \int _{Q(L)}\,\text{d}\boldsymbol{r}\;|f(\boldsymbol{r})|}.\end{eqnarray}$$

Performing a coarse graining of the whole domain, the sign singularity of the measure can be quantitatively estimated through the scaling exponent of the cancellation function (or cancellation exponent) $\unicode[STIX]{x1D705}$ , defined as

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D712}(l)=\mathop{\sum }_{Q_{i}(l)}|\unicode[STIX]{x1D707}_{i}(l)|\sim l^{-\unicode[STIX]{x1D705}},\end{eqnarray}$$

the sum being over all disjoint subsets  $Q_{i}(l)$ covering the domain  $Q(L)$ . In a turbulent field, cancellations between positive and negative fluctuations occur when integrating over large-size subsets, thus providing a small contribution to the signed measure. Conversely, for a subset of the typical size of the structures, the increasing presence of sign-defined patches of field reduces the cancellations. The cancellation exponent can thus provide an effective measure of the way the field cancellations change through the scale. For example, for a smooth field the cancellation function does not depend on the scale, and $\unicode[STIX]{x1D705}=0$ ; on the other hand, a homogeneous field with random discontinuities (like a Brownian noise) has $\unicode[STIX]{x1D705}=d/2$ . Exponents between those two values indicate the presence of smooth structures embedded in random fluctuations. Moreover, their values can be related to the geometrical properties of structures, through the phenomenological relationship $\unicode[STIX]{x1D705}=(d-D)/2$ , where $D$ is the fractal dimension of the typical structures of the field (Sorriso-Valvo et al. Reference Sorriso-Valvo, Carbone, Noullez, Politano, Pouquet and Veltri2002). Conversely, exponents $\unicode[STIX]{x1D705}>d/2$ indicate cancellations that are more efficient than for a random field, suggesting the presence of pairs or groups of structures that efficiently cancel each other. In this case, the phenomenological relationship for the fractal dimension $D$ does not hold.

Figure 7. (a) The cancellation functions of the $j_{z}$ component, $\unicode[STIX]{x1D712}(l)$ , at three different times, for the AP run. Power-law fits are indicated in two ranges of scale above and below the break at $l_{\star }$ . (b) The time evolution of the fractal dimension $D_{small}$ (upper window, thin lines) and $D_{large}$ (lower window), for all three runs.

Cancellation analysis has recently been used to track the time evolution of turbulence in a two-dimensional numerical simulation of the hybrid Vlasov–Maxwell equations (De Vita et al. Reference De Vita, Sorriso-Valvo, Valentini, Servidio, Primavera, Carbone and Veltri2014), showing its ability to highlight the transition to turbulence. Moreover, could accurately describe the formation and the geometrical characteristics of the main structures generated by the nonlinear interaction, independent of the type of turbulence. Here we adopt a similar approach, and use the current components $j_{i}$ (with $i=x,y,z$ ) to estimate the cancellation function (3.2) at each time in the three simulations, obtaining a time-dependent estimate of the fractal dimension $D$ of the structures that form in the system.

The figure 7(a) shows three example of scaling of the cancellation function $\unicode[STIX]{x1D712}(l)$ for the component $j_{z}$ in the AP run, at three different times. Cancellation functions computed for the other current density components show similar features for all runs (not shown). At early times, an apparent break in the slope suggests that a double power-law scaling range may be identified, roughly covering the inertial range. Interestingly, the break between the two scaling laws is located near $l_{\star }\simeq 7d_{i}$ , i.e. close to the initial current sheets separation. This feature is confirmed by the analysis of a run with only four current sheets (not shown), for which the cancellation function has a break again near the initial current sheets separation scale $15d_{i}$ . For the four current sheet case, this break much more clearly persists at late times. Given the apparent break, the cancellation functions have been fitted with two power laws, in the small-scale and large-scale ranges, respectively below and above $l_{\star }$ . The exponents $\unicode[STIX]{x1D705}$ obtained from the fit have been used to estimate the corresponding fractal dimensions $D_{small}$ and $D_{large}$ , whose time evolution is shown in the figure 7(b) (thin and thick lines, respectively), for the three runs (different colours and line style). Observing the time evolution of $D$ , the following features can be highlighted in the large- and small-scale ranges.

(i) In the small-scale range (thin lines), the field is initially smooth ( $D_{small}=d=3$ ). Successively, the parameter $D_{small}$ gently decreases during the onset of turbulence, indicating the formation of small-scale structures. In the AP run, for $t\unicode[STIX]{x1D6FA}_{i}\gtrsim 100$ structures reach a steady, quasi-2-D geometry ( $D_{small}\simeq 2.2$ ), corresponding to the formation of current sheets. A similar behaviour was observed in the small-scale evolution of two-dimensional Vlasov–Maxwell turbulence (De Vita et al. Reference De Vita, Sorriso-Valvo, Valentini, Servidio, Primavera, Carbone and Veltri2014), where a 3-D-extrapolated dimension $D_{small}\simeq 1.5$ was measured. In the SH run, the presence of velocity shears accelerates the nonlinear transfer, favouring an early appearance of small-scale structures, which eventually reach a steady geometry around the same time as for the AP run. On the contrary, in the GF run saturation is reached at a much later time, $t\unicode[STIX]{x1D6FA}_{i}\simeq 250$ , showing that the guide field effectively slows down the small-scale nonlinear interactions. However, the fractal dimension settles at $D_{small}\simeq 1.8$ , so that slightly disrupted current sheets are formed.

(ii) In the large-scale range (thick lines), the early stage of the simulation is characterized by the presence of the initial current sheets. Because of the alternation in the current sign, cancellations are enhanced ( $\unicode[STIX]{x1D705}>d/2$ ) on those scales, and the fractal dimension cannot be estimated. For the AP run, at $t\unicode[STIX]{x1D6FA}_{i}\gtrsim 100$ nonlinear interactions and magnetic reconnection progressively disrupt the initial current layers, eventually forming quasi-2-D large-scale structures ( $D_{large}\simeq 1.9$ at $t\unicode[STIX]{x1D6FA}_{i}\simeq 300$ ). This is also visible in figure 7(a), where the difference between small-scale and large-scale scaling exponents decrease with time. This is again the signature of the transition to turbulence. However, at the final time of the simulation the dimension seems to be still increasing, and has not reached a clear saturation. Moreover, structures in the two ranges of scales reach different dimension, $D_{small}\neq D_{large}$ , so that a scale separation due to the initial conditions is still present. The GF and SH runs show similar features. However, in these cases the formation of structures is even slower than for the AP run, so that at the end of the simulation $D_{large}\simeq 1$ and the difference between large- and small-scale structures is more evident. This is in agreement with the slower evolution of magnetic islands in the presence of super-Alfvénic shears (as described previously) and the two-dimensional nature of the turbulence in the GF case.

The above observations suggest that, although the spectral indexes reach steady values around $t\unicode[STIX]{x1D6FA}_{i}\simeq 250$ , the turbulence has not fully developed yet, and the large-scale signature of the initial current sheets is still affecting the formation of the structures. This is even more evident when guide field or velocity shears slow down the nonlinear interactions and large-scale magnetic reconnection.

Figure 8. Distribution of fluctuations as a function of residual energy $\unicode[STIX]{x1D70E}_{r}$ and cross-helicity $\unicode[STIX]{x1D70E}_{c}$ , for the cases AP (a,d), SH (b,e) and SH2d (c,f). Distribution functions have been calculated for the periods $t\unicode[STIX]{x1D6FA}_{i}=40$ –60 in (a–c) and $t\unicode[STIX]{x1D6FA}_{i}=280$ –300 in (d–f), with lag $l=10d_{i}$ .

Finally, we can characterize the magnetic fluctuations present in the simulation by examining the distribution of the cross-helicity and residual energy of fluctuations as described by Bavassano & Bruno (Reference Bavassano and Bruno2006), and utilized more recently for solar wind turbulence by Wicks et al. (Reference Wicks, Roberts, Mallet, Schekochihin, Horbury and Chen2013). We calculate the increments in magnetic field $\boldsymbol{B}$ and bulk velocity $\boldsymbol{V}$ from 1-D trajectories through the simulation domain with constant $y$ and $z$ such that

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{b}(x,l)=\frac{\boldsymbol{B}(x)-\boldsymbol{B}(x+l)}{\sqrt{\unicode[STIX]{x1D707}_{0}m_{i}n_{0}(x,l)}} & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}\boldsymbol{v}(x,l)=\boldsymbol{V}(x)-\boldsymbol{V}(x+l), & \displaystyle\end{eqnarray}$$

where $l$ is the spatial lag. The local mean field and density are given by

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{B}_{0}(x,l)=\frac{1}{l}\int _{x^{\prime }=x}^{x^{\prime }=x+l}B(x^{\prime })\,\text{d}x^{\prime }, & \displaystyle\end{eqnarray}$$

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{n}_{0}(x,l)=\frac{1}{l}\int _{x^{\prime }=x}^{x^{\prime }=x+l}n(x^{\prime })\,\text{d}x^{\prime }, & \displaystyle\end{eqnarray}$$

and the Alfvénic part of the fluctuations are therefore given by:

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D739}\boldsymbol{b}_{\bot }(x,l)=\unicode[STIX]{x1D6FF}\boldsymbol{b}(x,l)\cdot (1-\hat{\boldsymbol{b}}_{0}(x,l)\hat{\boldsymbol{b}}_{0}(x,l)), & \displaystyle\end{eqnarray}$$

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D739}\boldsymbol{v}_{\bot }(x,l)=\unicode[STIX]{x1D6FF}\boldsymbol{v}(x,l)\cdot (1-\hat{\boldsymbol{b}}_{0}(x,l)\hat{\boldsymbol{b}}_{0}(x,l)), & \displaystyle\end{eqnarray}$$

where $\hat{\boldsymbol{b}_{0}}=\boldsymbol{B}_{0}/B_{0}$ . We can therefore calculate the scale-dependence, normalized cross-helicity and residual energy respectively as follows:

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{c}(x,l)=\frac{2\unicode[STIX]{x1D6FF}\boldsymbol{v}_{\bot }(x,l)\cdot \unicode[STIX]{x1D6FF}\boldsymbol{b}_{\bot }(x,l)}{|\unicode[STIX]{x1D6FF}\boldsymbol{v}_{\bot }(x,l)|^{2}+|\unicode[STIX]{x1D6FF}\boldsymbol{b}_{\bot }(x,l)|^{2}}, & \displaystyle\end{eqnarray}$$

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{r}(x,l)=\frac{\unicode[STIX]{x1D6FF}\boldsymbol{v}_{\bot }(x,l)|^{2}-|\unicode[STIX]{x1D6FF}\boldsymbol{b}_{\bot }(x,l)|^{2}}{|\unicode[STIX]{x1D6FF}\boldsymbol{v}_{\bot }(x,l)|^{2}+|\unicode[STIX]{x1D6FF}\boldsymbol{b}_{\bot }(x,l)|^{2}}. & \displaystyle\end{eqnarray}$$

The combined distribution functions of the cross-helicity and residual energy are shown for the anti-parallel and velocity shear simulations in figure 8. We have selected a scale $l=20$ , within the apparent inertial range of the magnetic spectra. The figure shows the combined distributions for every available time step during the period $t\unicode[STIX]{x1D6FA}_{i}=40-60$ , during the period of peak intermittency. For all simulations, we find that the power is distributed along the outer edge of the circle. Here $\unicode[STIX]{x1D70E}_{c}^{2}+\unicode[STIX]{x1D70E}_{r}^{2}=1$ , and therefore $|\unicode[STIX]{x1D6FF}\boldsymbol{v}_{\bot }||\unicode[STIX]{x1D6FF}\boldsymbol{b}_{\bot }|=|\unicode[STIX]{x1D6FF}\boldsymbol{v}_{\bot }\cdot \unicode[STIX]{x1D6FF}\boldsymbol{b}_{\bot }|$ . Hence, the magnetic field and velocity fluctuations are closely aligned. In the anti-parallel case, there is much greater power in the lower half of the circle, $\unicode[STIX]{x1D70E}_{r}<0$ and the fluctuations are therefore strongly magnetically dominated. Additionally, the fluctuations are found in the region $\unicode[STIX]{x1D70E}_{c}\approx 0$ , which implies that forward and backward propagating Elsasser fluctuations are balanced; unsurprising for a simulation which is initialized with zero cross-helicity. The fluctuations present in the guide field case at the same time have the same distribution as for the anti-parallel case in figure 8(a).

In the 3-D super-Alfvénic velocity shear case at early times, shown in figure 8(b), we find high power along the full circumference of the circle $\unicode[STIX]{x1D70E}_{c}^{2}+\unicode[STIX]{x1D70E}_{r}^{2}=1$ . This implies that simulation contains regions in which fluctuations are velocity dominated $\unicode[STIX]{x1D70E}_{r}>0$ , and strongly unbalanced $\unicode[STIX]{x1D70E}_{c}\approx \pm 1$ . This is consistent with the persistence of super-Alfvénic flows in the simulation, where velocity dominates the energy partition. As with the anti-parallel case, there is little power in the central region $\unicode[STIX]{x1D70E}_{c,r}\approx 0$ , for which magnetic field and velocity fluctuations are unaligned.

At the end of the simulation, $t\unicode[STIX]{x1D6FA}_{i}=300$ , the distribution of fluctuations in the 3-D SH case shown in figure 8(e) more closely resembles the anti-parallel case in figure 8(a,d), i.e. with very little power in the velocity-dominated $\unicode[STIX]{x1D70E}_{r}>0$ region compared to the magnetically dominated region. Hence, mixing of the super-Alfvénic bulk flows is more efficient than reconnection of magnetic flux. However, if the velocity shear simulation is performed in two dimensions, as for the SH2d case shown in figure 8(c,f), we find that the full ring of fluctuations seen in figure 8(b) persists even to $t\unicode[STIX]{x1D6FA}_{i}=300$ . The extra degree of freedom in the 3-D simulations therefore allows for more efficient mixing of the bulk flows, and structures which generate velocity-dominated fluctuations are less prevalent in the end state.

In the solar wind, analysis of turbulent magnetic fluctuations present in data from the Wind spacecraft have shown significant power in the region $\unicode[STIX]{x1D70E}_{c}\approx 1$ and $\unicode[STIX]{x1D70E}_{r}\lesssim 0$ (Wicks et al. Reference Wicks, Roberts, Mallet, Schekochihin, Horbury and Chen2013). Fluctuations in this region are strongly unbalanced, with very pure outward propagating Elsasser fluctuations, are Alfvénically equipartitioned ( $\unicode[STIX]{x1D6FF}\boldsymbol{b}_{\bot }\propto \unicode[STIX]{x1D6FF}\boldsymbol{v}_{\bot }$ ) and aligned. Hence, the character of the fluctuations generated by a system of reconnecting current sheets is significantly different from that observed in the near-Earth solar wind.