Finite Dissipation in Anisotropic Magnetohydrodynamic Turbulence

In presence of an externally supported, mean magnetic field a turbulent, conducting medium, such as plasma, becomes anisotropic. This mean magnetic field, which is separate from the fluctuating, turbulent part of the magnetic field, has considerable effects on the dynamics of the system. In this paper, we examine the dissipation rates for decaying incompressible magnetohydrodynamic (MHD) turbulence with increasing Reynolds number, and in the presence of a mean magnetic field of varying strength. Proceeding numerically, we find that as the Reynolds number increases, the dissipation rate asymptotes to a finite value for each magnetic field strength, confirming the K\'arm\'an-Howarth hypothesis as applied to MHD. The asymptotic value of the dimensionless dissipation rate is initially suppressed from the zero-mean-field value by the mean magnetic field but then approaches a constant value for higher values of the mean field strength. Additionally, for comparison, we perform a set of two-dimensional (2DMHD) and a set of reduced MHD (RMHD) simulations. We find that the RMHD results lie very close to the values corresponding to the high mean-field limit of the three-dimensional runs while the 2DMHD results admit distinct values far from both the zero mean field cases and the high mean field limit of the three-dimensional cases. These findings provide firm underpinnings for numerous applications in space and astrophysics wherein von K\'arm\'an decay of turbulence is assumed.


I. INTRODUCTION AND BACKGROUND
Turbulence is a ubiquitous, although incompletely understood phenomena. In turbulent astrophysical plasmas, velocity and magnetic fields are often equally important, and magnetohydrodynamic (MHD) turbulence, the case considered here, becomes an appropriate description. Energy supply for turbulence usually originates at large scales due to some type of stirring or driving mechanism, after which the energy cascades to smaller scales, finally reaching the dissipation scale where the energy is dissipated into heat by viscosity and resistivity. For incompressible hydrodynamics, the single scalar viscosity ν parameterizes microscopic nonideal effects, such that, for laminar flows the energy dissipation rate vanishes when ν → 0. However, for turbulent flows, following the hypothesis stated by Taylor [1] and von Kármán & Howarth [2] and later employed by Kolmogorov [3], in the zero viscosity, infinite Reynolds number Re ∼ 1/ν → ∞ limit, the energy dissipation rate approaches a constant, nonzero limit [4]. (The renowned Kolmogorov theory of universal scaling in the inertial range (K41) [5], in essence assumes this limiting behavior.) This so-called "dissipative anomaly" is well supported in experiments and computations [6,7].
For the case of incompressible MHD, there are two relevant dissipation coefficients, viscosity (ν) and resistivity (µ), associated with velocity and magnetic field, * whm@udel.edu respectively. The limit of zero viscosity and resistivity (ν, µ → 0) is again of interest with regard to fully developed turbulent systems. In practice one finds that, for essentially all astrophysical as well as terrestrial turbulent systems, the mechanical and magnetic Reynolds numbers (inverse of viscosity and resistivity in non-dimensional units) are of large, and sometimes colossal magnitude, justifying the success of turbulence phenomenologies like K41. For example, the effective mechanical and magnetic Reynolds number, in the solar wind, are of the order of 10 5 or more [8] (also see [9]), and much larger in the interstellar medium. Therefore, one might be tempted to assume ν, µ = 0 as an approximation. However, in analogy to the hydrodynamics case, the limit ν, µ → 0 is very different from the ν, µ = 0 case. This is manifested in the counterintuitive phenomena that as the limit ν, µ → 0 is approached, the total turbulent energy dissipation rate does not vanish, but remains finite. MHD energy dissipation is exactly zero for the ideal case (ν, µ = 0). Onsager [10] conjectured that this problem of energy-dissipation anomaly arises due to lack of smoothness in the velocity field increments in the context of three-dimensional hydrodynamic turbulence. Based on this idea of lack of smoothness in the velocity field, Duchon & Robert [11] derived a local form of dissipation which was generalized to incompressible Hall MHD by Galtier [12]. An apparent explanation of anomalous dissipation was given in a work by Cichowlas et al. [13] who showed that even for ideal Eulerian flows, the large-wavenumber modes play the role of an effective eddy viscosity leading to an approximate Kolmogorov scaling at the inertial scales. This idea was further extended to two-dimensional MHD by Krstulovic et al. [14]. Here we will examine this apparently singular behavior (see e.g., [15][16][17][18]) not from a rigorous mathematical perspective, but from an empirical perspective based on accurate spectral method computations.
Taylor first suggested, based on empirical arguments, that the dissipation rate in a fully developed turbulent system becomes independent of viscosity. Kármán and Howarth [2] established Taylor's results more rigorously, assuming that the shape of the two-point correlation function remains unchanged. The first experimental verification of Kármán and Howarth's result came from an experiment performed by Batchelor and Townsend [19][20][21][22] using wind tunnel measurements. Politano and Pouquet [23] generalized the von Kármán-Howarth equations for isotropic MHD. Wan et al. [24] investigated nonuniversality of Kármán-Howarth like decay in MHD in the presence of factors like mean magnetic field, helicity, etc.
Recently, Wu et al. [25] and Parashar et al. [26] have shown using particle in cell (PIC) simulations that even in the case of weakly collisional plasmas, the von Kármán decay law remains valid. Most of these studies were interested in the temporal behavior of the fluctuation amplitude and its relationship with the large-scale fluctuations and an energy-containing length scale. For neutral fluids, the variation of dimensionless dissipation rate, C ǫ with Reynolds number was studied experimentally in shear flows [27], numerically in homogeneous incompressible flows [28], and in weakly compressible flows [29]. In all cases, an asymptotic value of C ǫ ≈ 0.5 was found. We remark here that the dimensionless dissipation rate, C ǫ , as defined in Pearson et al. [29] is proportional to one of the constants in von Kármán decay phenomenology (See Appendix B of [30]). Mininni and Pouquet [18] carried out direct numerical simulations (DNS) of decaying isotropic MHD turbulence, demonstrating that the mean dissipation rate per unit mass, ǫ, remains finite and becomes independent of viscosity and resistivity as the Reynolds number (Taylor-scale Reynolds number, R λ , in this case) increases. Dallas and Alexakis [31] performed a similar analysis showing the variation of dimensionless dissipation rate with Taylor-scale Reynolds number with initial velocity and current density being highly correlated. A comparison with [18] revealed that the dimensionless dissipation rate saturates to a finite value but the level of saturation depends on the strength of initial cross-correlation. Recently, Linkmann et al. [32,33] have performed a series of investigations for similar analysis in isotropic MHD. By fitting a model equation, an asymptotic value of dimensionless dissipation rate, C ǫ,∞ = 0.265 ± 0.013, was found for nonhelical decaying MHD with no mean field. This value is considerably different from the fluid case. McComb et al. [34] derived a similar model equation for fluid turbulence.
The presence of a mean magnetic field suppresses the decay rate in MHD systems and renders the system anisotropic. Early studies by Hossain et al. [35,36] at low Reynolds number demonstrated that the mean field initially inhibits dissipation, but the effect of Alfvénic propagation soon saturates for sufficiently large values of the mean field. More recently, Bigot et al. [37][38][39] investigated the energy decay in the presence of a strong uniform magnetic field. Zhdankin et al. [40] also studied the effect of a mean field at dissipation sites in MHD turbulence. Naturally, one would expect a decrease of the dimensionless dissipation rate (or, the von Kármán decay constant) with increasing strength of mean field. We present here a further analysis of the dissipation rate for increasing values of mean field. We find that for each value of the mean field, the dissipation rate asymptotes to a finite value in the limit of large Reynolds number. These results are relevant for systems in which the anisotropy due to presence of a mean magnetic field cannot be neglected. Such situations are often realized in astrophysical plasmas near stellar objects. Our results show that the value of the dimensionless dissipation rate is well separated from the isotropic value for a strong mean field. We first reestablish the isotropic case by comparing with results from [33] and then proceed in similar fashion to anisotropic cases.

II. EQUATIONS AND APPROACH
The equations of three dimensional incompressible MHD, in the absence of external forcing, are written as, where u is the fluctuating velocity field, B is the total magnetic field which we assume can be decomposed into a spatially uniform mean and a fluctuating part, B(r, t) = B 0 + b(r, t). Without any loss of generality, we choose B 0 = B 0ẑ . P is the total (thermal+magnetic) pressure field, j = (∇ × B)/µ 0 is the current density, ν is the kinematic viscosity, and µ is the magnetic diffusivity. Eq.
(3) enforces incompressibility. For simplicity, we assume the fluid density ρ is constant and set it to unity. For this study, we consider unit magnetic Prandtl number (Pm = ν/µ), i.e., equal viscosity and resistivity, With this constraint, the MHD equations Eqs. (1)-(4) can be written in terms of the Elsasser variables, z ± = u ± b, as, Here, v A is the Alfvén velocity defined as v A = B 0 / √ 4π. The Elsasser energies are given by Following the argument associated with preservation of the functional form of the two-point correlation functions, one can generalize the von Kármán-Howarth [2] result to isotropic MHD [23] or MHD isotropic in a plane perpendicular to a strong mean field [24], where α ± and β ± are constants, and L ± are the energycontaining scales corresponding to the two Elsasser variables. Precise definitions of L ± depend on the phenomenology used. We shall come back to the definition of energy-containing scale later in the paper. For small cross helicity (H c ≃ 0), the " + " and " − " variables remain almost equal, so that one can write, where, Z ≃ Z + ≃ Z − . It is apparent from these equations that the sequence of assumptions (especially, selfpreservation of the correlations) leading to this point implies that the energy dissipation rate becomes independent of the dissipation coefficients ν = µ. It is one of the main goals of this paper to measure the finite energy dissipation rate in the asymptotic limit of large Reynolds number for MHD in presence of a mean field at zero magnetic and cross helicity. The mean turbulent energy dissipation rate per unit mass in incompressible MHD, for unit magnetic Prandtl number (Pm), can be written as, where ω = ∇×u is the vorticity and · · · denotes spatial averaging.
Linkmann et al. [32,33] proposed a definition of the non-dimensional energy dissipation rate in MHD turbulence as, where where where z ± (k, t) denotes the Fourier transform of the respective Elsasser variables. Assuming isotropy, the rms values of the Cartesian components of Z 2 ± are given by, The W ± are the Elsasser variable analogs of the u ′ in E u = 3(u ′ ) 2 /2, which is the classical definition often employed in isotropic hydrodynamic work [19]. Here we use where E(t) is the total energy and H c (t) is the cross helicity For H c ≃ 0, one expects C + ǫ ≃ C − ǫ . Further, Linkmann et al. [32,33] also suggested a generalized Reynolds number defined as In this paper we have studied the variation of dimensionless dissipation rate, as defined by Eqs. (13)-(14) with generalized Reynolds number defined in Eq. (20). For comparison, we also report the integral scale Reynolds number defined as where and u ′ denotes the rms speed with E u = 3(u ′ ) 2 /2. Here E u (k) denotes the omnidirectional kinetic energy spectrum and E u = E u (k)dk is the total kinetic energy. Similarly, the Taylor microscale Reynolds number is defined as where The effect of a mean magnetic field on magnetohydrodynamic turbulence was first investigated and quantified by Shebalin et al. [41] in two dimensions, and generalized to three dimensions by Oughton et al. [42]. Spectral anisotropy in three dimensions is quantified by defining Shebalin angles, θ Q as where k is the wave vector component along the direction of the mean magnetic field, k ⊥ is the wave vector in the plane perpendicular to the mean field, and the summations extend over all values of wave vector k. In the present notation, k = k z and k 2 ⊥ = k 2 x + k 2 y . In general Q can be any vector field, although here we consider only velocity field u, and magnetic field b.
We define where E b = δb 2 /2 is the total magnetic field fluctuation energy. Therefore, √ r is the ratio of fluctuation amplitude and mean magnetic field strength.

III. NUMERICAL METHOD
We solve the MHD equations in a periodic box of dimension (2π) 3 using a pseudo-spectral method in a Fourier basis. We employ second-order Runge-Kutta (RK2) scheme for time advancement and the 2/3 rule for dealiasing. Instead of using an automatically adjustable time step, in this study the time step was held constant for a set of runs, but halved when instabilities occurred. All the runs have been initialized with a modal spectrum proportional to 1/[1 + (k/k 0 ) 11/3 ], where the "knee" of the spectrum is k 0 = 4, and only Fourier modes within the band 1 ≤ k ≤ 15 are excited. The total kinetic and magnetic energy are both normalized to 0.5 initially. All the measurements have been made at the time of highest dissipation. In all runs, the kinetic, magnetic, and cross helicities are small initially and remain so during the simulation times considered. Further, since velocity and magnetic fields are generated independently, we may assume that cross correlations involving higher derivatives [31] are suitably small throughout.
In all runs, the maximum resolved wavenumber, k max , is greater than the dissipation wavenumber (reciprocal of the Kolmogorov length scale η) defined as Wan et al. [43] showed that the ratio k max /k diss = k max η should be at least three for sufficient numerical accuracy of fourth-order moments. However, for studies of the present kind, as was earlier reported by Pearson et al. [29], we find that k max η ≥ 1 suffices and increasing resolution further does not substantially change the results. Recall that the more strict requirement proposed by Wan et al [43] pertains mainly to higher-order statistics and coherent structures, while accurate portrayal of lower-order quantities such as energy spectra, which control instantaneous dissipation, have less stringent requirements. This can also be seen from Table II of [43] where k diss becomes accurate in the range 1 ≤ k max η ≤ 1.9. Since k 4 diss = ǫ/ν 3 , this implies that dissipation rate also becomes accurate in this range. The details of the simulations used in this study are found in Table I. Measurements are made at instant of peak dissipation, and shortly after this time the simulations are stopped. Ideally, one should consider performing an ensemble of runs and then calculate the statistics. However, the computational cost being prohibitively expensive, particularly for the high mean field cases, we defer such refinements to a later time.

IV. RESULTS
The dimensionless dissipation rates, obtained from the runs in Table I, are plotted as function of the generalized Reynolds number in Fig. 1. Linkmann et al. [32,33] showed that for isotropic MHD one can fit a simple model equation to C ǫ , provided R − ≥ 80, where C ǫ,∞ is the asymptotic value of C ǫ as R − tends to ∞, A is a time dependent coefficient, and O(R −2 − ) represents terms of higher order, which we are neglecting here. For lower Reynolds numbers, second order and higher terms cannot be neglected and may be important to derive model equations such as those described in [32,33]. Since the main goal of this paper is not to build such model equations, but to measure the asymptotic value of the dimensionless dissipation rate and investigate its variation with mean magnetic field, we consider only high Reynolds numbers (R − ≥ 80) and work only with the leading-order term.
We fit Eq. (29) to each set of runs; see Table I. The solid lines in Fig. 1 are obtained from this fitting. We use a least-square fitting technique to fit the polynomial TABLE I. Parameters for three dimensional spectral simulations: Mean magnetic field strength B0, grid resolution N 3 , viscosity and resistivity ν and µ (set equal), Taylor scale Reyolds number R λ , Integral scale Reynolds number RL, Generalized Reynolds number R−, Dimensionless dissipation rate Cǫ, Shebalin angles for velocity field θu, and for magnetic field θ b , Square of ratio of fluctuation amplitude and mean magnetic field r = (δb/B0) 2 (Not Applicable for B0 = 0 case), time step dt, Ratio of maximum resolved wavenumber and dissipation wavenumber kmaxη. See text for definitions of the quantities. All measurements have been made near the time of maximum dissipation. Eq. (29) to the data for each value of the mean magnetic field, which ranges over values 0, 1, 3, 8, and 15. In all cases we see that the fits agree reasonably well with the the sets of individual data points, and in some cases the agreement is excellent. The asymptotic values, i.e., the values of C ǫ,∞ , obtained by this procedure, are reported in Table II. The value for the isotropic case is in agreement with the one reported in [33], within uncertainty.
The values of C ǫ,∞ for B 0 > 0 decrease first with increasing strength of the mean field, but then saturate. This can be seen by comparing the two sets of data obtained for B 0 = 8.0 and B 0 = 15.0, which almost lie on top of each other. Quantitatively, it can be seen from Table II that the two values of C ǫ for the B 0 = 8.0 and B 0 = 15.0 case lie within their limits of simulation uncertainty. This result is reminiscent of the Hossain et al. [36] study, that observed a suppression of dissipation rate due to a mean field of moderate strength compared to the dissipation rate in the isotropic case. This effect, however, saturates at higher mean field strengths. This is explained by noticing that the mean magnetic field suppresses spectral transfer parallel to the mean field, so that the spec-  Table I, and sets of RMHD and 2DMHD runs. Additionally, we report the value of the constant A in Eq. (29).
To compare the convergence of the non-dimensional dissipation rate in three-dimensional MHD in the presence of a strong mean field with two-dimensional MHD (2DMHD), we perform a set of 2DMHD simulations with decreasing viscosity and resistivity (See [41] for governing equations and other details). The two-dimensional simulations are performed in almost identical conditions as the three dimensional ones. Here, the box dimension is (2π) 2 and the initial spectrum is excited in 1 ≤ k ≤ 15 band, with the modal spectrum proportional to 1/[1 + (k/4) 8/3 ] so that the omni-directional spectrum is ∼ k −5/3 at large k. As can be seen in Fig. 1, the values corresponding to 2DMHD are close to neither the B 0 = 0 case nor the B 0 = 8, 15 limit. This result demonstrates the unique nature of two-dimensional decay of magnetofluids. Although spectral transfer becomes progressively two-dimensional as the mean magnetic field is increased, the limit is not identical with the perfectly two-dimensional system. This discontinuity is presumably due to the additional invariant in 2DMHD, namely the mean squared magnetic potential which is not conserved for 3DMHD, even in the presence of a strong mean field. This extra invariant puts additional constraints on the dynamics of decay of energy-containing eddies in 2DMHD. In particular, with regard to the presence of a second ideal invariant in 2D, we note that the mean-square potential is expected to be back-transfered to longer wavelengths during the cascade. This induces the growth of large magnetic islands, which contribute at most very weakly to the direct cascade. Therefore, one might view the standard estimate of the similarity scale (i.e., energy-containing scale), given in Eq. (15), as somewhat of an overestimate. This may point the way to understanding the overestimate given by Eq. (15) when applied to the 2D case.
To investigate the issue of reduction of dimensionality further, we also employ a set of reduced MHD (RMHD) simulations. RMHD is often used as a simpler model of plasma in the presence of an energetically strong mean magnetic field (see [45] for governing equations and other discussions). Again, the RMHD runs have almost identical conditions as the 3DMHD ones except that here, the measurements have been made near the instant of maximum dissipation associated with the k = 1 modes. Here, we set B 0 = 1, noting that this is in rescaled "code" units. In unscaled units this corresponds to physical situations with δb ≪ B 0 , as must be the case for RMHD [e.g., 45]. Fig. 1 shows that, contrary to the 2DMHD results, the RMHD values almost superimpose on the 3DMHD values for high mean field cases. These results show that, at least for the problem of global dissipation, RMHD may be a more realistic approximation to the full three-dimensional system in the presence of a strong mean magnetic field than 2DMHD is.
An additional view of the approach to asymptotic values of the dissipation rate is provided in Fig. 2. Panel (a) shows the empirical estimate of the asymptotic decay rate C ǫ,∞ plotted against the mean field strength B 0 , and panel (b) shows C ǫ,∞ plotted versus 1/ √ r = B 0 /δb. Still another view is provided by examining how the asymptotic decay rate varies with spectral anisotropy measured by the value of the Shebalin angle θ b at the time of peak dissipation; this is shown in Figure 2, panel (c). Here we see that the anisotropy saturates at large mean field values, so that the subsequent values of dissipation rate cluster near the associated saturated Shebalin angle. Since the values of the Shebalin angles also loosely depend on the Reynolds number, the values of θ b corresponding to ν = 0.0002 are considered here.
In Fig. 2 we also show the values of C ǫ,∞ , extrapolated from the 2DMHD and RMHD runs. The values of δb/B 0 and the ratio of perpendicular to parallel characteristic length scales ℓ ⊥ /ℓ are small but somewhat arbitrary for RMHD and there are an infinite choice of values. The 2DMHD values are only for comparison. Therefore, we have plotted the C ǫ,∞ corresponding to the 2DMHD and RMHD simulations at the right edge of the horizontal axes of the three subplots, arbitrarily.

V. SUMMARY
We have performed a quantitative measurement of variation of the non-dimensional dissipation rate in MHD turbulence in presence of a mean magnetic field using Fourier pseudo-spectral simulations. We find that the dissipation rate approaches a nonzero asymptotic value for increasing Reynolds numbers (mechanical and magnetic), and for increasing mean DC magnetic field strength. This confirms the generalizations of the von Kármán-Howarth theory of hydrodynamics to the case of magnetohydrodynamics [35,36]. This conclusion provides essential confirmation of the underlying theory, normally assumed, upon which research is based in several key areas of space and astrophysics. Two prominent examples are the derivation and use of so-called "third-order laws" or Yaglom relations, such as the well studied Politano-Pouquet relations [23], and the estimation of decay rates in turbulence transport theories [49]. Generalized Yaglom laws have been widely employed in studies of interplanetary turbulence [50][51][52], while turbulence-based global models have found important applications (as well as agreement with observations) in simulation of the outer heliospheric plasma [49,53,54], and in the inner solar wind and corona [30].
We note that the saturation of dimensionless dissipation rate is obtained for low values of δb/B 0 where weak turbulence could become important in certain circumstances [46,47]. For weak turbulence, the leading order behavior is that of waves, and energy transfer across scales is achieved through resonant interactions among wave modes, mediated by interaction with 2D (nonpropagating) modes or quasi-2D modes. Although the present simulation results, including the RMHD cases, fall within the strong-turbulence regime (due to nonpropagating fluctuations, see [48]), one might expect a similar dissipation anomaly in weak turbulence. However, we have not examined the weak turbulence regime here.
Clearly there is scope for application of the present findings to numerous astrophysical plasmas. As discussed before, the Reynolds number in these systems are very large and often these systems have a strong mean mag-netic field along with the turbulent component. For example, the ratio of the rms fluctuation to the mean magnetic field in the solar wind at 1 AU is ∼ 0.5. Other astrophysical systems can be even more anisotropic. The solar corona has an even stronger mean magnetic field, which makes the plasma highly anisotropic in those regions. The Parker Solar Probe (PSP) Mission, which was launched on August 12, 2018, will make in-situ measurements close to the Sun, in the solar corona. The results presented in this paper will be helpful to investigate PSP data and for modeling the solar corona.
The analysis presented here also enables one to directly probe the validity of the dissipative anomaly, that is, the assumption of finite dissipation at (approaching) infinite Reynolds number for anisotropic MHD. The same may not be true for pure hydrodynamic turbulence [17], since two dimensional hydrodynamic turbulence, unlike its MHD counterpart, does not maintain finite dissipation in the limit of vanishing viscosity.
Beyond the explanation based on anisotropic spectral transfer given above, the saturation of dissipation rate may also be related to the reconnection rate in the system, especially since magnetic reconnection provides an efficient mechanism of dissipation. We note that as the mean magnetic field becomes stronger, the system becomes quasi-two dimensional and the islands may begin to form thin current sheets, facilitating reconnection in the classical sense. It would be interesting to investigate further how the statistics of the reconnection rates change as the strength of the mean field increases (see [55]).
It may be interesting to perform a study similar to the current one for other systems that become anisotropic. One example is rotating turbulence. It has been shown that under rotation, a turbulent fluid system becomes anisotropic and approaches a two dimensional state, sim-ilar to anisotropic MHD [56].
Another important direction of further investigation is eddy viscosity which is often used in global and homogeneous turbulence simulations as an approximation to represent effects of unresolved scales. In order for the results of an eddy viscosity or any other such models to be physical, it is important that the results presented in this paper are maintained in such approximations, at least qualitatively. We plan to address these problems using Large-Eddy Simulations (LES) in future.