{'fraction_non_alphanumeric': 0.033503277494537506, 'fraction_numerical': 0.003941562058180884, 'mean_word_length': 4.443563432835821, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Controlling the level of the ideal invariant fluxes for MHD turbulence using TURBO spectral solver

The ideal invariants present in the formalism of magnetohydrodynamics (MHD), i.e. global quantities that are conserved in the absence of sources and dissipative effects, play an important role in various theoretical and numerical studies of MHD turbulence. The fluxes of these ideal invariants represent separate channels that transfer the information across different scales in a turbulent system. Once a statistically stationary state of turbulence is reached, the amount of any ideal invariant quantity introduced in the system by a forcing mechanism equals the amount of the same quantity removed by the dissipative effects from the system. For highly developed turbulence, these two mechanisms act predominantly at different scales that are largely separated. Since the ideal invariant quantities cascade between scales, a constant flux is generated with great implication on the state of the system. Numerically, controlling the ideal invariant fluxes levels for a turbulent MHD system is important for the analysis of fundamental MHD turbulence properties. We propose a forcing mechanism that controls the three ideal invariants of MHD turbulence: the total energy, the cross-helicity and the magnetic helicity. This forcing is implemented in the freely available TURBO solver, that is also briefly presented.


I. INTRODUCTION
The motion of a fluid is described mathematically by a nonlinear evolution equation, for which no general analytical solution exists. Except for laminar flows in very simple geometries, even the use of a perturbative approach is impossible. This known problem, caused by the coupling of different scales in a flow by the nonlinear terms has triggered the use of numerical methods in the study of fluid flows. Usually, solutions of a flow are obtained using a numerical solver, in an incremental manner, starting from a set of initial conditions and a prescribed set of boundary conditions. For a typical engineering problem, the correct description of the physical properties and geometry of the boundary conditions is crucial, as they are responsible for an entire class of instabilities introduced in the flow which in turn characterizes the subsequent motion. However, for the study of fully developed turbulence, where a huge range of scales are present and a clear separation exists between the geometry dependent large scales and the universal, dissipative small scales, periodic boundary conditions represent the preferred choice. In this situation, the driving instabilities are introduced in the form of an external force that acts only at a particular scale, usually a large one. The use of periodic boundaries conditions for a flow has the advantage of enabling us to easily translate the motion problem in terms of Fourier modes which can be solved numerically using spectral solvers. For this purpose we make use of the TURBO code, which we will briefly introduce in Section II of this paper.
When the fluid is electrically conducting like in the case of a plasma or of a liquid metal, the momentum balance equation is influenced by the Lorentz force and the number of non-linear terms increases compared to the case of fluid turbulence. This coupling of the Navier-Stokes equation with the Maxwell reduced equations, which leads to the magnetohydrodynamic (MHD) equations, increases the complexity of the interaction between scales. Compared to hydrodynamical flows, the ideal invariants change in the case of MHD, as the interplay between the velocity and the magnetic fields needs to be considered. The importance of the MHD ideal invariants will be discussed in Section III.
MHD turbulence, considered as an initial value problem with periodic boundary conditions, tends to decay in absence of driving mechanism. To maintain a stationary turbulent state, it is needed to add sources in the evolution equations that mimic the effects of the various instabilities that may appear in the large, geometry dependent scales of a realistic system. One way of achieving this would be through a constrain on the large scale value of one or both fields. This is often implemented in spectral space by freezing the value of a limited set of Fourier modes characterizing the largest scales in the system [10,12]. Alternatively, an external artificial force can be used [1,15]. We propose such a force in Section IV which controls the injection rate of the three ideal invariants of MHD turbulence: the total energy, the cross-helicity and the magnetic helicity.

II. TURBO SOLVER
The TURBO code 1 is designed to solve the equations for an incompressible fluid in a three dimensional slab geometry with periodic boundary conditions in the three directions. The real space representation of the incompressible MHD equations under the influence of an external force and in the presence of a constant magnetic field B 0 are written as, where u = u(x, t) is the fluid velocity field, b = b(x, t) is the magnetic field expressed in Alfvèn units and p = p(x, t) is the total, hydrodynamic + magnetic, pressure field divided by the constant mass density. Due to the incompressibility condition, the pressure p is not an independent variable and can be formally eliminated by solving the Poisson equation, The kinematic fluid viscosity is ν and the magnetic diffusivity is η. The divergence free external force fields f u = f u (x, t) and f b = f b (x, t) act on the velocity and magnetic fields, respectively. The two forces are part of a forcing mechanism that imposes the injection rates of the MHD ideal invariant quantities. A kinetic only forcing method (f b ≡ 0), used previously in the literature for similar studies [7,15], can also be employed. In the velocity evolution equation, f u can be considered as divergence free since the pressure will enforce the incompressibility of the velocity field by eliminating any ∇ · f u contribution of the force. On the other hand, f b must always be divergence free as a consistency condition for the magnetic field. If desired, a Coriolis force f C = Ω × u acting on the flow can also be considered. The force appears as result of a reference system rotation, with the angular velocity Ω. In that case, the centrifugal acceleration that depends explicitly on the distance to the rotation axis can be lumped into the pressure term due to the incompressibility condition and the use of periodic boundary conditions is still appropriate. In addition to the MHD equation, the TURBO code can also solve the evolution equations for a set passive scalars fields: where c α = c α (x, t) are passive scalar(s), each of which is characterized by a diffusion coefficient κ α . There is the possibility to include source or sink terms or even chemistry terms in the scalar equations through the function σ α . In the TURBO code, the space discretization is based on a Fourier representation of the quantities of interest. For a given quantity, the physical Q and the spectralQ representations are related using the direct and the inverse discrete Fourier transforms 2 : where N is the total number of modes in a given direction. In practice, the code allows the use of different numbers of modes in the three directions N x , N y and N z and different box sizes L x , L y and L z . Knowing the number of modes and the box size for each direction allows us to define the wavenumber space. Assuming a cubic box for simplicity of notations, the wavenumbers are defined as: where n ∈ [−N/2+1, N/2]. From the above definition, we see that the smallest nonzero wavenumber is k 0 = 2π/L, which for the typical box length choice L = 2π, becomes unity. The largest wavenumber accounted in a simulation (k max = πN/L) and the spacing between two neighbour grid points (∆ = L/N ) are related as ∆ = π/k max . The main advantages of spectral methods are definitely the accuracy and the simplicity of the representation of the differentiation operator ∇, which reduces to a multiplication in spectral space, ∇ u(x) → i kû(k). The weak point of this method is the mandatory use of periodic boundary conditions, a choice considered as hard coded in the TURBO solver. Spectral methods are thus not adequate for exploring very complex geometries but are extremely useful for investigating the fundamental properties of turbulence.
The time evolution is based on a modified Williamson, four-step, third-order low storage Runge-Kutta method [18]. Since the equations are discretized in space, it is desirable that the transfer of information be limited to neighbor grid point in one time step ∆t. This implies the CLF criterium which simply states that the time step ∆t has to be smaller than the time necessary for the wave with the largest propagation speed (determined on the velocity u max or on the Alfvén velocity b max ) to propagate between the smallest distance present between two grid points (∆). The linear terms are solved in an analytical manner for each mode k by performing an appropriate change of variableŝ u →û (û, k, ν, Ω) 3 andb →b (b, k, η) before formally performing the Taylor series expansion. Although this will still give us a third-order accuracy in a global sense for the solution, the linear terms do not affect the value of the time step.
Since the evolution equations are solved is spectral space, the value of the nonlinear terms for a mode k has to be computed. Determining the nonlinear terms directly in Fourier space, although possible, is prohibitive from a computational cost point of view, as they are represented by convolutions of modes. Instead, inverse Fourier transforms are used for the velocity and the magnetic field and the nonlinear terms are computed in real space before being Fourier transformed back to the spectral representation. In practice, the nonlinear terms appearing in the evolution equation for a velocity mode correspond to the divergence of a symmetric tensor u i u j − b i b j , while the nonlinear term appearing in the evolution equation for a magnetic field mode appears as the divergence of an anti-symmetric tensor This has an impact on the number of Fourier transforms that need to be perform: six for the symmetric tensor and only three for the anti-symmetric tensor, instead of nine transforms for each term. Furthermore, the trace of the tensor u i u j − b i b j can be lumped into the pressure term, leading to another reduction of the computational efforts. Two orthogonal projections with respect to k of the resulting nonlinear terms ensure that the solenoidal conditions are enforced for both the velocity and magnetic fields.
Computing the nonlinear terms in real space reduces the numerical method type from a purely spectral method to a pseudo-spectral one. This approach of computing the nonlinear terms, although much faster then computing the convolution, adds another problem known as aliasing. The aliasing error originates in the fact that a plane wave with a wavenumber k takes exactly the same values on the grid as a plane wave with a wavenumber k + N k 0 . When the nonlinearities are computed in real space, the aliasing error becomes a serious issue. For a one dimensional quantity q represented by modes with As a consequence, the mode of q 2 that corresponds to k = k 0 (−N + 2) is undistinguishable from the mode k = 2k 0 . Two approaches, known as dealiasing methods, are considered to eliminate this difficulty [14].
The first method consists in assuming that only the modes of q with k = k 0 × [(−M/2 + 1), M/2] are non zero and imposing that the other modes remain zero after each nonlinear computation. The two-third method consists simply taking M = 2/3 N and keeping all the modes outside the range k = k 0 × [(−N/3 + 1), N/3] to zero. This method represents the simplest way of fully removing the aliasing error.
The second method is based on the property that shifting the grid by a distance d results in a modification for the modes by a phase e i k d . For one-dimensional systems, the aliasing error can be removed exactly by computing the nonlinearity twice on two grids shifted respectively by ∆/2 and −∆/2 and by summing the two computations. The contributions that do not lead to aliasing errors are unaffected by this procedure. In three dimensional systems, height evaluations with different shifts are needed to get an alias-free computation of the nonlinear terms. This is of course quite prohibitive and an approximation is made. The height shift computations of the nonlinear terms are done as part of the sub-step of a Runge-Kutta scheme. The phase-shift dealiasing method represents the reason behind the four step implementation for the third order Williamson Runge-Kutta time advancement method. This dealiasing method represents an approximative process of removing the aliasing errors.
The TURBO code is parallelized through the use of MPI. The numerical cube is split in the ydirection in spectral space, which limits the maximum number of processors used for one run to N y . To improve the use of available parallel computing resources, an additional parallelization direction is made over the number of instances performed during each run, i.e. number of "cubes" solved at the same time that can exchange data among themselves. This last feature of the TURBO code allows the direct computation of ensemble average quantities.

III. MHD IDEAL INVARIANTS
In MHD turbulence, the three quadratic 4 ideal invariants have an important role in the dynamics of turbulence and the resulting cascades that appear for fully developed turbulence regimes [4,9]. For ideal MHD, that is an inviscid flow with zero magnetic diffusivity, total energy, cross-helicity and magnetic-helicity are conserved in absence of forcing. Magnetic-helicity H m = a · b is a purely magnetic invariant, defined as the scalar product between the magnetic potential a and the magnetic field b = ∇ × a, where . . . denotes volume average. Since it possesses an inverse cascade (transfer form small to the large scales [8]), magnetic-helicity plays an important role in dynamo phenomena, see [2,3]. The other two quadratic invariants, total energy E = E u + E b and cross-helicity H c = u · b should be studied together as they affect each-other. Since kinetic E u = u · u /2 and magnetic E b = b · b /2 energies are not conserved individually, the use of Elsässer variables z ± = u ± b might be more appropriate. In fact, in the Elsässer representation, the cross-helicity and total energy information is contained in the definition of two ideal invariants E + = z + · z + /4 and E − = z − · z − /4, known as pseudo-energy. The E ± ideal invariants are positively defined, which represents an advantage in spectral study. Since cross-helicity is related to the degree of alignement between the velocity and the magnetic field, the cross-helicity level defined as, provides global information regarding the alignment present in a system (ρ c ∈ [−1, 1]). Point-wise, the presence of a large cross-helicity level ρ c (x) = H c (x)/E(x) gives rise to the phenomena of nonlinear depletion in the evolution equations. Nonlinear depletion weakens the ability of the nonlinear terms to mix the flow compared to hydrodynamic case, which results in a different turbulent mixing time for MHD turbulence [19]. Moreover, it was shown in [11] that MHD turbulence is composed of zones of highly aligned (ρ c (x) ∼ +1) and highly anti-aligned structures (ρ c (x) ∼ −1) which are thought to develop naturally through the process of dynamical alignment [6]. Therefore, a non-zero value for the global parameter ρ c denotes a preference in the generation of one type of aligned structures. This situation is known as imbalanced MHD turbulence, Figure 1. It is interesting to note that the presence of cross-helicity does not necessarily affect the equipartition of energy between the u and b fields. However, it does change the pseudo-energy levels of z + and z − and therefore the strength of the co-propragating and contra-propagating Alfvèn waves that scatter on each-other. Although the ideal invariants are not conserved in the presence of viscosity and magnetic diffusivity, they are still redistributed between different scales without loss or gain through their nonlinear fluxes. For X, Y , Z standing in for the fields, the flux through a spherical surface in the Fourier space defined by a radius k c is defined as: The average represents here the sum over all the Fourier modes. For a turbulent state, the value of the fluxes in the inertial-inductive range are expected to be constant. In the Elsässer formalism, the values Π + +,− and Π − −,+ reached in the inertial-inductive range, enter in the phenomenological definition of pseudo-energy scaling laws [16].

IV. THE FORCING MECHANISM
As it is known, a forcing mechanism can be used to reach different turbulent regimes. In this study, the injection rates of the nonlinear ideal invariants will be used as the control parameters. These injection rates can easily be computed from equations (1)(2). For instance, the injection rates of total energy and of cross-helicity (H c = {û(k) ·b(−k)} ) are given by where, thanks to the Parseval theorem, the volume average . . . can be identified as the average over the number of modes. The parameters ε u and ε b represent the power injected byf u andf b respectively. Since the cross-helicity is bounded by the total energy, the sum of σ u and σ b , which individually denote the cross-helicity injected byf u andf b respectively, has to fulfil the condition: −ε ≤ σ u + σ b ≤ +ε. As such, the cross-helicity parameter σ is bounded in the interval [−1, 1]. Selecting the force control parameters in such a way to fix ε and σ will enforce the dissipation level for the energy and crosshelicity, once the stationary regime is reached as shown on Figure 2. This behavior is true only for ideal invariant quantities. For example, selecting the injection level of kinetic helicity will not enforce the kinetic helicity dissipation for MHD turbulence as it would in a purely hydrodynamic flow. Numerically, we consider the forcesf u (k) andf b (k) to be local quantities in Fourier space, which act in the same manner on all the modes N f within a wavenumber shell defined by the interval s f = [k inf , k sup ]. Sinceû,b andf u,b are divergence free, we use a helical decomposition [17] for the definition of the force. The helical decomposition projects a vectorâ(k) on a complex basis h ± : where h ± = e 1 × e 2 ± ie 1 and e 1 = (λ × e 2 )/ λ × e 2 , e 2 = k/k. The wave-vector λ is taken to be arbitrary and non-parallel to k. We warn the reader that in this section, the ± lower indices refer to the helical basis h ± , the vector projections on this basis and their contributions to the different scalar quantities of interest, while the upper indices denote, as usual, quantities in the Elsässer representation. The helical decomposition is very useful in ensuring zero divergence and in computing the curl operator. Indeed, the vectors h ± are eigenmodes of the curl operator, ik × h ± = ±kh ± . As such the vorticity is now defined asω ± (k) = ±kû ± (k) and the electric current has the formĵ ± (k) = ±kb ± (k).
We choose the projection of the forces on h ± , to have the form: if |k| ∈ s f and zero otherwise. Since the α's and β's parameters that need to be determined are considered to be real, the forcing method presented here does not influence the phases of the fields, which ensures that no change is made in the type of turbulent structures present. Another way of injecting cross-helicity into the system would be achieved by imposing the alignment of u and b in the real space, which in turn would modify the phases of the fields and potentially the turbulence behavior.
For a mode k, the kinetic energy (E u To simplify the numerical implementation of the force, we make the following assumption: we consider σ u ± = σε u ± and σ b ± = σε b ± . In practice, these equalities impose that both forcingf u and f b are responsible for the same amount of cross helicity injection in the system. Moreover, they also impose that cross helicity is injected at the same rate in both the h + and h − components of the velocity and magnetic fields. These assumptions, which could be reconsidered in the future, allow us to fix the eight real parameters α u ± (k), β u ± (k), α b ± (k) and β b ± (k) by giving only five control parameters, namely the energy injection rates ε u ± , ε b ± and the cross-helicity parameter σ, where the respective injection rates are assumed to be the same for all the N f forced modes. Selecting a large number of forced modes ensures that no anisotropy effect is induced by the forcing mechanism.
Because of the Cauchy-Schwarz inequality, we have the condition H c . We see that for the equality case, we develop a pole for the parameters α u ± (k), β u ± (k),α b ± (k) and β b ± (k). In an effort to advert this, we consider the condition σε ± ≤ 4ε u ± ε b ± on the control parameters, where ε ± = ε u ± + ε b ± . Since the injection rates per mode for the kinetic-helicity and magnetic helicity depend on k, the global kinetic-helicity injection rate h and the global magnetic-helicity injection rate χ are found as, From the above expression of the forces, we find the forces that act in the Elsässer form of the MHD equations asf ± =f u ±f b . If we chose a purely mechanical forcing of the turbulence (f b = 0) for this forcing method, we obtain the relationf + =f − =f u . For this case ε b = 0, which numerically is found to be unstable unless no cross-helicity is injected, σ = 0.
As a note, we see that using the helical decomposition we can redefine the hydrodynamical force we used in previous studies [7,15] (f b ≡ 0). From conditions imposed on the energy and kinetic helicity injection levels, we take the projections on h ± of the force to be: The condition for the kinetic helicity injection rate is automatically fulfilled for any selection of ε ± .

V. CONCLUSIONS AND DISCUSION
A forcing mechanism which controls the injection level of the three ideal invariants of MHD turbulence: the total energy, the cross-helicity and the magnetic helicity has been developed and results obtained from the implementation of this force into the TURBO solver have been presented. This type of force represents a useful tool for the spectral study of turbulence, since it allows to control the level to which all three ideal quadratic MHD invariant fluxes relax to.
As seen from Figure 3, the presence of a non-zero cross-helicity injection level causes a change in the two pseudo-energy dissipation rates, which in turn causes a separation of the respective fluxes levels once a statistical stationary regime is reached. Looking at the spectra of E + and E − , we observe a 5/3 scaling. For imbalanced MHD turbulence, we also observe a difference in the pseudo-energy spectra levels. However, the scaling exponent tends to remain 5/3. By controlling the cross-helicity injection level, this type of force should provide a help in the study of solar wind turbulence, which represents a well known case of imbalanced MHD turbulence. Also, controlling the magnetic helicity injection level should be of great help in the field of galactic and solar dynamo physics.