Reynolds-number dependence of the dimensionless dissipation rate in homogeneous magnetohydrodynamic turbulence
RReynolds-number dependence of the dimensionless dissipation rate in homogeneousmagnetohydrodynamic turbulence
Moritz Linkmann,
1, 2, ∗ Arjun Berera, and Erin E. Goldstraw
3, 2 Department of Physics & INFN, University of Rome Tor Vergata,Via della Ricerca Scientifica 1, 00133 Rome, Italy SUPA, School of Physics and Astronomy, University of Edinburgh, Peter Guthrie Tait Road, EH9 3FD, UK School of Mathematics and Statistics, University of St. Andrews, KY16 9SS, UK (Dated: September 4, 2018)This paper examines the behavior of the dimensionless dissipation rate C ε for stationary and non-stationary magnetohydrodynamic (MHD) turbulence in presence of external forces. By combiningwith previous studies for freely decaying MHD turbulence, we obtain here both the most generalmodel equation for C ε applicable to homogeneous MHD turbulence and a comprehensive numeri-cal study of the Reynolds number dependence of the dimensionless total energy dissipation rate atunity magnetic Prandtl number. We carry out a series of medium to high resolution direct numericalsimulations of mechanically forced stationary MHD turbulence in order to verify the predictions ofthe model equation for the stationary case. Furthermore, questions of nonuniversality are discussedin terms of the effect of external forces as well as the level of cross- and magnetic helicity. Themeasured values of the asymptote C ε, ∞ lie between 0 . (cid:54) C ε, ∞ (cid:54) .
268 for free decay, wherethe value depends on the initial level of cross- and magnetic helicities. In the stationary case wemeasure C ε, ∞ = 0 . I. INTRODUCTION
The dynamics of conducting fluids is relevant to many areas in geo- and astrophysics as well as in engineering andindustrial applications. Often the flow is turbulent, and the interaction of the turbulent flow with the magnetic fieldleads to considerable complexity. Being a multi-parameter problem, techniques that have been successfully appliedto turbulence in nonconducting fluids sometimes fail to deliver unambiguous predictions in magnetohydrodynamic(MHD) turbulence. This concerns e.g. the prediction of inertial range scaling exponents by extension of Kolmogorov’sarguments [1] to MHD, and considerable effort has been put into the further understanding of inertial range cascade(s)in MHD turbulence [2–9]. The difficulties are partly due to the many different configurations that can arise in MHDturbulence because of e.g. anisotropy, different levels of vector field correlations, different values of the dissipationcoefficients and different types of external forces, and as such are connected to the question of universality in MHDturbulence [10–21]. The behavior of the (dimensionless) dissipation rate is representative of this problem, in the sensethat the aforementioned properties of MHD turbulence influence the energy transfer across the scales, i.e. the cascadedynamics [11, 22–26], and thus the amount of energy that is eventually dissipated at the small scales.The behavior of the total dissipation rate in a turbulent non-conducting fluid is a well-studied problem. As such ithas been known for a long time that the total dissipation rate in both stationary and freely decaying homogeneousisotropic turbulence tends to a constant value with increasing Reynolds number following a well-known characteristiccurve [27–32]. For statistically steady isotropic turbulence this curve can be approximated by the real-space stationaryenergy balance equation, where the asymptote is connected to the maximal inertial flux of kinetic energy [30]. Thecorresponding problem in MHD has received much less attention, however, recent numerical results for freely decayingMHD turbulence at unity magnetic Prandtl number report similar behavior. Mininni and Pouquet [33] carried outdirect numerical simulations (DNSs) of freely decaying homogeneous MHD turbulence without a mean magnetic field,showing that the temporal maximum of the total dissipation rate ε ( t ) became independent of Reynolds number at aTaylor-scale Reynolds number R λ (measured at the peak of ε ( t )) of about 200. Dallas and Alexakis [34] measured thedimensionless dissipation rate C ε also from DNS data for free decay for random initial fields with strong correlationsbetween the velocity field and the current density. Again, it was found that C ε → const. with increasing Reynoldsnumber. Interestingly, a comparison with the data of Ref. [33] showed that the approach to the asymptote wasslower than for the data of Ref. [33], suggesting an influence of the level of certain vector field correlations on theapproach to the asymptote. A theoretical model for dissipation rate scaling in freely decaying MHD turbulencewas put forward recently [35] based on the von K´arm´an-Howarth energy balance equations (vKHE) in terms of ∗ [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] J a n Els¨asser fields [36]. For unity magnetic Prandtl number it predicts the dependence of C ε on a generalized Reynoldsnumber R − ≡ z − L + / ( ν + µ ), with z − denoting the root-mean-square value of one Els¨asser field, L + the integralscale corresponding to the other Els¨asser field, while ν and µ are the kinematic viscosity and the magnetic resistivity,respectively. The model equation has the following form C ε = C ε, ∞ + CR − + DR − + O ( R − − ) , (1)where C and D are time-dependent coefficients depending on several parameters, which themselves depend on themagnetic, cross- and kinetic helicities. The predictions of this equation were subsequently tested against data obtainedfrom medium to high resolution DNSs of freely decaying homogeneous MHD turbulence leading to a very goodagreement between theory and data.In summary, there is compelling numerical and theoretical evidence for finite dissipation in freely decaying MHDturbulence at least for unity magnetic Prandtl number P m = ν/µ , while so far no systematic results for the stationarycase have been reported. In this paper we extend the derivation carried out in Ref. [35] to include the effects of externalforces and we present the first systematic study of dissipation rate scaling for stationary MHD turbulence. In orderto be able to test the model equation against DNS data for a large range of generalized Reynolds numbers, weconcentrate on the case P m = 1. The most general form of Eq. (1) for nonstationary flows with large-scale externalforcing is derived, which can be applied to freely decaying and stationary flows by setting the corresponding termsto zero. This generalization of Eq. (1) is the first main result of the paper, it is applicable to both freely decayingand stationary MHD turbulence. It implies that the dissipation rate of total energy is finite in the limit R − → ∞ inanalogy to hydrodynamics, and highlights the dependence of the coefficients C and D on the external forces. As such,Eq. (1) predicts nonuniversal values of the asymptotic value C ε, ∞ of the dimensionless dissipation rate in the infiniteReynolds number limit and of the approach to the asymptote for a variety of MHD flows. The resulting theoreticalpredictions for the stationary case are compared to DNS data for stationary MHD turbulence for three different typesof mechanical forcing while the results for the freely decaying case [35] are reviewed for completeness. The DNS datashows good agreement with Eq. (1) and the different forcing schemes have no measurable effect on the values of thecoefficients in Eq. (1). The measured values of C ε, ∞ lie between 0 . (cid:54) C ε, ∞ (cid:54) .
268 for free decay, where the valuedepends on the initial level of cross- and magnetic helicities. In the stationary case we measure C ε, ∞ = 0 . II. THE TOTAL DISSIPATION IN TERMS OF ELS ¨ASSER FIELDS
In this paper we consider statistically homogeneous MHD turbulence in the absence of a background magnetic field.The flow is taken to be incompressible, leading to the following set of coupled partial differential equations ∂ t u = − ρ ∇ P − ( u · ∇ ) u + 1 ρ ( ∇ × b ) × b + ν ∆ u + f u , (2) ∂ t b = ( b · ∇ ) u − ( u · ∇ ) b + µ ∆ b + f b , (3) ∇ · u = 0 and ∇ · b = 0 , (4)where u denotes the velocity field, b the magnetic induction expressed in Alfv´en units, ν the kinematic viscosity, µ themagnetic resistivity, P the thermodynamic pressure, f u and f b are external mechanical and electromagnetic forces,which may be present, and ρ denotes the density which is set to unity for convenience. Equations (2)-(4) are consideredon a three-dimensional domain Ω, which due to homogeneity can either be the full space R or a subdomain [0 , L box ) with periodic boundary conditions. The MHD equations (2)-(4) can be formulated more symmetrically using Els¨asservariables z ± = u ± b [37] ∂ t z ± = − ρ ∇ ˜ P − ( z ∓ · ∇ ) z ± + ( ν + µ )∆ z ± + ( ν − µ )∆ z ∓ + f ± , (5) ∇ · z ± = 0 , (6)where f ± = f u ± f b and the pressure ˜ P consists of the sum of the thermodynamic pressure P and the magneticpressure ρ | b | /
2. Which formulation of the MHD equations is chosen often depends on the physical problem, for someproblems the Els¨asser formalism is technically convenient, while the formulation using the primary fields u and b facilitates physical understanding. The ideal invariants total energy E ( t ), cross-helicity H c ( t ) and magnetic helicity H m ( t ) are given in the respective formulations of the MHD equation by E ( t ) = 12 (cid:90) Ω d k (cid:104)| ˆ u ( k , t ) | + | ˆ b ( k , t ) | (cid:105) = 14 (cid:90) Ω d k (cid:104)| ˆ z + ( k , t ) | + | ˆ z − ( k , t ) | (cid:105) , (7) H c ( t ) = (cid:90) Ω d k (cid:104) ˆ u ( k , t ) · ˆ b ( − k , t ) (cid:105) = 14 (cid:90) Ω d k (cid:104)| ˆ z + ( k , t ) | − | ˆ z − ( k , t ) | (cid:105) , (8) H m ( t ) = (cid:90) Ω d k (cid:104) ˆ a ( k , t ) · ˆ b ( − k , t ) (cid:105) = 14 (cid:90) Ω d k (cid:28)(cid:20) i k k × ( ˆ z + ( k , t ) − ˆ z − ( k , t )) (cid:21) · ( ˆ z + ( − k , t ) − ˆ z − ( − k , t )) (cid:29) , (9)with ˆ b , ˆ u and ˆ z ± denoting the respective Fourier transforms of the magnetic, velocity and Els¨asser fields, whileˆ a is the Fourier transform of the magnetic vector potential a . The angled brackets indicate an ensemble average.Equation (9) is gauge-independent as shown in Appendix A.We now motivate the use of the Els¨asser formulation for the study of the dimensionless dissipation coefficient inMHD. In hydrodynamics, the dimensionless dissipation coefficient C ε,u is defined in terms of the Taylor surrogateexpression for the total dissipation rate, U /L u , where U denotes the root-mean-square (rms) value of the velocityfield and L u the integral scale defined with respect to the velocity field, as C ε,u ≡ ε kin L u U . (10)However, in MHD there are several quantities that may be used to define an MHD analogue to the Taylor surrogateexpression, such as the rms value B of the magnetic field, one of the different length scales defined with respect toeither b or u , or the total energy.Since the total dissipation in MHD turbulence should be related to the flux of total energy through different scales,one may think of defining a dimensionless dissipation coefficient for MHD in terms of the total energy. However, thiswould lead to a nondimensionalization of the hydrodynamic transfer term u · ( u · ∇ ) u with a magnetic quantity. Thiscan be seen by considering the analog of the von K´arm´an-Howarth energy balance equation in real space [38] statedhere for the case of free decay − d t E ( t ) = ε ( t ) = − ∂ t ( B uuLL ( r, t ) + B bbLL ( r, t )) + 32 r ∂ r (cid:18) r B uuuLLL ( r, t ) + r C bbuLLL ( r, t ) (cid:19) + 6 r C bub ( r, t ) + 1 r ∂ r (cid:0) r ∂ r ( νB uuLL ( r, t ) + µB bbLL ( r, t )) (cid:1) , (11)where B uuLL , B bbLL and B uuuLLL are the longitudinal structure functions, C bbuLLL the longitudinal correlation function and C bub another correlation function. The longitudinal structure and correlation functions are given by B uuLL ( r, t ) = (cid:104) ( δu L ( r , t )) (cid:105) , (12) B bbLL ( r, t ) = (cid:104) ( δb L ( r , t )) (cid:105) , (13) B uuuLLL ( r, t ) = (cid:104) ( δu L ( r , t )) (cid:105) , (14) C bbuLLL ( r, t ) = (cid:104) u L ( x , t ) b L ( x , t ) b L ( x + r , t ) (cid:105) , (15)where r = | r | and v L = v · r /r denotes the longitudinal component of a vector field v , that is its component parallelto the displacement vector r , and δv L ( r ) = [ v ( x + r ) − v ( x )] · r r , (16)its longitudinal increment. The function C bub is defined through the third-order correlation tensor C bubij,k ( r , t ) = (cid:104) ( u i ( x ) b j ( x ) − b i ( x ) u j ( x )) b k ( x + r ) (cid:105) = C bub ( r, t ) (cid:16) r j r δ ik − r i r δ jk (cid:17) . (17)As can be seen from their respective definitions, the functions C bbuLLL and C bub scale with B U while the function B uuuLLL scales with U . If Eq. (11) were to be nondimensionalized with respect to the total energy then the purelyhydrodynamic term B uuuLLL would be scaled partially by a magnetic quantity.This problem of inconsistent nondimensionalization can be avoided by working with Els¨asser fields, which requiresan expression for the total dissipation rate ε ( t ) in terms of Els¨asser fields. The total rate of energy dissipation inMHD turbulence is given by the sum of Ohmic and viscous dissipation ε ( t ) = ε mag ( t ) + ε kin ( t ) , (18)where ε mag ( t ) = µ (cid:90) Ω d k k (cid:104)| ˆ b ( k , t ) | (cid:105) , (19) ε kin ( t ) = ν (cid:90) Ω d k k (cid:104)| ˆ u ( k , t ) | (cid:105) . (20)Similarly, the total dissipation rate can be decomposed into its respective contributions from the Els¨asser dissipationrates ε ( t ) = 12 (cid:0) ε + ( t ) + ε − ( t ) (cid:1) , (21)where the Els¨asser dissipation rates are defined as ε ± ( t ) = ν + (cid:90) Ω d k k (cid:104)| ˆ z ± ( k , t ) | (cid:105) + ν − (cid:90) Ω d k k (cid:104) ˆ z ± ( k , t ) · ˆ z ∓ ( − k , t ) (cid:105) , (22)with ν ± = ( ν ± µ ). The total dissipation rate relates to the sum of the Els¨asser dissipation rates ε + ( t ) + ε − ( t ) = ε ( t ) + ε H c ( t ) + ε ( t ) − ε H c ( t ) = 2 ε ( t ) , (23)where the cross-helicity dissipation rate ε H c is given by ε H c ( t ) = 12 (cid:0) ε + ( t ) − ε − ( t ) (cid:1) . (24)Since this paper is concerned with both stationary and nonstationary flows, the total energy input rate ι must also beconsidered. Similar to the dissipation rate, the input rate can be split up into either kinetic and magnetic contributionsor the Els¨asser contributions ι ± ( t ) ι ( t ) = ι mag ( t ) + ι kin ( t ) (25) ι ( t ) = 12 (cid:0) ι + ( t ) + ι − ( t ) (cid:1) . (26)The latter equation can be rewritten as ι + ( t ) = ι ( t ) + 12 (cid:0) ι + ( t ) − ι − ( t ) (cid:1) = ι ( t ) + ι H c ( t ) , (27)where ι H c denotes the input rate of the cross-helicity. III. DERIVATION OF THE EQUATION
Since the total dissipation rate can be expressed either in terms of the Els¨asser fields or the primary fields u and b ,it should be possible to describe it also by the vKHE for z ± [36]. For the freely decaying case no further complicationarises as the rate of change of total energy, which figures on the left-hand side of the energy balance, equals the totaldissipation rate. However, in the more general case the rate of change of the total energy is given by the difference ofenergy input and dissipation. That is, in the most general case the total energy dissipation rate is given by ε ( t ) = ι ( t ) − d t E ( t ) . (28)For the stationary case d t E ( t ) = 0 and one obtains ε ( t ) = ι ( t ). For the freely decaying case ι ( t ) = 0 and the change intotal energy is due to dissipation only, that is − d t E ( t ) = ε ( t ). In terms of Els¨asser variables ε ( t ) can also be expressedas ε ( t ) = ι ( t ) − d t E ( t ) = ι ( t ) − d t E ± ( t ) ∓ d t H c ( t ) , (29)where E ± ( t ) denote the Els¨asser energies. Since we have related the total dissipation rate to the rate of change of theEls¨asser energies, we are now in a position to consider the energy balance equations for z ± , which are stated here forthe most general case of homogeneous forced nonstationary MHD flows without a mean magnetic field − ∂ t E ± ( t ) + I ± ( r, t ) = − ∂ t B ±± LL ( r, t ) − ∂ r r (cid:18) r C ±∓± LL,L ( r, t ) (cid:19) + 34 r ∂ r (cid:0) r ∂ r ( ν + µ ) B ± LL ( r, t ) (cid:1) + 34 r ∂ r (cid:0) r ∂ r ( ν − µ ) B ∓ LL ( r, t ) (cid:1) , (30)where I ± ( r, t ) are (scale-dependent) energy input terms and C ±∓∓ LL,L ( r, t ) = (cid:104) z ± L ( x , t ) z ∓ L ( x , t ) z ± L ( x + r , t ) (cid:105) , (31) B ±± LL ( r, t ) = (cid:104) ( δz ± L ( r , t )) (cid:105) , (32) B ±∓ LL ( r, t ) = (cid:104) δz ± L ( r , t ) δz ∓ L ( r , t ) (cid:105) , (33)are the third-order longitudinal correlation function and the second-order structure functions of the Els¨asser fields,respectively. As can be seen from the definition, the third-order correlation function scales with ( z ± ) z ∓ , where z ± denote the respective rms values of the Els¨asser fields. This permits a consistent nondimensionalization of the Els¨asservKHE using the appropriate quantities defined in terms of Els¨asser variables. As such the complication that aroseif the energy balance was written in terms of b and u can be circumvented. This motivates the definition of thedimensionless Els¨asser dissipation rates as C ± ε ( t ) ≡ ε ( t ) L ± ( t ) z ± ( t ) z ∓ ( t ) , (34)where L ± ( t ) = 3 π E ± ( t ) (cid:90) Ω d k k − (cid:104)| z ± ( k , t ) | (cid:105) , (35)are the integral scales defined with respect to z ± [39] . For balanced MHD turbulence, i.e. H c = 0, one should expect C + ε ( t ) = C − ε ( t ), since E ± ( t ) = 2 E ( t ) ± H c ( t ) = 2 E ( t ) . (36)Therefore all quantities defined with respect to the rms fields z + and z − should be the same in this case. Finally, thedimensionless dissipation rate C ε ( t ) is defined as C ε ( t ) = C + ε ( t ) + C − ε ( t ) ≡ ε ( t ) L + ( t ) z + ( t ) z − ( t ) + ε ( t ) L − ( t ) z − ( t ) z + ( t ) . (37)Using the definition given in Eq. (34), the Els¨asser energy balance equations (30) can now be consistently nondi-mensionalized. For conciseness the explicit time and spatial dependences are from now on omitted, unless there is aparticular point to make. A. Dimensionless von K´arm´an-Howarth equations
By introducing the nondimensional variables σ ± = r/L ± [12] and non-dimensionalising Eq. (30) as proposed in thedefinitions of C ± ε given in Eq. (34) one obtains − (cid:0) d t E ± − I ± (cid:1) L ± z ± z ∓ = − σ ± ∂ σ ± (cid:32) σ ± C ±∓± LL,L z ± z ∓ (cid:33) − L z ± z ± z ∓ ∂ t B ±± LL µ + νL ± z ∓ σ ± (cid:18) σ ± ∂ σ ± B ±± LL z ± (cid:19) + ν − µL ± z ± σ ± (cid:18) σ ± ∂ σ ± B ±∓ LL z ± z ∓ (cid:19) . (38)Before proceeding further, the scale-dependent forcing term on the left-hand side of this equation needs to be analyzedin some detail in order to clarify its relation to the energy input rates ι and ι ± . The Els¨asser energy input I ± is givenby I ± ( r ) = 3 r (cid:90) r dr (cid:48) r (cid:48) (cid:104) z ± ( x + r (cid:48) ) · f ± ( x ) (cid:105) . (39)Since the energy input rate is given by ι ± = (cid:104) z ± ( x ) · f ± ( x ) (cid:105) , the correlation function can be expressed as (cid:104) z ± ( x + r ) · f ± ( x ) (cid:105) = ι ± φ ± ( r/L f ) , (40)where φ ± are dimensionless even functions of r/L f satisfying φ ± (0) = 1 and L f the characteristic scale of the forcing.At scales much smaller than the forcing scale, i.e. for r/L f <<
1, for suitable types of forces φ ± ( r/L f ) can beexpanded in a Taylor series [40], leading to the following expression for the energy input I ± ( r ) = 3 r (cid:90) r dr (cid:48) r (cid:48) ι ± (cid:34) (cid:18) rL f (cid:19) ∂ φ ± ∂ ( r/L f ) (cid:12)(cid:12)(cid:12) r/L f =0 + O (cid:32)(cid:18) rL f (cid:19) (cid:33)(cid:35) . (41)In the limit of infinite Reynolds number the inertial range extends through all wavenumbers, formally implying that L f → ∞ , where Eq. (41) implies I ± ( r ) → ι ± . Therefore it should be possible to split the term I ± ( r ) into a constant, ι ± , and a scale-dependent term J ± ( r ), which encodes the additional scale dependence introduced by realistic, finiteReynolds number forcing. For consistency, this scale-dependent term must vanish in the formal limit Re → ∞ . Thiscan be achieved by writing I ± ( r ) in terms of the correlation of the force and Els¨asser field increments I ± ( r ) = ι ± − r (cid:90) r dr (cid:48) r (cid:48) (cid:104) δ z ± · δ f ± (cid:105) . (42)Therefore we define J ± ( r ) = − r (cid:90) r dr (cid:48) r (cid:48) (cid:104) δ z ± · δ f ± (cid:105) , (43)where lim Re →∞ J ± ( r ) = 0. Hence the energy input I ± ( r ) can be expressed as the sum of the scale-independentenergy input rate ι ± and a scale-dependent term which vanishes in the formal limit Re → ∞ I ± ( r ) = ι ± + J ± ( r ) , (44)with lim Re →∞ J ± ( r ) = 0. Substitution of Eq. (44) into the nondimensionalized energy balance Eq. (38) leads tothe dimensionless version of the Els¨asser vKHE for homogeneous MHD turbulence in the most general case fornonstationary flows at any magnetic Prandtl number C ± ε = − ∂ σ ± σ ± (cid:32) σ ± C ±∓± LL,L z ± z ∓ (cid:33) + L ± z ± z ∓ (cid:18) ± d t H c − ∂ t B ±± LL − J ± ∓ ι H c (cid:19) + 1 R ∓ ∂ σ ± σ ± (cid:18) σ ± ∂ σ ± B ±± LL z ± (cid:19) + 1 R (cid:48)± ∂ σ ± σ ± (cid:18) σ ± ∂ σ ± B ±∓ LL z ± z ∓ (cid:19) , (45)where R ∓ and R (cid:48)± denote generalized large-scale Reynolds numbers given by R ∓ = z ∓ L ± / ( ν + µ ) and R (cid:48)± = z ± L ± / ( ν − µ ) . (46)In order to express Eq. (45) more concisely, the following dimensionless functions are defined g ±∓± = C ±∓± LL,L z ± z ∓ , (47) h ±± = B ±± LL z ± , (48) h ±∓ = B ±∓ LL z ± z ∓ , (49) H ±± = L ± z ± z ∓ ∂ t B ±± LL , (50) F ± = L ± z ± z ∓ J ± , (51) G ± = L ± z ± z ∓ d t H c , (52) Q ± = L ± z ± z ∓ ι H c , (53)such that Eq. (45) can be written as C ± ε = − ∂ σ ± σ ± (cid:18) σ ± g ±∓± (cid:19) ± G ± − H ±± − F ± ∓ Q ± + 3 R ∓ ∂ σ ± σ ± (cid:0) σ ± ∂ σ ± h ±± (cid:1) + 3 R (cid:48)∓ ∂ σ ± σ ± (cid:0) σ ± ∂ σ ± h ±∓ (cid:1) . (54)This equation can be applied to the two simpler cases of freely decaying and stationary MHD turbulence by settingthe corresponding terms to zero. For the case of free decay there are no external forces therefore F ± = 0, while forthe stationary case the terms G ± and H ± vanish. A further simplification concerns the case P m = 1, that is ν = µ ,where the inverse of the generalized Reynolds numbers R (cid:48)± vanish. In this case the evolution of C ± ε depends only on R ∓ , and an approximate analysis using asymptotic series is possible. Most numerical results are concerned with thiscase due to computational constraints, hence it would be very difficult to test an approximate equation against DNSdata if not only Re but also P m needs to be varied. From now on the magnetic Prandtl number is therefore set tounity, keeping in mind that the analysis could be extended to
P m (cid:54) = 1 provided the approximate equation derived inthe following section is consistent with DNS data.
B. Asymptotic analysis for the case
P m = 1
Equation (54) suggests a dependence of C ± ε on 1 /R ∓ , however, the structure and correlation functions also havea dependence on Reynolds number, which describes their deviation from their respective inertial-range forms. Thehighest derivative in Eq. (54) is multiplied by the small parameter 1 /R ∓ , which suggests that this equation may beviewed as singular perturbation problem amenable to asymptotic analysis [41]. The Els¨asser vKHE was rescaled bythe rms values of the Els¨asser fields and the corresponding integral length scales, where the integral scales are bydefinition the large-scale quantities, the interpretation in hydrodynamics usually being that they represent the size ofthe largest eddies. As such, the nondimensionalization was carried out with respect to quantities describing the largescales, that is, with respect to ‘outer’ variables. Hence outer asymptotic expansions of the nondimensional structureand correlation functions are considered with respect to the inverse of the (large-scale) generalized Reynolds numbers1 /R ∓ . We point out that the case P m (cid:54) = 1 would require expansions in two parameters, where the cases
P m >
P m < R (cid:48)± between the two cases.The formal asymptotic series of a generic function f [used for conciseness in place of the functions on the right-handside of Eq. (54)] up to second order in 1 /R ∓ reads f = f + 1 R ∓ f + 1 R ∓ f + O ( R − ∓ ) . (55)After substitution of the expansions into Eq. (54), collecting terms of the same order in 1 /R ∓ , one arrives at equationsdescribing the behavior of C + ε and C − ε C ± ε = C ± ε, ∞ + C ± R ∓ + D ± R ∓ + O ( R − ∓ ) , (56)up to second order in 1 /R ∓ , using the coefficients C ± ε, ∞ , C ± and D ± defined as C ± ε, ∞ = − ∂ σ ± σ ± (cid:18) σ ± g ±∓± (cid:19) ± G ± − H ±± ∓ Q ± , (57) C ± = 3 ∂ σ ± σ ± (cid:20) σ ± (cid:18) ∂ σ ± h ±± − g ±∓± (cid:19)(cid:21) ∓ F ± − H ±± , (58) D ± = 3 ∂ σ ± σ ± (cid:20) σ ± (cid:18) ∂ σ ± h ±± − g ±∓± (cid:19)(cid:21) ∓ F ± − H ±± , (59)in order to write Eq. (54) in a more concise way. The zero-order term in the expansion of the function F ± vanishes,since F ± corresponds to the scale-dependent part J ± of the energy input which vanishes in the limit R ∓ → ∞ , hence F ± = 0. According to the definition of C ε in Eq. (37), the asymptote C ε, ∞ is given by C ε, ∞ = C + ε, ∞ + C − ε, ∞ , (60)and using the definition of the generalized Reynolds numbers, which implies R + = ( L − /L + )( z + /z − ) R − one can define C = C + + L − L + z + z − C − , (61)( D is defined analogously), resulting in the following expression for the dimensionless dissipation rate C ε = C ε, ∞ + CR − + DR − + O ( R − − ) . (62)Since the time dependence of the various quantities in this problem has been suppressed for conciseness, it has to beemphasized that Eq. (62) is time dependent, including the Reynolds number R − . Equation (62) in conjunction witheqs. (57)-(59) is the most general asymptotic expression for the Reynolds number dependence of C ε developed so far.It is applicable for freely decaying, stationary and non-stationary MHD turbulence in the presence of external forces,and it may be applied to the corresponding problem in non-conducting fluids by setting b = 0. As such it extendsprevious results for freely decaying MHD turbulence [35], as well as for the stationary case in homogeneous isotropicturbulence of non-conducting fluids [30].For nonstationary MHD turbulence at the peak of dissipation the term H ±± in Eq. (57) vanishes for constant fluxof cross-helicity (that is, d t H c = 0), since in the infinite Reynolds number limit the second-order structure functionwill have its inertial range form at all scales. By self-similarity the spatial and temporal dependences of e.g. B ++ LL should be separable in the inertial range, that is B ++ LL ( r, t ) ∼ ( ε + ( t ) r ) α , (63)for some value α , and ∂ t B ++ LL ∼ αε + ( t ) α − d t ε + r α . (64)At the peak of dissipation d t ε + | t peak = d t ε | t peak − d t H c = d t ε | t peak = 0 , (65)which implies H ++0 ( t peak ) = 0. Equation (57) taken for nonstationary flows at the peak of dissipation is thus identicalto Eq. (57) for stationary flows, which suggests that at this point in time a nonstationary flow may behave similarlyto a stationary flow. We will come back to this point in Sec. IV. Due to selective decay, that is the faster decay ofthe total energy compared to H c and H m [25], in most situations one could expect d t H c to be small compared to ε in the infinite Reynolds number limit. In this case G ± (cid:39) C ± ε, ∞ ( t peak ) = − ∂ σ ± σ ± (cid:18) σ ± g ±∓± (cid:19) , (66)which recovers the inertial-range scaling results of Ref. [36] and reduces to Kolmogorov’s four-fifth law for b = 0. C. Relation of C ε, ∞ to energy and cross-helicity fluxes In analogy to hydrodynamics, the asymptotes C ± ε, ∞ should describe the total energy flux, that is the contributionof the cross-helicity flux to the Els¨asser flux should be canceled by the respective terms G ± and Q ± in Eq. (57).However, since this is not immediately obvious from the derivation, further details are given here. For nonstationaryturbulence at the peak of dissipation, Eq. (57) for the asymptotes C ± ε, ∞ reduces to C ± ε, ∞ = − ∂ σ ± σ ± (cid:18) σ ± g ±∓± (cid:19) ± G ± ∓ Q ± . (67)The dimensional version of this equation is ε = − ∂ r r (cid:18) r C ±∓± LL,L (cid:19) ± d t H c ∓ ι H c , (68)where it is assumed that the function C ±∓± LL,L has its inertial range form corresponding to g ±∓± . The function C ±∓± LL,L can also be expressed through the Els¨asser increments [36] C ±∓± LL,L = 14 (cid:0) (cid:104) ( δz ± L ( r )) δz ∓ L ( r ) (cid:105) − (cid:104) z ± L ( x ) z ± L ( x ) z ∓ L ( x + r ) (cid:105) (cid:1) , (69)which can be written in terms of the primary fields u and b as C ±∓± LL,L = 14 23 (cid:104) ( δu L ( r )) − b L ( x ) u L ( x + r ) (cid:105)∓
14 23 (cid:104) ( δb L ( r )) − u L ( x ) b L ( x + r ) (cid:105) , (70)(see e.g. Ref. [36]). The two terms on the first line of Eq. (70) are the flux terms in the evolution equation of the totalenergy, while the two terms on last line correspond to the flux terms in the evolution equation of the cross-helicity[36]. Now Eq. (68) can be expressed in terms of the primary fields ε = − ∂ r r (cid:18) r C ±∓± LL,L (cid:19) ± d t H c ∓ ι H c = − ∂ r r (cid:18) r (cid:104) ( δu L ( r )) − b L ( x ) u L ( x + r ) (cid:105) (cid:19) ± ∂ r r (cid:18) r (cid:104) ( δb L ( r )) − u L ( x ) b L ( x + r ) (cid:105) (cid:19) ± d t H c ∓ ι H c = ε T ± ε H c ± d t H c ∓ ι H c = ε T , (71)where ε T is the flux of total energy and ε H c the cross-helicity flux, which must equal − d t H c + ι H c for nonstationaryMHD turbulence. Thus the contribution from the third-order correlator C ±∓± LL,L resulting in ε H c is canceled by d t H c − ι H c , or, after nondimensionalization, the cross-helicity flux ε H c L ± / [( z ± ) z ∓ ] is canceled by G ± − Q ± . The two simplercases of freely decaying and stationary MHD turbulence are recovered by setting either Q ± = 0 (free decay) or G ± = 0(stationary case). D. Nonuniversality
Since C ε, ∞ is a measure of the flux of total energy across different scales in the inertial range, differences for thevalue of this asymptote should be expected for systems with different initial values for the ideal invariants H m and H c . The flux of total energy and thus the asymptote C ε, ∞ is an averaged quantity. This implies that cancellationsbetween forward and inverse fluxes may take place leading on average to a positive value of the flux, that is, forwardtransfer from the large scales to the small scales. In case of H m (cid:54) = 0, the value of C ε, ∞ should therefore be less thanfor H m = 0 due to a more pronounced inverse energy transfer in the helical case, the result of which is less averageforward transfer and thus a smaller value of the (average) flux of total energy. For H c (cid:54) = 0 the asymptote C ε, ∞ isexpected to be smaller than for H c = 0, since alignment of u and b weakens the coupling of the two fields in the0induction equation, leading to less transfer of magnetic energy across different scales and presumably also less transferof kinetic to magnetic energy.Furthermore, from an analysis of helical triadic interactions in ideal MHD carried out in Ref. [42] it may be expectedthat high values of cross-helicity have a different effect on the asymptote C ε, ∞ , depending on the level of magnetichelicity. The analytical results suggested that the cross-helicity may have an asymmetric effect on the nonlineartransfers in the sense that the self-ordering inverse triadic transfers are less quenched by high levels of H c comparedto the forward transfers. The triads contributing to inverse transfers were mainly those where magnetic field modesof like-signed helicity interact, and so for simulations with maximal initial magnetic helicity the dynamics will bedominated by these triads. If the inverse fluxes are less affected by the cross-helicity than the forward fluxes, thenthe expectation is that for a comparison of the value of C ε, ∞ between systems with (i) high H m and H c , (ii) high H m and H c = 0, (iii) H m = 0 and high H c and finally (iv) H m = 0 and H c = 0, the value of C ε, ∞ should diminish morebetween cases (i) and (ii) compared to between cases (iii) and (iv). Such a comparison is carried out in Sec. IV usingDNS data.As can be seen from Eqs. (57)-(59), the force does not explicitly enter in the asymptote C ε, ∞ but does so in thecoefficients C and D . Therefore a dependence of C and D , and hence of the approach to the asymptote, on the forcemay be expected, while C ε, ∞ appears to be unaffected by the external force. However, different external forces willlead to different energy transfer scenarios, e.g. mainly dynamo and inertial transfer or mainly conversion of magneticto kinetic energy due to a strong Lorentz force, therefore the asymptote will be implicitly influenced by that. In short,nonuniversal values of C ε, ∞ , C and D are expected depending on the level of the ideal invariants and the type ofexternal force. We will address this point in further detail in Secs. IV and V. IV. COMPARISON TO DNS DATA
Before comparing Eq. (62) with DNS data the numerical method is briefly outlined. Equations (2)-(4) are solvednumerically in a three-dimensional periodic domain of length L box = 2 π using a fully de-aliased pseudospectral MHDcode [43, 44]. Both the initial magnetic and velocity fields are random Gaussian with zero mean with energy spectragiven by E mag,kin ( k ) = Ak exp( − k / (2 k ) ) , (72)where A (cid:62) k which locates the peak of the initial spectrum is taken to be k = 5 unless otherwise stated. No background magneticfield is imposed.Several series of simulations have been carried out for stationary and freely decaying MHD turbulence. In thecase of free decay the dependence of the asymptote on the initial level of the ideal invariants is studied. For thestationary simulations all helicities are initially negligible while the influence of different forcing methods is assessedby applying three different external mechanical forces labeled f , f and f to maintain the simulations in stationarystate, resulting in three different series of stationary DNSs. The forces always act at wavenumbers k (cid:54) k f = 2 . i.e. at the large scales. The first type of mechanical force f corresponds to the DNS series ND in Tbl. III and isgiven by ˆ f ( k , t ) = ( ι kin / E f ) ˆ u ( k , t ) for 0 < | k | < k f ;= 0 otherwise , (73)where ˆ f ( k , t ) is the Fourier transform of the forcing and E f is the total energy contained in the forcing band. Thesecond type of mechanical force f , which corresponds to the DNS series HF in Tbl. III is a random δ ( t )-correlatedprocess. It is based on a decomposition of the Fourier transform of the force into helical modes and has the advantagethat the helicity of the force can be adjusted at each wavevector [45], which gives optimal control over the helicityinjection. For all simulations using this type of forcing the relative helicity of the force was set to zero. The thirdtype of mechanical force [46] corresponds to the DNS series SF in Tbl. III and is given by f = f (cid:88) k f sin k f z + sin k f y sin k f x + sin k f z sin k f y + sin k f x , (74)where f is an adjustable constant. This type of force is nonhelical by construction.All three forces have been used in several simulations of stationary homogeneous MHD turbulence. The schemelabeled f was shown by Sahoo et al. [47] to keep the helicities at negligible levels even though zero helicity injection1 Run id N k max η mag k max η kin R − R L R λ ε mag /ε ε kin /ε µ = ν k C ε σ ρ c (0)H1 128 .
009 5 10 0.756 0.008 0H2 256 .
008 5 10 0.704 0.007 0H3 512 .
002 15 10 0.608 0.001 0H4 256 . .
003 5 10 0.450 0.004 0H8 512 . .
002 5 10 0.384 0.003 0H10 512 . .
001 5 10 0.310 0.004 0H12 1024 . . . . .
002 5 1 0.380 - 0.6CH06H2 512 . . . . R L denotes the integral-scaleReynolds number, R λ the Taylor-scale Reynolds number, R − the generalized Reynolds number, µ the magnetic resistivity, k the peak wavenumber of the initial energy spectra, k max the largest resolved wavenumber, η mag = ( µ /ε mag ) / and η kin = ( ν /ε kin ) / the magnetic and kinetic Kolmogorov microscales, respectively, at the peak of total dissipation, C ε the dimensionless total dissipation rate, σ the standard error on C ε and ρ c (0) the initial relative crosshelicity. All Reynolds numbers are measured at the peak of total dissipation. cannot be guaranteed with this forcing scheme. At very low Reynolds number this conservation of helicities appearsto be broken and induces peculiar self-ordering effects [48, 49]. The adjustable helicity forcing f has been extensivelyused in the literature [45, 50–52], mainly when nonzero levels of kinetic [45] or magnetic [50–52] helicity injection arerequired. The third forcing scheme f has been employed in the simulations by Dallas and Alexakis [46], where itwas shown that despite zero injection of all helicities, the system self-organized into large-scale fully helical states assoon as electromagnetic forces were applied. A. Decaying MHD turbulence
In this section we review the numerical results of Ref. [35]. In order to compare data for nonstationary systems atdifferent generalized Reynolds numbers, we measure all quantities at the peak of dissipation [10, 33]. Ensembles of upto 10 runs per data point were used in order to calculate statistics. Evidently, a larger ensemble would be desirable,however, especially at high resolution the computational cost of a single run is already substantial. Therefore wecompromised on the ensemble size in favor of running larger simulations, which is essential for the present study.Four series of simulations were carried out which differ between each other in the initial values of the ideal invariants H m and H c . Series H refers to a series with maximal initial H m while H c = 0, series CH06H was initialized withmaximal H m and relative cross-helicity ρ c = H c / ( U B ) = 0 .
6. Series NH and CH06NH label simulations initializedwith H m = 0 and differ in the initial level of ρ c , with ρ c = 0 for series NH and ρ c = 0 . η mag = ( µ /ε mag ) / and η kin = ( ν /ε kin ) / , thatis k max η mag,kin (cid:62)
1. Further details of series H and CH06H are given in Tbl. I while details corresponding to seriesNH and CH06NH are shown in Tbl. II.Figure 1(a) shows fits of Eq. (62) to DNS data for datasets that differ in the initial value of H m and H c . As canbe seen, eq. (62) fits the data very well. For the series H runs and for R − >
70 it is sufficient to consider terms offirst order in R − , while for the series NH the first-order approximation is valid for R − > Run id N k max η mag k max η kin R − R L R λ ε mag /ε ε kin /ε µ = ν k C ε σ ρ c (0)NH1 256 .
004 5 10 0.587 0.005 0NH2 256 .
003 5 10 0.530 0.004 0NH3 512 .
002 5 10 0.468 0.004 0NH4 512 . . .
001 5 5 0.358 0.002 0NH7 1024 . . . .
002 5 1 0.482 - 0.6CH06NH2 512 . . . . R L denotesthe integral-scale Reynolds number, R λ the Taylor-scale Reynolds number, R − the generalized Reynolds number, µ the magneticresistivity, k the peak wavenumber of the initial energy spectra, k max the largest resolved wavenumber, η mag = ( µ /ε mag ) / and η kin = ( ν /ε kin ) / the magnetic and kinetic Kolmogorov microscales, respectively, at the peak of total dissipation, C ε the dimensionless total dissipation rate, σ the standard error on C ε and ρ c (0) the initial relative crosshelicity. All Reynolds numbers are measured at the peak of total dissipation. Fig 1(b), where a function of the form
C/R n − was fitted to data from series H for R − >
70 after subtraction of theasymptote C ε, ∞ . The fit resulted in the value n = 1 . ± .
01 for the exponent, and the data from series NH, CH06Hand CH06NH are consistent with this result. Furthermore, Fig. 1(a) shows that the cross-helical CH06H runs gaveconsistently lower values of C ε compared to the series H runs, while little difference was observed between seriesCH06NH and NH. The asymptotes were C ε, ∞ = 0 . ± .
008 for the H series, C ε, ∞ = 0 . ± .
013 for the NHseries, C ε, ∞ = 0 . ± .
006 for the CH06H series and C ε, ∞ = 0 . ± .
005 for the CH06NH series.As predicted by the qualitative theoretical arguments outlined previously, the measurements show that the asymp-tote calculated from the nonhelical runs is larger than for the helical case, as can be seen in Fig. 1. The asymptotes ofthe series H and NH do not lie within one standard error of one another. Simulations carried out with H c (cid:54) = 0 suggestlittle difference in C ε for magnetic fields with initially zero magnetic helicity. For initially helical magnetic fields C ε is further quenched if H c (cid:54) = 0. In view of nonuniversality, an even larger variance of C ε, ∞ can be expected onceother parameters such as external forcing, plasma β , P m , etc., are taken into account. Here, attention is restrictedto nonuniversality related to different levels of cross- and magnetic helicity. The effect of external forcing will beanalyzed in the next section.
B. Stationary MHD turbulence
In this case measurements are taken after the simulations have reached a stationary state. The value of C ε andthe corresponding statistics for each data point are calculated from time-series obtained by evolving the stationarysimulations for a minimum of 9 large-eddy turnover times, as specified in Tbl. III. All runs of the series ND satisfy k max η mag,kin (cid:62) .
37 and as such are sufficiently resolved. The runs SF4, HF3 and HF4 are marginally resolved, andwe point out that the fitting procedure only involved data obtained from the ND series. The data points obtainedfrom the series SF and HF are included for comparison purposes. Further details of the stationary simulations aregiven in Tbl. III. All helicities are initially negligible and remain so during the evolution of the simulations.Figure 1 (a) shows error-weighted fits of Eq. (62) to DNS data obtained from series ND. As can be seen, Eq. (62)fits the data well, provided terms of second order in R − are included to accommodate the data points at low R − .For R − >
80, it is sufficient to consider terms of first order in R − only. The power law scaling to first order in R − is shown in further detail in Fig. 1 (b), where a function of the form C/R n − was fitted to data from series ND for R − >
80 after subtraction of the asymptote C ε, ∞ . The fit resulted in the value n = 1 . ± .
02 for the exponent,thus confirming that Eq. (62) describes the variation of C ε at moderate to high R − well already at first order in R − . This is similar to results in isotropic hydrodynamic turbulence, where the corresponding hydrodynamic equation3 Run id N k max η mag k max η kin R − R L R λ ε mag /ε ε kin /ε µ = ν C ε σ C ε t/T ND1 128 3.22 2.53 27.93 78.68 42.65 0.28 0.72 0.01 0.280 0.0028 30ND2 128 2.59 2.52 35.00 90.97 48.24 0.47 0.53 0.009 0.317 0.0065 26ND3 128 2.48 2.25 38.31 106.34 54.22 0.40 0.60 0.008 0.269 0.008 14ND4 128 2.12 2.12 45.63 119.92 59.59 0.50 0.50 0.007 0.285 0.0021 27ND5 128 1.85 1.93 53.12 137.38 65.32 0.54 0.46 0.006 0.290 0.0057 21ND6 128 1.46 1.58 66.21 173.33 75.83 0.58 0.42 0.0045 0.285 10 − f , f and f , respectively. R L denotes the integral-scale Reynolds number, R λ the Taylor-scale Reynolds number, R − the generalized Reynolds numbergiven in (46), µ the magnetic resistivity, k max the largest resolved wavenumber, η mag = ( µ /ε mag ) / and η kin = ( ν /ε kin ) / the magnetic and kinetic Kolmogorov microscales, respectively, C ε the dimensionless total dissipation rate defined in (37), σ C ε the standard error on C ε and t/T the run time in stationary state in units of large-eddy turnover time T = L u /U . All Reynoldsnumbers are time averages. C εu = C ε, ∞ u + C u /R L agreed well with the data for R L >
80 [30]. In this context we point out that the lowest valueof R L in Ref. [30] was R L = 81 .
5. For R − (cid:54)
80 we find the data to be consistent with Eq. (62) once terms of secondorder in R − are taken into account. However, it can be very difficult to extract two power laws clearly from numericaldata, especially if the leading and the subleading coefficients are of opposite sign [54, 55]. This is the case here, thesubleading coefficient D is always negative while the leading coefficient C is always positive.Figure 1 also shows that the result is independent of the forcing scheme, as the datasets obtained from simulationsusing the three different forcing functions are consistent with each other. This is likely to change if the strategy ofenergy input is fundamentally changed, for example if the helical content and/or the characteristic scale of the forceare altered, or if an electromagnetic force is used. We will come back to this point in Sec. V. The independence of C ε of the forcing scheme established here only shows independence of the specific implementation of the forcing. Theasymptote has been calculated to be C ε, ∞ = 0 . ± . C ε, ∞ obtained from the different simulation series is provided in Tbl. IV.Comparing the results from stationary and decaying simulations with the same level of helicities could shed somelight on the effect of external forces in the context of nonuniversality. As such, we compare the measured value of C ε, ∞ = 0 .
223 to the value calculated for the series of decaying simulations NH, which results in a difference of about12%. Since C ε, ∞ does not depend explicitly on the external forces, the difference between the measured values mayoriginate from dynamical effects. In relation to the effect of initial cross- and magntic helicities on the value of C ε, ∞ discussed in Sec. IV A, we expect further variance in the measured value of C ε, ∞ depending on the level of cross-and magnetic helicities of the external forces. As can be seen from a comparison of the curves shown in Fig. 1, thecoefficient C in Eq. (62) is not the same for the stationary and decaying cases. This is expected since C ± and hence4 C e R − (a) NDHFSFNHHCH06HCH06NH C e − C e , ¥ R − (b) NDHFSFNHHCH06HCH06NH
FIG. 1. (a) Eq. (62) fitted to the different datasets H, NH, CH06H, CH06NH (free decay) and ND (stationary). Data fromthe stationary series SF and HF are shown for comparison purposes. The black lines refer to fits using Eq. (62) up to firstorder, while the gray lines use Eq. (62) up to second order in 1 /R − . The error bars show one standard error. As can be seen,the respective asymptotes differ between the data sets H, NH, CH06H, CH06NH and ND, while data from series SF and HFis compatible with data from series ND. (b) Fit of the expression C/R n − corresponding to eq. (62) after subtraction of theasymptote C ε, ∞ on a logarithmic scale to the datasets H (free decay, top line) and ND (stationary, bottom line). Data fromthe series NH, CH06H, CH06NH, SF and HF are shown for comparison purposes. The resulting exponents are n = 1 . ± . n = 1 . ± .
02 for the stationary case. C depend explicitly on the energy input, as can be seen from Eq. (58). Series id C ε, ∞ σ C ε, ∞ descriptionH 0.241 0.008 decaying, helicalNH 0.265 0.013 decaying, nonhelicalCH06H 0.193 0.006 decaying, helical, cross-helicalCH06NH 0.268 0.005 decaying, nonhelical, cross-helicalND 0.223 0.003 stationary, nonhelicalRef. [30] 0.468 0.006 non-conducting, stationaryTABLE IV. Summary of measured values for the asymptotic dimensionless dissipation rate C ε, ∞ from the different simulationseries. For comparison the measured values from a series of DNSs of stationary isotropic hydrodynamic turbulence [30] areincluded. The error on C ε, ∞ is denoted by σ C ε, ∞ and refers to the error obtained from the error-weighted fitting procedure. Some further observations can be made from a comparison of the datasets concerning the variance of the magneticand kinetic contributions, ε mag and ε kin , to the total energy dissipation rate. The fractions ε mag /ε and ε kin /ε aregiven columns 8 and 9, respectively, of Tbls. I-III for the different datasets. For the series H and CH06H it can beseen that the kinetic dissipation fraction ε kin /ε grows with increasing R − (or R L ), of course at the expense of themagnetic dissipation fraction ε mag /ε . The reason for this behaviour could be connected with the inverse cascade ofmagnetic helicity becoming more efficient at higher R − resulting in a higher residual inverse transfer of magneticenergy and thus slightly less magnetic energy to be dissipated at the small scales. This interpretation is supportedby the observation that no such variation is present for the nonhelical series NH and CH06NH, where we measure ε mag /ε (cid:39) . ε mag /ε (cid:39) . ε mag /ε > ε kin /ε .The stationary datasets show yet a different variation of ε mag /ε and ε kin /ε with R − . For R − <
45 we find ε mag /ε < ε kin /ε , while for R − >
46 the results are similar to free decay with ε mag /ε > ε kin /ε . The magneticdissipation fraction increases with R − from ε mag /ε (cid:39) . ε mag /ε (cid:39) .
7, where it appears to reach a plateau.The fluctuating magnetic field is maintained by the velocity field fluctuations through a nonlinear dynamo processin the present DNSs of stationary MHD turbulence. Hence in the statistically stationary state ε mag is in balancewith the dynamo term (cid:104) b · ( b · ∇ ) u (cid:105) , and the measured values of ε mag /ε and ε kin /ε imply that at lower Reynolds5number the nonlinear dynamo is less efficient in maintaining small-scale magnetic field fluctuations than at higherReynolds number. Similar conclusions have been reached in a study of the magnetic Prandtl number dependence ofthe ratio ε kin /ε mag [56], where the efficiency of the dynamo at different values of P m was linked to measured valuesof ε kin /ε mag . V. CONCLUSIONS
The behavior of the dimensionless dissipation coefficient C ε in homogeneous MHD turbulence with P m = 1 and nobackground magnetic field is well described by C ε = C ε, ∞ + CR − + DR − + O ( R − − ) . (75)This equation was derived from the energy balance equation for z ± in real space (the vKHE) by outer asymptoticexpansions in powers of 1 /R ∓ , leading necessarily to a large-scale description of the behavior of the dimensionlessdissipation rate . The approximative equation (75) has been shown to agree well with data obtained from mediumto high resolution DNSs of both decaying MHD turbulence at the peak of dissipation and statistically steady MHDturbulence sustained by large-scale forcing. The measurements for C ε, ∞ ranged from 0 . (cid:54) C ε, ∞ (cid:54) .
268 betweenthe different series of simulations. Interestingly, the measured values of C ε, ∞ for MHD are smaller than the measuredvalue of C ε, ∞ (cid:39) . R − → ∞ originates from the sum of the nonlinear terms in the momentum andinduction equations, that is, it measures the total transfer flux, which is expected to depend on the values of the idealinvariants. As predicted, the values of the respective asymptotes from the datasets differ, suggesting a dependence of C ε, ∞ on different values of the helicities, and thus a connection to questions of universality in MHD turbulence. Formaximally helical magnetic fields C ε, ∞ is smaller than for nonhelical fields. This is expected from the inverse cascadeof magnetic helicity. The dependence of C ε, ∞ on the remaining ideal invariant, the cross-helicity, is more complex.Since C ε, ∞ describes the flux of total energy across the scales, this flux is expected to diminish for increasing cross-helicity. This is indeed the case for helical magnetic fields, where C ε, ∞ depends on the cross-helicity in the expectedway. Surprisingly, for nonhelical magnetic fields C ε, ∞ does not depend on the cross-helicity. This is consistent withthe asymmetric effect of the cross-helicity on forward and inverse fluxes of total energy suggested by the analysis oftriad interactions in Ref. [42], where high levels of cross-helicity were found to quench forward transfer more thaninverse transfer. A similar effect can be inferred from predictions obtained from statistical mechanics [22], wherethe simultaneous presence of cross- and magnetic helicities resulted in inverse transfers of both magnetic and kineticenergy. In this case, the forward flux of total energy should be lower than for all other cases, which is consistent withthe numerical results presented here. Concerning stationary nonhelical MHD turbulence, we found that C ε, ∞ differedby about 12% from the value measured for nonhelical decaying turbulence (series NH) at the peak of dissipation.According to the results from the asymptotic analysis, there is no explicit dependence of C ε, ∞ on the external force.As such, the difference in the measured value of C ε, ∞ between the stationary and the decaying systems may be dueto dynamical effects, which may be interpreted as further support for nonuniversal values of C ε . The approach tothe asymptote is predicted to differ between decaying and stationary systems due to the explicit dependence of thecoefficient C in Eq. (62) on the forcing. This is indeed observed in the simulations.The numerical results showed that C ε, ∞ is universal with respect to different forcing schemes applied to the samefield in the same wavenumber range, thus confirming that the particular functional form and stochasticity of a large-scale force is irrelevant to the small-scale turbulent dynamics as long as the ideal invariants remain the same for thedifferent forcing schemes. However, as mentioned in Sec. IV, this is expected to change if the strategy of energy inputis changed. The effect of large-scale magnetic forcing on the scaling of ε with different rms quantities was investigatedrecently [26]. Numerical results showed that ε ∼ U /L f u for a large region of parameter space even in the presenceof electromagnetic forces. Only once the large-scale magnetic field became very strong a different scaling related tothe dominance of magnetic shear over mechanical shear was found: ε ∼ BU /L f b . Differences may also be expectedfor forces applied at smaller scales. The analysis presented here relies on taking outer asymptotic expansions of allscale-dependent functions in the vKHE, including the energy input from the forcing. Here it was crucial to assumethat the system was forced at the large scales, as the limit of infinite Reynolds number was defined as energy input atthe lowest wavenumbers k → k → ∞ . This clearly precludes theapplication of the present analysis to situations where the system is forced at intermediate or small scales. Therefore,it can be expected that systems forced at intermediate scales deviate from the 1 /R − -scaling of C ε . For hydrodynamics,this is the case [54, 67, 68]. Furthermore, experimental and numerical results for nonstationary flows [69, 70] suggest6even further variance possibly due to the influence of the time-derivative of the second-order structure function in thevKHE.The results presented here were restricted to homogeneous MHD turbulence at P m = 1 without a mean magneticfield. In general, further variance in the measured value for C ε, ∞ is expected depending on e.g. the magnetic Prandtlnumber [56] or the influence of a background magnetic field. The presence of a background magnetic field, whichleads to spectral anisotropy and the breakdown of the conservation of magnetic helicity [71], will introduce severaldifficulties to be overcome when generalizing the analytical approach. The spectral flux will then depend on thedirection of the mean field [12, 72] and a more generalized description and role for the magnetic helicity would beneeded [73, 74]. Other questions concern the generalization of this approach to MHD flows with magnetic Prandtlnumbers P m (cid:54) = 1, the effect of compressive fluctuations or the influence of other vector field correlations on thedissipation rate and/or the approach to the asymptote as observed by Dallas and Alexakis [10].
VI. ACKNOWLEDGMENTS
We thank V. Dallas for useful discussions. A. B. acknowledges support from the UK Science and Technol-ogy Facilities Council, M. L. received support from the UK Engineering and Physical Sciences Research Coun-cil (EP/K503034/1) and E. E. G. was supported by a University of Edinburgh Physics and Astronomy Fellow-ship. This work has made use of the resources provided by the UK National Supercomputing facility ARCHER,( ), made available through the Edinburgh Compute and Data Facility (ECDF, ). The research leading to these results has received funding from the European Union’sSeventh Framework Programme (FP7/2007-2013) under grant agreement No. 339032.
Appendix A: Gauge independence of Equation (9)
In order to prove that Eq. (9) is correct for an arbitrary choice of gauge, we first express the current density j interms of the vector potential a j = ∇ × ( ∇ × a ) = − ∆ a + ∇ ( ∇ · a ) , (A1)which holds in any gauge. In Fourier space this relation becomesˆ j = k ˆ a + i k ( i k · ˆ a ) , (A2)hence one obtains (cid:90) Ω d k (cid:28) k ˆ j · ˆ b ∗ (cid:29) = (cid:90) Ω d k (cid:28)(cid:18) ˆ a + i k k ( i k · ˆ a ) (cid:19) · ˆ b ∗ (cid:29) = (cid:90) Ω d k (cid:68) ˆ a · ˆ b ∗ (cid:69) = H m , (A3)since b is a solenoidal vector field. Equation (9) now follows by writing the Fourier coefficients of the magnetic fieldand the current density in the Els¨asser formulationˆ b = ˆ z + − ˆ z − j = i k × ˆ b = i k × ˆ z + − ˆ z − . (A4) [1] A. N. Kolmogorov, C. R. Acad. Sci. URSS , 301 (1941).[2] P. S. Iroshnikov, Soviet Astronomy , 566 (1964).[3] R. H. Kraichnan, Phys. Fluids , 1385 (1965).[4] P. Goldreich and S. Sridhar, Astrophys. J. , 763 (1995).[5] S. Boldyrev, ApJ , L37 (2005).[6] S. Boldyrev, Phys. Rev. Lett. , 115002 (2006).[7] A. Beresnyak and A. Lazarian, ApJL , L175 (2006).[8] J. Mason, F. Cattaneo, and S. Boldyrev, Phys. Rev. Lett. , 255002 (2006).[9] G. Gogoberidze, Phys. Plasmas , 022304 (2007).[10] V. Dallas and A. Alexakis, Phys. Fluids , 105106 (2013).[11] V. Dallas and A. Alexakis, Phys. Rev. E , 063017 (2013). [12] M. Wan, S. Oughton, S. Servidio, and W. H. Matthaeus, J. Fluid Mech. , 296 (2012).[13] A. A. Schekochihin, S. C. Cowley, and T. A. Yousef, in IUTAM Symposium on Computational Physics and New Perspectivesin Turbulence , edited by Y. Kaneda (Springer, Berlin, 2008) pp. 347–354.[14] P. D. Mininni, Annu. Rev. Fluid Mech. , 377 (2011).[15] R. Grappin, A. Pouquet, and J. L´eorat, Astron. Astrophys. , 51 (1983).[16] A. Pouquet and P. Mininni and D. Montgomery and A. Alexakis, in IUTAM Symposium on Computational Physics andNew Perspectives in Turbulence , edited by Y. Kaneda (Springer, 2008) pp. 305–312.[17] A. Beresnyak, Phys. Rev. Lett. , 075001 (2011).[18] S. Boldyrev, J. C. Perez, J. E. Borovsky, and J. J. Podesta, Astrophys. J. , L19 (2011).[19] R. Grappin and W.-C. M¨uller, Phys. Rev. E , 026406 (2010).[20] E. Lee, M. E. Brachet, A. Pouquet, P. D. Mininni, and D. Rosenberg, Phys. Rev. E , 016318 (2010).[21] S. Servidio, W. H. Matthaeus, and P. Dmitruk, Phys. Rev. Lett. , 095005 (2008).[22] U. Frisch, A. Pouquet, J. L´eorat, and A. Mazure, J. Fluid Mech. , 769 (1975).[23] A. Pouquet, U. Frisch, and J. L´eorat, J. Fluid Mech. , 321 (1976).[24] A. Pouquet and G. S. Patterson, J. Fluid Mech. , 305 (1978).[25] D. Biskamp, Nonlinear Magnetohydrodynamics. , 1st ed. (Cambridge University Press, 1993).[26] A. Alexakis, Phys. Rev. Lett. , 084502 (2013).[27] K. R. Sreenivasan, Phys. Fluids , 1048 (1984).[28] K. R. Sreenivasan, Phys. Fluids , 528 (1998).[29] W. D. McComb, Homogeneous, Isotropic Turbulence: Phenomenology, Renormalization and Statistical Closures (OxfordUniversity Press, 2014).[30] W. D. McComb, A. Berera, S. R. Yoffe, and M. F. Linkmann, Phys. Rev. E , 043013 (2015).[31] P. K. Yeung, X. M. Zhai, and K. R. Sreenivasan, PNAS , 1263312638 (2015).[32] S. Jagannathan and D. A. Donzis, J. Fluid Mech. , 669 (2016).[33] P. D. Mininni and A. G. Pouquet, Phys. Rev. E. , 025401 (2009).[34] V. Dallas and A. Alexakis, Astrophys. J. , L36 (2014).[35] M. F. Linkmann, A. Berera, W. D. McComb, and M. E. McKay, Phys. Rev. Lett. , 235001 (2015).[36] H. Politano and A. Pouquet, Phys. Rev. E , R21 (1998).[37] W. M. Els¨asser, Phys. Rev. , 183 (1950).[38] S. Chandrasekhar, Proc. Roy. Soc. London. Series A , 435 (1951).[39] The scaling is ill-defined for the (measure zero) cases u = ± b , which correspond to exact solutions to the MHD equationswhere the nonlinear terms vanish. Thus no turbulent transfer is possible, and these cases are not amenable to an analysiswhich assumes nonzero energy transfer [36].[40] E. A. Novikov, Soviet Physics JETP , 1290 (1965).[41] T. S. Lundgren, Phys. Fluids , 638 (2002).[42] M. F. Linkmann, A. Berera, M. E. McKay, and J. J¨ager, J. Fluid Mech. , 61 (2016).[43] A. Berera and M. F. Linkmann, Phys. Rev. E , 041003(R) (2014).[44] M. Linkmann, Self-organisation in (magneto)hydrodynamic turbulence , Ph.D. thesis, University of Edinburgh (2016).[45] A. Brandenburg, Astrophys. J. , 824 (2001).[46] V. Dallas and A. Alexakis, Phys. Fluids , 045105 (2015).[47] G. Sahoo, P. Perlekar, and R. Pandit, New Journal of Physics , 013036 (2011).[48] W. D. McComb, M. F. Linkmann, A. Berera, S. R. Yoffe, and B. Jankauskas, J. Phys. A: Math. Theor. , 25FT01(2015).[49] M. F. Linkmann and A. Morozov, Phys. Rev. Lett. , 134502 (2015).[50] W. C. M¨uller, S. K. Malapaka, and A. Busse, Phys. Rev. E , 015302 (2012).[51] S. K. Malapaka and W.-C. M¨uller, Astrophys. J. , 21 (2013).[52] M. Linkmann and V. Dallas, Phys. Rev. E , 053209 (2016).[53] The data is publicly available, http://dx.doi.org/10.7488/ds/247 .[54] L. Biferale, A. S. Lanotte, and F. Toschi, Phys. Rev. Lett. , 094503 (2004).[55] D. Mitra, J. Bec, R. Pandit, and U. Frisch, Phys. Rev. Lett. , 194501 (2005).[56] A. Brandenburg, ApJ , 12 (2014).[57] J. Jim´enez, A. A. Wray, P. G. Saffman, and R. S. Rogallo, J. Fluid Mech. , 65 (1993).[58] L.-P. Wang, S. Chen, J. G. Brasseur, and J. C. Wyngaard, J. Fluid Mech. , 113 (1996).[59] P. K. Yeung and Y. Zhou, Phys. Rev. E , 1746 (1997).[60] N. Cao, S. Chen, and G. D. Doolen, Phys. Fluids , 2235 (1999).[61] B. R. Pearson, P. A. Krogstad, and W. van de Water, Phys. Fluids , 1288 (2002).[62] Y. Kaneda, T. Ishihara, M. Yokokawa, K. Itakura, and A. Uno, Phys. Fluids , L21 (2003).[63] B. R. Pearson, T. A. Yousef, N. E. L. Haugen, A. Brandenburg, and P. A. Krogstad, Phys. Rev. E , 056301 (2004).[64] D. A. Donzis, K. R. Sreenivasan, and P. K. Yeung, J. Fluid Mech. , 199 (2005).[65] W. J. T. Bos, L. Shao, and J.-P. Bertoglio, Phys. Fluids , 045101 (2007).[66] P. K. Yeung, D. A. Donzis, and K. R. Sreenivasan, J. Fluid Mech. , 5 (2012).[67] C. R. Doering, Annu. Rev. Fl. Mech. , 109 (2009).[68] B. Mazzi and J. C. Vassilicos, J. Fluid Mech. , 65 (2004).[69] P. C. Valente, R. Onishi, and C. B. da Silva, Phys. Rev. E , 023003 (2014). [70] J. C. Vassilicos, Annu. Rev. Fluid Mech. , 95 (2015).[71] W. H. Matthaeus and M. L. Goldstein, J. Geophys. Res. , 6011 (1982).[72] M. Wan, S. Servidio, S. Oughton, and W. H. Matthaeus, Phys. Plasmas , 090703 (2009).[73] M. A. Berger, J. Geophys. Res. , 2637 (1997).[74] M. A. Berger, Plasma Physics and Controlled Fusion41