Nonlinear Effects on the Precessional Instability in Magnetized Turbulence

: By means of direct numerical simulations (DNS), we study the impact of an imposed uniform magnetic ﬁeld on precessing magnetohydrodynamic homogeneous turbulence with a unit magnetic Prandtl number. The base ﬂow which can trigger the precessional instability consists of the superposition of a solid-body rotation around the vertical ( x 3 ) axis (with rate Ω ) and a plane shear (with rate S = 2 ε Ω ) viewed in a frame rotating (with rate Ω p = ε Ω ) about an axis normal to the plane of shear and to the solid-body rotation axis and under an imposed magnetic ﬁeld that aligns with the solid-body rotation axis ( B (cid:107) Ω ) . While rotation rate and Poincaré number are ﬁxed, Ω = 20 and ε = 0.17, the B intensity was varied, B = 0.1, 0.5, and 2.5, so that the Elsasser number is about Λ = 0.1, 2.5 and 62.5, respectively. At the ﬁnal computational dimensionless time, St = 2 ε Ω t = 67, the Rossby number Ro is about 0.1 characterizing rapidly rotating ﬂow. It is shown that the total (kinetic + magnetic) energy ( E ) , production rate ( P ) due the basic ﬂow and dissipation rate ( D ) occur in two main phases associated with different ﬂow topologies: (i) an exponential growth and (ii) nonlinear saturation during which these global quantities remain almost time independent with P ∼ D . The impact of a "strong" imposed magnetic ﬁeld ( B = 2.5 ) on large scale structures at the saturation stage is reﬂected by the formation of structures that look like ﬁlaments and there is no dominance of horizontal motion over the vertical (along the solid-rotation axis) one. The comparison between the spectra of kinetic energy E ( κ ) ( k ⊥ ) , E ( κ ) ( k ⊥ , k (cid:107) = 1, 2 ) and E κ ) ( k ⊥ , k (cid:107) = 0 ) at the saturation stage reveals that, at large horizontal scales, the major contribution to E ( κ ) ( k ⊥ ) does not come only from the mode k (cid:107) = 0 but also from the k (cid:107) = 1 mode which is the most energetic. Only at very large horizontal scales at which E ( κ ) ( k ⊥ ) ∼ E ( κ ) 2 D ( k ⊥ ) , the ﬂow is almost two-dimensional. In the wavenumbers range 10 ≤ k ⊥ ≤ 40, the spectra E ( κ ) ( k ⊥ ) and E ( κ ) ( k ⊥ , k (cid:107) = 0 ) respectively follow the scaling k − 2 ⊥ and k − 3 ⊥ . Unlike the velocity ﬁeld the magnetic ﬁeld remains strongly three-dimensional for all scales since E ( m ) 2 D ( k ⊥ ) (cid:28) E ( m ) ( k ⊥ ) . At the saturation stage, the Alfvén ratio between kinetic and magnetic energies behaves like k − 2 (cid:107) for Bk (cid:107) / ( 2 ε Ω ) < 1.


Introduction
Rotating flows are omnipresent in several geophysical and astrophysical systems (oceans, Earth's atmosphere [1], gaseous planets, stars [2]). The role of the Coriolis force, however, is not as dominant in the dynamics of the Earth's atmosphere as it is in the atmosphere of giant planets, regarding its impact on zonal flows (e.g., see Galperin [3]). This relative weakness of the Coriolis force in Earth's where x = (x 1 , x 2 , x 3 ) T denotes the physical coordinate. We indicate that the vorticity vector of the plane shear balances the gyroscopic torque created by the misalignment of Ω = Ωe 3 and Ω p = εΩe 1 . The additional plane shear, which results from the precession, was found in previous experimental studies (see, e.g., [25]). Mahalov [24] considered an infinite cylinder in which a tilted (sheared) streamline solution can exist under precession. From that study, Salhi and Cambon [23] have derived the solution (1) and called it as Mahalov's solution. A similar solution to (1) was derived by Kerswell [11] from the Poincaré's solution for a precessing oblate spheroid [13]: for the Mahalov's solution the velocity vector of plane shear U s = −2εΩx 2 e 3 aligns with the solid rotation axis while for the Kerswell's solution U s = −2εΩx 3 e 2 is normal to the solid-body rotation axis but for the two solutions the absolute vorticity is the same, ∇ × U + 2Ω p = 2Ωe 3 . The limit case ε = 0 corresponds to the case of rotating turbulence that was widely addressed in the literature to characterize the role of the Coriolis force (see, e.g., [8,26]).
Mahalov's solution corresponds to an extensional flow, characterized by space-uniform vector fields (vorticity field, magnetic field) and/or space-uniform velocity gradients (shear, strain). Disturbances (or fluctuations) to that base flow are expanded in terms of Fourier modes with a (periodic) time-dependent wave vector. This is for the viewpoint of conventional stability analysis: in this context, a single disturbance mode cannot interact with itself, in purely incompressible turbulence (e.g., see [27]), so that the whole base flow plus the disturbance flow can be seen as an exact solution of Navier-Stokes equation in an unbounded domain, provided that some admissibility conditions be fulfilled. Regarding more realistic flow cases, linear spectral theory (LST), also called rapid distortion theory (RDT), for homogeneous turbulence can have the same starting point as the previous stability analysis but initial data can consist of a superposition of many individual Fourier modes, with a dense spectrum. As an advantage with respect to stability analysis, very large Reynolds numbers can be reached only with the particular choice of initial dense spectra. In counterpart, the linear dynamics from given initial data are only valid in a "rapid" time of evolution.
The third aspect investigated in this article is the interplay with MHD, that is relevant for the study of some magnetized planetary atmospheres. In addition to details found in [11,23,24] without MHD, new instability tongues are presented in Salhi et al. [28] for magnetized precessing flow (see also [8]). The aim of the present study is to provide new DNS results for explicit nonlinear dynamics on magnetized precessing rotating flows including stability analysis and related LST.
For the Mahalov's solution (1), the wave-vector component along the solid-body rotation axis, k , remains time-independent while the other two components are periodic in time with period T = 2π/Ω. This makes it possible at any time t to compare DNS results characterizing the precessing turbulence dynamics at the planes k = constant with those in forced rotating turbulence (see later). In counterpart, for the Kerswell's solution, k is periodic in time with period 2π/Ω (see [11,23,29,30]).
In rapidly rotating turbulence (i.e., the Rossby number is small) the emergence of large-scale vortices characterizes a quasi-two-dimensional (2D) state perpendicular to the fixed rotation axis [26,31], and, at increasing rotation rate, energy accumulates in the vicinity of the so-called "spectral buffer layer" around k = 0 [26]. The Rossby number Ro = u l /(2lΩ), where l is a characteristic length scale and u l is the associated velocity, measures the ratio of the advection term on the Coriolis force in the Navier-Stokes equations, also a small value of Ro means a dynamics mainly driven by rotation. Typical large-scale planetary flows are characterized by Ro ∼ 0.1 [32]. At k = 0, the inertial wave frequency vanishes. The separation between the dynamics of velocity modes with non zero and zero inertial wave frequencies, or equivalently, with k = 0 (wave modes) and k = 0 (vortex modes) was discussed and observed in previous studies of rotating flows [33][34][35][36].
In the astrophysical context, Barker [30] used the Kerswell's solution and studied, by means of DNS, the effect of nonlinear interactions on the precessional instability for several values of the Poincaré number (0.01 ≤ ε ≤ 0.3) while rotation rate was fixed, Ω = 1. He found that the instability led to turbulence and the formation of columnar vortices similar to those observed in rapidly rotating turbulence. The similarity between the observed elongated structures in the two flow cases (precessing and rotating) was investigated in more detail in the study by Khlifi et al. [10]. In that study, the Mahalov's solution (Equation (1)) was considered and DNS was performed for several values of the couple (Ω, ε). It is shown that the time development of the kinetic energy occurs in two main phases associated with different flow topologies: (i) an exponential growth characterizing three-dimensional turbulence dynamics and (ii) nonlinear saturation during which the kinetic energy remains almost time independent, the flow becoming quasi-2D-dimensional. The latter stage, wherein the development of kinetic energy remains insensitive to the initial state, shares an important common feature with other quasi-two-dimensional rotating flows such as rotating Rayleigh-Bénard convection [37,38], or the large atmospheric scales [9,[39][40][41]: in the plane k = 0, the kinetic energy spectrum scales as k −3 ⊥ for k S < k < k Ω where k S and k Ω are the Zeman "precessional" and "rotational" wavenumbers, respectively (defined by Equation (29)). In addition, an inverse cascade develops for (k ⊥ , k) < k S , and the spherically averaged kinetic energy spectrum exhibits a k −2 inertial scaling for k S < k < k Ω .
Some DNS results for explicit nonlinear dynamics on magnetized precessing rotating flows [30] indicate that, in the presence of even a weak magnetic field, the properties of the precessing flow are significantly altered and the formation of the columnar structures (observed in precessing hydrodynamic turbulence, hereinafter, PHT) is inhibited. In the present DNS study, we analyze the impact of an imposed magnetic field on turbulence dynamics by investigating various spectra for kinetic and magnetic energies. For the Mahalov's solution the imposed magnetic field must align with the solid-body rotation axis to satisfy the admissibility conditions [27]. Denoting by B the associated Alfvén velocity vector. We perform DNS for which rotation rate and Poincaré's number are fixed, Ω = 20 and ε = 0.17 while the magnetic field intensity is varied, B = 0.1, 0.5 and B = 2.5, so that Λ = B 2 /(2Ωη) = 0.1, 2.5 and Λ = 62.5 respectively. The value B = 2.5 can be seen as a "strong" intensity. Initially (at t = 0) the Rossby number is about Ro ∼ 0.01 and reaches Ro ∼ 0.1 at the final computational time St = 2εΩt = 67. In addition, we perform DNS for the case of PHT with Ω = 20 and ε = 0.17 and for the case of rotating MHD turbulence (hereinafter, RMHDT) with Ω = 20 and B = 2.5.
Note that RMHDT with B Ω was investigated in recent studies [42][43][44][45][46]. For instance, it is found that, at large scales, the equipartition between magnetic and kinetic energies observed in the case of non-rotating MHD turbulence in a weak turbulence regime (u l B) is broken by the dynamics of inertial waves, so that, the Alfvén ratio between kinetic and magnetic energies behaves like k −2 for at L Ω = kB/(2Ω) < 1 [42,45]. A similar behavior is found for the staturation stage of precessing MHD turbulence (hereinafter, PMHDT), as shown in the present study.
The paper is organized as follows. Model equations are given in Section 2 in physical space, for base (mean) flow and disturbance (fluctuating) one, with additional relationship for statistics and energetics. Section 3 is devoted to their spectral counterpart, mainly in a continuous unbounded domain, with time-dependent wave-vector. Some results from the linear stability analysis of precessing MHD flows by Salhi et al. [28] are reported in this section. Numerical results from DNS are given and discussed in Section 4. Concluding remarks are offered in Section 5.

Incompressible MHD Equations
We consider a precessing homogeneous turbulence of an incompressible electrically conducting fluid in the presence of a uniform background magnetic field. In the following, mean and fluctuating magnetic fields are scaled using Alfvén velocity units i.e., they are divided by √ ρ 0 µ 0 where ρ 0 and µ 0 are the constant density and the magnetic permeability. Incompressible MHD equations provide the simplest theoretical framework for studying large-scale turbulent motion of an electrically conducting fluid. Accordingly, the equations for the instantaneous velocityũ and magnetic fieldb are [47] wherep is the pressure (including the magnetic part and the centrifugal potential), ν is the kinematic viscosity and η is the magnetic diffusivity. We choose ν = η, so that the magnetic Prandtl number is Pm = η/ν = 1. In the momentum conservation Equation (2a), (ũ·∇)ũ is the inertial term and −2Ω p ×ũ and b ·∇ b are the Coriolis and Lorentz forces per unit mass, respectively. In the induction Equation (2b), (ũ·∇)b and b ·∇ ũ respectively represent the advection and stretching of magnetic field. Instantaneous velocity and magnetic fields are decomposed into base (or mean) flow and disturbances (or fluctuations),ũ = U + u(x, t) andb = B + b(x, t). As indicated previously, the base flow considered in the present study (see Equation (1) and Figure 1) corresponds to a vertical x 3 solid-body rotation viewed in a rotating frame about the x 1 axis, (Ω p = εΩe 1 ), with additional plane shear, U = A·x, in an external uniform magnetic field B = (0, 0, B) T , where the matrix A is given by (1). Obviously, the base flow (1) is an admissible solution, i.e., a solution of Equation (2) [23,28]. This ensures consistency with statistical homogeneity [8].

Equations for the Kinetic and Magnetic Energies
The resulting MHD equations for the fluctuating velovity u(x, t) and magnetic field b(x, t) can be written as follows: in which p(x, t) denotes the fluctuating modified pressure. We assume that the turbulence is statistically homogeneous, so that the resulting equations for the kinetic and magnetic energies per unit mass E (κ) = 1 2 u i u i and E (m) = 1 2 b i b i , where · denotes ensemble (statistical) averaging, are in which P (κ) = 2εΩ u 2 u 3 denotes the turbulent kinetic energy production which is proportional to the Reynolds stress component u 2 u 3 , while P (m) = −2εΩ b 2 b 3 denotes the turbulent magnetic energy production which is proportional to the Maxwell stress component b 2 b 3 . Both P (κ) and P (m) are due to the background shear with rate S = 2εΩ. Note that rotation produces no energy.
Here, D (κ) = ν ∂ x j u i ∂ x j u i and D (m) = η ∂ x j b i ∂ x j b i denote respectively the viscous and magnetic (or Joule) dissipation rates per unit mass and represents the energy exchanges between magnetic and kinetic energies due to the imposed magnetic field. Note that the term F (κm) has no explicit contribution on the total (kinetic + magnetic, E = E (κ) + E (m) ) energy evolution, where P = P (κ) + P (m) and D = D (κ) + D (m) .

Time-Dependent Wave Vector
Fourier space representation is useful to quantify the energy contents at various scales. Therefore, we represent the fluctuating fields as the superposition of Fourier modes using the Fourier transforms where k(t) is a time-dependent wave vector. In order to remove the explicit dependence on x in the resulting equations for the Fourier amplitudesû(k, t) andb(k, t) one has to ensure that k(t) varies in time according tok = −A T ·k (i.e., the eikonal equation), This amounts to following characteristic lines of the mean flows, although expressed in spectral variables (details in Ref. [23,28]). Accordingly, the wave-vector components evolve as where k j0 with j = 1, 2, 3 is the initial wave vector component. Therefore, the wave-vector trajectories are circles with sheared centers, as their counterparts in physical space, since k·x = k 0 ·x 0 . We note that both the wavenumbers are time-independent.

Equations for the Fourier Amplitudes
The resulting equations for the Fourier coefficientsû(k, t) andb(k, t) arėû wherep ( ) is the spectral counterpart of the linear part of pressure fluctuations. Here,u is the time derivative. The two terms t (κ) and t (m) are bilinear operators for the kinetic and magnetic energy transfers, written as [48] t (κ) where P αβ = δ αβ − k −2 k α k β is the projection operator and δ αβ is the Kronecker delta.
In the direct numerical simulation approach used here, the fully nonlinear equations of motion (11) are solved using a modified version of the Fourier spectral code Snoopy by [49]. Note that in this entire paper, as for the basic nonlinearity seen as a convolution product in 3D Fourier space, according to previous equations, both the physical space of coordinates and the spectral space are considered as continuous. Consequently, integrals are used, and this description will prevail for the definition of statistical quantities, such as averaged energy spectra, single-point statistics, etc., obtained by integration of 3D basic spectra. On the other hand, the calculation of terms connected to nonlinearity in DNS follows the classical pseudo-spectral procedure inherited from Orszag and Patterson [50]: a regular mesh is used in a triple-periodic box, with related discretization and the convolution product is evaluated from its physical counterpart by inverse fast Fourier transform (FFT). More details can be found at the beginning of Section 4. In addition, we use here a less classical, modified, version of the Snoopy code [49] for DNS in the presence of mean distortion: as for SLT, the equations are solved using the comoving system of coordinates; this is equivalent to following mean (elliptical here) trajectories as characteristic lines in physical space and to using time-dependent wave vectors in Fourier space. Finally, the classical pseudo-spectral procedure for nonlinear terms is applied to an initial grid or an initial vector k 0 , as in Equation (9).

Linear Stability Analysis
In the inviscid (ν = 0) and ideal (η = 0) linear limit, the differential system (11) reduces to a seventh order Floquet system˙û The use of the constraints k iûi = 0 and k ibi = 0 allows us to eliminate the Fourier amplitude of fluctuating pressurep (14) and to reduce the latter seventh order differential system to a fourth order one. The rank of the differential system can be reduced using new variables, corresponding to solenoidal modes, that are individually divergence-free in physical space, or normal to the wavevector in Fourier space. The three Fourier components of the velocity field reduce to only two components, denotedĉ 1 andĉ 2 , related to a poloidal/toroidal divergence-free decomposition (see [8] for details). The same holds for the magnetic field, withĉ 3 andĉ 4 components. In terms of these variables [28] The resulting fourth order Floquet system takes the form where τ = Ωt. The elements D jm are not reported here for the sake of brevity (see [28]). We denote by Φ(τ) the fundamental matrix solution of the Floquet system (16), dΦ/dτ = D·Φ with Φ(0) = I 4 where I 4 is the unit matrix and by M = Φ(2π) the Floquet multiplier matrix [51,52]. Because it is found that Det Φ = (k 0 /k) 2 and hence Det M = 1 [28,53]. This implies that the product of the eigenvalues of M is unity, say [51]. A necessary condition for the onset of linear instability is a resonance where two Floquet multipliers coincide (if an eigenvalue is at the onset of instability, it must have multiplicity 2 or higher). In the stable case, all the eigenvalues lie on the unit circle [52].
The resonant cases for ε = 0 are characterized by the condition ω i − ω j = , where is an integer and ω j (j = 1, 2, 3, 4) denotes the frequency (normalized by Ω) associated with the eigenvalue λ j of the matrix solution M(2π) Here, L Ω is a local Lehnert number [54], The first two frequencies (ω 1,2 ) correspond to hydrodynamic modes or fast modes, while the second two (ω 3,4 ) refer to magnetic modes [52], slow modes [55], or magnetostrophic modes [44] since they vanish in the purely hydrodynamic case (L Ω = 0). Previous studies for magneto-inertial wave turbulence indicate that the equipartition between kinetic and magnetic energy due to Alfvén waves is broken by inertial waves for L Ω < 1 [43,45].
The stability analysis at first-order in ε by Salhi et al. [28] shows that there are three subharmonic instabilities associated to ω 1 − ω 2 = 1 (hydrodynamic modes), ω 1 − ω 3 = 1 (mixed modes), and ω 3 − ω 4 = 1 (magnetic modes). The analytical results for the maximal growth rate σ m of the three subharmonic instabilities normalized by ε versus the local Lehner number L Ω are displayed in the right panel of Figure 2. Note that the magnetic modes instability occurs only for L Ω > √ 5/4 ≈ 0.559. As it can be seen, for the magnetic modes instability σ m /ε increases with L Ω approaching 1 2 at L Ω 1, while for the hydrodynamic modes instability σ m /ε, which takes the value 5 √ 15/32 at L Ω = 0, decreases with L Ω approaching 1 2 at L Ω 1. In counterpart, for the mixed modes instability, σ m /ε approaches zero at L Ω 1. Nonlinear interactions act to saturate the linear precessional instability since the time evolution of the turbulent (magnetic) energy exhibits two important stages: an exponential growth stage during which the energy grows exponentially with time followed by a saturation stage where the energy remains almost constant (see Section 4, see also [10]). The present simulations indicate that, at L Ω < 1, the quasi-equipartition between kinetic and magnetic energies observed in non-rotating MHD turbulence is broken while, at L Ω > 1, the magnetic energy slightly exceeds the kinetic energy (see later). Variation of the maximal growth rate σ m /ε versus the local Lehner number L Ω = ω A /ω R = Bk/(2Ω) for the precessing magnetohydrodynamics (MHD) instability, where ω A = Bk and ω R = 2Ωk /k are the frequencies of Alfvén and inertial waves, respectively. Hydrodynamic modes: Salhi et al. [28]).

Spectral Density of Energy
The resulting equation for the spectral density of kinetic energy e (κ) = 1 2 u i u * i and that for the spectral magnetic energy density e (m) = 1 where denote the spectral density of kinetic and magnetic energy production, respectively, and again characterizes the energy exchange between the kinetic and magnetic components due to the background magnetic field B. Here, T (κκ) , T (κm) , T (mκ) , and T (mm) are transfer-like terms that involve triple velocity and magnetic fluctuation correlations and are such that The energy exchange between kinetic and magnetic component due to nonlinear interactions is characterized by both terms T (κm) and T (mκ) . At k = 0, the inviscid linear solution shows that there is no linear instability: the vertical modê u 3 remains identical to its initial value, as in the pure rotating case (ε = 0), while the horizontal mode (û 2 = −(k 1 /k 2 )û 1 ) performs a precessing motion (with period Ω) induced by the vertical mode, where k ⊥ = k 2 1 + k 2 2 . Obviously, at k = 0, the nonlinear interactions can generate a coupling between vertical and horizontal modes as shown by the equation for the horizontal and vertical spectral density of energies, (the same notation is used for the fluctuating magnetic fields), in which the subscript 2D refers to the k = 0 mode, for example, e

Initial Conditions
We initialize direct numerical simulations from homogeneous isotropic hydrodynamic turbulence generated by a separate simulation started from random initial conditions with the initial energy spectrum [56] dk is the initial kinetic energy. In order to let the turbulence develop a significant inertial zone as well as nonlinear transfer, an isotropic precomputation [10,57] is done before applying the mean flow (1). Only during this pre computation, a large-scale forcing is applied until a statistical steady state, with a Taylor-microscale-based Reynolds number Re λ ∼ 40. This state holds for a new initial one before starting the effect of the mean flow; the procedure is not sufficient to reach a developed Kolmogorov-like energy spectrum but the trend is satisfactory. In all the runs, the magnetic fluctuations are initially zero. Note that in the case of PHT, it is found that the saturation stage, during which the turbulent kinetic energy remains almost time independent, is insensitive to the initial Reynolds number independently of the initial conditions (precomputed initial isotropic or solenoidal random noise) [10].
The domain is a periodic cube with edge length L 0 = 2π with n 0 = 256 grid points to a side. Aliasing errors are removed using the 2/3 dealiasing rule and as a result the minimum and maximum wavenumbers are k min = 1 and k max = n 0 /3 respectively. We also limit the effects of L 0 −periodic boundary conditions, by computing all integral scales in the x ( = 1, 2, 3) directions, and imposing that max L ( ) jj < 0.5L 0 (no sum over subscript j ). By this, we control numerical confinement, even if one integral length scale grows faster due to precessional instability.
In our simulations, the initial value of Rossby number is defined as follows where Ro = 0.01 characterizes rather a case of rapidly rotating turbulence, which is known to mimic some of the features of the 2D turbulence. Here, l = u l E/(πD) is a characteristic length scale, u l = 2E (κ) /3 is the associated velocity, and E and D respectively denote the total (kinetic+magnetic) energy and dissipation rate. Typical large-scale planetary flows are characterized by Ro ∼ 0.1 [32], as already indicated, and the liquid metals in the Earth's outer core are strongly affected by rotation with Ro ∼ 10 −6 [58]. For a giant planet like Jupiter, it is believed that the Rossby number may even be smaller [59], whereas for the solar convective region (where the magnetic field is believed to be enhanced) Ro ∼ 1 [44]. The present simulations show that, at the saturation stage (see below), the Rossby number can reach the value Ro = 0.1 (see the inset of Figure 3). The value of the Poincaré parameter considered in the present study is about ε = 0.17, so that, the value of the shear rate S = 2εΩ is S = 6.8. In addition to the case of PHT, we consider three values for the strength of the imposed magnetic field in the case of MHD precessing turbulence, namely, B = 0.1, 0.5, 2.5. The respective values of the Elsasser number, which characterizes the relative strengths of the Lorentz and Coriolis forces, are Λ = 0.1, 2.5 and 62.5 (see Table 1). We also consider the case of RMHDT (ε = 0) with B = 2.5, so that Λ = 62.5, and B is parallel to the solid-body rotation axis. Note that the Elsasser number is commonly employed as a measure of the strong field dynamo, which is synonymous with the magnetostrophic balance (i.e., the pressure gradient, the Coriolis force, and the Lorentz force are balanced). It is believed that the Elsasser number of the Earth's core is about unity [59]. All the runs are run up to a maximum computational dimensionless time t + ≡ St = 67.

Time Evolution of Global Quantities
In this section, we study the temporal behavior of several global quantities to characterize the precessing MHD flow dynamics and the influence of the strength of the imposed magnetic field on it.
The present simulations indicate that the time evolution of total (kinetic+magnetic) energy, E = E (κ) + E (m) , exhibits two important stages: (i) an exponential growth of E(t) followed by (ii) a saturation stage at St > 20, during which E(t) remains almost constant. This is also clear when observing the evolution of the energy growth rate, which slightly fluctuates about zero during the saturation stage for all runs, as shown by Figure 3a.
Recall that P (κ) = S u 2 u 3 (respectively, P (m) = −S b 2 b 3 ) is the kinetic (magnetic) production rate and D (k) (respectively, D (m) ) denotes the kinetic (magnetic) dissipation rate. For the runs 2-4 (PMHDT cases), a similar behavior is found for the total dissipation rate, D = D (κ) + D (m) , and the total production rate P = P (κ) + P (m) where they become slowly balanced such that P ∼ D at the saturation stage (see Figure 3b). Indicate that the Rossby number Ro is about 0.05 for the runs 3 and 4 and about 0.1 for the run 2 (see the inset of Figure 3a). For the run 1 (PHT case), both D (κ) and P (κ) exhibit an exponential growth at St < 20 and decay with time at the saturation stage during which they become slowly balanced such that P (κ) ∼ D (κ) (see also Khlifi et al. [10]).
At the exponential growth stage, the kinetic quantities E (κ) , D (κ) , and P (κ) , are approximately insensitive to the B intensity as illustrated by Figure 4a and its inset displaying the time evolution of E (κ) (t)/E (κ) (0) and P (κ) . In counterpart, the generation of large and small scales magnetic fluctuations during the exponential growth stage is sensitive to the B intensity as shown by Figure 4b and its inset displaying E (m) (t)/E (κ) (0) and D (m) (t)/D (κ) (0), respectively. Recall that our simulations start from homogeneous isotropic hydrodynamic turbulence, so that, the initial magnetic energy and dissipation rate (or Joule dissipation) are zero. From Figure 4b, it clearly appears that for both the runs 2 and 3 (for which B = 0.1 and B = 0.5, respectively) the magnetic energy increases with time even at the saturation stage while for the run 4 (for which B = 2.5), it remains almost constant at the saturation stage. As for the behavior of the magnetic (or Joule) dissipation rate at the saturation stage, we observe that D (m) remains almost constant for the three runs (2-4, see the inset of Figure 4b). In the following, we analyze the dynamics of the saturation stage and comment about the exponential growth stage when relevant. For non-rotating MHD flows in the weak turbulence regime (u l B), the Alfvén ratio between kinetic and magnetic energies, Γ = E (κ) /E (m) , is equal to one, thus characterizing an equipartition between kinetic and magnetic energies. However, For non-rotating MHD flows in the strong turbulence regime (u l B), this ratio is slightly smaller than one, indicating the presence of non-Alfvénic fluctuations [60]. For RMHDT (ε = 0), inertial waves have an important impact on the dynamics of large scales since the equipartition between kinetic and magnetic energies is broken. At large scales such that L Ω = Bk/(2Ω) 1, the generation of magnetic fluctuations is reduced more and more by the inertial waves as the rotation rate increases (whereas the B intensity remains constant). This leads to an increase of the Alfvén ratio with rotation rate [42,43]. In addition, it was shown that, scale by scale, the Alfvén ratio behaves like k −2 for L Ω 1 and is about one for L Ω > 1 ( [45], see also the next section).
In Figure 5a we plot the time evolution of the Alfvén ratio for the three runs 2-4 (i.e., PMHDT cases). We observe that, at the saturation stage, the Alfvén ratio decreases as the B intensity increases (whereas rotation rate remains constant). For both the runs 2 (for which B = 0.1) and 3 (for which B = 0.5), Γ(t) decreases with time reaching the value Γ ∼ 10 at the final computational time. In counterpart, for the run 4 (for which B = 2.5), the Alfvén ratio Γ and also the ratio D (κ) /D (m) are almost constant (at the saturation stage), Γ ∼ 2.5 and D (κ) /D (m) ∼ 0.8 (see Figure 5a and its left inset), while the ratio P (κ) /P (m) oscillates around 5.5 (see the right inset of Figure 5a). Indicate that for the same value of the Elsasser number (Λ = 62.5), the contribution to total (kinetic + magnetic) energy coming from the magnetic component is more important in the case of PMHDT than in the case of RMHDT. Indeed, from Figure 5b displaying the behavior of Γ versus the dimensionless time Ωt/π for the runs 4 and 5, we can see that Γ in the later case (run 5) is approximately six times larger than in the former case (run 4). In order to see if the dominance of the horizontal motion (over the vertical one) observed in the case of PHT [45] or in the case of rapidly rotating turbulence [26,31] is altered or not in PMHDT we examine the behavior of the ratio of horizontal to vertical kinetic energies and the ratio of associated integral length scales. Here, the horizontal and vertical kinetic energies are defined as E We consider only the integral length scales in the direction parallel to the solid-body rotation axis, defined by Equation (24), or equivalently, 22 increases due to rotation (see e.g., [45,61]). Figure 6 shows the time development of v for the runs 1-4. Note that in isotropic turbulence, one has E With respect to the PHT case (run 1), the ratios E v are strongly reduced at the saturation stage only for the PMHDT case with B = 2.5 (run 4). Therefore, we may conclude that the dominance of the horizontal motion (over the vertical one) observed in the PHT case (run 1) is prevented in the presence of a "strong" imposed magnetic field. The visualizations of the vertical vorticity field for both the runs 1 and 4 corroborates such a conclusion (see Figure 7). Indeed, the left panel of Figure 7 related to the PHT case (run 1) indicates the presence of columnar structures along the solid-body rotation which are indicative of a quasi-2D state of turbulence (e.g., see [26]). The recent numerical study by Goto et al. [62] of an incompressible fluid confined in precessing containers (sphere, spheroid, and cylinder) indicates that a pair of twisted large-scale vortices appears in the developed turbulence (ε ≈ 0.1) and shear flow around these vortices stretches and amplifies the small-scale turbulent eddies (the mechanism of small-scale turbulent vortices is consistent with the experimental results for the precessing sphere [18]).  From the right panel of Figure 7 related to the run 4, we observe structures which look like filaments. Note that the inhibition of the formation of columnar structures in PMHDT was also observed in the DNS by Barker [30] as already indicated. On the other hand, in their review of interstellar turbulence (observations and processes), Elmegreen and Scalo [63] indicated that magnetism can slow collapsing cores and remove angular momentum on large and small scales, and it may contribute to filamentary structures that control the accretion rate of gas onto protostars. According to the DNS study by Bigot et al. [60], the formation of filaments is reported within current and vorticity sheets in strongly magnetized flows.

Spectra for Kinetic and Magnetic Energies
In this section, we analyze some spectra of kinetic and magnetic energies since they provide information about the distribution of energy at different scales. A first simple dimensional analysis suggests to delineate the range of anisotropic scales using a single threshold wavenumbers based on shear rate S = 2εΩ or on rotation rate Ω and on kinetic dissipation rate D (κ) , as proposed by Corrsin [64] and Zeman [65] for shear flows and rotating flows, respectively. The corresponding threshold length scales would then be l S = 2π/k S and l Ω = 2π/k Ω . For rotating flows, the wavenumber k Ω identifies a scale at which turbulence eddy turnover time is commensurate with the time scale of system rotation. It characterizes scalewise strength of the effect of rotation: it is strong on scales where k < k Ω and weak on scales where k > k Ω (isotropy is restored for wavenumbers larger than k Ω [66,67]). The present simulations indicate that, at the saturation stage, the precessional k S and rotational k Ω numbers are about 4 and 40, respectively (see also [10]).
In the light of the analysis presented in Section 4.2, we principally consider the runs 1 (PHT) and 4 (PMHDT with B = 2.5).

Radial Spectra
The radial spectrum of kinetic energy is defined as follows where (k, θ, ϕ) is a spherical coordinates system for the wave vector k such that k = k cos θ, k 1 = k ⊥ cos ϕ and k ⊥ = k 2 1 + k 2 2 . Similar definitions hold for other radial spectra. Theoretical arguments for forced rapidly rotating flows suggest that for the intermediate wave numbers where k f denotes the forcing wavenumber, E (κ) (k, t) behaves like E (κ) = C Ω ΩD (κ) 1/2 k −2 , where C Ω = 1.22 − 1.87 [65]. This scaling is also observed in a DNS study of a rotating turbulence under precessionlike perturbation generated by a change of the rotation axis at a fixed time instant [68]. Figure 8a shows the radial (spherically averaged) spectrum of total energy at different times, St = 30, 40, 50 and St = 60 as well as the average E(k, t) t over the time window St ∈ [25,65] for the run 4. The time development of the kinetic energy spectrum and its average E (κ) (k, t) t over the same time window is reported in the inset. As it can be seen, there are no noticeable differences between the spectra and their time average indicating a statistical steady state of turbulence at the saturation stage. Hereinafter, we use the notation E(k) instead of E(k, t) t for the sake of clarity. The same notation holds for the other spectra. Figure 8b displays E (κ) (k) for both the runs 1 (PHT) and 4 (PMHDT case with B = 2.5 or Λ = 62.5). On the basis of the case of PHT, the kinetic energy is strongly reduced at k < k S = 4 even if it also reduced at the other scales, and the wavenumbers range at which E (κ) (k) is compatible with the k −2 scaling is also reduced. For the PHT, E (κ) (k) well follows the k −2 scaling at k S = 4 ≤ k ≤ k Ω = 40. In Figure 8c we plot the radial spectrum of magnetic energy, E (m) (k) for the runs 2-4. It is clear that, at almost all scales, E (m) (k) increases as the intensity of the applied magnetic field is increased. Recall that our simulations start from isotropic hydrodynamic turbulence with zero magnetic fluctuations. Moreover, at 1 ≤ k ≤ 4, E (m) (k) increases like ∼ k α where α = 1.15 for run 3 and α = 1.72 for run 4. Otherwise, for the three runs (2)(3)(4), the amplitude of the magnetic energy is much smaller than that of the kinetic energy in the large scales but of the same order in the small scales as shown by Figure 8d. This is essential for the presence of the 2D-inverse cascade for the kinetic energy since if the magnetic field fluctuations were strong enough in the large scales the flow would behave as a 2D-MHD flow with a direct cascade [69]. The spectral energy flux which is reduced when compared to the PHT case (see the inset of Figure 8b) is positive for k < 5 indicating an inverse energy cascade. Negative energy flux corresponds to a forward cascade (e.g., see [68]). In rapidly rotating (or precessing) hydrodynamic turbulence, energy accumulates in the vicinity of the so-called "spectral buffer layer" around k = 0 [10,26]. This physical phenomenon corresponds to the emergence of large-scale vortices and characterizes a quasi-2D state perpendicular to the fixed rotation axis [31], as already indicated. In Figure 9a we plot the spectrum of the 2D kinetic energy at the plane k = 0, for the runs 1 and 4. It can be seen that E (κ) 2D k ⊥ , k = 0 for run 4 (PMHDT with B = 2.5) is reduced, especially at large scales, when compared to run 1 (PHT). It is rather the horizontal motion that is strongly affected in the presence of a "strong" imposed magnetic field (see Figure 9b) since the vertical kinetic energy E (κ) v2D (k ⊥ ) is almost the same as its counterpart in PHT (see the inset of Figure 9b). This may be explained as follows. The equation for E where P dϕ is the 2D horizontal production term (the 2D vertical one is zero). At k = 0, there is no linear precessional instability (see § III.4) and the coupling between horizontal and vertical kinetic energies are only due to the nonlinear transfers which are strongly reduced by rotation (e.g., see [8]). In counterpart, the 2D horizontal kinetic energy E  10 100 Even if the dynamics of the spectral buffer layer around k = 0 [26] may not be completely captured by simulations due to discretization between the planes k = 0 and k = 1, similarities between the dynamics of this spectral buffer layer and the dynamics of its near-neighbor 3D modes can be drawn. Figure 10a,b compares the following kinetic energy spectra E (κ) and for the runs 1 et 4. Note that E (κ) 3D (k ⊥ , k = ) characterizes the dynamics of the 3D modes. For run 1 (PHT), it is the k = 0 that is the most energetic one, at least for k ⊥ ≤ 4, since . We note that, at k S < k ⊥ < k Ω , the behavior of E (κ) (k ⊥ ) well follows the k −2 ⊥ scaling whereas the behavior of E scaling and the three modes (k = 0, 1, 2) have substantially the same contribution to kinetic energy (see Figure 11a). Comparing the saturation stage of PHT with other quasi-2D flows, the spectrum k −3 ⊥ was observed in "isotropically" forced rotating turbulence [36,[70][71][72][73] and in rapidly rotating Rayleigh-Bénard convection [37,38]. 10 100 3D (k ⊥ , k = 1, 2) for runs 1 and 4, respectively. The dashed straight lines, which indicate the k −2 ⊥ and k −3 ⊥ scaling, delimit the range k S = 4 ≤ k ⊥ ≤ k Ω = 40. 10 100 0.00 0.08 0.16 0.24 1 10 100 2D (k ⊥ ), P (m) (k ⊥ , k = 1, 2), and 1 10 P (κ) (k ⊥ ).
For the PMHDT case with a "strong" imposed magnetic field (run 4), the contribution to E (κ) (k ⊥ ) coming from E (κ) 2D (k ⊥ ) is dominant only in the very thin layer around k ⊥ = 1 while the contribution coming from E (κ) Figure 10b). For instance, 3D (k ⊥ , k = 1) for k ⊥ = 4 at which E (κ) (k ⊥ ) reaches its maximum. We note that the spectra E (κ) (k ⊥ ) and E (κ) 2D (k ⊥ ) are respectively compatible with the scaling k −2 ⊥ and k −3 ⊥ only for the range 10 < k ⊥ < 40. Therefore, we may conclude that for run 4 (PMHDT with B = 2.5) the major contribution to kinetic energy does not come only from the k = 0 mode but also from the k = 1 mode which is the most energetic one. Thus the velocity field in run 4 is quasi-2D only for a very thin layer around k = 0, otherwise, it is three-dimensional. Figure 11a compares the magnetic energy spectra E (m) (k ⊥ ), E 3D (k ⊥ , k = 1, 2) for all scales and the contribution to E (m) (k ⊥ ) coming from the mode k = 1 is more important than those of the modes k = 0 and k = 2. This is mainly due to the fact that P (m) 3D (k ⊥ , k = 1) (i.e., the spectrum of the magnetic production at the k = 1 plane, see Equation (20)) is significantly larger than P (m) 2D (k ⊥ ) and does not strongly differ from P (m) (k ⊥ ) (see Figure 12b).
We indicate that at large scales, the main contribution to the spectrum F (κm) (k ⊥ ), which characterizes the energy exchange between kinetic and magnetic energies (see Equation (21), comes from the mode k = 1 (not shown). The maximum value of F (κm) (k ⊥ ) is approximately the same as the one of P (m) (k ⊥ ) which is 10 times lower than the maximal value of the kinetic production P (κ) (k ⊥ ). These maximum values are reached at k ⊥ ≈ 4.

Alfvén Ratio between Kinetic and Magnetic Energies
As indicated previously, for the case of rotating MHD weak wave turbulence (u l B), linear theory shows that, scale by scale, the Alfvén ratio between kinetic and magnetic energies behaves like so that, Γ(k) ∼ k −2 for large scales L Ω 1 and Γ(k) ∼ 1 indicating equipartion between kinetic and magnetic energies for small scales L Ω 1 (see Salhi et al. [45]). The present DNS results for the case of freely decaying rotating MHD turbulence (run 5) are in agreement with LST predictions as shown by Figure 12a displaying Γ(k) versus L Ω for Ωt/π = 3, 19, and 32. A similar behavior is found for the saturation stage of precessing MHD turbulence (PMHDT) as it clearly appears from Figure 12b displaying Γ(k) versus L S = αBk/S, where the coefficient α depends on the B intensity: α = 0.50 for B = 0.5 (run 3) and α = 0.25 for B = 2.5 (run 4). Indeed, Γ(k) behaves like k −2 for L S < 1 and Γ(k) is sligtly less than one for L S > 1.
In Figure 12c we plot Γ(k ) versus L S for the runs 3 and 4. We observe that Γ(k ) is about one for L S ∼ 1 for both the flow cases. At L S < 1, Γ(k ) is greater than unity and behaves like k −2 while, at k > 1, it is slightly smaller than one. The variation of the Alfvén ratio Γ(k ⊥ , k = 0) at the plane k = 0 versus L S⊥ = αBk ⊥ /(2εΩ) is shown in Figure 12d for runs 3 and 4. It can be seen that, at k ⊥ < 1, Γ(k ⊥ , k = 0) decreases with slope slightly steeper than k −2 . This is mainly due to the fact that the smallest contribution of the modes k = constant to the magnetic energy corresponds to that of the mode k = 0 (see Figure 11a). Recall that the energy exchanges between kinetic and magnetic energies due to the imposed magnetic field vanish at k = 0, and the contribution of that mode to the magnetic production due to the background shear is small (see Figure 11b).

Concluding Remarks
By means of direct numerical simulations (DNS), we studied the impact of an imposed uniform magnetic field on precessing magnetohydrodynamic turbulence with a unit magnetic Prandtl number. Our simulations start from homogeneous isotropic hydrodynamic turbulence, so that the initial magnetic energy and dissipation rate (or Joule dissipation) are zero. Since our main goal is to study the impact of an imposed magnetic field on the dynamics of large scales and intermediate scales, we have limited ourselves to a resolution of 256 3 grid points. The mean flow, which can trigger precessional instability, corresponds to the so-called Mahalov solution (see Equation (1)). Disturbances to that base flow are expanded in terms of Fourier modes with a time-dependent wave vector. For instance, the component (k ) along the rotation axis is time-independent while the other two components of the wave vector are periodic in time. In the background of a linear stability analysis, the mode k = 0 at which both the dispersion frequency (ω A = k B) of Alfvén waves and the one (ω R = 2k Ω/k) of inertial waves vanish, is stable. For k = 0 there is a parametric (precessional) instability resulting from the resonant growth of two magneto-inertial waves when the difference in their frequencies coincides with one of the "distortion" frequencies of the basic flow [28]. At large local Lehnert's number L Ω = ω A /ω R = Bk/(2Ω) 1, the maximal growth rate of instability due to the resonance of slow waves, i.e., those characterized by the frequency ω s = 1 2 −ω R + ω 2 R + 4ω 2 A , approaches the one of instability due the resonance of fast waves, i.e., those characterized by the frequency The instability acts to enhance turbulence until nonlinear interactions (of that turbulence) saturate it. In other words, the energy production due to the background shear, P = 2εΩ u 2 u 3 − b 2 b 3 , feeds the nonlinear interactions and in turn they saturate the instability even in the presence of a "strong" imposed magnetic field, so at the saturation stage, the dissipation (D) and production rates become slowly balanced (P ∼ D). Therefore, at the saturation stage, the turbulence behaves as being forced at large scales by the energy production (rather by kinetic energy production since it is 10 times higher than the magnetic one, see Figure 11b). At the saturation stage, total (kinetic + magnetic) energy, dissipation, and production rates are almost time independent (see Figure 3).
For precessing hydrodynamic turbulence (PHT), similarities between large scale structures of the saturation stage and those occurring in rapidly rotating turbulence were drawn in previous studies [10,30]. The impact of a "strong" imposed magnetic field on large scale structures of the saturation stage is traduced by the formation of structures that look like filaments (see Figure 7). In strongly magnetized non-rotating flows the formation of filaments is reported within current and vorticity sheets [60]. The dominance of the horizontal motion (the motion in the plane perpendicular to the solid-body rotation axis) over the vertical one observed in PHT is altered in the presence of a "strong" imposed magnetic field. Indeed, for PMHDT, it is found that the ratio of horizontal to vertical kinetic energies is about one at the saturation stage (see Figure 6a). It is rather the horizontal motion that is most affected in the presence of a "strong" magnetic field especially at large scales (Figure 9b illustrates this).
At the exponential growth stage, the development of kinetic energy, dissipation, and production rates is not sensitive to the B intensity (B) while their magnetic counterparts increase as B is increased (see Figure 5a). At the saturation stage, the Alfvén ratio between kinetic and magnetic energies decreases as B is increased and for fixed B the contribution to total energy coming from the magnetic fluctuations are more important in the presence of precession than without precession (ε = 0, see Figure 5b).
The dynamics of the k = 0 mode and its nearest neighbors modes (k = 1, 2) were investigated in detail in addition to radial (spherically averaged) spectra of kinetic and magnetic energies for all the runs presented in the present study. In fact, the state of turbulence at the saturation stage of PHT is 2D-3C (C for componality) as in rapidly rotating nonmagnetized flows and is mainly due to the dynamics of the spectral buffer layer around k = 0 [26].
The radial spectrum of kinetic energy E (κ) (k) at the saturation stage of PHT follows the k −2 scaling for k ranging between the precessional k S = 4 and rotational k Ω = 40 Zeman's wavenumbers and there is an inverse cascade of energy at k < 8 (see Figure 8b). For the case of PMHDT, the kinetic energy is reduced at all scales especially at large ones, and the wavenumbers range at which E (κ) (k) is compatible with the k −2 scaling is also reduced when compared to the case of PHT. The amplitude of the magnetic energy is much smaller than the one of kinetic energy in the large scales but of the same order in the small scales (see Figure 8d).
The comparison between the spectra E (κ) (k ⊥ ), E (κ) (k ⊥ , k = 1, 2), and E κ) (k ⊥ , k = 0) at the saturation stage of PMHDT under a "strong" imposed magnetic field reveals that, at large horizontal scales, the major contribution to E (κ) (k ⊥ ) does not come only from the mode k = 0 but also from the k = 1 mode which is the most energetic. Only at very large horizontal scales at which E (κ) (k ⊥ ) ∼ E (κ) 2D (k ⊥ ), the flow is almost two-dimensional (see Figure 10b). Unlike the velocity field the magnetic field remains strongly three-dimensional for all scales since E (m) 2D (k ⊥ ) E (m) (k ⊥ ) (see Figure 11a). Note that the wavenumbers range at which the spectra E (κ) (k ⊥ ) and E (κ) (k ⊥ , k = 0) respectively follow the scaling k −2 ⊥ and k −3 ⊥ is reduced when compared to the PHT case. The behavior, scale by scale, of the Alfvén ratio between kinetic and magnetic energies is also investigated. For the case of RMHDT, the DNS results are in agreement with the linear theory solution: Γ(k) = E (κ) (k)/E (m) (k) = (1 + 2L 2 k )/(2L 2 k ) where L k = Bk/(2Ω). This is so that Γ(k) behaves like k −2 for large scales. A similar behavior is found for Γ(k) versus L S = Bk/(2εΩ) and also for Γ(k ) = E (κ) (k )/E (m) (k ) versus L S = Bk /(2εΩ).