Variable Energy Fluxes and Exact Relations in Magnetohydrodynamics Turbulence

: In magnetohydrodynamics (MHD), there is a transfer of energy from the velocity ﬁeld to the magnetic ﬁeld in the inertial range itself. As a result, the inertial-range energy ﬂuxes of velocity and magnetic ﬁelds exhibit signiﬁcant variations. Still, these variable energy ﬂuxes satisfy several exact relations due to conservation of energy. In this paper, using numerical simulations, we quantify the variable energy ﬂuxes of MHD turbulence, as well as verify several exact relations. We also study the energy ﬂuxes of Elsässer variables that are constant in the inertial range.


Introduction
Magnetohydrodynamics (MHD) provides a framework to study the dynamics of flows in interstellar medium, galaxies, accretion disks, stars and planet interiors, solar wind, Tokamak, etc. [1]. In these systems, typically, the kinetic and magnetic Reynolds numbers are significantly large; hence, they exhibit turbulent behavior. The three-dimensional (3D) MHD turbulence and the 3D hydrodynamics turbulence have several common features. For example, both these systems exhibit nonlinear interactions among multiple scales. In addition, the energy injected at large scales cascades to intermediate and then to small scales. The dissipative terms destroy the fluid energy at small scales [2]. This cascade is quantified by energy flux.
Since inertial range of hydrodynamic turbulence does not have any forcing (except very weak viscous force), the energy flux in the inertial range remains constant [2][3][4]. However, MHD turbulence has a magnetic field in addition to the velocity field, and there are energy exchanges among the two fields via nonlinear interactions. Consequently, there are many energy fluxes, e.g., velocity to velocity, magnetic to magnetic and velocity to magnetic. Dar et al. [5] and Verma [6] showed that MHD turbulence has six energy fluxes related to the velocity and magnetic fields. Due to the energy exchange between the velocity and magnetic fields, the energy flux for the velocity field is no longer constant. Verma et al. [7] showed that in a typical turbulent magnetofluid, the inertial-range kinetic energy flux is depleted due to the energy transfer from the velocity to the magnetic field. The magnetic energy flux also exhibits variability due to these transfers. MHD turbulence is also described in terms of Elsässer variables, which are z ± = u ± b. The fluxes associated with the Elsässer variables are constant in the inertial range [8,9]. This is because there is no energy exchange between the z + and z − variables. The mean magnetic field also affects the energy fluxes. In this paper, we study the energy fluxes in the presence of nonzero cross helicity, but no magnetic helicity and mean magnetic field.
There are several models for MHD turbulence that describe the associated energy spectra and fluxes. Kraichnam [10], Iroshnikov [11] and Dobrowolny et al. [12] assumed the turbulence to be homogeneous and isotropic and argued that the kinetic and magnetic energy spectra scale as E(k) ≈ ( B 0 ) 1/3 k −3/2 , where B 0 is the magnitude of the mean magnetic field or of the large-scale magnetic field. Marsch [13] argued that in the absence of a mean magnetic field, the energy spectra of Elsässer variables exhibit a k −5/3 spectra. Using renormalization group theory, Verma [6,14] argued that the effective mean magnetic field scales as k −1/3 , and hence the energy spectra of kinetic and magnetic fields scale as E(k) ≈ k −5/3 . Goldreich and Sridhar [15] argued that in the presence of a strong external magnetic field, the flow becomes anisotropic and that E(k ⊥ ) ≈ k −5/3 ⊥ , where k ⊥ is the wavenumber perpendicular to the mean magnetic field. In this paper, we do not discuss the energy spectrum of MHD turbulence in detail, but focus on the energy fluxes.
Some of the important past works related to energy fluxes of MHD turbulence are as follows. Dar et al. [5] formulated the energy fluxes for MHD turbulence in terms of mode-to-mode energy transfers and computed the fluxes for two-dimensional (2D) MHD turbulence. Debliquy et al. [16], Alexakis et al. [17] and Carati et al. [18] computed these fluxes, as well as shell-to-shell energy transfers, for 3D MHD turbulence. Note that the simulations of Debliquy et al. [16] are for decaying turbulence. Using numerical simulations, Verma et al. [19] computed the energy fluxes of z ± and showed consistency with the predictions of Marsch [13]. In the presence of mean magnetic field, Teaca et al. [20] and Sundar et al. [21] computed the corresponding energy fluxes. Verma [6,22,23,24] computed the energy fluxes using field-theoretic formalism. In addition, Verma [9] also described several exact relations among these fluxes using the analytical formalism of variable energy flux.
Solar wind provides an important platform to test the theories of MHD turbulence. Matthaeus and Goldstein [25] and Tu and Marsch [26] analyzed the solar wind data and observed a near k −5/3 energy spectra for u, b, z ± fields. Parashar et al. [27] studied the variations of spectral indices as a function of cross helicity. Verma et al. [28] estimated the energy flux in the solar wind using the energy spectrum and assuming Kolmogorov-like turbulence phenomenology for the MHD turbulence [13]. A more commonly-used strategy for the estimation of energy flux in the solar wind is to employ structure functions and the four-third rule [29]; Sorriso-Valvo et al. [30] employed this method. Following a similar procedure, Bandyopadhyay et al. [31] computed the energy flux near the Sun and argued that energy flux is enhanced here.
The effects of the magnetic field on the dynamics of peristaltic and nanofluid flows have been extensively studied in literature [32][33][34][35][36][37]. Eldesoky et al. [34] investigated the combined effects of the magnetic field and heat transport in peristaltic flow, and showed that the magnetic field enhances the thermal energy of the fluid. It has been argued that the magnetic field also affects the dynamics of biological nanofluids, which could have practical implications. Abdlesalam and Sohail [36] showed that the velocity distribution of the bioconvective flow of viscous fluid reduces with an enhancement of the magnetic field. An analytical investigation of Abdelsalam, Velasco-Hernández and Zaher [37] showed that the propulsive velocity of swimming sperms increases with the Hartmann number, a measure of the strength of the magnetic field. Apart from the effect of the magnetic field, the dynamics of peristaltic flows and nanofluid flow are also affected by the particle concentration in the flow, viscosity of the flow and temperature [38][39][40].
In MHD turbulent flow, the flow properties are usually studied using the Newtonian model, where the relation between stress and strain rate is linear [17,18,41]. However, there are flows in nature, where the relation of stress and strain rate are nonlinear. Such flows are studied using non-Newtonian models like Williamson flow, Casson Carreau and Jeffrey flow, etc. The dynamics of non-Newtonian flows and the flows through the porous medium like Darcy-Forchheimer flow are also affected by the magnetic field and have been broadly investigated in the literature [42][43][44][45][46][47][48]. A numerical investigation of MHD Williamson nanofluid by Rasool et al. showed that the drag force exerted by the medium on the flow increases with the magnetic parameter, a measure of the strength of magnetic field [43]. Rasool et al. [48] investigated the dynamics of convective MHD nanofluid flow bounded by non-linear stretching surface and reported a decrease in flow velocity with an increase in magnetic parameter. Ali et al. [45] reported a decrease of Nusselt number with an increase in magnetic parameter for Darcy-Forchheimer rotating flow of a Casson Carreau nanofluid.
In this paper, we compute the various energy fluxes of forced MHD turbulence and validate several exact relations with numerical results. We perform a direct numerical simulation of forced MHD turbulence on a 512 3 grid with kinetic hyperviscosity and magnetic hyperdiffusivity. We set random initial conditions for both the velocity and magnetic fields, and inject kinetic energy at wavenumber shell (2, 3). We observe that the numerical results satisfy several exact relations of the fluxes in MHD turbulence. We also find that the fluxes of the Elsässer variables are constant in the inertial range.
The outline of this paper is as follows. In Section 2, we present various energy fluxes of MHD turbulence and the exact relations relating them. In Sections 3 and 4, we present the details of our numerical simulation and verification of exact results. Finally, we conclude in Section 5.

Energy Fluxes and Exact Relations
The equations for incompressible MHD turbulence in the absence of a mean magnetic field are [49] ∂u ∂t where u and b are respectively the velocity and magnetic fields, p is the total pressure (a sum of kinetic and magnetic pressures); ν 0 and η 0 are respectively the viscosity and diffusivity. The random large-scale force term is F ext , and represent respectively the Lorentz force and the stretching of the magnetic field by the velocity field. Note that the magnetic field b is normalized using b = b cgs / 4πρ to convert it in the units of velocity; here, b cgs is the magnetic field in cgs units, and ρ is the fluid density. The scale-dependent properties of a system are customarily studied using Fourier transforms. In Fourier space, the evolution equations for the kinetic and magnetic modal energies are [8] where u(k) and b(k) are respectively the Fourier transforms of the velocity and magnetic fields, E u (k) = |u(k) 2 |/2 and E b (k) = |b(k) 2 |/2 are respectively the modal kinetic and magnetic energies, T uu (k, t) and T bb (k, t) are nonlinear energy transfers arising due to (u · ∇)u and (u · ∇)b, respectively, D u (k) and D b (k) are respectively the kinetic and magnetic energy dissipation rates, F ext (k, t) is the energy injection rate due the external force and F ub (k) and F bu (k) are the cross energy transfers, i.e., from magnetic to kinetic and vice versa. The formulas for the above transfers are as follows [8]: where , and * represent the imaginary and real parts, and the conjugate of a complex number, respectively. In addition, q = k − p.
The kinetic and magnetic energy spectra are computed using Apart from the kinetic and magnetic energy spectra, we also have energy spectra related to the Elsässer variables: where z ± (k) = u(k) ± b(k) are the Elsässer fields. In MHD turbulence, there are four nonlinear terms in the governing equations. These terms yield six energy fluxes which are illustrated in Figure 1b and are defined as follows [5,6,8]: where the superscript and subscript of Π represent respectively the giver and receiver modes, and < and > denote respectively the modes inside and outside the sphere of radius k 0 . For example, Π X < Y > (k 0 ) denotes the rate of energy transfer from the wavenumbers inside the sphere of radius k 0 of field X to the wavenumbers outside the sphere of field Y. It is Note that there are energy exchanges among the velocity field and the magnetic field. However, there is no such cross transfer between z + and z − . For these variables, the corresponding energy fluxes are The total energy flux in MHD turbulence, which is a sum of energy transfers from the velocity and magnetic modes inside the sphere to the modes outside the sphere, is It can be easily shown that 2Π total (k 0 ) = Π In 3D hydrodynamic turbulence, the kinetic energy flux Π u < u > remains constant in inertial range and it is equal to the kinetic energy dissipation rate. This flux, however, is not constant in MHD turbulence [7]. However, the total energy transferred to the inertial range of the velocity field is dissipated by the viscous force. This leads to where u is the kinetic energy dissipation rate. Refer to Figure 1b for an illustration. In the above expression, is the total energy transferred from the magnetic field to the velocity field of the inertial range. Note that each of the above fluxes, Π u < u > , Π b < u > and Π b > u > , vary with k, but the sum is an approximate constant, hence, we call them variable energy fluxes [8,9].
In a similar vein, we derive the following relation for the energy flowing into the magnetic-field channel: where b is the magnetic energy dissipation rate. In addition, the energy injected into the large-scale velocity field is transferred to the inertial-range velocity and magnetic fields. These fluxes are dissipated at small scales. Consequently, where is the total dissipation rate, which is equal to the total-energy injection rate. Similar balance of energy transfers in the magnetic channel leads to Schematic diagram of Figure 1b helps us understand the above relations. The above four equations represent exact relations. The equality holds statistically in spite of the fact that each of the fluxes exhibits significant variations with k. In the next section, we validate the above relations using numerical simulations.

Governing Equations and Simulation Method
We employ a pseudo-spectral code named TARANG [50,51] to solve the following equations in a 3D cubic periodic box of size (2π) 3 (as shown in Figure 1a): In the above equations, we employ the hyperviscous and hyperdiffusive terms to increase the inertial range and to suppress the dissipation range. In addition, the velocity, length and time are non-dimensionalized using characteristic velocity (U 0 ), box size (2π), and the eddy turn over time (2π/U 0 ). We use a fourth-order Runge-Kutta scheme for time marching and the Courant-Friedrich-Lewis (CFL) condition for optimizing the time step ∆t. We employ 2/3 rule for idealizing.
We performed a numerical simulation on a 512 3 grid with the hyperviscous and hyperdiffusion parameters ν = η = 3 × 10 −7 . We use the Craya-Herring basis [52,53] to generate the random initial conditions for both velocity and magnetic field. In a Craya-Herring basis, the three basis vectors are, wheren is chosen as any arbitrary direction, andk is the unit vector for wavenumber k.
The velocity field for 3D incompressbile flow in the Craya-Herring basis is written as We start our simulation with the following initial condition: u 2 (k) = (E u /2N 3 ) (exp(iφ 1 (k)) + exp(iφ 2 (k))), (37) where N 3 is the total number of modes, E u = 0.5 is the total kinetic energy and the phases φ 1 (k) and φ 2 (k) are chosen randomly in the band [0, 2π]. Then we transform the velocity field from Craya-Herring basis to Cartesian basis. We use a similar approach to generate the random initial condition for magnetic field with total initial magnetic energy set as 0.25.
Note that the mean magnetic field in our simulation is zero, i.e., B 0 = 0. We define Reynolds number Re = u rms L 3 /ν, and magnetic Reynolds number Rm = u rms L 3 /η, where L = 2π is the box size [54]. The L 3 factor arises due to the ∇ 4 factor of the hyperdiffusion term. For the steady state, both these numbers are 8.2 × 10 8 . However, we caution that they are not reliable estimates of Reynolds number and magnetic Reynolds number, and they should not to be compared with earlier simulations.
We apply random force to all the velocity modes in a wavenumber shell (2, 3), which is denoted by k f = 2. We follow the procedure outlined in Alvelius [55] and Maffioli [56]. Note that we force only the velocity modes. The kinetic energy injection rate is 0.4. The simulation parameters are tabulated in Table 1. In addition, the simulation workload is summarized in Table S1 and Figure S1 of the Supplementary Material. In Figure 2a, we exhibit the time series of the total kinetic energy (solid red curve), the total magnetic energy (solid green curve) and the total energy (solid blue curve). Figure 2b exhibits the corresponding dissipation rates. Clearly, the system reaches a steady state in 2 to 3 eddy turnover time. During the steady state, the magnetic energy dissipation rate is larger than the kinetic energy dissipation rate and the total dissipation rate matches with the energy injection rate (Figure 2b). Interestingly, the kinetic energy and the magnetic energy are equipartitioned during the steady state.
In the next section we will describe the numerical results related to the energy fluxes.

Numerical Results on the Energy Fluxes
Before embarking on energy flux studies, we compute the spectra of the velocity and magnetic fields. We observe that the spectra show significant fluctuations; hence, to obtain smooth curves in the inertial range, we compute the cumulative spectra, which is Note that if the inertial-range energy spectrum scales as k α , then the cumulative spectrum would scale as k α+1 . In Figure 3, we exhibit the cumulative spectra of the velocity, magnetic and z ± variables averaged over time interval 9 to 12 (in non-dimensional time units). Figure 3. Time-averaged cumulative spectra of (a) the kinetic energy (solid red curve) and magnetic energy (solid blue curve), and (b) z + (solid red curve) and z − (solid blue curve). These spectra are close to k −2/3 indicating consistency with the Kolmogorov-like spectra for MHD turbulence. Figure 3, in the wavenumber range of (3,20), the cumulative spectra follow power laws to a reasonable accuracy. The exponents of the power law for u, b, z + , z − variables are −0.73 ± 0.01, −0.68 ± 0.01, −0.71 ± 0.03 and −0.66 ± 0.01 respectively, which are close to −2/3. Hence, we conclude that energy spectra for these fields are reasonably close to Kolmogorov's k −5/3 power law. The deviations of the exponents from −5/3 may be related to the variable energy flux and intermittency [2,8]. Note that the kinetic energy spectrum, ∼k −0.73 , is steeper than the magnetic energy spectrum, ∼k −0.68 . This differential can be attributed to the energy transfer from the velocity field to the magnetic field [7]. We have compared our results with those from the past, e.g., observational works on solar wind [25,26] and numerical works [57][58][59][60]. We observed general consistency between these results. Note, however, that several authors report k −3/2 energy spectra, but a detailed discussion on this topic will take us beyond the scope of the paper. We also remark that higher-resolution simulations, which are planned in future, would provide a better handle on the spectral exponents.

As shown in
Following the main objectives of the paper, we compute various energy fluxes of MHD turbulence for 80 concentric spheres. The radii of the first 16 spheres are linearly binned as [1, 2, 3,. . . , 16], and those of the last two spheres are 128 and 256. The radii of the intermediate spheres are binned as follows: where r 16 = 16, s = log 2 (r max /32)/(n − 5), with r max being the radius of the largest wavenumber sphere, and n as the total number of spheres.
In Figure 4, we exhibit the plots of various energy fluxes, which are averaged over a time interval of t = 9 to 12. The black solid curve represents the total-energy injection rate. These plots reveal the following interesting properties of the energy fluxes of MHD turbulence: 1.
The energy fluxes corresponding to the total energy and z ± are nearly constant in the wavenumber band (3,20), consistent with the power-law regime of the energy spectra discussed earlier. Note that the inertial-range energy flux of the total energy matches with [Π z + + Π z − ]/2; in addition, these fluxes are equal to the energy supply rate and the total-energy dissipation rate, consistent with the conservation of energy.

2.
As shown in Figure 4a, in the wavenumber band k = (3, 6), the kinetic energy flux dips sharply, while the magnetic energy fluxes, Π u < b > and Π u < b< , grow rapidly. This observation indicates energy transfer from u to b. Note that Π b < b > picks up significantly after this band (k = (3, 6)).

3.
The energy fluxes Π u > b < and Π u > b > are negative and become significant beyond wavenumber band (3,6). These fluxes indicate energy transfers from the magnetic field to the intermediate-scale velocity field. Consequently, Π u < u > grows and becomes significant beyond k = 10.

4.
The energy fluxes corresponding to the velocity and magnetic fields exhibit significant variability due to cross energy transfers. However, the fluxes of z ± are nearly constant in the inertial range due to lack of such transfers. We also compute the flux of cross helicity, which is Π Hc (k) = (Π z + (k) − Π z − (k))/4, and exhibit this flux in Figure 4b. We need to further explore the evolution of cross helicity flux in MHD turbulence.
and Π total , and (b) Π z + , Π z − and Π Hc averaged over time frame t = 9 − 12. Note, the solid black curve shows the total-energy injection rate.
We compare the results of energy fluxes of MHD turbulence with the earlier works reported in literature [17,18,41]. Our results for the transfer of energy between magnetic field, and the transfer of energy from the velocity field to the magnetic field are consistent with those of Alexakis, Mininni and Pouquet [17], Carati et al. [18] and Alexakis, Mininni and Pouquet [41].
Finally, we come to the exact relations for the energy fluxes of MHD turbulence. The four subfigures of Figure 5 provide numerical verification of Equations (26)-(29) as follows: 1. Figure 5a demonstrates that in the inertial range, the sum Π u < u > + Π b < u > + Π b > u > matches with the kinetic energy dissipation rate u . Note that the sum represents the total energy transfer to the inertial-range velocity modes that gets dissipated in the dissipation range; this is the reason for the equality of Equation (26). 2. Figure 5b illustrates a similar balance between the energy transfer to the inertial-range magnetic modes and the magnetic-energy dissipation rate b (see Equation (27)).

3.
The energy supplied to the large-scale velocity modes get transferred to the inertialrange velocity and magnetic modes. It leads to the exact relation of Equation (28). This relation is verified in Figure 5c.

4.
The magnetic field is not forced externally. Instead, the large-scale magnetic modes (b < ) receive energy from the velocity modes as The energy received by the large-scale magnetic modes cascades to the inertial range of the magnetic field as This relation is verified for the wavenumber range k = (3,18). See Figure 5d for an illustration.
Thus, we verify several important exact relations regarding the energy fluxes of MHD turbulence.

Conclusions
In MHD turbulence, there are complex energy transfers among the velocity and magnetic fields. The energy fluxes of MHD turbulence provide a measure for these transfers. In the past, these fluxes have been computed using numerical simulations [5,[16][17][18]. In this paper, we describe the subtle variations of these fluxes in the framework of variable energy flux. The energy fluxes related to the velocity and magnetic fields vary with wavenumber. However, several combinations of these fluxes are constant; these are the exact relations related to the energy fluxes.
In this paper, we describe variable energy fluxes and exact relations of MHD turbulence using numerical simulations. We summarize our results as follows: 1.
Our work is focused on the energy fluxes of forced MHD turbulence, in contrast to those of decaying MHD turbulence, studied earlier by Debliquy et al. [16]. A close comparison between the two sets of energy fluxes shows that the decaying and the forced MHD turbulence have several critical differences. For example, we observe positive Π u < b < , while Debliquy et al. [16] reported negative Π u < b < .

2.
We employed hyperviscous and hyperdiffusive terms in our simulation to increase the extent of the inertial range. For our simulation, the flux for the total energy is nearly constant in the inertial range, which is k = (3,20). The extent of our inertial range is larger than that of Debliquy et al. [16], who do not employ hyperdiffusion.

3.
For our numerical simulations, the spectral indices for u, b and z ± are close to −5/3, rather than −3/2 [10]. A word of caution, however, is that the inertial range is quite narrow due to the moderate resolution (512 3 ) of our simulation. For a better understanding of the spectral indices and the energy fluxes, we need a broader inertial range that is possible with high-resolution simulations; we plan for such simulations in the near future.

4.
It has been recently reported that the energy fluxes of MHD turbulence satisfy several exact relations [8,9]. These relations are based on energy conservation principles. In this paper, we validate four such exact relations, which are Equations (26)- (29).
In summary, we employ numerical simulations to understand the variations of inertialrange energy fluxes of MHD turbulence. We also demonstrate several exact relations related to energy fluxes. This study provides valuable insights into the dynamics of MHD turbulence.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/fluids6060225/s1. Table S1: Workload of the simulation: the grid-resolution N, the total number of processor p , data division along x direction p col , data division along y direction p row , total RAM used in GB M RAM , simulation time t (eddy turnover time) and simulation ran time in hours T. Figure S1: Schematic diagram for the pencil decomposition of a 3D array. From Chatterjee et al. [51]. Reprinted with permission from Elsevier.  Institutional Review Board Statement: This paper has not been reviewed by any institutional board.
Informed Consent Statement: The paper does not employ any clinical research or personal data, hence no consent is required.

Data Availability Statement:
The numerical data will be made available up on request.

Acknowledgments: The authors thank Franck Plunian and Rodion
Stepanov for valuable discussions. Parts of the simulations were run on IIT Kanpur's HPC clusters and Shaheen II of KAUST supercomputing laboratory, Saudi Arabia (Project No. k1416).

Conflicts of Interest:
The authors have no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

MHD Magnetohydrodynamics 3D
Three dimensions 2D Two dimensions Nomenclature k, p, q Fourier wavenumbers u Velocity field b Magnetic field z ± Elsässer variables t Time p Pressure F ext Random large-scale force F u Lorentz force F b Stretching of magnetic field by velocity field E u (k) Modal kinetic energy E b (k) Modal magnetic energy E u (k) Kinetic energy spectrum E b (k) Magnetic energy spectrum E z ± (k) Elsässer energy spectra T uu (k), T bb (k) Nonlinear modal energy transfers D u (k) Modal kinetic energy dissipation rate D b (k) Modal magnetic energy dissipation rate F ext (k) External modal energy injection rate F ub (k), F bu (k) Cross energy transfers among velocity and magnetic modes u Kinetic energy dissipation rate b Magnetic energy dissipation rate Total dissipation rate inj Kinetic energy injection rate Π X < Y > Energy flux from the wavenumbers inside the sphere of radius k 0 of field X to outside the sphere of field of field Y, e.g., X = Y = u, b, z ± (ê 1 ,ê 2 ,ê 3