Scaling and anisotropy of magnetohydrodynamic turbulence in a strong mean magnetic field
aa r X i v : . [ phy s i c s . p l a s m - ph ] A ug Scaling and anisotropy of magnetohydrodynamic turbulence in a strong meanmagnetic field
Roland Grappin ∗ LUTH, Observatoire de Paris
Wolf-Christian M¨uller † Max-Planck Institut f¨ur Plasmaphysik, 85748 Garching, Germany (Dated: November 2, 2018)We present a new analysis of the anisotropic spectral energy distribution in incompressible mag-netohydrodynamic (MHD) turbulence permeated by a strong mean magnetic field. The turbulentflow is generated by high-resolution pseudo-spectral direct numerical simulations with large-scaleisotropic forcing. Examining the radial energy distribution for various angles θ with respect to B reveals a specific structure which remains hidden when not taking axial symmetry with respect to B into account. For each direction, starting at the forced large-scales, the spectrum first exhibits anamplitude drop around a wavenumber k which marks the start of a scaling range and goes on up toa dissipative wavenumber k d ( θ ). The 3D spectrum for k ≥ k is described by a single θ -independentfunctional form F ( k/k d ), the scaling law being the same in every direction. The previous propertiesstill hold when increasing the mean field from B = 5 up to B = 10 b rms , as well as when passingfrom resistive to ideal flows. We conjecture that at fixed B the direction-independent scaling regimeis reached when increasing the Reynolds number above a threshold which raises with increasing B .Below that threshold critically balanced turbulence is expected. It is known that in presence of a mean magnetic fieldassumed here to point in the z -direction, B = B ˆ e z ,nonlinear interactions in incompressible magnetohydro-dynamics (MHD) are weakened in the field-parallel di-rection. The MHD approximation allows to describethe large-scale dynamics of astrophysical plasmas, i.e.ionized gases, like the interstellar medium or the solarcorona. Due to the above-mentioned anisotropy, thenonlinear energy transfer in MHD turbulence proceedspreferably to larger perpendicular spatial wavenumbers[1–4]. Direct numerical simulations (DNS) show thatthe field-perpendicular energy spectrum exhibits self-similar inertial-range scaling in wavenumber ∼ k − m with m = 5 / B , or m = 3 / B .Iroshnikov [11] and Kraichnan [12] proposed the firsttheory of the effect of a mean magnetic field on incom-pressible MHD turbulence. They remarked that any flowcan be decomposed into a sum of weakly interactingwaves with different wavevectors k , the term weak inter-action meaning that the characteristic time of deforma-tion of the waves is much longer than their periods. Thisled to the prediction of a slow cascade, with a spectralslope m = 3 / m = 5 / < ω > − = < ( k . B ) − > ≃ kB − (1)ignoring deliberately the waves with wave vectors per-pendicular (or almost perpendicular) to the mean field,for which the deformation time should clearly be smallerthan their period, and hence the interaction strong. This was criticized by Goldreich and Sridhar [13], who deniedthe possibility of the previous weak cascade to occur, andargued that the perpendicular strong cascade leading to k − / ⊥ should be the only one present. For a large enoughmean field, the perpendicular cascade should thus be re-stricted to a thin subset around the k ⊥ axis in Fourierspace, the subset becoming thinner when the mean fieldincreases.More precisely, the subdomain in ( k k , k ⊥ ) space wherethe perpendicular cascade is believed to occur is definedby a critical balance [13] between the characteristic timeof nonlinear interaction, τ NL ≃ ( k ⊥ u λ ) − , and the Alfv´entime τ A ≃ ( k k B ) − , the fluctuations becoming corre-lated along the guide field up to a distance ≃ B τ NL ,where u λ is the typical magnitude of fluctuations at thescale λ ≃ /k ⊥ . Assuming a scaling law in the perpendic-ular direction, and spectral transfer dominated by strongcoupling, i.e., χ = τ A /τ NL &
1, one obtains the 3D spec-trum: E ( k k , k ⊥ ) = k − m − q − ⊥ f ( χ ) (2)where χ = k − q k q ⊥ k − k b rms /B , with m = 5 / q =2 / f ( χ ) ≃ | χ | ≥ f ( χ ) negligible for | χ | ≪ B is significantly larger than the magneticfluctuation rms value. First, the spectral slope is ob-served to become m = 3 / m = 5 /
3; second, the time-scales ratio χ = τ A /τ NL becomes significantly smaller than unity[14] in the excited part of the spectrum, showing thatthe strict critical balance condition χ = 1 is too restric-tive to describe the anisotropy of the cascade, or, in otherwords, that the cascade is more extended in the obliquedirections than predicted by the critical balance condi-tion. This has led several authors to suggest modifica-tions which either still assume that the anisotropy is dic-tated by the critical balance condition [15, 16], or proposethat the time-scales ratio χ decreases with increasing B [17]. All these phenomenologies predict a spectral formdifferent from Eq. (2), with different spectral slopes inthe perpendicular and parallel directions.In the present paper, we analyze the angular spectrumand find by taking slices along the radial directions thata unique spectral form (and slope) holds in all directions,with the radial power-law range extent depending on theangle with the mean field. As a result, a significant por-tion of this spectrum lies in a domain where the time-scales ratio χ is sub-critical, that is, much smaller thanunity.This study focuses on representative states of fully-developed turbulence permeated by a strong mean mag-netic field with B = 5 b rms from high-resolution di-rect numerical simulations of quasi-stationary MHD tur-bulence forced at large scales. The forcing is realizedby freezing all modes (velocity and magnetic field) with k ≤ ω = ∇ × v and themagnetic field b are given by ∂ t ω = ∇ × [ v × ω − b × ( ∇ × b )] + µ ∆ ω ,∂ t b = ∇ × ( v × b ) + η ∆ b , ∇ · v = ∇ · b = 0 . The equations are solved by a standard pseudospectralmethod with spherical mode truncation to alleviate alias-ing errors. The numerical resolution is 1024 ×
256 col-locations points with reduced resolution in the directionof B [8] and with kinematic viscosity µ and resistiv-ity η set to µ = η = 9 · − . The analyzed datais the temporal average of five snapshots of the three-dimensional Fourier energy distribution taken equidis-tantly within about four to five field-perpendicular large-scale turnover times, T , ⊥ , of quasi-stationary turbulencewhere T ,γ = L /v γ rms = π/ h v γ i / R d k ′ δ ( k ′ γ ) | v γ ( k ′ ) | , γ ∈ { x, y, z } , cf. [18], and v rms ≈ b rms = 1 with T , ⊥ ≈ . T , k ≈ .
7. The normalized cross helicity ρ = h v · b i / ( h v i / h b i / ) and the Alfv´en ratio h v i / h b i fluctuate around 23% and 93%.The three-dimensional (3D) energy spec-trum, E ( k x , k y , k z ), relates to the total energy E tot = R d x ( v + b ) / R d kE ( k x , k y , k z ).Fig. 1 shows contour levels of E in two mutually FIG. 1: Energy contour levels of three-dimensional spectralenergy density E ( k x , k y , k z ) : (a)plane k z = 0, (b) plane k y = 0, (c) average of energy density over all planes containing B . (d) anisotropic scaling law between wave numbers k ⊥ and k k (see text), with two compensated scalings: k k ( k ⊥ ) k − ⊥ (diamonds) and k k ( k ⊥ ) k − / ⊥ (crosses). orthogonal planes containing the origin. The field-perpendicular k x - k y plane (Fig. 1,a) displays an isotropicenergy distribution, as expected. Anisotropy inducedby B appears in planes containing the B direction(Fig. 1,b).Spectral anisotropy is traditionally diagnosed by 1Dspectra, e.g., E ⊥ ( k x ) = R d k ′ E δ ( | k x | − k ′ x ) and E k ( k z ) = R d k ′ E δ ( | k z | − k ′ z ). These are shown com-pensated by k / (a) and by k / (b) in Fig. 2. The per-pendicular spectrum exhibits a power-law with m = 3 / k ⊥ , k k ) sharingthe same 1D energy density [14]. The anisotropy expo-nent q is seen to lie between q = 1 and q = 2 /
3, which isalso obtained in [14] for B = 5.While planar integration yields some informationabout anisotropy, it mixes all wavenumbers perpendic-ular to the chosen direction and thus blurs the separa-tion between inertial and dissipative scales if the dissi-pative scale is not constant over the planar domain ofintegration. More importantly, no information on in-termediate directions between parallel and perpendic-ular is available. Thus, spherical coordinates ( k, θ, φ )with respect to the mean field axis along ˆ e z are con-sidered. As E is isotropic in the azimuthal plane,the φ -dependence of E ( k, θ, φ ) is eliminated by aver-aging over φ ∈ [0 , π ] which strongly decreases sta-tistical noise and yields the φ -averaged 3D spectrum E ( k, θ ) = 1 / (2 π ) R dφE ( k, θ, φ ) whose isocontours in( k k , k ⊥ ) are shown in Fig. 1,c. We define the cor-responding one-dimensional (1D) spectrum E ( k, θ ) as E ( k, θ ) = k E ( k, θ ). The total energy is thus E tot =2 π R k dk R π/ E ( k, θ ) sin( θ ) dθ . FIG. 2: Plane-integrated one-dimensional perpendicular andparallel energy spectra compensated by (a) k / and (b) k / ,respectively. The symbol k stands for the respective field-perpendicular and field-parallel wavenumber The properties of E ( k, θ ) are shown in Fig. 3. A scal-ing range is seen starting at k ≃ k d ( θ ). The 1D scaling exponent mE ( k, θ ) = A ( θ ) k − m − (3)is independent of θ . The anisotropy appears at fixed k as a θ -dependence of the spectral amplitude A ( θ ) andof the dissipative wavenumber k d ( θ ). Normalizing thewavenumber by k d shows that the spectrum follows asingle functional form F whatever θ : E ( k, θ ) = F ( k/k d ) = F ( k/k d ( θ )) − m − (4)where F is a constant (the amplitude of the spectrum atthe dissipative scale). Note that the first equality holdsalso beyond the dissipative range, the second being validfor k ≤ k d . A corollary is that : A ( θ ) ∝ k d ( θ ) m +2 . (5)Fig. 3,a shows the self-similar wavenumber intervals ofall E ( k, θ ) spectra, starting in the range k ≃ θ . The dissipative wavenumber k d ( θ ) is estimated by locating the maximum of k E ( k ) ineach direction θ : varying θ from π/ k d ( θ ) from about 100 to 14 while the spectral energyat fixed k decreases in the inertial range from 1 to 10 − .Fig. 3,b indicates that the spectral θ -dependence can benearly eliminated by normalizing with k d (Eq. (4)). Therelation between spectral amplitude A ( θ ) and k d ( θ ) as given by Eq. (5) is confirmed by Fig. 3,c which shows that k / d (dotted curve) closely follows A ( θ ) = E ( k, θ ) k / (Eq. (3)), with k ≤
20 to eliminate the dissipative rangein the parallel direction.A simple model of anisotropic spectrum with a spec-tral exponent being the same in all directions has beenproposed previously in the context of shell models of tur-bulence [19] as well as solar wind turbulence [20, 21]. Itreads: A ( θ ) = ( cos θ/ε + sin θ ) − (1+ m/ (6)We tried to use this model to adjust the energy con-tours of our numerical simulations, and found that theglobal anisotropy between the perpendicular and parallelamplitude requires ε = 0 . sine and cosine as: A ( θ ) = ( cos θ/ε + sin θ ) − (1+ m/ (7)leads to a reasonable good fit to the simulation resultsin all directions, as seen both in the dotted-dashed curvein Fig. 3,c for the amplitude variation vs θ and in thedotted contours in the ( k k , k ⊥ ) plane in Fig. 3,d.The energy contours of the critical balance spectrum(Eq. 2, with B = 5 and k = 5) are also represented bydashed lines in Fig. 3,d. They isolate a small cone aboutthe k k axis in the whole plane, so excluding a large partof the angular structure of the true angular spectrum.Note that the dissipative wavenumber is determined upto an error of about a factor two (due to errors in interpo-lating the spectrum), which leads to the noisy appearanceof the curve of k d in Fig. 3,c, to the finite thickness of thenormalized spectra in Fig. 3,b and to variations of abouta factor 4 in the constant F in Eq. (4).To determine the scaling exponent m of E ( k, θ ) ∼ k − m , the 1D spectra averaged over four θ -intervals andcompensated by k / and k / are shown in Fig. 4. Thespectrum with θ = π/ m = 3 / m = 5 / / k d , and on theleft by an intermediate range which separates the inertialfrom the forcing range. The start of the inertial range isthus growing from k ≃ θ → π/ B = 5 √ FIG. 3: Details of spectral properties: (a) E ( k, θ ) for θ ranging from 0 to π/
2, (b) E ( k/k d , θ ) (c) k / d (dotted), k / E ( k, θ ) for 8 ≤ k ≤
20 (solid), Eq.(7) with m = 3 / E ( k k , k ⊥ ): sim-ulation data (solid); Eq. 2 (dashed line, the oblique line trac-ing the boundary x=1 with k = 5, B = 5); Eqs. 3,7 with β = 3, ε = 0 .
158 (dotted)FIG. 4: Spectra averaged in four subsets of directions: A:14 ≤ θ ≤ ; B: 42 ≤ θ ≤ ; C: 67 ≤ θ ≤ ; D:79 ≤ θ ≤ ; direction θ = 90 as thick curve. (a): spectracompensated by k / , the two oblique dashed lines indicatingroughly the inertial range; (b): spectra compensated by k / range decreasing substantially, so that the ratio of bothranges varies with B as k d ( π/ /k d (0) = ( A ( π/ /A (0)) m +2 ≃ B (8)Note that the fit by Eq. (7) remains as good as in Fig. 3,dafter decreasing ε by a factor √
2, as expected sinceEq. (7) implies k d ( π/ /k d (0) = 1 /ε hence from Eq. (8) ε ≃ /B . Increasing again B up to 10 confirmed thistrend, but the parallel range becomes too small in thatcase to allow a good determination of the associated dis-sipative wavenumber. This precludes checking that thespectrum normalized by k d is angle-independent for small θ . This difficulty could be alleviated by increasing thenumerical resolution which would allow increasing theReynolds number. We choose instead below to comparewith ideal MHD simulations with 512 resolution andwith B = 5 and 10. FIG. 5: Anisotropy in ideal runs; (a-b) B = 5 (c-d) B = 10.(a)-(c): k / d (dotted), k / E ( k, θ ) (solid) (a) for 6 ≤ k ≤ ≤ k ≤
8, Eq.(7) (dashed-dotted) with m = 5 / ε = 0 .
158 (c) ε = 0 . /
2. (b)-(d): energy contourlevels of E ( k k , k ⊥ ): simulation data (solid); Eq. 2 (dashedline, the oblique line tracing the boundary χ = 1 with k = 5, B = 5 (b) and B = 10 (d)); Eqs. 3,7 with β = 3, ε = 0 . ε = 0 . / In the ideal case, Fourier space is divided in a largescale range presenting spectral properties close to thoseof a standard turbulent spectrum with dissipation, anda small scale range where the spectral slope increases,the latter scales playing the role of a dissipative range(cf. [22] in the hydrodynamic case). The boundary be-tween the two domains slowly shifts with time to everlarger scales. It can be identified with the dissipativewavenumber, thus determined here as the minimum ofthe 1D spectrum. The simulations are initialized with aquasi-stationary state of the resistive run and are con-tinued without any dissipation and with the chosen B in the same numerical setup until the energetically risingsmall scales begin to pollute the scaling region. Choosingthe appropriate time, power-law ranges in all directionscan be properly identified even with a large field B = 10.The resulting power-law is found to be now m = 5 / B = 5 and 10, contrary to the resistive runs.This difference in scaling between the resistive and idealruns might be attributed to a different role played by thebottlenkeck effect [23] in these two setups. All reportedfindings, cf. Eqs.(3-8), are however confirmed when set-ting m = 5 /
3, as seen in Fig. 5 which shows (cf. Fig. 3,c)the anisotropy in two ideal runs with B = 5 (a) and B = 10 (b).Let us come back to the resistive case. As alreadyfound in [14], we find that the excited part of the k k , k ⊥ space is not restricted to regions where χ ≥
1. Theimportant point is that the 3D-energy contours, and aswell the boundary of the power-law range, ignore the iso-contours of χ as shown in Fig. 3 (see also Fig 5b,c): theform of the 3D-energy contours, as well as the angle-independent spectral slope, actually suggest an isotropiccascade, that is, a cascade along radial directions. Totest this idea, we define a θ -dependent effective Reynolds Re ∝ A ( θ ) / based on the energy density A ( θ ) k − m − at k . For θ = 0 → π/ Re increases by a factor 30,while k d grows by a factor 10, as k d ∝ Re α (9)where α = 2 / ( m + 2) ≃ / m being the 1D slope.The exponent in Eq. (9) is substantially smaller than thevalue α = 3 / m = 5 / / m =3 /
2) obtained by equating the input flux at k and thedissipative flux ǫ ≃ νk d u ( k d ) ). This means that if θ increases from zero to π/
2, the inertial range increasesmore slowly than it would if the dissipative loss at smallscales would balance the input energy rate at the ( k ≃ θ → π/ q = 2 / q = 2 / B /b rms = 5 or 10. Be-sides, it is remarkable that the fit by our anisotropy func-tion in Fig. 5 is as good when B = 5 b rms as when B = 10 b rms which is an indication that the spectralproperties of the turbulence are correctly revealed by us-ing our method.A possible way to reconcile both pictures is the fol-lowing. It takes into account the fact that there is asecond parameter, the Reynolds number. Indeed, asthe anisotropy increases with B , we find that the scal-ing range in the parallel direction decreases accordingly(Eq. (8)), which implies that, at a fixed viscosity, the par-allel power-law range disappears at even moderate B .In our case, k d ⊥ ≃
100 and k d k ≃
10 when B = 5 b rms ,while the start of the scaling range is k ≃ −
8, pre-venting to increase significantly B . The regime at high B will thus depend on the Reynolds number ( Re ). In-creasing B at fixed perpendicular Re depletes the field-parallel cascade and could ultimately lead to criticallybalanced turbulence. If B and Re are however largeenough, e.g. under astrophysical conditions, allowing forscaling in all directions then the properties describedin this work are expected to hold. We thus proposedirection-independent scaling for high Re and criticallybalanced turbulence at low Re . The crossover Re -value isexpected to increase with B . A strong indication in thissense is found in [14] where the anisotropy scaling expo-nent q is seen to be between 1 and 2 / B = 5 b rms ashere (Fig. 1,d), while it begins to cluster around 2 / B /b rms ≥
10. Such a good agreement with q = 2 / B /b rms ≃ A ( θ ) f ( k ) (ii) the anisotropy function A ( θ )is best expressed (Eq. 8) as the ratio of the perpendicu-lar over the parallel power-law range extent, which scaleslinearly with B in the B /b rms = 5 ,
10 interval consid-ered here. These results offer important tests for futuretheories of anisotropic turbulence.We thank G. Belmont, J. L´eorat, and A. Busse forseveral fruitful discussions. ∗ [email protected] † [email protected][1] H. R. Strauss, Phys. Fluids , 134 (1976). [2] D. Montgomery and L. Turner, Phys. Fluids , 825(1981).[3] J. V. Shebalin, W. H. Matthaeus, and D. Montgomery,Journal of Plasma Physics , 525 (1983).[4] R. Grappin, Physics of Fluids , 2433 (1986).[5] J. Cho and E. T. Vishniac, The Astrophysical Journal , 273 (2000).[6] J. Cho, A. Lazarian, and E. T. Vishniac, The Astrophys-ical Journal , 291 (2002).[7] J. Maron and P. Goldreich, The Astrophysical Journal , 1175 (2001).[8] W.-C. M¨uller and R. Grappin, Phys. Rev. Lett. ,114502 (2005).[9] P. D. Mininni and A. Pouquet, Phys. Rev. Lett. ,254502 (2007).[10] J. Mason, F. Cattaneo, and S. Boldyrev, Phys. Rev. E , 36403 (2008).[11] P. S. Iroshnikov, Astronomicheskii Zhurnal , 742(1963).[12] R. H. Kraichnan, Physics of Fluids , 1385 (1965).[13] P. Goldreich and S. Sridhar, Astrophysical Journal ,763 (1995).[14] B. Bigot, S. Galtier, and H. Politano, Phys. Rev. E , 66301 (2008).[15] S. Boldyrev, Phys. Rev. Lett. , 115002 (2006).[16] G. Gogoberidze, Physics of Plasmas , 2304 (2007).[17] S. Galtier, A. Pouquet, and A. Mangeney, Physics ofPlasmas , 2310 (2005), (c) 2005: American Instituteof Physics.[18] O. Zikanov and A. Thess, Applied Mathematical Mod-elling , 1 (2004).[19] V. Carbone and P. Veltri, Geophysical and AstrophysicalFluid Dynamics , 153 (1990).[20] J. A. Tessein, C. W. Smith, B. T. MacBride, W. H.Matthaeus, M. A. Forman, and J. E. Borovsky, The As-trophysical Journal , 684 (2009).[21] V. Carbone, F. Malara, and P. Veltri, Journal of Geo-physical Research , 1763 (1995).[22] C. Cichowlas, P. Bona¨ıti, F. Debbasch, and M. Brachet,Phys. Rev. Lett. , 264502 (2005).[23] A. Beresnyak and A. Lazarian, arXiv astro-ph.GA (2010), 1002.2428v2.[24] T. S. Horbury, M. Forman, and S. Oughton, Phys. Rev.Lett.101