Scaling and anisotropy of magnetohydrodynamic turbulence in a strong mean magnetic field
We present a new analysis of the anisotropic spectral energy distribution in incompressible magnetohydrodynamic (MHD) turbulence permeated by a strong mean magnetic field. The turbulent flow is generated by high-resolution pseudo-spectral direct numerical simulations with large-scale isotropic forcing. Examining the radial energy distribution for various angles with respect to reveals a specific structure which remains hidden when not taking axial symmetry with respect to into account. For each direction, starting at the forced large-scales, the spectrum first exhibits an amplitude drop around a wavenumber which marks the start of a scaling range and goes on up to a dissipative wavenumber . The 3D spectrum for is described by a single -independent functional form , the scaling law being the same in every direction. The previous properties still hold when increasing the mean field from up to , as well as when passing from resistive to ideal flows. We conjecture that at fixed the direction-independent scaling regime is reached when increasing the Reynolds number above a threshold which raises with increasing . Below that threshold critically balanced turbulence is expected.
pacs:47.65.+a, 47.27.Eq, 47.27.Gs
It is known that in presence of a mean magnetic field assumed here to point in the -direction, , nonlinear interactions in incompressible magnetohydrodynamics (MHD) are weakened in the field-parallel direction. The MHD approximation allows to describe the large-scale dynamics of astrophysical plasmas, i.e. ionized gases, like the interstellar medium or the solar corona. Due to the above-mentioned anisotropy, the nonlinear energy transfer in MHD turbulence proceeds preferably to larger perpendicular spatial wavenumbers Strauss (1976); Montgomery and Turner (1981); Shebalin et al. (1983); Grappin (1986). Direct numerical simulations (DNS) show that the field-perpendicular energy spectrum exhibits self-similar inertial-range scaling in wavenumber with (Cho and Vishniac, 2000; Cho et al., 2002) for weak to moderate , or (Maron and Goldreich, 2001; Müller and Grappin, 2005; Mininni and Pouquet, 2007; Mason et al., 2008) for strong .
Iroshnikov (1963) and Kraichnan (1965) proposed the first theory of the effect of a mean magnetic field on incompressible MHD turbulence. They remarked that any flow can be decomposed into a sum of weakly interacting waves with different wavevectors , the term weak interaction meaning that the characteristic time of deformation of the waves is much longer than their periods. This led to the prediction of a slow cascade, with a spectral slope different from the Kolmogorov prediction .
This theory used an isotropic measure of the propagation time based on the modulus of the wave vector,
ignoring deliberately the waves with wave vectors perpendicular (or almost perpendicular) to the mean field, for which the deformation time should clearly be smaller than their period, and hence the interaction strong. This was criticized by Goldreich and Sridhar (1995), who denied the possibility of the previous weak cascade to occur, and argued that the perpendicular strong cascade leading to should be the only one present. For a large enough mean field, the perpendicular cascade should thus be restricted to a thin subset around the axis in Fourier space, the subset becoming thinner when the mean field increases.
More precisely, the subdomain in (, ) space where the perpendicular cascade is believed to occur is defined by a critical balance (Goldreich and Sridhar, 1995) between the characteristic time of nonlinear interaction, , and the Alfvén time , the fluctuations becoming correlated along the guide field up to a distance , where is the typical magnitude of fluctuations at the scale . Assuming a scaling law in the perpendicular direction, and spectral transfer dominated by strong coupling, i.e., , one obtains the 3D spectrum:
where , with and , for and negligible for .
Eq. (2) actually suffers from two limitations when the mean field is significantly larger than the magnetic fluctuation rms value. First, the spectral slope is observed to become (Maron and Goldreich, 2001; Müller and Grappin, 2005; Mininni and Pouquet, 2007; Mason et al., 2008) instead of the strong cascade value ; second, the time-scales ratio becomes significantly smaller than unity (Bigot et al., 2008) in the excited part of the spectrum, showing that the strict critical balance condition is too restrictive to describe the anisotropy of the cascade, or, in other words, that the cascade is more extended in the oblique directions than predicted by the critical balance condition. This has led several authors to suggest modifications which either still assume that the anisotropy is dictated by the critical balance condition (Boldyrev, 2006; Gogoberidze, 2007), or propose that the time-scales ratio decreases with increasing (Galtier et al., 2005). All these phenomenologies predict a spectral form different from Eq. (2), with different spectral slopes in the perpendicular and parallel directions.
In the present paper, we analyze the angular spectrum and find by taking slices along the radial directions that a unique spectral form (and slope) holds in all directions, with the radial power-law range extent depending on the angle with the mean field. As a result, a significant portion of this spectrum lies in a domain where the time-scales ratio is sub-critical, that is, much smaller than unity.
This study focuses on representative states of fully-developed turbulence permeated by a strong mean magnetic field with from high-resolution direct numerical simulations of quasi-stationary MHD turbulence forced at large scales. The forcing is realized by freezing all modes (velocity and magnetic field) with in an energetically roughly isotropic and equipartitioned state. The driving of magnetic energy could be realized physically by large-scale fluctuations of electrical current although here it is mainly applied to achieve a state of approximate equipartion of kinetic and magnetic energy. Decaying test simulations have confirmed that the forcing does not modify the results presented in the following. The dimensionless equations of resistive MHD formulated with the vorticity and the magnetic field are given by
The equations are solved by a standard pseudospectral method with spherical mode truncation to alleviate aliasing errors. The numerical resolution is collocations points with reduced resolution in the direction of (Müller and Grappin, 2005) and with kinematic viscosity and resistivity set to . The analyzed data is the temporal average of five snapshots of the three-dimensional Fourier energy distribution taken equidistantly within about four to five field-perpendicular large-scale turnover times, , of quasi-stationary turbulence where , , cf. (Zikanov and Thess, 2004), and with , . The normalized cross helicity and the Alfvén ratio fluctuate around and .
The three-dimensional (3D) energy spectrum, , relates to the total energy .
Fig. 1 shows contour levels of in two mutually orthogonal planes containing the origin. The field-perpendicular - plane (Fig. 1,a) displays an isotropic energy distribution, as expected. Anisotropy induced by appears in planes containing the direction (Fig. 1,b).
Spectral anisotropy is traditionally diagnosed by 1D spectra, e.g., and . These are shown compensated by (a) and by (b) in Fig. 2. The perpendicular spectrum exhibits a power-law with , in agreement with previous works (Maron and Goldreich, 2001; Müller and Grappin, 2005; Mininni and Pouquet, 2007; Mason et al., 2008), while the parallel spectrum does not show any convincing scaling range. An anisotropic scaling law is build from the previous 1D spectra by plotting (Fig. 1,d) the modes sharing the same 1D energy density (Bigot et al., 2008). The anisotropy exponent is seen to lie between and , which is also obtained in (Bigot et al., 2008) for .
While planar integration yields some information about anisotropy, it mixes all wavenumbers perpendicular to the chosen direction and thus blurs the separation between inertial and dissipative scales if the dissipative scale is not constant over the planar domain of integration. More importantly, no information on intermediate directions between parallel and perpendicular is available. Thus, spherical coordinates with respect to the mean field axis along are considered. As is isotropic in the azimuthal plane, the -dependence of is eliminated by averaging over which strongly decreases statistical noise and yields the -averaged 3D spectrum whose isocontours in () are shown in Fig. 1,c. We define the corresponding one-dimensional (1D) spectrum as . The total energy is thus .
The properties of are shown in Fig. 3. A scaling range is seen starting at – down to a dissipative wavenumber . The 1D scaling exponent
is independent of . The anisotropy appears at fixed as a -dependence of the spectral amplitude and of the dissipative wavenumber . Normalizing the wavenumber by shows that the spectrum follows a single functional form whatever :
where is a constant (the amplitude of the spectrum at the dissipative scale). Note that the first equality holds also beyond the dissipative range, the second being valid for . A corollary is that :
Fig. 3,a shows the self-similar wavenumber intervals of all spectra, starting in the range – with a slight dependence on . The dissipative wavenumber is estimated by locating the maximum of in each direction : varying from to leads to a drop of from about to while the spectral energy at fixed decreases in the inertial range from to . Fig. 3,b indicates that the spectral -dependence can be nearly eliminated by normalizing with (Eq. (4)). The relation between spectral amplitude and as given by Eq. (5) is confirmed by Fig. 3,c which shows that (dotted curve) closely follows (Eq. (3)), with to eliminate the dissipative range in the parallel direction.
A simple model of anisotropic spectrum with a spectral exponent being the same in all directions has been proposed previously in the context of shell models of turbulence Carbone and Veltri (1990) as well as solar wind turbulence Tessein et al. (2009); Carbone et al. (1995). It reads:
We tried to use this model to adjust the energy contours of our numerical simulations, and found that the global anisotropy between the perpendicular and parallel amplitude requires . However, the model fails to reproduce correctly the detailed angular anisotropy, that is, the contours in oblique directions, because our energy contours differ much from ellipsoids, as is seen in Fig. 3,d which shows a zoom of the energy contours as solid lines (the model of Eq. (6) would produce circular contours in this figure). However we found that switching from the to for the exponents of the and as:
leads to a reasonable good fit to the simulation results in all directions, as seen both in the dotted-dashed curve in Fig. 3,c for the amplitude variation vs and in the dotted contours in the plane in Fig. 3,d.
The energy contours of the critical balance spectrum (Eq. 2, with and ) are also represented by dashed lines in Fig. 3,d. They isolate a small cone about the axis in the whole plane, so excluding a large part of the angular structure of the true angular spectrum.
Note that the dissipative wavenumber is determined up to an error of about a factor two (due to errors in interpolating the spectrum), which leads to the noisy appearance of the curve of in Fig. 3,c, to the finite thickness of the normalized spectra in Fig. 3,b and to variations of about a factor 4 in the constant in Eq. (4).
To determine the scaling exponent of , the 1D spectra averaged over four -intervals and compensated by and are shown in Fig. 4. The spectrum with is represented by the bold line. The -scaling is seen to be dominant, except possibly for group B which follows (Fig. 4,a). The extent of the inertial range is shown by oblique dashed lines; it is bounded on the right by , and on the left by an intermediate range which separates the inertial from the forcing range. The start of the inertial range is thus growing from to for .
Increasing the mean field up to in test simulations (not shown) leads to further decrease of the power-law range in the parallel direction, with the perpendicular range increasing slightly, and the parallel range decreasing substantially, so that the ratio of both ranges varies with as
Note that the fit by Eq. (7) remains as good as in Fig. 3,d after decreasing by a factor , as expected since Eq. (7) implies hence from Eq. (8) . Increasing again up to confirmed this trend, but the parallel range becomes too small in that case to allow a good determination of the associated dissipative wavenumber. This precludes checking that the spectrum normalized by is angle-independent for small . This difficulty could be alleviated by increasing the numerical resolution which would allow increasing the Reynolds number. We choose instead below to compare with ideal MHD simulations with resolution and with and .
In the ideal case, Fourier space is divided in a large scale range presenting spectral properties close to those of a standard turbulent spectrum with dissipation, and a small scale range where the spectral slope increases, the latter scales playing the role of a dissipative range (cf. (Cichowlas et al., 2005) in the hydrodynamic case). The boundary between the two domains slowly shifts with time to ever larger scales. It can be identified with the dissipative wavenumber, thus determined here as the minimum of the 1D spectrum. The simulations are initialized with a quasi-stationary state of the resistive run and are continued without any dissipation and with the chosen in the same numerical setup until the energetically rising small scales begin to pollute the scaling region. Choosing the appropriate time, power-law ranges in all directions can be properly identified even with a large field . The resulting power-law is found to be now both with and , contrary to the resistive runs. This difference in scaling between the resistive and ideal runs might be attributed to a different role played by the bottlenkeck effect Beresnyak and Lazarian (2010) in these two setups. All reported findings, cf. Eqs.(3-8), are however confirmed when setting , as seen in Fig. 5 which shows (cf. Fig. 3,c) the anisotropy in two ideal runs with (a) and (b).
Let us come back to the resistive case. As already found in (Bigot et al., 2008), we find that the excited part of the space is not restricted to regions where . The important point is that the 3D-energy contours, and as well the boundary of the power-law range, ignore the iso-contours of as shown in Fig. 3 (see also Fig 5b,c): the form of the 3D-energy contours, as well as the angle-independent spectral slope, actually suggest an isotropic cascade, that is, a cascade along radial directions. To test this idea, we define a -dependent effective Reynolds based on the energy density at . For , increases by a factor 30, while grows by a factor 10, as
where , being the 1D slope. The exponent in Eq. (9) is substantially smaller than the value for (or in the case ) obtained by equating the input flux at and the dissipative flux ). This means that if increases from zero to , the inertial range increases more slowly than it would if the dissipative loss at small scales would balance the input energy rate at the () wavenumber which marks the large-scale boundary of the inertial range. Hence, for the nonlinear radial energy flux must be depleted while the contrary is true for the parallel directions.
By examining the solar wind turbulence, it has been shown Tessein et al. (2009) that the power-law index of the fluctuation spectra is independent of the angle between the wave vector and the mean interplanetary magnetic field, which is fully compatible with the results reported here from direct simulations. Other studies of solar wind turbulence have however reached a different conclusion (Horbury et al., 2008). The latter study used wavelet transforms, which allows to define parallel and perpendicular directions with respect to local averages of the magnetic field. Indeed, it has been argued by (Cho and Vishniac, 2000; Cho et al., 2002) that the critical balance phenomenon, and the associated spectral laws, (in particular the anisotropy index relating the perpendicular and parallel wave numbers) emerge only when considering, instead the mean field, the local average field. In the work by (Bigot et al., 2008) who consider as we do only the global mean field, the law appears clearly only when the mean field is large enough, which can be explained by the fact that in this limit the local and mean field approach coincide. According to this viewpoint, the results reported here should be a simple artifact of the fact that our frame is not attached to the local mean field, but to the global mean field.
However, while the effect could indeed appear for interplanetary turbulence where the fluctuation level is high, it is hardly the case here, since or . Besides, it is remarkable that the fit by our anisotropy function in Fig. 5 is as good when as when which is an indication that the spectral properties of the turbulence are correctly revealed by using our method.
A possible way to reconcile both pictures is the following. It takes into account the fact that there is a second parameter, the Reynolds number. Indeed, as the anisotropy increases with , we find that the scaling range in the parallel direction decreases accordingly (Eq. (8)), which implies that, at a fixed viscosity, the parallel power-law range disappears at even moderate . In our case, and when , while the start of the scaling range is , preventing to increase significantly . The regime at high will thus depend on the Reynolds number (). Increasing at fixed perpendicular Re depletes the field-parallel cascade and could ultimately lead to critically balanced turbulence. If and are however large enough, e.g. under astrophysical conditions, allowing for scaling in all directions then the properties described in this work are expected to hold. We thus propose direction-independent scaling for high and critically balanced turbulence at low . The crossover -value is expected to increase with . A strong indication in this sense is found in (Bigot et al., 2008) where the anisotropy scaling exponent is seen to be between 1 and at as here (Fig. 1,d), while it begins to cluster around at . Such a good agreement with is found already at in (Cho and Vishniac, 2000) because of a moderate Reynolds number.
We have reported here two previously unknown properties of the MHD angular energy spectra: (i) its functional form is (ii) the anisotropy function is best expressed (Eq. 8) as the ratio of the perpendicular over the parallel power-law range extent, which scales linearly with in the interval considered here. These results offer important tests for future theories of anisotropic turbulence.
Acknowledgements.We thank G. Belmont, J. Léorat, and A. Busse for several fruitful discussions.
- Strauss (1976) H. R. Strauss, Phys. Fluids 19, 134 (1976).
- Montgomery and Turner (1981) D. Montgomery and L. Turner, Phys. Fluids 24, 825 (1981).
- Shebalin et al. (1983) J. V. Shebalin, W. H. Matthaeus, and D. Montgomery, Journal of Plasma Physics 29, 525 (1983).
- Grappin (1986) R. Grappin, Physics of Fluids 29, 2433 (1986).
- Cho and Vishniac (2000) J. Cho and E. T. Vishniac, The Astrophysical Journal 539, 273 (2000).
- Cho et al. (2002) J. Cho, A. Lazarian, and E. T. Vishniac, The Astrophysical Journal 564, 291 (2002).
- Maron and Goldreich (2001) J. Maron and P. Goldreich, The Astrophysical Journal 554, 1175 (2001).
- Müller and Grappin (2005) W.-C. Müller and R. Grappin, Phys. Rev. Lett. 95, 114502 (2005).
- Mininni and Pouquet (2007) P. D. Mininni and A. Pouquet, Phys. Rev. Lett. 99, 254502 (2007).
- Mason et al. (2008) J. Mason, F. Cattaneo, and S. Boldyrev, Phys. Rev. E 77, 36403 (2008).
- Iroshnikov (1963) P. S. Iroshnikov, Astronomicheskii Zhurnal 40, 742 (1963).
- Kraichnan (1965) R. H. Kraichnan, Physics of Fluids 8, 1385 (1965).
- Goldreich and Sridhar (1995) P. Goldreich and S. Sridhar, Astrophysical Journal 438, 763 (1995).
- Bigot et al. (2008) B. Bigot, S. Galtier, and H. Politano, Phys. Rev. E 78, 66301 (2008).
- Boldyrev (2006) S. Boldyrev, Phys. Rev. Lett. 96, 115002 (2006).
- Gogoberidze (2007) G. Gogoberidze, Physics of Plasmas 14, 2304 (2007).
- Galtier et al. (2005) S. Galtier, A. Pouquet, and A. Mangeney, Physics of Plasmas 12, 2310 (2005), (c) 2005: American Institute of Physics.
- Zikanov and Thess (2004) O. Zikanov and A. Thess, Applied Mathematical Modelling 28, 1 (2004).
- Carbone and Veltri (1990) V. Carbone and P. Veltri, Geophysical and Astrophysical Fluid Dynamics 52, 153 (1990).
- Tessein et al. (2009) J. A. Tessein, C. W. Smith, B. T. MacBride, W. H. Matthaeus, M. A. Forman, and J. E. Borovsky, The Astrophysical Journal 692, 684 (2009).
- Carbone et al. (1995) V. Carbone, F. Malara, and P. Veltri, Journal of Geophysical Research 100, 1763 (1995).
- Cichowlas et al. (2005) C. Cichowlas, P. Bonaïti, F. Debbasch, and M. Brachet, Phys. Rev. Lett. 95, 264502 (2005).
- Beresnyak and Lazarian (2010) A. Beresnyak and A. Lazarian, arXiv astro-ph.GA (2010), eprint 1002.2428v2.
- Horbury et al. (2008) T. S. Horbury, M. Forman, and S. Oughton, Phys. Rev. Lett. 101, 175005 (2008).