Extended Scaling Laws in Numerical Simulations of MHD Turbulence
Magnetised turbulence is ubiquitous in astrophysical systems, where it notoriously spans a broad range of spatial scales. Phenomenological theories of MHD turbulence describe the self-similar dynamics of turbulent fluctuations in the inertial range of scales. Numerical simulations serve to guide and test these theories. However, the computational power that is currently available restricts the simulations to Reynolds numbers that are significantly smaller than those in astrophysical settings. In order to increase computational efficiency and, therefore, probe a larger range of scales, one often takes into account the fundamental anisotropy of field-guided MHD turbulence, with gradients being much slower in the field-parallel direction. The simulations are then optimised by employing the reduced MHD equations and relaxing the field-parallel numerical resolution. In this work we explore a different possibility. We propose that there exist certain quantities that are remarkably stable with respect to the Reynolds number. As an illustration, we study the alignment angle between the magnetic and velocity fluctuations in MHD turbulence, measured as the ratio of two specially constructed structure functions. We find that the scaling of this ratio can be extended surprisingly well into the regime of relatively low Reynolds number. However, the extended scaling becomes easily spoiled when the dissipation range in the simulations is under-resolved. Thus, taking the numerical optimisation methods too far can lead to spurious numerical effects and erroneous representation of the physics of MHD turbulence, which in turn can affect our ability to correctly identify the physical mechanisms that are operating astrophysical systems.
Magnetised turbulence pervades the universe. It is likely to play an important role in the transport of energy, momentum and charged particles in a diverse range of astrophysical plasmas. It is studied with regards to its influence on the generation of magnetic fields in stellar and planetary interiors, small-scale structure and heating of stellar winds, the transport of angular momentum in accretion discs, gravitational collapse and star formation in molecular clouds, the propagation and acceleration of cosmic rays, and interstellar scintillation (e.g., Biskamp, 2003; Kulsrud, 2004; McKee & Ostriker, 2007; Goldstein, Roberts & Matthaeus, 1995; Brandenburg & Nordlund, 2011; Schekochihin & Cowley, 2007). The effects of magnetised turbulence need to be taken into account when analysing astrophysical observations and also when modelling astrophysical processes.
The simplest theoretical framework that describes magnetised plasma turbulence is that of incompressible magnetohydrodynamics (MHD),
where the Elsässer variables are defined as , is the fluctuating plasma velocity, is the fluctuating magnetic field normalized by , is the Alfvén velocity based upon the uniform background magnetic field , , is the plasma pressure, is the background plasma density, represents forces that drive the turbulence at large scales and for simplicity we have taken the case in which the fluid viscosity is equal to the magnetic resistivity. Energy is transferred to smaller scales by the nonlinear interactions of oppositely propagating Alfvén wavepackets (Kraichnan, 1965). This can be inferred directly from equation (1) by noting that in the absence of forcing and dissipation, if then any function is an exact nonlinear solution that propagates parallel and anti-parallel to with the Alfvén speed. The efficiency of the nonlinear interactions splits MHD turbulence into two regimes. The regime in which the linear terms dominate over the nonlinear terms is known as ‘weak’ MHD turbulence, otherwise the turbulence is ‘strong’. In fact, it has been demonstrated both analytically and numerically that the MHD energy cascade occurs predominantly in the plane perpendicular to the guiding magnetic field. This ensures that even if the turbulence is weak at large scales it encounters the strong regime as the cascade proceeds to smaller scales. MHD turbulence in astrophysical systems is therefore typically strong.
For strong MHD turbulence, Goldreich & Sridhar (1995) argued that the linear and nonlinear terms in equations (1) should be approximately balanced at all scales, known as the critical balance condition. Consequently, Goldreich & Sridhar (1995) postulated that the wave packets get progressively elongated in the direction of the guide field as their scale decreases (with the field-parallel lengthscale and field-perpendicular scale related by ) and that the field-perpendicular energy spectrum takes the Kolmogorov form .
Recent high resolution direct numerical simulations with a strong guide field () do indeed verify the strong anisotropy of the turbulent fluctuations, however, the field-perpendicular energy spectrum appears to be closer to (e.g., Maron & Goldreich, 2001; Müller & Grappin, 2005; Mason et al., 2008; Perez & Boldyrev, 2009, 2010). A resolution to this contradiction was proposed in Boldyrev (2006). Therein it was suggested that in addition to the elongation of the eddies in the direction of the guiding field, the fluctuating velocity and magnetic fields at a scale are aligned within a small scale-dependent angle in the field perpendicular plane, . In this model the wavepackets are three-dimensionally anisotropic. Scale-dependent dynamic alignment reduces the strength of the nonlinear interactions and leads to the field-perpendicular energy spectrum .
Although the two spectral exponents and are close together in numerical value, the physics of the energy cascade in each model is different. The difference between the two exponents is especially important for inferring the behaviour of processes in astrophysical systems with extended inertial intervals. For example, the two exponents can lead to noticeably different predictions for the rate of turbulent heating in coronal holes and the solar wind (e.g., Chandran et al., 2010; Chandran, 2010; Podesta & Borovsky, 2010). Thus, there is much interest in accurately determining the spectral slope from numerical simulations. Unfortunately, the Reynolds numbers that are currently accessible by most direct numerical simulations do not exceed a few thousand, which complicates the precise identification of scaling exponents. Techniques for careful optimisation of the numerical setup and alternative ways of differentiating between the competing theories are therefore much sought after.
Maximising the extent of the inertial range is often achieved by implementing physically motivated simplifying assumptions. For example, since the turbulent cascade proceeds predominantly in the field-perpendicular plane it is thought that the shear-Alfvén waves control the dynamics while the pseudo-Alfvén waves play a passive role (see, e.g., Maron & Goldreich (2001)). If one neglects the pseudo-Alfvén waves (i.e. removes the fluctuations parallel to the strong guide field) one obtains a system that is equivalent to the reduced MHD system (RMHD) that was originally derived in the context of fusion devices by Kadomtsev & Pogutse (1974) and Strauss (1976) (see also Biskamp (2003); Oughton, Dmitruk & Matthaeus (2004)). Incompressibility then enables the system to be further reduced to a set of two scalar equations for the Elsässer potentials, resulting in a saving of approximately a factor of two in computational costs. Further computational savings can be made by making use of the fact that the wavepackets are elongated. Hence variations in the field-parallel direction are slower than in the field-perpendicular plane and a reduction in the field-parallel resolution would seem possible. Indeed, this is widely used as an optimisation tool in numerical simulations of the inertial range of field-guided MHD turbulence (e.g., Maron & Goldreich, 2001; Mason et al., 2008; Grappin & Müller, 2010; Perez & Boldyrev, 2010). The accumulated computational savings can then be re-invested in reaching larger Reynolds numbers for the field-perpendicular dynamics.
Additionally, it is advantageous to seek other ways of probing the universal scaling of MHD turbulence. In this work we point out a rather powerful method, which is based on the fact that there may exist certain quantities in MHD turbulence that exhibit very good scaling laws even for turbulence with relatively low Reynolds numbers. The situation here is reminiscent of the well known phenomenon of extended self-similarity in hydrodynamic turbulence (Benzi et al., 1993). We propose that one such “stable” object is the alignment angle between the velocity and magnetic fluctuations, which we measure as the ratio of two specially constructed structure functions. This ratio has been recently measured in numerical simulations in an attempt to differentiate among various theoretical predictions (Beresnyak & Lazarian (2006); Mason et al. (2006)). Also, it has recently been shown by Podesta et al. (2009) that the same measurement is accessible through direct observations of solar wind turbulence. Scale-dependent alignment therefore has practical value: its measurement may provide an additional way of extracting information about the physics of the turbulent cascade from astrophysical observations. In the present work we conduct a series of numerical simulations with varying resolutions and Reynolds numbers. We find that as long as the simulations are well resolved, the alignment angle exhibits a universal scaling behavior that is virtually independent of the Reynolds number of the turbulence. Moreover, we find that the length of scaling range for this quantity extends to the smallest resolved scale, independently of the Reynolds number. This means that although the dissipation spoils the power-law scaling behaviour of each of the structure functions, the dissipation effects cancel when the ratio of the two functions is computed and the universal inertial-range scaling extends deep in the dissipation region.
The described method allows the inference of valuable scaling laws from numerical simulations, experiments, or observations of MHD turbulence with limited Reynolds number. However, one can ask how well the extended-scaling method can be combined with the previously mentioned optimisation methods relying on the reduced MHD equations and a decreased parallel resolution. We check that reduced MHD does not alter the result. However, when the dissipation region becomes under-resolved (as can happen, for example, when the field-parallel resolution is decreased), the extended scaling of the alignment angle deteriorates significantly. Thus the optimisation technique that works well for viewing the inertial range of the energy spectra should not be used in conjunction with the extended-scaling measurements that probe deep into the dissipation region. The remainder of this paper will report the findings of a series of numerical measurements of the alignment angle in simulations with different Reynolds numbers and different field-parallel resolutions in both the MHD and RMHD regimes. The aim is to address the need to find an optimal numerical setting for studying strong MHD turbulence and to raise caution with regards to the effects that implementing simplifying assumptions in the numerics can have on the solution and its physical interpretation.
2 Numerical results
We simulate driven incompressible magnetohydrodynamic turbulence in the presence of a strong uniform background magnetic field, . The MHD code solves equations (1,2) on a periodic, rectangular domain with aspect ratio , where the subscripts denote the directions perpendicular and parallel to and we take . A fully dealiased 3D pseudospectral algorithm is used to perform the spatial discretisation on a grid with a resolution of mesh points. The RMHD code solves the reduced MHD counterpart to equations (1,2) in which (see Perez & Boldyrev (2008)).
The domain is elongated in the direction of the guide field in order to accommodate the elongated wavepackets and to enable us to drive the turbulence in the strong regime while maintaining an inertial range that is as extended as possible (see Perez & Boldyrev (2010)). The random forces are applied in Fourier space at wavenumbers , , where we shall take or . The forces have no component along and are solenoidal in the -plane. All of the Fourier coefficients outside the above range of wavenumbers are zero and inside that range are Gaussian random numbers with amplitudes chosen so that . The individual random values are refreshed independently on average every , i.e. the force is updated approximately times per turnover of the large-scale eddies. The variances control the average rates of energy injection into the and fields. The results reported in this paper are for the balanced case . In all of the simulations performed in this work we will set the background magnetic field in velocity units. Time is normalised to the large scale eddy turnover time . The field-perpendicular Reynolds number is defined as .
The system is evolved until a stationary state is reached, which is confirmed by observing the time evolution of the total energy of the fluctuations, and the data are then sampled in intervals of the order of the eddy turnover time. All results presented correspond to averages over approximately 30 samples. We conduct a number of MHD and RMHD simulations with different resolutions, Reynolds numbers and field-parallel box sizes. The parameters for each of the simulations are shown in Table 1.
For each simulation we calculate the scale-dependent alignment angle between the shear-Alfvén velocity and magnetic field fluctuations. We therefore define velocity and magnetic differences as and , where is a point-separation vector in the plane perpendicular to . In the MHD case the pseudo-Alfvén fluctuations are removed by subtracting the component that is parallel to the local guide field, i.e. we construct (and similarly for ) where . In the RMHD case fluctuations parallel to are not permitted and hence the projection is not necessary. We then measure the ratio of the second order structure functions
where the average is taken over different positions of the point in a given field-perpendicular plane, over all such planes in the data cube, and then over all data cubes. By definition of the cross product
where is the angle between and and the last approximation is valid for small angles. We recall that the theoretical prediction is (Boldyrev, 2006).
Figure 1 illustrates the ratio (4) as a function of the separation for two MHD simulations (M1 and M2) corresponding to a doubling of the resolution from to mesh points with the Reynolds number increased from to . Excellent agreement with the theoretical prediction is seen in both cases. As the resolution and Reynolds number increase, the scale-dependence of the alignment angle persists to smaller scales. Indeed, we believe that the point at which the alignment saturates can be identified as the dealiasing scale, corresponding in configuration space to for the simulations, respectively. This is verified in Figure 2 that shows that alignment is largely insensitive to the Reynolds number (provided that the system is turbulent) and Figure 3 that shows that the saturation point decreases by a factor of approximately 2 as the resolution doubles at fixed Reynolds number. Thus as computational power increases, allowing higher resolution simulations to be conducted, we expect to find that scale-dependent alignment persists to smaller and smaller scales.
The fact that even in the lower Reynolds number cases scale-dependent alignment is clearly seen over quite a wide range of scales is particularly interesting, as in those cases only a very short inertial range can be identified in the field-perpendicular energy spectrum, making the identification of spectral exponents difficult (see Figure 1 in Mason et al. (2006)). In the larger cases, we can estimate the inertial range of scales in configuration space to be the range of over which the energy spectrum displays a power law dependence. The field-perpendicular energy spectrum for the Case M2 is shown in Figure 1 in Mason et al. (2008), with the inertial range corresponding to approximately , i.e. . Comparison with Figure 1 shows that a significant fraction of the region over which the scaling is observed corresponds to the dissipative region, i.e. that ratios of structure functions appear to probe deeper than the inertial range that is suggested by the energy spectra.
We now consider the effect on the alignment ratio of decreasing the field-parallel resolution. Figure 4 shows the results from three MHD simulations (M3, M4 & M5) for which the field-parallel resolution decreases by a factor of two, twice. As the resolution decreases the extent of the self-similar region diminishes and the scale-dependence of the alignment angle becomes shallower. If one were to calculate the slope for the lowest field-parallel resolution case (M5) one would find a scale-dependence that is shallower than the predicted power law exponent of . This may lead one to conclude (incorrectly) that scale-dependent alignment is not a universal phenomenon in MHD turbulence. However, the effect is obviously a result of the poor resolution rather than being an attribute of the alignment mechanism itself.
Finally, we mention that for the three cases illustrated in Figure 4, the field-perpendicular energy spectra (not shown) display no appreciable difference. Since the Reynolds number is moderate the inertial range in -space is quite short. However, when the spectra are compensated with and the former results in a better fit in all cases. This happens for two reasons. First, the stronger deviation from the alignment scaling occurs deeper in the dissipation region, that is, further from the inertial interval where the energy spectrum is measured. Second, according to the relationship between the scaling of the alignment angle and the energy spectrum, a noticeable change in the scaling of the alignment angle leads to a relatively small change in the scaling of the field-perpendicular energy spectrum.
There are two main conclusions that can be drawn from our results. The first is that the measurement of the alignment angle, which is composed of the ratio of two structure functions, appears to display a self-similar region of significant extent, even in the moderate Reynolds number case which requires only a moderate resolution. We have checked that plotting the numerator and the denominator of the alignment ratio separately as functions of the increment displays only a very limited self-similar region, from which scaling laws cannot be determined. A clear scaling behaviour is also not found when one plots the numerator versus the denominator as is the case in extended self-similarity (Benzi et al., 1993). The result is interesting in its own right. It also has important practical value as it allows us to differentiate effectively between competing phenomenological theories through numerical simulations conducted in much less extreme parameter regimes than would otherwise be necessary. The result could be especially useful if it extends to ratios of structure functions for which an exact relation, such as the Politano & Pouquet (1998) relations, is known for one part, as it would then allow the inference of the scaling of the other structure function. Reaching a consensus on the theoretical description of magnetised fluctuations in the idealised incompressible MHD system represents the first step towards the ultimate goal of building a theoretical foundation for astrophysical turbulence.
The second main result that can be drawn from our work is that the measurement of the alignment angle appears to probe deep into the dissipation region and hence it is necessary to adequately resolve the small scale physics. As the field-parallel resolution is decreased, numerical errors contaminate the physics of the dissipative range and affect measurement of the alignment angle. As the decrease in resolution is taken to the extreme, the errors propagate to larger scales and may ultimately spoil an inertial range of limited extent. We propose that similar contamination effects should also arise through any mechanism that has detrimental effects on the dissipative physics. Mechanisms could include pushing the Reynolds number to the extreme or using hyperdiffusive effects. For example, our results may provide an explanation for the numerical findings by Beresnyak (2011) who noticed a flattening of the alignment angle in simulations of MHD turbulence with a reduced parallel resolution and strong hyperdiffusivity.
We also point out that the result recalls the phenomenon of extended self-similarity in isotropic hydrodynamic turbulence (Benzi et al., 1993), which refers to the extended self-similar region that is found when one plots one structure function versus another, rather than as a function of the increment. Our finding is fundamentally different however, in the sense that the self-similar region only becomes apparent when one plots ratios of structure functions versus the increment, rather than structure functions versus other structure functions. Our result appears to be due to a non-universal features in the amplitudes of the functions, rather than their arguments, cancelling when the ratios are plotted. Whether such a property holds for other structure functions in MHD turbulence is an open and intriguing question. This is a subject for our future work.
- Benzi et al. (1993) Benzi, R., Ciliberto, S., Tripiccione, R., Baudet, C., Massaioli, F. & Succi, S. 1993, Phys. Rev. E, 48, R29
- Beresnyak (2011) Beresnyak, A. 2011, Phys. Rev. Lett., 106, 075001
- Beresnyak & Lazarian (2006) Beresnyak, A. & Lazarian, A. 2006, ApJ, 640, L175
- Biskamp (2003) Biskamp, D. 2003, Magnetohydrodynamic Turbulence (Cambridge University Press, Cambridge)
- Boldyrev (2006) Boldyrev, S. 2006, Phys. Rev. Lett., 96, 115002
- Boldyrev et al. (2009) Boldyrev, S., Mason, J. & Cattaneo, F. 2009, ApJ, 699, L39.
- Brandenburg & Nordlund (2011) Brandenburg, A. & Nordlund A. 2011, Rep. Prog. Phys., 74, 046901
- Chandran et al. (2010) Chandran, B. D. G., Li, B., Rogers, B. N., Quataert, E., & Germaschewski, K. 2010, ApJ, 720, 503
- Chandran (2010) Chandran, B. D. G. 2010, ApJ, 720, 548
- Goldreich & Sridhar (1995) Goldreich, P. & Sridhar, S. 1995, ApJ, 438, 763
- Goldstein, Roberts & Matthaeus (1995) Goldstein, M. L., Roberts, D. A., & Matthaeus, W. H. 1995, Ann. Rev. Astron. Astrophys., 33, 283
- Grappin & Müller (2010) Grappin, R. & Müller, W.-C. 2010, Phys. Rev. E, 82, 026406
- Kadomtsev & Pogutse (1974) Kadomtsev, B. B. & Pogutse, O. P. 1974, Sov. Phys. JETP, 38, 283
- Kolmogorov (1941) Kolmogorov, A.N. 1941, Dokl. Akad. Nauk SSSR, 32, 16 (reprinted in Proc. R. Soc. Lond. A, 434, (1991) 15)
- Kraichnan (1965) Kraichnan, R. H. 1965, Phys. Fluids, 8, 1385
- Kulsrud (2004) Kulsrud, R. M. 2004, Plasma physics for astrophysics (Princeton University Press)
- Maron & Goldreich (2001) Maron, J., & Goldreich, P. 2001, ApJ, 554, 1175
- Mason et al. (2006) Mason, J., Cattaneo, F. & Boldyrev, S. 2006, Phys. Rev. Lett., 97, 255002
- Mason et al. (2008) Mason, J., Cattaneo, F. & Boldyrev, S. 2008, Phys. Rev. E, 77, 036403
- McKee & Ostriker (2007) McKee, C. F. & Ostriker, E. C., 2007, Ann. Rev. Astron. Astrophys., 45, 565
- Müller & Grappin (2005) Müller, W.-C. & Grappin, R. 2005, Phys. Rev. Lett., 95, 114502
- Oughton, Dmitruk & Matthaeus (2004) Oughton, S., Dmitruk, P. & Matthaeus, W. H. 2004, Physics of Plasmas, 11, 2214
- Perez & Boldyrev (2008) Perez, J. C. & Boldyrev, S. 2008, Astrophys. J., 672, L61
- Perez & Boldyrev (2009) Perez, J. C. & Boldyrev, S. 2009, Phys. Rev. Lett., 102, 025003
- Perez & Boldyrev (2010) Perez, J. C. & Boldyrev, S. 2010, Phys. Plasmas, 17, 055903
- Perez et al. (2011) Perez, J.C., Mason, J., Boldyrev, S. & Cattaneo, F. 2011, in preparation
- Podesta & Borovsky (2010) Podesta, J. J. & Borovsky, J. E. 2010, Phys. Plasmas, 17, 112905
- Podesta et al. (2009) Podesta, J. J., Chandran, B. D. G., Bhattacharjee, A., Roberts, D. A., & Goldstein, M. L. 2009, J. Geophys. Res., 114, A01107
- Politano & Pouquet (1998) Politano, H. & Pouquet, A. 1998, Geophys. Res. Lett, 25, 273
- Schekochihin & Cowley (2007) Schekochihin, A. A. & Cowley, S. C. 2007, in Magnetohydrodynamics: Historical Evolution and Trends, (Springer, Dordrecht), 85
- Strauss (1976) Strauss, H. 1976, Phys. Fluids, 19, 134