2.1 Development and saturation of the turbulent magnetic power spectrum

As modes driven at the scale of the simulation domain rise in amplitude, they begin to interact nonlinearly, leading to a cascade of their energy to smaller scales. After sufficient time has passed, the turbulent energy dissipated by kinetic mechanisms rises to match the energy input into the turbulence by the driving, causing the turbulent energy spectrum to reach a statistically steady state. In figure 2, we display the one-dimensional perpendicular magnetic power spectrum, $E_{\bot }(k_{\bot })=\int \text{d}k_{z}\int \text{d}\unicode[STIX]{x1D703}k_{\bot }|\unicode[STIX]{x1D6FF}\boldsymbol{B}_{\bot }(\boldsymbol{k})|^{2}/8\unicode[STIX]{x03C0}$ , where the wavevector components are expressed in cylindrical coordinates $(k_{\bot },\unicode[STIX]{x1D703},k_{z})$ . The spectra shown in figure 2 range from over-strong turbulence with $\unicode[STIX]{x1D712}>1$ (Run 1, blue), critically balanced, strong turbulence with $\unicode[STIX]{x1D712}\simeq 1$ (Run 4, red) to weak turbulence with $\unicode[STIX]{x1D712}<1$ (Run 6, yellow). The $E_{\bot }(k_{\bot })$ spectra plotted here are obtained through averaging over the full time interval from the saturation of the turbulent cascade (at time $t_{2}$ , as defined below) until the end of the simulation at $t_{f}$ . These time averaged $E_{\bot }(k_{\bot })$ spectra are fit to the function $C(k_{\bot }\unicode[STIX]{x1D70C}_{i})^{\unicode[STIX]{x1D6FC}}\exp (-k_{\bot }\unicode[STIX]{x1D70C}_{e})$ (Alexandrova et al. Reference Alexandrova, Lacombe, Mangeney, Grappin and Maksimovic2012; TenBarge et al. Reference TenBarge, Howes and Dorland2013), where $C$ and $\unicode[STIX]{x1D6FC}$ are fitting parameters and $\unicode[STIX]{x1D70C}_{e}=v_{th_{e}}/\unicode[STIX]{x1D6FA}_{e}$ is the electron Larmor gyroradius.

Figure 2. The time-averaged one-dimensional perpendicular magnetic energy spectrum $E_{\bot }(k_{\bot })$ computed for three kinetic Alfvén wave turbulence simulations: Run 1 (blue) with $\overline{\unicode[STIX]{x1D712}}=4.9$ , Run 4 (red) with $\overline{\unicode[STIX]{x1D712}}=1.3$ and Run 6 (yellow) with $\overline{\unicode[STIX]{x1D712}}=0.44$ . Fits to $C(k_{\bot }\unicode[STIX]{x1D70C}_{i})^{\unicode[STIX]{x1D6FC}}\exp (-k_{\bot }\unicode[STIX]{x1D70C}_{e})$ are overplotted.

Figure 3. The time evolution of (a) the perpendicular magnetic energy spectrum $E_{\bot }(k_{\bot }\unicode[STIX]{x1D70C}_{i}=20.6)$ , (b) the spectral index $\unicode[STIX]{x1D6FC}$ from the fit to the function $A(k_{\bot }\unicode[STIX]{x1D70C}_{i})^{\unicode[STIX]{x1D6FC}}\exp (-k_{\bot }\unicode[STIX]{x1D70C}_{e})$ , (c) the normalized total heating rate, ${\hat{H}}=H/P_{in}$ , the normalized rate of change of energy of turbulent fluctuations ${\hat{W}}={\dot{W}}/P_{in}$ ( ${\dot{W}}=-\text{d}W/\text{d}t$ ) and the sum ${\hat{H}}+{\hat{W}}$ , where $P_{in}$ is the power input into the plasma by the Langevin antenna driving. Columns from left to right are for Run 1 ( $\overline{\unicode[STIX]{x1D712}}=4.9$ ), Run 2 ( $\overline{\unicode[STIX]{x1D712}}=3.4$ ) and Run 3 ( $\overline{\unicode[STIX]{x1D712}}=1.6$ ).

We evaluate the time scale of the saturation of the turbulent magnetic energy spectrum by following the evolution of the amplitude of the spectrum at $k_{\bot }\unicode[STIX]{x1D70C}_{i}=20.6$ , a position in the middle of the dynamic range of the simulation, as a function of the normalized time $t/t_{l}$ . Figure 3(a) displays the time evolution of the one-dimensional perpendicular magnetic spectrum $E_{\bot }(k_{\bot }\unicode[STIX]{x1D70C}_{i}=20.6)$ for Run 1 ( $\overline{\unicode[STIX]{x1D712}}=4.9$ ), Run 2 ( $\overline{\unicode[STIX]{x1D712}}=3.4$ ) and Run 3 ( $\overline{\unicode[STIX]{x1D712}}=1.6$ ), left to right. These runs correspond to a range of turbulent amplitudes, from over-strong turbulence (above critical balance) with $\unicode[STIX]{x1D712}>1$ to approximately critically balanced, strong turbulence with $\unicode[STIX]{x1D712}\simeq 1$ . These plots demonstrate the existence of two time scales that characterize the evolution of the turbulent magnetic field when driving by counterpropagating kinetic Alfvén waves from uniform magnetic field conditions. First, the energy in Fourier modes with $k_{\bot }\unicode[STIX]{x1D70C}_{i}=20.6$ rises as a steep power law in time $\propto t^{b_{1}}$ for times $t<t_{1}$ , denoted by the first vertical blue line in figure 3. Next, the rate of increase of energy in mode $k_{\bot }\unicode[STIX]{x1D70C}_{i}=20.6$ decreases to less steep approximate power law $\propto t^{b_{2}}$ , with $b_{2}<b_{1}$ , over the time interval $t_{1}<t<t_{2}$ . For times $t>t_{2}$ , the energy in modes with $k_{\bot }\unicode[STIX]{x1D70C}_{i}=20.6$ remains approximately constant, indicating that the region of the turbulent cascade has reached a statistically steady state.

We interpret these two time scales in the following way. At times $t<t_{1}$ , the energy of the modes with $k_{\bot }\unicode[STIX]{x1D70C}_{i}=20.6$ is dominated by the energy input from nonlinear interactions among lower wavenumber modes, leading to simple power-law increase in amplitude. At $t=t_{1}$ , the amplitude of the modes with $k_{\bot }\unicode[STIX]{x1D70C}_{i}=20.6$ has risen to a sufficient level that these modes are able to interact nonlinearly with counterpropagating fluctuations to transfer their energy to yet higher wavenumber modes; therefore, the rate of energy increase for the $k_{\bot }\unicode[STIX]{x1D70C}_{i}=20.6$ modes diminishes for $t_{1}<t<t_{2}$ . The second time scale $t_{2}$ signifies the time when this nonlinear energy loss to higher wavenumber modes increases sufficiently to balance the nonlinear energy input from lower wavenumber modes, effectively marking the saturation of a fully developed turbulent cascade.

To further illuminate the physical significance of these two time scales, we present two more measure of the turbulent cascade as a function of time. First, in figure 3(b), we present the time evolution of spectral index $\unicode[STIX]{x1D6FC}$ of the one-dimensional perpendicular magnetic spectrum when fit to the functional form $C(k_{\bot }\unicode[STIX]{x1D70C}_{i})^{\unicode[STIX]{x1D6FC}}\exp (-k_{\bot }\unicode[STIX]{x1D70C}_{e})$ with free adjustable parameters $\unicode[STIX]{x1D6FC}$ for the spectral exponent and $C$ as an arbitrary amplitude.

Figure 3(c), we plot diagnostics of the global power balance in each of the simulations. Further details of the power balance diagnostics in AstroGK are described in the Appendix of TenBarge et al. (Reference TenBarge, Howes and Dorland2013), so here we present only a brief description. Energy conservation dictates that the rate of change total fluctuating power (Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009) in the turbulent plasma in AstroGK, ${\dot{W}}=-\text{d}W/\text{d}t$ , is given by ${\dot{W}}=P_{in}-H$ , where $P_{in}$ is the power input into the plasma by the Langevin antenna driving and $H$ is the energy lost from the fluctuations due to collisional (and numerical) dissipation. Note that the heating $H\geqslant 0$ because collisional dissipation can only heat the plasma, removing energy irreversibly from turbulent field and plasma fluctuations and leading to thermal heating of the plasma species (which is not contained in $W$ ). In a statistically steady state, one has a time averaged $\langle {\dot{W}}\rangle \simeq 0$ , so therefore the time-averaged antenna power input into the plasma is $\langle P_{in}\rangle \simeq \langle H\rangle \geqslant 0$ , meaning that, on average, all of the input power is collisionally dissipated by the plasma. For plotting convenience, we choose to normalize the power balance by dividing by $P_{in}$ , leading to a normalized expression for instantaneous energy conservation, ${\hat{W}}+{\hat{H}}$ , where ${\hat{W}}={\dot{W}}/P_{in}$ and ${\hat{H}}=H/P_{in}$ . Energy is not conserved exactly by the AstroGK code due to a small amount of numerical dissipation, but the energy is typically conserved within a few per cent (TenBarge et al. Reference TenBarge, Howes and Dorland2013), as seen by the small departures of this sum from 1 in figure 3(c). Note that although the slope of the magnetic energy spectrum (figure 3 b) has reached a steady value at $t_{2}$ , the normalized heating rate ${\hat{H}}$ (figure 3 c) does continue to increase after this time. Although the definition of $t_{2}$ as the saturation time of the turbulent cascade agrees qualitatively with the expected scaling with turbulent amplitude, the slight increase of energy (manifesting an increasing heating rate) after $t_{2}$ suggests that the turbulence does experience a slow evolution beyond $t_{2}$ . Investigating how this slow evolution alters the properties of the turbulence will be explored in future works.

Figure 3(b), it is clear that second time scale $t_{2}$ is the saturation time for the development of the turbulent energy spectrum. This is the time required for the turbulent energy input at large scales to be transferred nonlinearly down to sufficiently scales that the turbulent energy is dissipated by kinetic mechanisms, irreversibly converted into plasma heat. It is interesting to note that, when normalized to the linear period $t_{l}$ of the kinetic Alfvén waves that are driven at the domain scale, we generally find $t_{2}/t_{l}<1$ , in contrast to the typical time scale assumed for strong turbulence that the saturation time $t_{2}\sim t_{l}$ . The likely explanation is the dispersive nature of kinetic Alfvén waves is that in the limit $k_{\bot }\unicode[STIX]{x1D70C}_{i}\gg 1$ , the kinetic Alfvén wave frequency scales as $\unicode[STIX]{x1D714}\propto k_{\bot }\unicode[STIX]{x1D70C}_{i}k_{\Vert }v_{A}$ , meaning that smaller-scale kinetic Alfvén waves have increasingly short period when normalized to the parallel wavenumber. Therefore, the full turbulent cascade develops more rapidly than the usual expectation from the non-dispersive inertial range. This result also demonstrates that driving turbulence simulations by counterpropagating kinetic Alfvén waves is a very efficient means of generating a turbulent spectrum.

Looking at the third row of figure 3, we see that the point at which the normalized rate of energy increase in the plasma due to driving ${\dot{W}}/P_{in}$ is overtaken by the rate of heating due to dissipation $H/P_{in}$ falls within the time interval $t_{1}<t<t_{2}$ . By the time $t_{2}$ is reached, the rate of heating $H/P_{in}$ dominates over the rate of energy increase in the turbulent plasma. At $t>t_{2}$ , eventually the time-averaged change of the turbulent plasma energy $\langle {\dot{W}}/P_{in}\rangle \rightarrow 0$ , leaving $H/P_{in}\simeq 1$ , indicating that all of the input power is collisionally dissipated, heating the plasma irreversibly.

The same qualitative behaviour of the development and saturation of the turbulent cascade is also found for cases of weak turbulence, $\unicode[STIX]{x1D712}<1$ , as shown in figure 4. Here, we plot the same rows of information as in figure 3 for Run 5 ( $\overline{\unicode[STIX]{x1D712}}=0.83$ ) and Run 6 ( $\overline{\unicode[STIX]{x1D712}}=0.44$ ). Again, two time scales are observed, a first time scale $t_{1}$ where the power-law increase of the $k_{\bot }\unicode[STIX]{x1D70C}_{i}=20.6$ decreases to a slower rate and the saturation time scale $t_{2}$ where spectral index of the perpendicular magnetic energy spectrum reaches a statistically steady value and the rate of heating $H/P_{in}$ dominates over the rate of energy increase in the turbulent plasma ${\dot{W}}/P_{in}$ .

For all of the Runs 1–7, the times $t_{1}$ and $t_{2}$ , normalized to the linear kinetic Alfvén wave period at the domain scale (effectively the outer time scale of the turbulent cascade), are presented in table 1. For the cases of over-strong turbulence $\unicode[STIX]{x1D712}>1$ , the saturation time $t_{2}$ is relatively constant, with a value between one-quarter and one-third of the domain-scale kinetic Alfvén wave period. As the turbulence weakens in the regime $\unicode[STIX]{x1D712}<1$ , the saturation time $t_{2}$ increases monotonically as the strength of the turbulence $\unicode[STIX]{x1D712}$ decreases, as expected theoretically (Sridhar & Goldreich Reference Sridhar and Goldreich1994; Howes, TenBarge & Dorland Reference Howes, TenBarge and Dorland2011a ).

2.2 Stochasticity of magnetic field lines

Here we describe the development of the stochasticity of the magnetic field lines in gyrokinetic turbulence as a function of the time-averaged nonlinearity parameter $\overline{\unicode[STIX]{x1D712}}$ . To study the magnetic field stochasticity, we employ Poincaré recurrence plots derived from snapshots of the magnetic field at different chosen times in the evolution of the turbulence simulations listed in table 1.

Figure 4. The time evolution of: (a) the perpendicular magnetic energy spectrum $E_{\bot }$ at $k_{\bot }\unicode[STIX]{x1D70C}_{i}=20$ , $E_{\bot }$ (a,b) the spectral index $\unicode[STIX]{x1D6FC}$ that corresponds to function $C(k_{\bot }\unicode[STIX]{x1D70C}_{i})^{\unicode[STIX]{x1D6FC}}\exp (-k_{\bot }\unicode[STIX]{x1D70C}_{e})$ fit to $E_{\bot }(k_{\bot })$ , (c) the normalized total heating rate, ${\hat{H}}=H/P_{in}$ , the normalized rate of change of energy of turbulent fluctuations ${\hat{W}}={\dot{W}}/P_{in}$ ( ${\dot{W}}=-\text{d}W/\text{d}t$ ) and the sum ${\hat{H}}+{\hat{W}}$ , where $P_{in}$ is the power input into the plasma by the Langevin antenna driving. Columns from left to right are for Run 5 ( $\overline{\unicode[STIX]{x1D712}}=0.83$ ) and Run 6 ( $\overline{\unicode[STIX]{x1D712}}=0.44$ ).

To construct the Poincaré plot, we begin with the magnetic field $\boldsymbol{B}(\boldsymbol{x})$ at some time $t$ . On the perpendicular plane at one end of the simulation domain at $z=0$ , we specify a sparse pattern of points with the colour of each point creating a bullseye pattern, as shown in figure 5. The magnetic field line passing through each point is traced through the domain to the far end of the simulation domain at $z=L_{\Vert }$ , and a point is plotted there, with colour matching that of the original field line position. That field line is periodically wrapped to $z=0$ , and the process is continued, with a coloured point plotted at each crossing at $z=L_{\Vert }$ . We trace through the box $20$ times for each field line, thereby plotting 20 coloured points on the plane for each field line. A fourth-order Runge–Kutta method with adaptive step size is used to trace each field line by integrating the ordinary differential equation $\text{d}\boldsymbol{r}/\text{d}l=\hat{\boldsymbol{b}}(\boldsymbol{x})$ , where $\hat{\boldsymbol{b}}(\boldsymbol{x})=\boldsymbol{B}(\boldsymbol{x})/|\boldsymbol{B}(\boldsymbol{x})|$ . If the field line passes through the boundaries of the simulation domain in the $x$ or $y$ directions, it is periodically wrapped to the opposite boundary.

Figure 5. Coloured points represent the coordinates of the starting points in the $z=0$ plane for tracing field lines to generate Poincaré recurrence plots.

Figure 6. Poincaré recurrence plots for over-strong and strong turbulence in Runs 1–4 at the times indicated. The corresponding time-averaged nonlinearity parameters are $\overline{\unicode[STIX]{x1D712}}=4.9$ (a), $\overline{\unicode[STIX]{x1D712}}=3.4$ (b), $\overline{\unicode[STIX]{x1D712}}=1.6$ (c) and $\overline{\unicode[STIX]{x1D712}}=1.3$ (d).

Figure 7. Poincaré recurrence plots for weak turbulence in Runs 5–7 at the times indicated. The corresponding time-averaged nonlinearity parameters are $\overline{\unicode[STIX]{x1D712}}=0.83$ (a), $\overline{\unicode[STIX]{x1D712}}=0.44$ (b) and $\overline{\unicode[STIX]{x1D712}}=0.14$ (c).

The Poincaré plots for Runs 1–4 are presented in figure 6 and for Runs 5–7 in figure 7. Stochasticity of the magnetic field is qualitatively assessed by the appearance of the Poincaré recurrence plots: a regular pattern of colours indicates a non-stochastic field, whereas a scrambling of the different colours indicates a stochastic field. In some cases, regions of the domain display a stochastic field while other regions appear more well ordered.

From the Poincaré plots presented in figures 6 and 7, we make the following observations. Early in the evolution at $t\ll t_{2}$ , when the turbulence is still not well developed and little energy has cascaded to small scales to generate a broadband turbulent spectrum, the Poincaré plots for all cases show a well-ordered pattern with no indication of stochasticity. After the turbulent cascade has saturated, at times $t>t_{2}$ , the over-strong and strong turbulence cases in Runs 1–4 with $\overline{\unicode[STIX]{x1D712}}\gtrsim 1$ all show that the magnetic field lines have become completely stochastic throughout the domain. For the weaker turbulence cases of Runs 5–7 with $\overline{\unicode[STIX]{x1D712}}<1$ , a remnant of order persists in the Poincaré plots even at $t\sim t_{2}$ when the turbulent cascade has saturated. At sufficiently long times $t\gg t_{2}$ , one sees the eventual development of a fully stochastic magnetic field throughout the domain for Run 5 with $\overline{\unicode[STIX]{x1D712}}=0.83$ and Run 6 with $\overline{\unicode[STIX]{x1D712}}=0.44$ . In contrast, the very weakest turbulence case with $\overline{\unicode[STIX]{x1D712}}=0.14$ in Run 7 shows persistent regions of order, even at long times $t>10t_{2}$ . Therefore, it seems that there is indeed a minimum amplitude threshold $\unicode[STIX]{x1D712}_{thresh}$ for the development of a completely stochastic magnetic field within the saturation time interval $t_{2}<t<13t_{2}$ , with order persisting for cases with nonlinear parameters $\unicode[STIX]{x1D712}\lesssim \unicode[STIX]{x1D712}_{thresh}\sim 0.1$ .

2.3 Separation of the magnetic field lines

The stochasticity is not the only feature that characterizes the magnetic field line wander in plasma turbulence. Two adjacent magnetic field lines may also separate from one another as we follow along one of the field lines. This magnetic field line separation may significantly impact the propagation of energetic particles in a turbulent plasma. Here we investigate the separation of stochastic magnetic field lines caused by plasma turbulence as a function of the nonlinearity parameter $\unicode[STIX]{x1D712}$ using gyrokinetic simulations (listed in table 1). We consider only Runs 1–6, excluding Run 7 as the magnetic field lines are not fully stochastic at the saturation phase of that simulation run.

Figure 8. Sketch describing how field lines pass through an extended domain built of six simulation boxes, showing the tracing of several wandering magnetic field lines up to distance $z=2L_{\Vert }$ .

To explore the separation of magnetic field lines, we assume an extended domain built up from the periodic extension of the single gyrokinetic simulation domain of size $L_{\bot }\times L_{\bot }\times L_{\Vert }$ . The magnetic field is initially defined by simulation over the ranges $0\leqslant z\leqslant L_{\Vert }$ and $0\leqslant x,y\leqslant L_{\bot }$ . To track the magnetic field lines in the extended spatial domain beyond these limits, we begin at a point in the $z=0$ plane and trace along the magnetic field line passing through that point (using the fourth-order Runge–Kutta method with adaptive step size) in the $+\hat{\boldsymbol{z}}$ direction. When the field line reaches the end of the single simulation domain at $z=L_{\Vert }$ , it enters an identical domain periodically extended to $z>L_{\Vert }$ , as depicted in figure 8. If the magnetic field line exits through the side of the domain at $x,y<0$ or $x,y>L_{\bot }$ , it is not periodically wrapped to the opposite side but instead extended into a periodically continued extended domain in the same direction, as shown in figure 8. To be clear, the vector magnetic field in each of the six rectangular black boxes in figure 8 is identical, derived from a single time snapshot of the gyrokinetic simulation.

To investigate the separation of magnetic field lines, we initially select a large number of field lines uniformly distributed within a circular region in $(x,y)$ plane at $z=0$ , (their uniform distribution is depicted in side projection in figure 9 c). These field lines are each followed through the extended simulation box (figure 9 b), after which their distribution appears to be well fit by a Gaussian distribution (figure 9 a, blue line) at $z=3L_{\Vert }$ . The full width at half-maximum (FWHM) of the Gaussian function represents the mean separation $\langle (\unicode[STIX]{x1D6FF}r)\rangle$ of the field line at given distance $z=n\times L_{\Vert }$ . For example, in figure 10 we plot the distribution of the field lines (stars) in $(x,y)$ plane at $z=12L_{\Vert }$ for Run 3 (at time $t/t_{l}=3.6$ ), showing that the distribution of the field lines is well fit by a Gaussian function (blue solid line). Here it is worth mentioning that the centre of the ensemble of the field lines in the $x$ – $y$ plane is changing with distance $z$ . The centre position of the field lines at a given distance $z$ can be calculated as the average of all field line positions in $x$ – $y$ plane at distance $z$ . The distributions of the field lines (see, e.g. figures 9 and 10) are estimated with respect to this moving centre. In our analysis we focus only on the spread of the magnetic field lines.

Figure 9. An initially uniform distribution of field lines within a small circle in the plane at $z=0$ (c) wanders through the simulation domain in Run 3 (b) to generate a Gaussian distribution at $z=3L_{\Vert }$ (a).

Figure 10. For Run 3 with $\overline{\unicode[STIX]{x1D712}}=1.6$ , the binned distribution of the field lines (asterisks) at $z=12L_{\Vert }$ versus the normalized separation length $r/\unicode[STIX]{x1D70C}_{i}$ , where $r=\sqrt{x^{2}+y^{2}}$ . The separation distance is divided into 11 bins. The distribution is fit to a Gaussian function with a FWHM $\langle \unicode[STIX]{x1D6FF}r\rangle$ .

Using this method to trace the magnetic field line separation for Runs 1–6 in table 1, we compute directly the mean square displacement $\langle (\unicode[STIX]{x1D6FF}r)^{2}\rangle$ as a function of $z/L_{\Vert }$ at different evolution times within the saturation phase of each simulation run. Note the we do not fit the FWHM of the best fit Gaussian, but rather calculate directly the mean square displacement $\langle (\unicode[STIX]{x1D6FF}r)^{2}\rangle$ of the magnetic field line distribution in order to account for any skewness in the resulting distribution that may prevent a good Gaussian fit. We fit the computed mean square displacement $\langle (\unicode[STIX]{x1D6FF}r)^{2}\rangle$ by the power-law function

where we expect the amplitude factor $A(\unicode[STIX]{x1D712})$ and the power-law exponent $p(\unicode[STIX]{x1D712})$ both to be a function of the nonlinearity parameter $\unicode[STIX]{x1D712}$ . Note that an exponent of $p=1$ corresponds to a diffusive behaviour of the field line separation, while $p>1$ indicates a superdiffusive process and $p<1$ indicates a subdiffusive process.

Figure 11. For Runs 1–6, plots of the mean square separation $\langle \unicode[STIX]{x1D6FF}r^{2}\rangle$ versus distance $z/L_{\Vert }$ at four different evolution times. The solid curves are the corresponding power-law fitting functions, with coefficients denoted.

Figure 12. (a) The mean value of the amplitude coefficient $\overline{A}$ versus the time-averaged nonlinearity parameter $\overline{\unicode[STIX]{x1D712}}$ . The blue line is the corresponding power-law exponent fitting curve $\overline{A}(\overline{\unicode[STIX]{x1D712}})=0.3\exp (1.1\overline{\unicode[STIX]{x1D712}})$ . (b) The average value of exponent $\overline{p}$ versus the time-averaged nonlinearity parameter $\overline{\unicode[STIX]{x1D712}}$ . The red line is the corresponding fitting curve $\overline{p}=1.46-0.82\exp (-0.80\overline{\unicode[STIX]{x1D712}})$ . The green horizontal line corresponds to value $\overline{p}=1$ which is associated with a standard diffusive process. Vertical and horizontal error bars correspond to standard deviations.

In figure 11, we plot the results of this analysis of the magnetic field line separation, including the values of the fitting parameters $A$ and $p$ . For each simulation run, we plot the results of $\langle (\unicode[STIX]{x1D6FF}r)^{2}\rangle$ as a function of the distance along the equilibrium magnetic field $z/L_{\Vert }$ at four different representative times after the turbulent cascade has saturated. Note that, in addition to the time-averaged nonlinearity parameter $\overline{\unicode[STIX]{x1D712}}$ given for each run in table 1, we also compute the instantaneous nonlinearity parameter $\unicode[STIX]{x1D712}$ at each of the chosen times in figure 11.

From the results in figure 11, we find the following general properties about the separation of magnetic fields lines as a function of the turbulence amplitudes $\unicode[STIX]{x1D712}$ . First, the amplitude of the field line separation, parameterized by the coefficient $A$ , generally increases monotonically with the time-averaged nonlinearity parameter $\overline{\unicode[STIX]{x1D712}}$ . However, it does not always follow that a larger instantaneous value of the nonlinearity parameter $\unicode[STIX]{x1D712}$ will yield a larger value of the coefficient $A$ . This finding suggests that the amplitude of the field line separation may indeed depend not only on the nonlinearity parameter, but also on the history of the turbulent tangling of the magnetic field. Therefore, to determine the functional forms of the amplitude $A(\unicode[STIX]{x1D712})$ and the power-law exponent $p(\unicode[STIX]{x1D712})$ requires a statistical approach, effectively leading us to determine averaged values of these parameters as a function of the averaged nonlinearity parameter, $\overline{A}(\overline{\unicode[STIX]{x1D712}})$ and $\overline{p}(\overline{\unicode[STIX]{x1D712}})$ .

Another important conclusion regards the power-law scaling exponent $p$ of the field line separation as a function of the nonlinearity parameter $\unicode[STIX]{x1D712}$ . Again, in general, we find a larger exponent $p$ for larger values of the time-averaged nonlinearity parameter $\overline{\unicode[STIX]{x1D712}}$ . For the cases of over-strong and strong turbulence in Runs 1–4 with $\overline{\unicode[STIX]{x1D712}}>1$ , the field line separation appears to follow a superdiffusive scaling, $p>1$ , as a function of $z$ ; in contrast, for the weak turbulence cases Runs 5 and 6, it follows a subdiffusive scaling, $p<1$ .

To quantitatively describe the field line separation as a function of the time-averaged nonlinearity parameter $\overline{\unicode[STIX]{x1D712}}$ , we compute the mean values of the amplitude coefficient $\overline{A}$ and the exponent $\overline{p}$ for each of the simulation Runs 1–6. We compute the mean and standard deviation for $\overline{A}$ and $\overline{p}$ by averaging over the four values of $p$ and $A$ shown in figure 11. We plot in figure 12 the average values $\overline{A}$ (a) and $\overline{p}$ (b) as a function of the time-averaged nonlinearity parameter $\overline{\unicode[STIX]{x1D712}}$ . Uncertainties in figure 12 are given by the standard deviation of each of these averaged quantities.

Figure 12(a) clearly shows that $\overline{A}$ nonlinearly increases with increasing $\overline{\unicode[STIX]{x1D712}}$ indicating more spreading of magnetic field in strong turbulence than in weak turbulence. The data are fit to an exponential form

with best fit values $\overline{a}_{0}=0.3$ and $\overline{b}_{0}=1.1$ .

Figure 12(b) displays the dependence of the mean exponent $\overline{p}$ on the time-averaged nonlinearity parameter $\overline{\unicode[STIX]{x1D712}}$ . We fit this behaviour to the functional form

where the best fit yields values $\overline{p}_{0}=1.46$ , $\overline{p}_{1}=0.82$ and $C_{1}=-0.80$ . Based on this statistically averaged functional form, we conclude that the diffusion approximation is a reasonable model for turbulence with a nonlinearity parameter $\overline{\unicode[STIX]{x1D712}}=1.\pm 0.15$ .