Energy-flux vector in anisotropic turbulence: application to rotating turbulence

Energy flux plays a key role in the analyses of energy-cascading turbulence. In isotropic turbulence, the flux is given by a scalar as a function of the magnitude of the wavenumber. On the other hand, the flux in anisotropic turbulence should be a geometric vector that has a direction as well as the magnitude, and depends not only on the magnitude of the wavenumber but also on its direction. The energy-flux vector in the anisotropic turbulence cannot be uniquely determined in a way used for the isotropic flux. In this work, introducing two ansatzes, net locality and efficiency of the nonlinear energy transfer, we propose a way to determine the energy-flux vector in anisotropic turbulence by using the Moore--Penrose inverse. The energy-flux vector in strongly rotating turbulence is demonstrated based on the energy transfer rate obtained by direct numerical simulations. It is found that the direction of the energy-flux vector is consistent with the prediction of the weak turbulence theory in the wavenumber range dominated by the inertial waves. However, the energy flux along the critical wavenumbers predicted by the critical balance in the buffer range between in the weak turbulence range and the isotropic Kolmogorov turbulence range is not observed in the present simulations. This discrepancy between the critical balance and the present numerical results is discussed and the dissipation is found to play an important role in the energy flux in the buffer range.


Introduction
It is one of the most important subjects in turbulence research how energy, which is provided by external force and dissipated by viscosity, is transferred among scales via nonlinear interactions. In the Kolmogorov turbulence, the energy cascades from large-scale eddies to the small-scale ones via nonlinear interactions. In the weak-wave turbulence, the energy is transferred via nonlinear resonant interactions among waves.
The assumption of the weak nonlinearity that the linear time scale is much smaller than the nonlinear time scale is violated at small or large wavenumbers in almost all the wave turbulence systems (Biven et al. , 2003Newell et al. 2001). In this case, the weak-wave turbulence and the strong turbulence coexist (Meyrand et al. 2018;Vinen & Niemela 2002;Yokoyama & Takaoka 2014). According to the conjecture of the critical balance (Goldreich & Sridhar 1995;Nazarenko 2011;Nazarenko & Schekochihin 2011), the energy is considered to be transferred along the wavenumbers at which the linear wave period is comparable with the eddy turnover time of the isotropic Kolmogorov turbulence in the buffer range between the wavenumber ranges of the weakwave turbulence and the isotropic Kolmogorov turbulence. The conjecture is being eagerly 2 N. Yokoyama and M. Takaoka Schematic energy flux in rotating turbulence. The wavenumbers of the two-dimensional vortex are indicated by the thick purple line on the k ⊥ axis. The wavenumber ranges of the weak-wave turbulence, the buffer range and the isotropic Kolmogorov turbulence are respectively coloured red, green and blue.
The scale-by-scale energy cascade is often investigated by energy flux. The constancy of the energy flux in the wavenumber space is intensively examined as a corollary of the energy cascade in the homogeneous isotropic turbulence (HIT). While the Kolmogorov theory predicts scaling properties, the quantitative feature of the flux for sufficiently large Reynolds numbers has been examined theoretically and numerically.
However, the energy flux in anisotropic turbulence is less elucidated because the analytical expression for the flux is not known as contrasted in HIT. While the flux can be treated as a scalar in HIT, it should be treated as a vector in anisotropic turbulence even if it is homogeneous in the real space. In this paper, the definition of the energy-flux vector is proposed, and it is applied to rotating turbulence.
Let us consider a homogeneous anisotropic turbulence system which has one distinguishing direction, say the z direction, and is statistically isotropic in the x and y directions perpendicular to the z direction. It is convenient to investigate the energy transfers in the k ⊥ -k plane, where k ⊥ = |k ⊥ | = (k 2 x + k 2 y ) 1/2 and k = |k z |. In figure 1, the energy flux expected in the rotating turbulence is schematically drawn. The rotating turbulence is a typical turbulence system where the weak-wave turbulence of inertial waves and the isotropic Kolmogorov turbulence of eddies as well as the two-dimensional columnar vortex coexist (Clark di Leoni et al. 2014;Yokoyama & Takaoka 2017). In the wavenumber range where the linear period of the inertial wave is shorter than the eddy turnover time, the weak turbulence of the inertial waves is dominant (WT in figure 1). According to the weak turbulence theory, the resonant interactions among the inertial waves transfer energy to waves which have small scales perpendicular to the rotational axis without changing its scales parallel to the rotational axis (Bellet et al. 2006;Galtier 2003;Yarom & Sharon 2014). On the other hand, the isotropic energy transfer due to the isotropic Kolmogorov turbulence appears in the wavenumber range where the Coriolis period is longer than the eddy turnover time (Mininni et al. 2012), i.e., the wavenumber is larger than the Zeman wavenumber k Ω (KT in figure 1). However, there is no concrete theory that quantitatively gives the energy transfer in the buffer range (green in figure 1). The critical balance predicts the isotropisation due to the redistribution of energy that is transferred anisotropically to the buffer range. If the local energy-flux vectors can be obtained, the arrows of the energy flux in the buffer range are added to figure 1.

3
The energy flux in HIT is defined as the flux going through the sphere with a radius |k| = k. To evaluate the energy flux in the anisotropic turbulence, the fluxes going through the cylindrical surface with a radius |k ⊥ | = k ⊥ and through the planes with |k z | = k have been used as a simple tool in anisotropic turbulence (Deusebio et al. 2014;Lindborg 2006;Scott & Arbic 2007). The conical energy flux, which evaluates the flux going through the surface of a cone in the wavenumber space, was also proposed by Sharma et al. (2019). The perpendicular energy flux, for example, is obtained by integration over all the parallel wavenumbers, and it corresponds to the flux going through a line parallel to the parallel-wavenumber axis in figure 1. These integrated energy fluxes inevitably consist of contributions from multiple turbulence ranges, the weak-wave turbulence, the buffer range, and the isotropic Kolmogorov turbulence. Thus, the local energy flux is expected to be identified to see the paths of the energy transfer.
The locality of the net energy transfer in the wavenumber space is required for the concept of the flux to be reasonable. The locality is one of the major concept naturally assumed in the cascading theory, originating from Richardson's poetical note (Richardson & Lynch 2007) and sophisticated by Kolmogorov theory (K41) (Kolmogorov 1941). On the other hand, the locality as well as the correction of some aspects of K41 has been widely studied in the literature, starting with the pioneering work of Kraichnan (1959Kraichnan ( , 1971. It should be noted that the nonlocal energy transfer should be distinguished from the energy transfer in a nonlocal triad that is a flat triad (Waleffe 1992(Waleffe , 1993. The sweeping effect, in which small-scale eddies are advected by a large-scale eddy, does not make the net nonlocal energy transfer (Ohkitani & Kida 1992). In fact, Sharma et al. (2019) reported that the net energy transfer is mainly local in the rotating turbulence when the triad interactions are considered as one-to-one interactions. According to the net locality, diffusion models given by partial differential equations in the wavenumber space are often employed in the weak turbulence (Dyachenko et al. 1992;Galtier et al. 2019;Hasselmann et al. 1985;Zakharov et al. 1992) and in the Kolmogorov turbulence (Connaughton & Nazarenko 2004;Leith 1967). The diffusion models are also applied to anisotropic turbulence systems (Galtier & Buchlin 2010;Matthaeus et al. 2009). In this paper, the energy-flux vector is obtained by using the Moore-Penrose inverse, to identify the direction of the energy flux in the anisotropic turbulence. The organisation of this paper is as follows. The procedure to obtain the energy-flux vector is proposed, and it is verified by examining the energy-flux vector in HIT in §2. The energy-flux vectors in the rotating turbulence are presented in §3. The validity of the proposed vector in anisotropic turbulence is examined by comparing it with the weak turbulence theory. The energy flux for a more commonly-used external forcing is shown in detail. It is consistent with the weak turbulence theory in the wavenumber range where the energy spectrum agrees with the weak turbulence theory, but is observed to be inconsistent with that predicted by the critical balance in the buffer range in this simulation. In §4, the reason for the inconsistency between the fluxes in the present simulation and in the critical balance is clarified. The last section is devoted to the summary.

Formulation of energy-flux vector
In this section, the energy-flux vector is proposed by introducing two ansatzes. The procedure to numerically obtain the energy-flux vector is described with a concrete example, which is the rotating turbulence. The rotating turbulence is a typical homogeneous anisotropic turbulence system that has been extensively investigated, and can be an appropriate testbed to examine the proposed idea, since it contains different kinds of 4 N. Yokoyama and M. Takaoka turbulence (see figure 1). In addition, it is easy to examine the direction of the energyflux vectors in the range of the inertial-wave turbulence.

Integrated energy fluxes
The governing equations for the velocity u of rotating turbulence in incompressible fluid are the Navier-Stokes equation with the Coriolis term and the divergence-free condition: where the centrifugal force is included in the pressure p. The rotation vector Ω = Ωe z is assumed to be constant. The governing equations of the isotropic turbulence are the same as (2.1) but Ω = 0. The kinematic viscosity is expressed by ν. The external force f is Gaussian white. Under the periodic boundary condition, the governing equation (2.1) is rewritten in the wavenumber space as Energy is transferred among wavenumbers via the nonlinear interactions which come from the advection term in (2.2a). The energy transfer rate for a wavenumber k via the nonlinear interactions among three-wavenumber modes is symbolically written as where u k and E k = |u k | 2 /2 are respectively the velocity and the energy of the wavenumber k, and · represents the ensemble average. Note that the energy transfer rate T k is statistically equal to the difference between energy input by the external force and the energy dissipation rate in the statistically steady states. The triad-interaction function T (u k ; u k1 , u k2 ) quantifies the energy transfer from or to the wavenumber k via the triad k + k 1 + k 2 = 0, and in both isotropic turbulence and rotating turbulence studied in this paper. Here, (1 ↔ 2) represents the terms with the suffices 1 and 2 interchanged in the preceding ones. The triad-interaction function T (u k ; u k 1 , u k 2 ) is symmetric under the interchange of k 1 and k 2 , and the energy detailed balance T (u k ; u k 1 , u k 2 ) + T (u k 1 ; u k 2 , u k ) + T (u k2 ; u k , u k1 ) = 0 holds. The triad-interaction function T (u k ; u k1 , u k2 ) is the sum of the energy transfer between k and k 1 and that between k and k 2 .
To confirm the cascade theory in the homogeneous isotropic turbulence (HIT), the energy flux is usually examined. The isotropic energy transfer rate T (k) is obtained from T k by integration over the solid angle of k, and is assumed to be a continuous function of k = |k|.
Energy-flux vector in anisotropic turbulence

5
Anisotropic turbulence generally has a distinguishing direction, which is, for example, the direction of the rotational axis in the rotating turbulence. The z direction is set to be such distinguishing direction in this paper. The x and y directions are perpendicular to the distinguishing direction. Suppose that the system statistically has the azimuthal symmetry with respect to the z direction, i.e., the azimuthal isotropy in the x-y plane. Then, the statistical quantities can be described in the k ⊥ -k plane, where k ⊥ = (k 2 x + k 2 y ) 1/2 and k = |k z |.
As a natural extension of the isotropic energy flux (2.5) to azimuthally symmetric turbulence, the energy fluxes perpendicular and parallel to the system's distinguishing direction have been used (e.g., Alexakis et al. 2007). The anisotropic energy transfer rate T (k ⊥ , k ) is obtained from T k by integration over the azimuthal angle of k ⊥ and the sign of k z . Because P ⊥ (k ⊥ ) and P (k ) are respectively obtained by the integration over k ⊥ and k , detailed local structures such as critical balance in the wavenumber space cannot be captured directly by these integrated energy fluxes. These energy fluxes, (2.5) and (2.6), are referred to as integrated fluxes in this paper.

Minimal-norm energy-flux vector
To quantitatively investigate the energy-transfer mechanism in anisotropic turbulence, the detailed structure local in the wavenumber space of the energy flux is needed to be investigated. In this paper, the scalar-valued energy flux defined in HIT is extended to a vector-valued energy flux in the anisotropic turbulence. Because of the energy cascade, the definition of the energy flux in HIT (2.5) implicitly assumes the net locality of the nonlinear interactions and the local energy conservation in the wavenumber space: As an extension to the energy flux in the anisotropic turbulence, we assume the local energy conservation: where div k is the divergence operator in the wavenumber space, and P k = (P xk , P yk , P zk ) is the energy-flux vector. In general, the nonlinear interactions due to the advection term contain both local and nonlocal interactions in the wavenumber space. However, when we try to draw the flux as a vector field, we must implicitly consider the flux to represent the local interactions. Therefore, the ansatz of the local energy conservation is a natural consequence from the present purpose to find the energy-flux vector. The local energy conservation (2.8) can be interpreted as an alternative expression of the net local interactions to the diffusion models where the energy transfer in the wavenumber space is approximated by partial differential equations (Galtier & Buchlin 2010;Matthaeus et al. 2009).
The origin of the anisotropy of the energy flux should be described here based on (2.8). The anisotropy of the energy-flux vector comes from that of the energy transfer rate. In the statistically steady state, the anisotropy of the energy transfer rate can be 6 N. Yokoyama and M. Takaoka explained in terms of the energy input due to the external force and the dissipation rate due to the viscous term. The velocity and hence the energy can be statistically anisotropic owing to the system's anisotropy. In the rotating turbulence, the anisotropy of the energy is remarkable at the small wavenumbers. Thus, the energy input given by the inner product of the velocity and the external force is anisotropic owing to the anisotropy of the velocity even if the external force is isotropic. The dissipation rate at small wave numbers is also anisotropic owing to the anisotropy of the energy even if the viscous term and the dissipation rate at the large wavenumbers outside of the isotropic Kolmogorov turbulence range are isotropic. The energy-flux vector is determined by the energy transfer rate all over the wavenumber domain. It is similar to the pressure obtained by solving the Poisson equation in the real space. Therefore, the energy-flux vector is anisotropic in the anisotropic turbulence even if the external force and the viscous term are isotropic. The local energy conservation (2.8) is extended to a discrete formulation that is convenient in the numerical analysis. Let us consider the energy balance in a coarsegrained cell that has the side lengths ∆k x , ∆k y and ∆k z in the wavenumber space as shown in figure 2. The energy balance in the (i, j, k) cell is obtained from the local energy conservation (2.8) as The energy fluxes incoming to and outgoing from the cell, P i±1/2,j,k x , P i,j±1/2,k y and P i,j,k±1/2 z , are defined on the cell faces, and the energy transfer rate of the cell T i,j,k is given as the sum of the energy transfer rates of the wavenumbers in the cell. It is represented in a matrix-vector form as where D ∈ R NxNyNz×(3NxNyNz−NxNy−NyNz−NzNx) corresponds to the difference in (2.9) derived from the divergence operator in (2.8), P ∈ R (3NxNyNz−NxNy−NyNz−NzNx)×1 is a solution column vector of the flux, and T ∈ R NxNyNz×1 is the column vector of the transfer rates. The numbers of the coarse-grained cells in the x, y and z directions are respectively N x , N y and N z . The number of the components of P is approximately three times larger than that of T, since P k consists of three components in the threedimensional (3D) turbulence while T k is a scalar. Note that the number of the components Energy-flux vector in anisotropic turbulence of P is smaller than three times that of T, because the energy flux to or from the outer range of the computational domain does not exist. Obviously, the divergence matrix D has linearly independent rows. Thus, the solution of (2.10) is not unique. The nature often adopts the most efficient way under constraints. For example, the minimal energy state where the Euclidean norm of the velocity vector is minimal is realised under the condition that the vorticity invariants are conserved in Euler flows (Vallis et al. 1989). The minimal-norm flux proposed here has the minimal Euclidean norm under the condition that the energy transfer rate is provided by the nonlinear interactions. The selection of the minimal-norm vector is the least-action principle in this system where the Euclidean norm of the energy-flux vector is considered as an action. We here introduce this principle to uniquely determine energy-flux vectors, the solution of (2.10).
It is the Moore-Penrose inverse, which is a generalised inverse, that can find an appropriate flux vector among the infinite number of the solutions. For D having linearly independent rows, the Moore-Penrose inverse is defined as D + = D T (DD T ) −1 . The Moore-Penrose inverse selects the solution as P * = −D + T = argmin P 2 such that (2.10) holds. Namely, P * selected by the Moore-Penrose inverse has the minimal Euclidean norm among the infinite number of the solutions of (2.10).
In addition, the minimal-norm solution P * is irrotational. Because the divergence of the rotation is 0 and because every null-space component of the matrix having linearly independent rows is orthogonal to any of the row-space component, P * does not have the null-space solenoidal component. Thus, the use of the Moore-Penrose inverse to (2.10) is equivalent to the assumption that the energy-flux vector is irrotational. The selection of the minimal-norm flux can be interpreted that the energy transfer is "efficient" in the sense that the minimal-norm flux excludes local circulations of the energy transfer.
The minimal-norm solution P * is obtained by numerically solving (2.10) as follows. The energy transfer rate for each wavenumber mode T k in the statistically steady state is obtained in the direct numerical simulation (DNS). The column vector of the transfer rates T is composed by coarse-graining of T k . By using the generalised minimal residual (GMRES) method, (DD T ) −1 T is obtained. For the convergence criterion, the relative error is below 10 −10 . The minimal-norm vector P * is obtained by applying −D T to the vector obtained in the previous step. Once P * is obtained, the energy-flux vector of the wavenumber k located at the center of the (i, j, k) cell is given as The energy-flux vector in the system with one distinguishing direction and the statis-8 N. Yokoyama and M. Takaoka tical isotropy in the directions perpendicular to the distinguished direction is reduced to the two-dimensional (2D) energy-flux vector in the k ⊥ -k plane. The perpendicular and parallel components of the energy-flux vector by averaging over the azimuthal angles and the signs of k z . Let us evaluate the perpendicular component of the energy-flux vector going through an arc in a coarse-grained cell represented by the red curve in figure 3. In this example, the incoming flux P ⊥in is the sum of the energy fluxes through the sides inside the arc, P i−1/2,j,k xin and P i,j−1/2,k yin . Here, the flux on the cut-cell edge P i−1/2,j,k x , for example, is divided into P i−1/2,j,k xin and P i−1/2,j,k xout according to the divided lengths. The outgoing flux P ⊥out is similarly obtained. The energy flux through the arc is obtained by a weighted average of the incoming and outgoing fluxes as where S in and S out respectively denote the areas inside and outside the arc, and S in + S out = ∆k x ∆k y . The perpendicular component of the energy-flux vector is obtained by averaging over these arcs as well as the sign of k z . The parallel component of the energyflux vector is obtained from the z component of the energy flux by averaging over the azimuthal angles and the signs of k z .

Energy-flux vector in homogeneous isotropic turbulence
To confirm the consistency of the minimal-norm energy-flux vector with the integrated energy fluxes (2.5) and (2.6), the minimal-norm energy-flux vector is obtained from the DNS of the well-known HIT. The DNS is performed with 512 3 grid points in a cubic box whose volume is (2π) 3 . The pseudo-spectral method with aliasing removal by the phase shift is employed to evaluate the nonlinear term, and hence the maximal wavenumber is approximately 512 √ 2/3 ≈ 240. The Runge-Kutta-Gill method is used for the time integration. The external force generated by the white noise is added in the wavenumber space to the wavenumber mode in k f −1/2 |k| < k f +1/2, where the forced wavenumber k f is set to 4 in this simulation.
In the statistically steady state, the nonlinear energy transfer rate for each wavenumber mode T k is obtained. The column vector of the energy transfer rates in the coarse-grained cell with the side ∆k x = ∆k y = ∆k z = ∆k = 3 is composed, and the solution vector of the energy flux is obtained as the minimal-norm vector. The energy-flux vector is converted to the 2D vector in the k ⊥ -k plane by using the azimuthal average. The vector in HIT is examined to have only the radial component in the k ⊥ -k plane. Once the isotropy of the minimal-norm energy-flux vector in HIT is verified, the direction of the energy flux in anisotropic turbulence can be examined by the minimal-norm vector below.
The 2D energy-flux vectors in HIT are drawn in figure 4. The 2D energy-flux vectors are obtained from the 3D energy-flux vectors by averaging over the azimuthal angles and the signs of k z as written in §2.2. Figure 4(a) is the enlarged view in the small wavenumber range, while the vectors in a whole computational domain are drawn by eliminating some vectors for visibility in figure 4(b). The energy-flux vectors radiate outward at all the wavenumbers as expected by the forward cascade of energy and its isotropy. The magnitudes of the energy-flux vectors are large near the origin, and become small as the magnitudes of the wavenumbers k = (k 2 ⊥ + k 2 ) 1/2 become large, because the areas of the spheres on which the flux is evaluated are proportional to k 2 and the energy conservation holds under (2.8).
The direction of the energy flux is quantitatively evaluated to justify the validity of   PSfrag replacements Figure 5. (a) Angle between P k and k measured in radian. The mean and the standard deviation are represented by the solid curve and the filled region, respectively. (b) Integrated energy fluxes in homogeneous isotropic turbulence. The integrated fluxes with the superscript v, P v ⊥ , P v , and P v , are obtained from the three-dimensional energy-flux vectors, while those without the superscript, P ⊥ , P , and P , are obtained according to (2.5) and (2.6).
the minimal-norm energy-flux vector. The angle between P k and k is defined as 13) and the angle is expected to be 0 in HIT. The mean as well as the standard deviation of the angle in the spherical shell k − ∆k/2 |k| < k + ∆k/2 is drawn in figure 5(a). The angles are small, that is, the energy-flux vectors are radial at all the wavenumbers. In particular, the energy-flux vectors in the inertial subrange are almost completely radial. The small but relatively large angles at the small wavenumbers are mainly due to the fluctuation of the random external force. The angles increase near the largest wavenumber owing to the boundary condition of the energy-flux vectors, but the magnitudes of the energyflux vectors are negligibly small. The isotropy of the energy flux in HIT is successfully validated.
To confirm the consistency of the minimal-norm energy-flux vector with the integrated energy fluxes, integrated energy fluxes are constituted from the energy-flux vectors. The integrated perpendicular flux P v ⊥ (k ⊥ ) is obtained by the summation of P x and P y through the arc (2.12) as well as the summation over k z . The integrated parallel flux P v (k ) is obtained by the summation of P z over k x , k y and the signs of k z . The isotropic energy flux P v (k) is obtained by averaging the energy fluxes through the inner and outer surfaces of the cells which cover the sphere with the radius k. These integrated energy fluxes constituted from the energy-flux vectors are compared with the integrated energy fluxes (2.5) and (2.6) in figure 5(b).
The perpendicular energy fluxes, P v ⊥ and P ⊥ , are almost equal to each other. The difference at the small perpendicular wavenumbers k ⊥ < 6 comes from the discreteness of the coarse-grained cell during the conversion from P x and P y to P v ⊥ . The conversion does not make the difference at the large wavenumbers because the difference between the incoming and outgoing fluxes due to discreteness of the cells is not so large there. Moreover, if the outgoing flux of the cell instead of the weighted average (2.12) is used as the flux to evaluate P v ⊥ , then P v ⊥ and P ⊥ are almost equal to each other at the small wavenumbers, though P v ⊥ is then shifted slightly to the smaller wavenumber at the large wavenumbers. Similarly, the radial energy fluxes, P v and P , are almost equal to each other, though the difference at the small wavenumbers also emerges owing to the discreteness of the cells. Because P v is evaluated exactly on the cell faces, and is not affected by the discreteness of the cells, P v and P are equal to each other all over the wavenumbers within the convergence criterion during the calculation by the GMRES method. Note that P v has its value at k i = 3(i − 1/2) owing to the coarse-graining with ∆k = 3 as written above, while P does at k j = j − 1/2, where i and j here denote positive integers. Therefore, these energy fluxes obtained from the energy-flux vectors are equal to the integrated energy fluxes except for the difference due to the discreteness of the cells. In this way, the magnitudes of the energy-flux vectors obtained by the Moore-Penrose inverse are quantitatively consistent with the integrated energy fluxes in HIT.

Application to strongly rotating turbulence
In this section, the energy flux as the minimal-norm solution of the continuity equation of energy is examined in rotating turbulence by comparing the vector with the theoretical predictions. Direct numerical simulations of rotating turbulence are performed according to (2.2) by using the pseudo-spectral method. In the following simulations of the rotating turbulence, the periodic box has dimensions of 2π × 2π × 8π so that k x , k y ∈ Z and k z ∈ Z/4. Here, the periodic box long in the z direction is used because of the anisotropy at the large scales. The non-dimensional numbers which characterise the rotating turbulence are the turbulent Reynolds number Re t = ε 1/3 /(νk 4/3 f ) = (k η /k f ) 4/3 and the rotational Reynolds number Re Ω = ε/(νΩ 2 ) = (k η /k Ω ) 4/3 , where ε is the energy dissipation rate and k η is the Kolmogorov wavenumber.

Theoretical prediction of energy flux in rotating turbulence
When the nonlinear term, the viscosity and the external force in the governing equation (2.2) are neglected, the linear inviscid equation can be written as where the complex amplitude a s k k is defined as a s k k = u k · h −s k k according to the helicalmode decomposition (Alexakis 2017;Galtier 2003;Smith & Waleffe 1999;Waleffe 1993), and s k = ±1 denotes the sign of the helicity of the inertial wave. The basis is expressed Energy-flux vector in anisotropic turbulence 11 as h s k k = (e 1 + is k e 2 )/ √ 2, where (e 1 , e 2 ) = (e z × k/|e z × k|, k × (e z × k)/|k × (e z × k)|) for k ⊥ = 0, and (e 1 , e 2 ) = (e x , e y ) for k ⊥ = 0. The linear dispersion relation is given by σ k = 2Ωk z /k. The linear inviscid equation (3.1) has the wave solutions u(x) ∝ h s k k e i(k·x−s k σ k t) + c.c. called inertial waves.
In the wave-dominant range, the period of the inertial wave is considered to be much shorter than the eddy turnover time, that is, 1/σ k ≪ (k 2 ε) −1/3 (Clark di Leoni et al. 2014). On the premise of the local nonlinear interaction, Galtier (2003) applied the weak turbulence theory to the inertial waves in the strongly rotating turbulence, and he found that only a small energy transfer along Ω is allowed by the resonance condition: (3.2) He also discussed the nonlocality of the nonlinear interactions, which generates the anisotropy. Waleffe (1993) applied his idea of the instability assumption on the nonlinear energy transfers to predict the anisotropic energy transfer among the wavenumber modes and obtained similar results. The energy flux parallel to the perpendicular wavenumber axis is theoretically expected in the wave-dominant range. At the wavenumbers larger than the Zeman wavenumber, where the Coriolis period is larger than the eddy turnover time, the rotation is negligible at such large wavenumbers. Thus, the isotropic Kolmogorov turbulence appears at the larger wavenumbers, and the energy flux is isotropic like that in HIT shown in figure 4.
The energy flux in the buffer range between in the weak turbulence range and the isotropic Kolmogorov turbulence range is expected to connect the energy flux parallel to the perpendicular wavenumber axis in the inertial-wave turbulence and the isotropic energy flux in the isotropic Kolmogorov turbulence. Nazarenko & Schekochihin (2011) predicted that the energy is transferred along the wavenumbers at which the period of the inertial wave is comparable with the eddy turnover time, after the energy is carried to such wavenumbers by the resonant interactions among inertial waves. That is, the critical balance predicts that the energy is transferred to large k 's in the buffer range in figure 1.

Comparison with the weak turbulence theory
In order to validate the minimal-norm energy-flux vector in anisotropic turbulence, it is compared with the energy flux in the weak turbulence theory. A direct numerical simulation of rotating turbulence where the external force is applied to the wavenumber modes k ⊥ ≈ 0 is performed to compare the energy-flux vectors directly with the prediction of the weak turbulence theory of the inertial waves. Here, the random external force is applied to the small perpendicular wavenumber modes k x , k y = 0, ±1 and |k z | 50 in a DNS with 256 × 256 × 1024 grid points. In this simulation, the rotational Reynolds number is evaluated as Re Ω ≈ 0.4. Although the turbulent Reynolds number is not well defined in this simulation because the forced wavenumbers are widely distributed, the turbulent Reynolds number is considered to be small to investigate the energy flux in the weak turbulence. Time averaging to obtain the energy transfer rate T k required in the local energy conservation (2.8) as well as the energy spectra is performed in the statistically steady state.
In the weak turbulence theory, the energy spectrum of the inertial waves is predicted as E(k ⊥ , k ) ∝ k −5/2 ⊥ k −1/2 (Galtier 2003). To observe the anisotropic spectra, the kinetic energy spectra for each k as function of k ⊥ and for each k ⊥ as function of k are drawn in figure 6. The energy spectrum for each k as a function of k ⊥ , for example, is defined as where the summations k ′ ⊥ ′ and k ′ ′ are respectively taken over k ⊥ − ∆k ⊥ /2 |k ′ ⊥ | < k ⊥ + ∆k ⊥ /2 and k − ∆k /2 |k ′ z | < k + ∆k /2, and ∆k ⊥ and ∆k are the bin widths to obtain the spectrum. The corresponding integrated spectra as a function of k ⊥ are also drawn in figure 6. The perpendicular-wavenumber spectra of the energy are close to k −5/2 ⊥ in the wavenumber range of 4 k ⊥ 30 and 10 k 30. Similarly, the parallel-wavenumber spectra of the energy are close to k −1/2 in the same wavenumber range. Moreover, the parallel-wavenumber spectra have abrupt drops at k ≈ 50, and the energy injected by the external force is rarely transferred to the large parallel wavenumbers k > 50. It is consistent with the prediction of the weak turbulence theory of the inertial waves, in which the resonant interactions transfer the energy only to the wavenumbers having the same k . These spectra demonstrate that the inertial-wave turbulence appears in the wavenumber range of 4 k ⊥ 30 and 5 k 30. The energy spectra show that this DNS is appropriate to compare the energy-flux vector with the prediction of the weak turbulence theory.
The energy-flux vector in the rotating turbulence whose spectra are drawn in figure 6 is obtained by the procedure written in §2.2. The perpendicular and parallel components of the energy-flux vectors are respectively obtained from the x and y components and the z components by averaging over the azimuthal angles and the signs of k z . The energy-flux vectors in the k ⊥ -k plane are drawn in figure 7.
The weak-wave turbulence is expected to exist at the wavenumbers, where the linear wave period τ w = (2Ωk /k) −1 is much shorter than the eddy turnover time τ e = (k 2 ε) −1/3 . The critical wavenumber is evaluated to appear at τ w = τ e /3 in magnetohydrodynamic turbulence (Meyrand et al. 2016) and stratified turbulence (Yokoyama & Takaoka 2019). On the other hand, the Coriolis force affects little and the isotropic Kolmogorov turbulence appears at the wavenumber range where the Coriolis period τ Ω = (2Ω) −1 is larger than the eddy turnover time τ e . The separation wavenumber of the isotropic Kolmogorov turbulence is known as the Zeman wavenumber k Ω . The buffer range should exist between the critical wavenumber and the Zeman wavenumber. The curves which represent the critical wavenumber (red) and the Zeman wavenumber (cyan) are drawn in figure 7. The energy-flux vectors at the wavenumber modes k 50 in the weak turbulence range are almost completely parallel to the k ⊥ axis, and the energy provided by the external force rarely goes to the modes k > 50. Thus, the wavenumber modes at k > 50 in the weak turbulence range have little energy. The minimal-norm vector can well reproduce the anisotropic energy flux in accordance with the resonant interactions among the inertial waves locally in the wavenumber space.
The energy flux in the range where τ w < τ e /3 and the weak-wave turbulence is expected to exist does not necessarily demonstrate the perpendicular flux. It results from the fact that the wavenumber modes at k > 50 have little energy and the modes in the range are subordinate to the modes at k < 50 having much larger energy.
The energy-flux vectors turn the direction near the wavenumber modes having τ w = τ e /3. In this case, the energy flux along the critical wavenumbers that gives the isotropisation of energy in the buffer range is not observed. In fact, the energy transferred via the resonant interactions moves on to the 2D modes k = 0.
The integrated energy fluxes obtained from the 3D energy-flux vectors in the simulation where the external force is applied to the small perpendicular wavenumbers are almost equal to the corresponding integrated energy fluxes as drawn in figure 8. The difference between P v ⊥ and P ⊥ at the small perpendicular wavenumbers k ⊥ < 6 comes from the discreteness of the coarse-grained cell during the conversion from P x and P y to P v ⊥ as seen in the isotropic turbulence (figure 5). The difference between P v and P at k < 50 where the external force directly affects comes from the discreteness of the coarse-grained cell during the averaging the energy fluxes over the azimuthal angles or through the inner and outer surfaces of the cells. The negative parallel flux, P < 0, appearing at k < 50 where the external force affects shows that most of the energy provided by the external force goes to the small parallel wavenumbers, and is dissipated there.   PSfrag replacements Figure 9. Energy spectrum, that given only by two-dimensional modes (k = 0), and that given by three-dimensional modes (k = 0). The vertical line indicates the Zeman wavenumber.

Energy flux in direct numerical simulation with isotropic forcing
Another DNS where more commonly-used external force is employed is performed with 512 × 512 × 2048 grid points. One may think that the resolution is not so high as that in recent high-resolution simulations. However, such high-resolution simulation is not easily performed because a long-time integration is required to obtain the weak inertial-wave turbulence where the energy is transferred by the resonant interactions.
The 3D three-component random force is added isotropically to the small wavenumbers in k f − 1/2 |k| < k f + 1/2, where the forced wavenumber k f is set to 4. In this DNS, the turbulent Reynolds number and the rotational Reynolds number are respectively evaluated as Re t ≈ 1.6 × 10 2 and Re Ω ≈ 2.8.
The energy spectrum as a function of the norms of wavenumbers obtained in the DNS is drawn in figure 9. The energy spectrum at the wavenumbers smaller than those of the external force k f ≈ 4 is as steep as k −3 . The 2D flows, which are uniform in the direction parallel to the rotational axis, i.e., k = 0, account for a large fraction of the energy at the small wavenumbers, though the number of the 2D modes is much smaller than that of the 3D modes which have k = 0. In fact, a large-scale columnar vortex is formed. Note that a statistically steady state can be achieved without small-wavenumber drag (Chan et al. 2012). At the wavenumbers larger than k f , the energy spectrum is dominated by the 3D flows, and is less steep than at the small wavenumbers. In this energy spectrum, the isotropic Kolmogorov turbulence cannot be clearly observed, because the Zeman wavenumber is close to the Kolmogorov wavenumber as well as the cutoff wavenumber. The energy spectra are drawn in figure 10. The spectra are close to the prediction of the weak turbulence theory, E(k ⊥ , k ) ∝ k −5/2 ⊥ k −1/2 , in the wavenumber range of 4 k ⊥ 30 and 1/4 k 4. The abrupt drops at k ≈ 4 = k f in the parallel-wavenumber spectra indicates that little energy is transferred to the large parallel wavenumbers k > k f consistently with the weak turbulence theory. These spectra demonstrate that the inertial-wave turbulence appears in the wavenumber range of 4 k ⊥ 30 and 1/4 k 4. The energy-flux vectors in the k ⊥ -k plane in the isotropically forced rotating turbulence are drawn in figure 11. The energy-flux vectors roughly in the range of k ⊥ 30 and k 4 are almost parallel to the k ⊥ axis. The perpendicular flux is consistent with the resonant interactions of the inertial waves in the weak turbulence theory. In fact, the  range where the perpendicular flux is observed corresponds to the range where the weak inertial-wave spectra are observed in figure 10.
On the other hand, the vector field of the energy flux has a sink at (k ⊥ , k ) ≈ (100, 0). In this simulation, the energy flux along the critical wavenumbers predicted by the critical balance is not observed. This discrepancy is discussed in the next section.
To confirm the consistency of the energy-flux vectors with the commonly-used integrated energy fluxes, the integrated fluxes are obtained from the 3D flux vectors and compared with those obtained by (2.5) and (2.6) in figure 12. The isotropic and perpendicular fluxes obtained from the 3D flux vectors are almost equal to those obtained by (2.5) and (2.6) except for the small wavenumbers k, k ⊥ < 8. The difference at the small wavenumbers comes from the discreteness of the cells and the averaging for k ⊥ and k. The parallel flux obtained from the 3D flux vectors and the corresponding integrated flux are the same within the convergence criterion adopted in the GMRES method. These agreement and disagreement have been seen in the energy flux in HIT ( figure 5(b)). Therefore, the energy-flux vectors are consistent with the integrated energy fluxes, which have been usually used in the anisotropic turbulence.

Discussion
The minimal-norm energy-flux vector successfully shows the local structures of the energy flux both in HIT corresponding to the strong turbulence (figure 4) and in the anisotropic weak-wave turbulence ( figure 7). Thus, one can expect that the minimalnorm vector can be applied also in the buffer range between the weak turbulence range and the strong turbulence range. Contrasting to the prediction by the critical balance, however, the energy flux along the wavenumbers at which the linear wave periods and the nonlinear eddy turnover times are comparable is not observed in the present simulations. In fact, the energy flux appears to have a sink in the buffer range.
It is important to clarify the reason for the discrepancy between the numerical result shown in the present paper and the prediction of the weak turbulence theory. To find the reason, the 2D spectrum of energy dissipation rate in the DNS with the isotropic external force is drawn in figure 13. Note again that the dissipation rate and the energy transfer rate statistically balance at all the wavenumbers except for the forced wavenumbers in the ideal statistically steady states. The 2D spectrum is defined as where 1/(2πk ⊥ ) is introduced to easily compare the 2D spectrum with the isotropic one. The dissipation rate is large at the forced wavenumbers and at the wavenumbers of the large-scale columnar vortex which has k ≈ 0. In addition, the dissipation is remarkable in the buffer range between the wave-dominant range surrounded by the critical wavenumber and the range of the isotropic turbulence where the wavenumber is larger than the Zeman wavenumber. The energy provided by the external force is transferred to this buffer range by the resonant interactions as seen in figure 11. The critical balance predicts that the energy which reaches the buffer range must start to be isotropised by changing the direction of transfer in this range. However, the sink due to the strong dissipation is observed in this range instead of the isotropisation. In this sense, the dissipation plays a dominant role in the buffer range in the present simulations, and the large dissipation there is not considered in the theory of the critical balance. A large range of the isotropic turbulence at the large wavenumbers is required to determine whether or not the energy is transferred along the critical wavenumbers. The wave-dominant range cannot be small because the resonant interactions are relatively sparse in the wavenumber space. Then, the highresolution simulations where both of the wavenumber ranges of the weak turbulence and the isotropic turbulence are large enough are required.
Let us consider the asymptotic behaviour of the large dissipation which appears in the buffer range in the present simulations at the large Reynolds numbers. Suppose that the higher-resolution simulation has a smaller kinematic viscosity without changing the magnitude of the external force and that of the anisotropy, i.e., the system's rotation rate. The buffer range in the higher-resolution simulation remains at almost the same wavenumber range as that in the lower-resolution simulation, because the linear time scales are unchanged and the total energy dissipation rate and hence eddy turnover time are almost the same. On the other hand, the relatively large dissipation in the buffer range should become smaller. Then, some of the energy pass through the dissipation in the buffer range, and the energy flux along the critical wavenumbers that brings the isotropisation could be possible.
The critical balance obviously assumes that the dissipation affects only at the large wavenumbers outside of the isotropic Kolmogorov turbulence range. In fact, however, the dissipation is not localised; the dissipation spectrum in the inertial subrange of the isotropic Kolmogorov turbulence slowly increases as k 1/3 , for example. Thus, much higher-resolution simulation is required to realise the critical balance. Note that the largeeddy simulations and the hyper-viscosity should be carefully used because the change of the dissipation can alter the direction of the energy flux.
The critical balance predicts the isotropisation in the buffer range between in the weak turbulence range and the isotropic Kolmogorov turbulence range. Such ranges are characterised by the time scales of dominant physical mechanisms: the linear wave period, the characteristic time of the anisotropic system and the eddy turnover time. When the external force is applied to the small perpendicular wavenumbers, τ w < τ e /3 in the inertial-wave turbulence range as shown in figures 6 and 7. On the other hand, when the external force is applied isotropically to the small wavenumbers, the inertial-wave turbulence range goes into the range where τ w ∼ τ e as shown in figures 10 and 11. It is mainly due to the high directionality of the energy transfer in the inertial-wave turbulence. The wavenumber range of the inertial-wave turbulence depends on the wavenumbers at which the external force is applied as pointed out by Nazarenko & Schekochihin (2011). The misalignment of the inertial-wave turbulence might indicate the existence of key mechanisms other than the time scales in the energy flux in the anisotropic turbulence.
When the energy due to the external force is injected at the middle wavenumbers, the energy is expected to cascade to the smaller wavenumbers as 2D turbulence and to the larger wavenumbers as 3D turbulence (Alexakis & Biferale 2018). If the high-resolution simulations could be performed, such split energy cascades should be observed. In fact, the energy flux toward small k ⊥ at k ⊥ < 10 and k ≈ 0 can be observed in figure 7(a). This 2D inverse cascade suggests that the split energy cascade can be observed by using the energy-flux vectors proposed here.
The two ansatzes, the net locality of the energy balance and the efficiency of the energy transfer, are introduced to obtain the minimal-norm energy-flux vector. It is assumed by the former that the local energy transfer is dominant in the wavenumber space, and the local energy conservation holds. Such locality is essential when the detailed local structures of the energy flux is discussed. It is assumed by the latter that the minimalnorm vector among other possible solutions of the local energy conservation is irrotational and excludes local circulations of the energy transfer. This efficiency of the energy transfer is natural because the local circulations make closed loops of the energy fluxes, and do not intrinsically contribute the net energy flux. Although these ansatzes are introduced, the energy-flux vectors are consistent with the integrated energy fluxes, which have been conventionally used.

Summary
In anisotropic turbulence, the energy flux in the wavenumber space should be considered as a vector at each wavenumber mode, because different kinds of turbulence coexist inhomogeneously at each scale. The integrated energy fluxes which have been used to exhibit the anisotropy in the energy transfer cannot reveal detailed local structures of the energy flux, because these integrated fluxes inevitably consist of multiple turbulence ranges. In this work, the minimal-norm energy-flux vector obtained by using the Moore-Penrose inverse was proposed to uniquely determine the energy-flux vectors based on the two ansatzes: the net locality of the energy balance and the efficiency of the energy transfer. The latter is equivalent to the irrotationality of the energy-flux vectors.
The minimal-norm energy-flux vector is tested in both strongly nonlinear isotropic turbulence and anisotropic weak-wave turbulence. The minimal-norm vectors successfully show the isotropic energy flux that radiates outward in homogeneous isotropic turbulence. In the rotating turbulence, the minimal-norm flux in the weak turbulence range demonstrates the energy flux parallel to the perpendicular wavenumber axis. Namely, the minimal-norm vectors successfully display the energy flux consistent with the energy transfer due to the resonant interactions among the inertial waves according to the weak turbulence theory.
On the other hand, the energy flux along the critical wavenumbers predicted by the critical balance was not observed in the buffer range between in the weak turbulence range and the isotropic Kolmogorov turbulence range in the rotating turbulence. The inconsistency of the energy-flux vector proposed here with the critical balance does not indicate the failure of the proposed method but results from the appearance of the large dissipation in the buffer range. This inconsistency will be resolved when much higherresolution simulations are performed. In addition, a two-dimensional turbulence such as the quasi-geostrophic turbulence is more suitable to obtain the coexistence of the large wavenumber ranges of the weak-wave turbulence and the strong turbulence as well as the buffer range. The energy-flux vector in the two-dimensional turbulence will be reported elsewhere.
Numerical computation in this work was carried out at the Yukawa Institute Computer Facility, Kyoto University and Research Institute for Information Technology, Kyushu University. This work was partially supported by JSPS KAKENHI Grant No. 15K17971, No. 16K05490, No. 17H02860, No. 18K03927, and No. 19K03677.
The authors report no conflict of interest.