Three-dimensional simulations of sheared current sheets: transition to turbulence?

Abstract

Systems of multiple current sheets arise in various situations in natural plasmas, such as at the heliospheric current sheet in the solar wind and in the outer heliosphere in the heliosheath. Previous three-dimensional simulations have shown that such systems can develop turbulent-like fluctuations resulting from a forward and inverse cascade in wave vector space. We present a study of the transition to turbulence of such multiple current sheet systems, including the effects of adding a magnetic guide field and velocity shears across the current sheets. Three-dimensional hybrid simulations are performed of systems with eight narrow current sheets in a triply periodic geometry. We carry out a number of different analyses of the evolution of the fluctuations as the initially highly ordered state relaxes to one which resembles turbulence. Despite the evidence of a forward and inverse cascade in the fluctuation power spectra, we find that none of the simulated cases have evidence of intermittency after the initial period of fast reconnection associated with the ion tearing instability at the current sheets. Cancellation analysis confirms that the simulations have not evolved to a state which can be identified as fully developed turbulence. The addition of velocity shears across the current sheets slows the evolution in the properties of the fluctuations, but by the end of the simulation they are broadly similar. However, if the simulation is constrained to be two-dimensional, differences are found, indicating that fully three-dimensional simulations are important when studying the evolution of an ordered equilibrium towards a turbulent-like state.

1 Introduction

In this paper we consider the evolution of systems of multiple current sheets, using hybrid simulations with particle-in-cell (PIC) ions to correctly model key aspects of the physical processes involved, such as reconnection. Such systems, as well as being an interesting case to study in their own right, occur in several situations. For example, in the heliosphere the magnetic field in the solar wind has a polarity connecting towards or away from the Sun, and regions of outward and inward polarity are separated by a thin region known as the heliospheric current sheet (HCS) (Owens & Forsyth 2013). Misalignment of the Sun’s rotation and magnetic axes leads to a warped HCS which extends at low heliolatitudes throughout the heliosphere. Combined with the spiral pattern of the heliospheric magnetic field, this leads to multiple, embedded current sheets with field reversals in the outer heliosphere, as observed by the Voyager spacecraft up to the heliospheric termination shock (HTS) and beyond into the heliosheath (Burlaga & Ness 2011).

Recent interest in the effects of these current sheets, particularly multiple current sheet (MCS) systems, has been stimulated by discrepancies between simple theories of shock acceleration and Voyager observations of the anomalous cosmic ray (ACR) component (Stone et al. 2005, 2008). Using PIC simulations Drake et al. (2010) investigated particle energization from energy release in a system of closely spaced current sheets. A similar system was studied by Burgess, Gingell & Matteini (2016) using a three-dimensional hybrid simulation (particle ions with fluid electron response) including pick up ions as found in the outer heliosphere. A major result of this study was that the highly symmetric and ordered initial state of a MCS system evolves to a complex three-dimensional state, with a long-time end state which appears to become close to well-developed turbulence. A power-law form of the power spectrum of magnetic fluctuations develops due to both inverse and forward cascades and fluctuations are anisotropic for the case when a guide field is present, as found for solar wind turbulence.

In this paper we address specifically the transition to turbulence of an MCS system, in order to fully characterize the resultant fluctuations. The concept of universality is often invoked when considering plasma turbulence, and it is pertinent to consider the properties of apparently turbulent fluctuations generated by different initial conditions in simulations. For this reason the study of the MCS system is interesting in that it might represent an alternative method for initializing simulations of turbulence. One advantage of the MCS system is that velocity shear can be included in the initial conditions, so one has the possibility of producing turbulence with various levels of mean velocity shear fluctuations. Furthermore, the characterization of turbulent or nonlinear fluctuations and intermittency in complex and kinetic systems, where dissipation mechanisms are not fully understood, is important for the ongoing study of structure formation in nonlinear systems (Matthaeus et al. 2015). Comprehensive analysis of intermittency in one such complex system, here the ion kinetic MCS system, can therefore provide a valuable point of comparison for future study.

The paper is organized as follows: we first describe the simulations we have performed and their initialization, then the results of a number of different analysis techniques (power spectra, intermittency metrics, cancellation analysis and cross-helicity against residual energy) are presented, and finally we conclude with some discussion of the results of the simulations and implications for the use of multiple current sheet systems to study plasma turbulence.

2 Simulations

We investigate the evolution of turbulence from a multiple current sheet system using the three-dimensional hybrid simulation code HYPSI, as used previously to study the evolution of heliospheric ion-scale current sheets (Gingell, Burgess & Matteini 2015; Burgess et al. 2016). The hybrid method combines a fully kinetic, PIC treatment of the ion species with a charge-neutralizing, massless and adiabatic election fluid. Maxwell’s equations are solved in the low-frequency Darwin limit, with zero resistivity, using the current advance method and cyclic leapfrog (CAM-CL) algorithm described by Matthews (1994).

The three-dimensional simulations use a grid of $(N_{x},N_{y},N_{z})=(120,120,120)$ cells with resolution $\unicode[STIX]{x0394}x=0.5d_{i}$ , where $d_{i}=v_{A}/\unicode[STIX]{x1D6FA}_{i}$ is the ion inertial length. The size of the full simulation domain is therefore $(60d_{i})^{3}$ . All boundaries are chosen to be periodic. Distance and time are normalized to units of the ion inertial length $d_{i}$ and inverse ion gyrofrequency $t_{\unicode[STIX]{x1D6FA}}=\unicode[STIX]{x1D6FA}_{i}^{-1}$ , respectively; velocity is normalized to the Alfvén speed $v_{A}$ . In order to reduce noise as far as possible given the constraints of the available computational resources, the ion phase space has been sampled with 200 pseudo-particles per computational cell.

The simulations are initialized with eight current sheets parallel to the $y$ – $z$ plane, with spacing $7.5d_{i}$ . We use a ‘force-free’ equilibrium which consists of a rotation in the magnetic field with uniform density. Each current sheet therefore has the following magnetic structure:

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle B_{x}(x)=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle B_{y}(x)=\pm B_{0}\tanh (x/L), & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle B_{z}(x)=\pm B_{0}\text{sech}(x/L), & \displaystyle\end{eqnarray}$$

where the $B_{0}$ is the background magnetic field strength. The magnetic field reverses over a length scale $L=d_{i}$ . In contrast to the Harris current sheet equilibrium often used in other studies (e.g. Drake et al. 2010), both the magnetic field strength and density are uniform over the current sheet. Ions are initialized with a Maxwellian distribution function, with plasma beta $\unicode[STIX]{x1D6FD}_{i}=0.5$ . For the plasma conditions chosen, we have verified that evolutionary differences which arise from the use of non-Maxwellian distribution functions required at kinetic scales (see Neukirch, Wilson & Harrison 2009; Wilson & Neukirch 2011) are negligible.

In addition to the initial conditions described above, referred to hereafter as the ‘anti-parallel’ or AP case, two variations are also described in this study. In the guide field (GF) case a uniform $B_{z}$ component is added, such that initially $B_{z}=B_{y}$ and the background field strength $B_{g}^{2}=B_{y}^{2}+B_{z}^{2}=B_{0}^{2}$ . This corresponds to a rotation of the background magnetic field without a change in the background field strength. In the super-Alfvénic velocity shear case (SH), a variable bulk velocity $\boldsymbol{u}_{SH}=(0,u_{y}(x),0)$ is included between each current sheet, such that each of the current sheets is subject to a different velocity shear. The initial bulk velocities between the current sheets in order of increasing $x$ coordinate are as follows: $u_{y}/v_{A}=[6,-8,2,1,-2,4,-3]$ . Each velocity shear layer is initialized with the same width as the magnetic shear layer $L$ , such that the initial bulk velocity is given by:

(2.4) $$\begin{eqnarray}u_{y,total}(x)=\mathop{\sum }_{i=1}^{7}\frac{1}{2}u_{y,i}\left(\tanh \left(\frac{x-x_{c,i}}{L}\right)-\tanh \left(\frac{x-x_{c,i+1}}{L}\right)\right),\end{eqnarray}$$

where $x_{c,i}$ is the centre of a given current sheet. Note that the net bulk velocity in the simulation is zero. In this case, each current sheet is subject to a super-Alfvénic velocity shear in addition to the magnetic shear. We also discuss a ‘2.5-dimensional’ (2.5-D) version of these initial conditions (SH2d), for which all three components of the fields and moments do not vary in the $z$ -direction, e.g. $B_{xyz}(x,y,t)$ .

The full set of simulations discussed in this paper, including relevant abbreviations, is given in table 1.

Figure 1. Line integral convolution visualization of the magnetic structure at $t\unicode[STIX]{x1D6FA}_{i}=300$ for (a) anti-parallel, (b) guide field and (c) super-Alfvénic shear simulations. (The plots are produced in MATLAB with routines from Toolbox image.)

Table 1. Summary of simulations discussed in this paper.

3 Results

The evolution of the simulations in time can be summarized as follows (Burgess et al. 2016): (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. (2016). 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. (2015).

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 1990; Cassak & Otto 2011; Landi & Bettarini 2012; Doss et al. 2015), 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. 2016). 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. 2009). 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 1995), 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 1995). 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 1995). 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. 1993) 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 1987), 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. 1992). Subsequently, it has been used to describe the scaling properties of MHD structures in two-dimensional simulations (Sorriso-Valvo et al. 2002; Graham, Mininni & Pouquet 2005; Martin et al. 2013) and for solar photospheric active regions (Abramenko, Yurchishin & Carbone 1998a , b ; Sorriso-Valvo et al. 2004; Sorriso-Valvo et al. 2015; De Vita et al. 2015). 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. 2002). 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. 1992). 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. 2002). 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. 2014), 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. 2014), 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 (2006), and utilized more recently for solar wind turbulence by Wicks et al. (2013). 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. 2013). 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.

4 Conclusions

In this study we have presented, using several methods of analysis, the first comprehensive analysis of an apparently turbulent or disordered system which has evolved by the relaxation of an ordered system with magnetic and velocity shear. We have presented results for a set of three-dimensional hybrid simulations of a multiple current sheet system, relevant to Kelvin–Helmholtz driven turbulence in the flanks of planetary magnetospheres (Stawarz et al. 2016), and the bunching of the heliospheric current sheet in the outer heliosphere (Burlaga & Ness 2011). Previous studies (Burgess et al. 2016) have shown that the MCS system evolves from an ordered equilibrium towards a disordered state of interacting magnetic islands (see figure 1). The evolution in all the simulations presented here is driven principally by the ion kinetic tearing instability, leading to reconnection across several current sheets as magnetic islands grow to scales larger than the sheet separation. We find that the rate of reconnection is reduced in the guide field and super-Alfvénic shear cases (compared to the anti-parallel AP case) due to stabilization of the tearing instability.

Indeed, the trace power spectrum of magnetic fluctuations, shown for the anti-parallel case in figure 2, demonstrates that power is distributed from the initial sheet separation scales in an apparent cascade to both smaller and larger scales. The inverse cascade is consequence of the rapid merging of magnetic islands during the nonlinear phase of the tearing instability, and slower coalescence of magnetic islands between adjacent current sheets in the later phase of the evolution. This second phase, in particular, is similar to the inverse energy cascade observed for 2-D fluid turbulence (Biskamp 2003). For all the simulations discussed here, the spectrum saturates by $t\unicode[STIX]{x1D6FA}_{i}\sim 150$ . The power-law exponent $\unicode[STIX]{x1D6FC}=2.6$ in the AP case, falling between typical solar wind values of $\unicode[STIX]{x1D6FC}=5/3$ in the fluid range and $\unicode[STIX]{x1D6FC}\sim 2.8$ at ion kinetic scales.

However, despite the apparent forward and inverse cascade visible in the evolution of the magnetic topology and spectra, none of the simulations discussed in this study have evidence of intermittency after the period of fast reconnection, $t\unicode[STIX]{x1D6FA}_{i}\gtrsim 100$ . For example, at late time the kurtosis of the PDF of magnetic fluctuations for all simulations $K(l)\lesssim 3$ over all scales $l$ , indicating a Gaussian or slightly sub-Gaussian PDF. The lack of intermittency is also supported by analysis of the scaling of the structure functions in figure 4. Due to this lack of intermittency, we can conclude, that the chaotic state observed at late times is not strictly a turbulent state.

During the period of active reconnection, $t\unicode[STIX]{x1D6FA}_{i}<150$ , Burgess et al. (2016) found that the tearing instability is able to generate regions of local temperature anisotropy, in agreement with recent observations (Hietala et al. 2015). These regions of local temperature anisotropy may lead to the localized growth of micro-instabilities (Gingell et al. 2015), and therefore act as a source of intermittency. Although the results presented in this paper demonstrate that the system is intermittent during the early reconnecting phase, we can conclude that any intermittency associated with ion kinetic micro-instabilities is transient, i.e. does not persist once the source of significant temperature anisotropy has become inactive.

Cancellation analysis reveals a clear break in the behaviour of the system at scales above and below the initial sheet separation scale, particularly in the SH and GF cases. In the GF case, the fractal dimension for scales below the separation scale $D_{small}\simeq 1.8$ , indicative of a slightly disrupted current sheets. Above separation scales $D_{large}\simeq 1$ , indicative of quasi-1-D structures (i.e. flux ropes within the magnetic islands). That $D_{small}\neq D_{large}$ in all cases at $t\unicode[STIX]{x1D6FA}_{i}=300$ suggests that signatures of the ordered initial conditions are still present in the simulation, and hence we have not yet reached a saturated state. This leaves the possibility that a truly turbulent (i.e. intermittent) state may still develop at much later times. However, the significant computational resources necessary to run a hybrid simulation over these time scales leads to the conclusion that the MCS configuration may not be as useful for generation of turbulence as those based upon decay of a superposition of wave modes. However, the system remains important to studies of turbulence arising from an initially ordered state, or those focused on the comparison of these systems to decaying turbulence.

The set of simulations presented in this paper also allows a direct comparison to be made of the turbulent or chaotic state in a system which includes super-Alfvénic velocity shear to the zero shear case. In the SH case, the later saturation of the magnetic spectrum and the slower growth of magnetic islands are consistent with a reduction in the reconnection rate. However, during the period of fast reconnection $t\unicode[STIX]{x1D6FA}_{i}\sim 40-50$ , the anomalous scaling of the structure functions, the higher kurtosis of magnetic fluctuations at small scales, and the lower fractal dimension at small scales all demonstrate that the super-Alfvénic shear drives fast, nonlinear structure formation below the sheet separation scale. During this period, we also find that the PDF of Alfvénic fluctuations in the SH case shows significant power in the velocity-dominated and strongly unbalanced regions of $(\unicode[STIX]{x1D70E}_{c},\unicode[STIX]{x1D70E}_{r})$ space (see figure 8). However, we can conclude that this clear difference between the character of the turbulence in the SH and AP cases is not persistent; by $t\unicode[STIX]{x1D6FA}_{i}=300$ the analyses presented in this paper show similar results for the AP and SH cases, or a trend in that direction. We further note that the differences observed in cases which include velocity shear layers are only significant if the shear is super-Alfvénic. Variations upon the SH initial conditions with reduced intra-sheet bulk velocities show neither reducing in the reconnection rate, nor velocity-dominated Alfvénic fluctuations in $(\unicode[STIX]{x1D70E}_{c},\unicode[STIX]{x1D70E}_{r})$ space.

In combination, these results are a clear demonstration of both the current limitations of simulations of turbulence, and the caution required in their analysis. We reiterate that a power law in the magnetic spectrum is not a sufficient measure of turbulence, and particular attention must be paid to quantitative measures of intermittency. The difference between the power-law exponent $\unicode[STIX]{x1D6FC}$ of the magnetic spectra measured here and previous measurements in the solar wind observation or more traditional turbulence simulations may be due to the limitation of the box size and grid resolution of our hybrid simulations. Likewise, we may find a more intermittent state over much longer time scales than the simulations discussed here. Finally, clear differences between otherwise identical 2-D and 3-D runs, such as those shown in figure 8, underline the importance of the application of fully three-dimensional simulations when studying the transition to turbulence from an ordered equilibrium.