Anisotropic Wave Turbulence for Reduced Hydrodynamics with Rotationally Constrained Slow Inertial Waves

Kinetic equations for rapidly rotating flows are developed in this paper using multiple scales perturbation theory. The governing equations are an asymptotically reduced set of equations that are derived from the incompressible Navier-Stokes equations. These equations are applicable for rapidly rotating flow regimes and are best suited to describe anisotropic dynamics of rotating flows. The independent variables of these equations inherently reside in a helical wave basis that is the most suitable basis for inertial waves. A coupled system of equations for the two global invariants: energy and helicity, is derived by extending a simpler symmetrical system to the more general non-symmetrical helical case. This approach of deriving the kinetic equations for helicity follows naturally by exploiting the symmetries in the system and is different from the derivations presented in an earlier weak wave turbulence approach that uses multiple correlation functions to account for the asymmetry due to helicity. Stationary solutions, including Kolmogorov solutions, for the flow invariants are obtained as a scaling law of the anisotropic wave numbers. The scaling law solutions compare affirmatively with results from recent experimental and simulation data. Thus, anisotropic wave turbulence of the reduced hydrodynamic system is a weak turbulence model for strong anisotropy with a dominant k⊥ cascade where the waves aid the turbulent cascade along the perpendicular modes. The waves also enable an appropriate closure of the kinetic equation through averaging of their phases.


Introduction
The theory of weak wave dynamics (i.e., the stochastic theory of nonlinear wave interactions) has been extensively studied since the seminal works of Kadomstev [1], Galeev et al. [2], Zakharov and Filonenko [3] and more recently reviewed by Zakharov [4], Balk [5], Choi et al. [6] and Nazarenko [7].This theory remains one of the few areas where a mathematical framework exists with predictive capabilities for studying the energetics and dynamics associated with fluid turbulence.In this regard, the theory of weak wave turbulence was further explored in recent papers [8][9][10][11] in the context of unstratified rotating turbulent flows.The theory utilizes the incompressible Euler equation for inviscid dynamics in an infinite space, in dimensionless form: Here u = (u, v, w) is the three-dimensional velocity field in the Cartesian geometry x = (x, y, z), p is the pressure field, and ∇ = (∂ x , ∂ y , ∂ z ) is the gradient operator.Equation ( 1) is characterized by system rotation 2Ω z along with length, advective velocity, time and pressure scales denoted by L, U, L U and 2ΩLU, respectively.The Rossby number Ro = U 2ΩL , describing the relative importance of nonlinear advection to the Coriolis acceleration force, is the sole non-dimensional parameter.Of particular interest is the regime Ro 1 for rotationally constrained flows along with the concomitant viewpoint that the dynamics can be partitioned into fast inertial waves evolving on the rotational timescale of O((2Ω) −1 ), and eddies evolving on the advection timescale of O( L U ).In nondimensional units these timescales are denoted by O(Ro) and O(1) respectively.In the Cartesian coordinate system the linear inertial-wave dispersion relation, obtained from Equation (1) for Fourier plane waves with wavevector k = (k x , k y , k z ) and of the form exp{i(k • x − ω k t)}, is given by Slow dynamics, to leading order, are geostrophically balanced, i.e., ẑ × u ≈ −∇p.
Underlying hypotheses for the theory of rotating wave turbulence are: (i) the separation of inertial and advective timescales, i.e., Ro 1 and (ii) the non-interaction between geostrophically balanced and inertial waves dynamics [16].A necessary criterion for this to occur in Equation ( 1) is (|u|, p) = o(Ro −1 ) which ensures that the nonlinear terms remain small compared to linear terms.It follows from this bound that wave amplitudes can be large outside the slow manifold.Within the Ro 1 regime (slow manifold), laboratory experiments [17,18] and sufficiently spatially resolved simulations [11,19,20] have clearly demonstrated the tendency for the inertial wave spectra to evolve anisotropically towards a slow manifold associated with axially invariant geostrophic dynamics (k z = 0).Using wave turbulence theory on Equation (1), Galtier [8] predicts an anisotropic energy spectrum E(k ⊥ , k z ) ∼ k −5/2 ⊥ k −1/2 z and a helicity spectrum ⊥ solution for E k ⊥ observed in simulations [21,22].A critical concern for this discrepancy is that the uniformity of the asymptotic approach of the weak-wave turbulence theory is lost as the slow manifold is approached.Notably, as k z → O(Ro), it is found that inertial waves still exist within the slow manifold and are slow.Such slow waves are not accounted for in weak wave turbulence literature [8,9].
Spectral power laws for anisotropic wave turbulence in the limit of rapid rotation are obtained in this manuscript.Specifically, it is found that e k ∼ k −3 ⊥ (equivalently, E k ⊥ ∼ k −2 ⊥ and H k ⊥ ∼ k −1 ⊥ ) and are in agreement with power law solutions obtained numerically and experimentally [17,[21][22][23] (see Figure 1).These solutions pertain to the slow manifold and hence differ from the power laws for the isotropic flow regime [24] and those obtained by Galtier [8,9] for the flow state outside the slow manifold.A schematic picture of the various regimes of rotating turbulence is also provided in this paper along with the respective power law solutions for each regime.This will serve in broad understanding of rotating turbulent dynamics and specifically, undergird the importance of anisotropy in the slow manifold and its influence on spectral power law solutions.It is also demonstrated that there is, in principle, a simpler approach in obtaining the coupled energy-helicity kinetic equations using the method of Hamiltonian reduction [25] and symmetries than the more algebraically tedious correlation calculations in [8,9].For clarity and for the sake of averting any confusion, it is important to note that the derivations of the kinetic equations presented here are not the same as the derivations in [8,9].This point is further elaborated in Sections 6.4 and 6.5.
⊥ comport with the spectra obtained analytically in the current article as is explained in later sections.The inertial range in the above plots span roughly over an order of magnitude in the log scale, i.e., roughly 70 to 90 wave numbers.The line showing the −5/3 slope is just for reference as is clearly mentioned in [21].

Reduced-Rotating Hydro-Dynamic Equations, R-RHD
We apply a multi-scale perturbation method directly within the slow manifold to a recently derived set of asymptotically reduced equations for rotationally constrained flows [14,15,26,27].In what follows, we adopt the nomenclature reduced rotating hydrodynamic (R-RHD) for the reduced equations describing an unstratified, non-buoyant rotating fluid.We note that R-RHD capture a slow manifold that is more precisely identified as k z ∼ O(Ro) and contains not only geostrophic columnar eddies but also anisotropic inertial waves characterized by scales: k z /k ⊥ Ro 1.The amplitude of these slow waves evolve on slower advective timescale.
At the outset, for convenience to aid our investigation, we centrally locate some notations and definitions involving position vector x and wavenumber vector k, in Appendix A.1.The detailed derivation of R-RHD is provided in [14,15].To summarize, the asymptotic framework for Equation ( 1) is established by assuming the small expansion parameter Ro and a multiple-scale expansion in the axial direction ∂ z = ∂ z + Ro∂ Z with the isotropic scale z = z ∼ (x, y) and the anisotropic columnar length scale Z = Roz.Fluid variables v = (u, p) T , where T denotes tranpose, are now written as an asymptotic series in terms of the small parameter, Ro: To leading order in Equation ( 1), we observe a point wise geostrophic balance: z × u 0 = −∇p 0 .It follows that fluid motions are horizontally non-divergent, i.e., ∇ ⊥ • u 0⊥ = 0, with u 0⊥ = ∇ ⊥ ψ where p 0 = ψ is the geostrophic stream function as in the classical theory of quasigeostrophy.The Taylor-Proudman constraint [28] associated with the geostrophic balance further requires vertical variations to be negligible on O(1) vertical scales, i.e., ∂ z v 0 ≡ 0. Importantly, in compliance with the Taylor-Proudman constraint, the multiple scales approach permits variations of O(Ro −1 ) on the Z-scale [14,15].Hereafter, in the following we set k z = Rok Z .
At the next order in Equation ( 1), the requisite solvability conditions, ensuring non-secular behavior of v 1 , lead to the R-RHD equations describing the evolution of unstratified slow geostrophically balanced motions: Here u ⊥ = ∇ ⊥ ψ, and ζ := ∇ 2 ⊥ ψ and W are the z components of the vorticity and velocity fields.Akin to classical quasigeostrophic theory, nonlinear vertical advection, W∂ Z , is an asymptotically subdominant process and does not appear.However, unlike quasigeostrophic theory the velocity field is isotropic in magnitude with |u ⊥ | ∼ |W|, hence the appearance of a prognostic equation for W. Physically, the R-RHD ( 5) and ( 6) state that unbalanced vertical pressure gradients drive vertical motions that are materially advected in the horizontal, in turn, vortical stretching due to vertical gradients in W produce vortical motions.The vertical velocity, W, also generates an ageostrophic velocity fleld, . Consistent with the Euler Equation (1), the R-RHD also conserve, in time, the volume-averaged (over a domain B) kinetic energy E V and helicity H V : In what follows, an investigation of wave turbulence will be applied to the R-RHD.

Geostrophic Inertial Waves and Eddies
We observe that upon linearization, the R-RHD support slow geostrophically balanced inertial waves of the form, with planar phase function Φ(k, 1) with dimension of 1 T and can be set to 1 without loss of generality.More explicitly, the phase function can be split into two terms: e iΦ(k,s k ω k )) = e −is k ω k t e i(k ⊥ •x ⊥ +k Z Z) , the temporal phase factor and the spatial Fourier basis.The region of validity of R-RHD is k Z ≤ k ⊥ , of special interest is the limit k Z k ⊥ .When k Z = 0 this expression represents vertically invariant modes with ω k = 0 (i.e., the 2D modes of turbulent eddies).The circularly polarized wave amplitude vector is where s k = ± denotes the handedness, '+' for right-handed circularly polarized waves (with positive helicity) and '−' for left-handed circularly polarized waves (with negative helicity).Here, The superscriptˆrepresents Fourier coefficients from the Fourier series expansion of a function.In subsequent sections, we invoke the infinite box limit (L → ∞) and switch to Fourier transform variables thusly ĉk → L 2π 2 ĉk ≡ c k and δ(0) = L 2π 2 (see ch. 5 and 6 in [7]).The power 2 arises from the fact that the analysis is performed in 2D, i.e., x = (x ⊥ , ẑ), unless otherwise specified.

Helical Basis for Circularly Polarized Inertial Waves
We note that in wavenumber space the unit vectors ( k ⊥ , z, k ⊥ ), with k ⊥ = k ⊥ /k ⊥ and k ⊥ = k ⊥ /k ⊥ , form a right-handed orthogonal basis.Within the slow manifold, we have k ⊥ ↔ k as the direction of wave vector and henceforth we will simply call this vector k (see Figure 2).The leading order velocity field associated with an inertial wave is given by where Ûs k k := ik ⊥ ψs k k .From (10), Here h represents the complex helical wave basis (see ref. [29] for details) within the slow manifold incorporating the leading order incompressibility criteria as with its counterpart that exists outside the slow manifold (ref. Figure 2), this wave basis exhibits the following property that enables switching across different handedness by a conjugation operation, These findings illustrate that the R-RHD are naturally set up in the helical wave coordinate basis.As stated above in the paragraph following Equation (10), we switch from Fourier series to Fourier transform in the next section, i.e., e.g., k unless explicitly mentioned.This switch demands a multiplication by a factor ( L 2π ) d , d = 2 which will be explicitly written in the subsequent sections.

⊥
. The wave propagation direction is given by the wave vector, k (which we simply call k in the body of the manuscript).Within the slow manifold where

Dynamical Limits of Turbulence and Multiple Scales Order Parameters
Before we present the derivation of the amplitude and kinetic equations, a note on the time-scales of the different variables is necessary along with their order of magnitude.

Time Period of Oscillations of Waves and Amplitudes
First we introduce the scalings typically involved in weak wave turbulence [7] as applied to rotating turbulence.In terms of the time periods of oscillations, we have where T f is the fast inertial wave time period proportional to 1 [8] (The subscript f refers to fast waves).Θ is the rotation rate of the system.T τ sets the time scale for weak amplitude oscillations.Typically, With these scalings, the wave time period is taken as Ro ∂ ∂T τ .The wave kinetic equations are subsequently derived in an infinite box limit and by taking a large time limit integration as is explained in detail in [7,8,30].
The R-RHD equations introduced in the previous section eliminate the fast inertial waves with frequency ω f := 2π T f but retain slow inertial waves with frequency ω s := 2π T s which have the dispersion relation 9) (The subscript s refers to slow R-RHD waves).This leads to a slight modification of Equation ( 14) as T f T s T τ T τ NL which we simply write below as, where now we may have We have ω s ω f indicating the slowness of the R-RHD waves compared with the fast inertial waves in [8].However, for the multiple scales analysis carried forth in the subsequent sections, we adopt a new set of slow (or strained) time variables with appropriate order of magnitude as stated below.

Notations and New Slow Time Variables for the Multiple Scales Derivation Presented in This Paper
For reasons that will soon become clear, we adopt a suitable convention in our multiple scales approach by introducing new strained time variables τ, t s and t f [31][32][33] (As an analogy, one may think of t f as the seconds hand on a clock, t s as the minutes hand on a clock and τ as the hour hand on a clock.The multiple scales variables τ, t s and t f must not be conflated with the time period variables T τ , T s and T f as they have a distinct interpretation as stated here).The fast wave time scale is t f , the R-RHD slow waves have time scale of t s , while the wave amplitudes introduced in the next section have the slowest time scale denoted by τ.We fix our analysis by setting the most stretched time scale to order one and designate it as the slow time.In other words, we set our frame of reference to this slow time universe.Thus we have the following scale separated temporal variables:  [31][32][33] and is articulated as follows.Note that the frame of reference is set to the most slowly elapsing dynamic denoted by τ.As an example, say = 0.01.This means for the dynamic in the frame of reference to elapse by unit time (τ = 1), the R-RHD waves would have elapsed by 100 units (i.e., t s = 100) and the fast inertial waves would have elapsed even more rapidly as seen by an observer in the set frame of reference.We conveniently have the case that within the scope of the amplitude dynamic elapsing at τ, the slow wave dynamic time scale is already infinitely large.However, recall that the R-RHD have filtered out the fast inertial waves by applying a method of multiple scales with Ro as the order parameter as explained in the previous section.This means that from within the frame of reference of the reduced universe of R-RHD, the fast inertial wave is no more accessible and mathematically this means we do not need to track the dynamical universe elapsing at the time scale t f .This means that in our multiple scales analysis we only need to track τ and t s (we simply write t for t s except where we need to distinguish between t s and t f ) with τ t, or equivalently, τ = t, Ro 1.
Relation ( 16) is equivalent to T τ = 1 T s for the same choice of .So in the following we provide order of magnitude estimates for terms involving the slow inertial R-RHD waves (ω, t) and the non-linear amplitude time scale (τ ∼ τ NL ) for non-linear advection in the k ⊥ modes.Consequently, we have ω ≡ ω k = k Z k ⊥ ω 1 ∼ O( ) representing slow (slow in comparison to the fast inertial waves that are filtered out by the R-RHD) frequency R-RHD waves with ω 1 ∼ O(1).This gives us the following limits

Definitions of Turbulence Regimes
In what follows, we will generally abbreviate reference to weak wave turbulence applied to full hydrodynamic equations by Galtier [8] as WT, anisotropic wave turbulence of the reduced hydrodynamic equations developed here, and that developed for reduced magnetohydrodynamics (MHD) [7,34], as AT, and critical balance [27] as CB.Moreover, unless otherwise stated, we will primarily refer to the works by Zakharov et al. [30] and Nazarenko [7] for the section on wave turbulence closure.The terminology anisotropic wave turbulence for reduced hydrodynamic equations is borrowed from the MHD literature [7,34] where the model was developed for reduced MHD equations [7].This will be explained later in Section 6.2.Based on the order of magnitude scalings prescribed above, we have the following definitions. WT: AT: In essence, one may envision the AT dynamical regime for hydrodynamics as a zoomed in version of the WT limit in Galtier [8].In this paper we investigate the subtle variations in the flow dynamic at this magnified scale concealed by the treatment in [8].Further, in this paper, we perform the multiple scales analysis in the anisotropic regime defined by the order parameter = k Z k ⊥ 1.So the region of validity of AT is also set by the limit where As will be shown subsequently, the dynamic and solutions are distinctly different for WT and AT.

Wave Amplitude Equations
In the framework of this paper we consider small amplitude dynamics, where smallness of the amplitude is measured in terms of the extent of anisotropy, i.e., the ratio k Z k ⊥ 1.However, rotating fluid motion includes significantly large wave amplitudes of o(Ro −1 ) outside the slow manifold.For connectivity and consistency with observed energy spectra as prescribed in this paper, one can therefore envision the scenario whereby the amplitudes of resonantly cascading inertial waves have been sufficiently attenuated once the dominant flow dynamics have reached the slow manifold.We therefore proceed with a multiple scales asymptotic approach in time, where where to ensure consistency of order in the expansion above.The order parameter 0 < 1 is a measure of the wavefield amplitude and τ = t denotes the slow advective timescale for weak amplitude modulations in comparison to the inertial wave propagation time t.We note that ∂ τ ∼ u 1 • ∇ ⊥ , indicating that τ defines the advective timescale.This ensures a necessary separation of temporal scales τ and t as explained in [16], thereby allowing for a multi-scale treatment of the system with distinct dynamics at every order.The choice of an order parameter, proportional to k Z /k ⊥ has been mentioned by Nazarenko and Scheckochihin [27] and serves as a guiding principle for the theory developed here.This way a systematic multi-scale treatment that captures the anisotropy in the rotating system is developed.A multi-scale approach exploits the presence of scale separation between slow and fast dynamics.This allows the relevance of reduced models for understanding the coupling between slow-fast dynamics.The interested reader is referred to a recent paper by Abramov [35] where a reduced model for multi-scale dynamics is presented that illustrates the power of multiple scales method in a general setting.
At O( ), the leading order solution can be interpreted as a complex field of waves undergoing resonant interactions on slow inertial timescale.At leading order, O( ), utilizing Equation ( 12), we have that planar inertial waves satisfy Here  is the Hamiltonian matrix and I 2 is the identity matrix.The solution to this system is the helical base vector U It is important to note that L H being non-Hermitian, an extended eigen basis, including both the left and right eigenvectors, is required.This extended basis is the complex helical basis h s k k introduced earlier in Equation (12).Also of importance to the application of a solvability condition at higher asymptotic orders is the solution h Given system ( 19), a complex wave field can now be expressed as a superposition of inertial waves and eddies (2D modes) represented in terms of the helical basis: At O( 2), the next asymptotic order, we have for each helical mode h where the subscript 1 has been dropped from the right hand side of the above equations.Application of the solvability condition, 1 2 (h 2k , gives the wave amplitude equation in terms of the Fourier transform variables over the full k space as follows: where V s k s p s q kpq := s k s p s q + s p is the interaction coefficient.Here φ(ω) = (s k ω k − s p ω p − s q ω q ).The intermediary steps in the derivation of Equation ( 22) are given in Appendix A.2.

Dimensional Consistency of Wave Amplitude Equation (22)
For dimensional analysis, we use M, L, T to represent mass, length and time.Moreover, [•] is used as the dimensional operator, e.g., [W] (Here l.h.s. and r.h.s.mean the left and right hand sides of the concerned equation respectively).
For s k = +, the sum is actually carried over the combinations given by the set (s p , s q ) = {(+, +), (+, −), (−, +)} as the terms with c − p c − q have null contribution due to the resonance condition.The converse holds for s k = −.Since Ψ is real valued in the physical space, the corresponding Fourier coefficients satisfy the symmetry condition ĉs k −k = ĉs k * k and the integration over p, q is performed on the positive spectral space.With this physical information, we re-write Equation ( 22) in symmetrical form, for s k = +, as follows: where the interaction coefficient L+++ is clearly symmetric in the second and third arguments, i.e., L+++ kpq = L+++ kqp by noting that Hereafter, the integration over p, q > 0 is implied.The symmetrization of the nonlinear form is discussed in detail in ch.6 of the textbook on wave turbulence by Nazarenko [7].

Spectral Tensor
Having deduced the wave amplitude equation ( 23), we now proceed with the primary objective of finding the functional forms for stationary energy and helicity spectra, often referred to as Kolmogorov solutions [36].The relation between energy and helicity in terms of the canonical complex variable, c s k k can be understood by carefully analyzing the average spectral tensor for a given spectral mode k.
We review some definitions of terms at the outset.Henceforth, we use the symbol Λ := L 2π to denote the dimensional constant of dimension L. Ensemble average of the squared amplitude function enables us to estimate the energy spectral density.For a given polarity s k , we have k δ(0), where the infinite box limit entails δ(0) = Λ 2 [7].Here δ(k − k ) ≡ δ k,k is the Dirac delta function and δ k k is the Kronecker delta function (In this paper, delta function written with arguments both in the superscript and subscript refers to the Kronecker delta whereas the delta function written with arguments either all in the subscript or all within parenthesis refers to the Dirac delta.In other words, δ i j is Kronecker delta whereas δ(x − y) or δ x,y refers to Dirac delta).RPA stands for random phase and amplitude [7].Thus we have the following relation that is valid when L → ∞, For reference, we list the dimensions of the relevant terms again: The average spectral tensor is defined as follows (see p. 179, Section 5.11 in [37]): The ensemble averaging of the phase factor results in the appearance of the Kronecker delta function above.The Kronecker and Dirac delta functions may be interpreted as selector functions, i.e., f (x i , x j )δ i j = f (x i , x i ) and f (x )δ(x − x)dx = f (x) respectively.The arguments and results that follow in the rest of the paper must be interpreted in the statistical sense whereby the ensemble phase averaging can be regarded as a time averaging operation upon assuming ergodic property of the flow variables.From the relations in (25), it is clear that the following is true, By solving the above equations, we get, These expressions are particularly useful in the derivation of the stationary energy and helicity spectra.

Zero Helicity Dynamics
In this section, we analyze the special case of a flow with zero helicity, i.e., h k ≡ 0 for all k such that e + k = e − k = 1 2 e k .From (25) this imposes the constraint on the complex amplitude functions c k where for sake of brevity we will write c This entails a reflection symmetry of the wave field and a reduction in the Hamiltonian description of the system [25] because a unique handedness, associated with one of the s variables, now describes the full system as the system described by s = + and s = − are mirror replicas of one another in the statistical sense.This point is explained in more detail in Sen [38].Newell and Rumpf [39] have exploited this reduction and have studied resonance wave dynamics for the nonlinear Schrödinger's system.Now, the inertial waves may be interpreted as occurring in helicity couplets involving h + k and h − k = h + * k with wavefield given by

The Hamiltonian
The interaction Hamiltonian for the R-RHD can be expressed as a power series of the complex amplitudes, C k := c k e −iω k t , as follows: The H (4) component denotes four-wave resonant interactions which are negligibly small and therefore not considered.Resonant self-interactions involving quadratic terms, captured by H (2) , evolve on the inertial timescale t and is already accounted for within the operator L H .The full Hamiltonian is H ≈ H (2) + H (3) in agreement with the weak wave expansion (18).The wave amplitude Equation ( 23) is obtained from the well known Hamilton equation where and the three wave interaction Hamiltonian H (3) is as follows [7,30]: At the order of the weak nonlinear interactions, Hamilton's Equation ( 31) is identical in form to Equation (23) where δH (3)   δC * k is the r.h.s. of Equation ( 23) after invoking the aforementioned Hamiltonian Thus the wave amplitude equation of the reduced system (c L+++ kpq c p c q e i(ω k −ω p −ω q )t δ k,pq − L+++ qpk c * p c q e i(ω k +ω p −ω q )t δ q,pk dpdq.(33)

Kinetic Equations
By multiplying Equation (33) with c * k and the complex conjugate of Equation ( 33) with c k , subsequently subtracting the latter from the former, followed by ensemble averaging and using the results from Section 4.2, we obtain an evolution equation for e k ≡ e + k = e − k : where Λ := L 2π and C k = c k e −iω k t as stated earlier.Henceforth, we simply write Lkpq for L+++ kpq .Here, (• • •) refers to imaginary part of the argument in parenthesis.The triple correlation can be defined as follows: where the superscript ∆ is a convenient way of writing the Dirac delta function δ k,pq .On evaluating the ensemble average and considering the infinite box limit, we have the following reduction: Here δ k pq is the Kronecker delta function that takes the value 1 when k = p + q and 0 otherwise.The relation (36) 34) can be written as Lkpq J kpq δ k,pq − 2 Lqpk J qpk δ q,pk dpdq = Lkpq J kpq δ k,pq − 2 Lpkq J pkq δ p,kq dpdq . ( Here we have simply written Lkpq ≡ L+++ kpq and hereafter the superscript +++ will be implied unless otherwise mentioned.We observe that Equation ( 37) is not closed in the sense that the left hand side of the equation is a second order correlation function that is expressed in terms of third order correlation functions on the right hand side.For Gaussian distributed fields, odd correlators are null.Hence, we must devise an estimate of J kpq using a property of the fourth order correlator for Gaussian statistics as described in the following section.

Closure
In order to close Equation (37), we need to express the triple correlator in terms of second order correlators.Note that Equation ( 37) is obtained at O( 2 ) as explained in Section 4.An estimate of the triple correlator in terms of the second order correlator can be obtained at the previous order of the multiple scales perturbation, i.e., at O( ), as explained in this section.On applying Wick's theorem to the Gaussian distributed wave field, quadruple correlation functions are defined as follows [7,30]: The applicability of Wick's theorem relies on the assumption that the field is Gaussian distributed and this point is further discussed in Section 6.5.Let us first recall some relevant order of magnitude for the terms involved in this section.Since φ(ω) := (ω k − ω p − ω q ) ∼ O( ), we take ).We use the definition of J ∆ kpq (τ) in (35) and Equation (33) to write a differential equation for J kpq .The solution of this differential equation provides us an estimate of J kpq at O( ) to close the Equation (37) that was derived at O( 2 ).
The ensemble average of the phase factor can be succinctly captured by the Kronecker delta (This can be readily checked by assuming ergodicity and computing the ensemble average as a large time average), which is equal to δ k,pq Λ 2 at the infinite box limit, Λ := L 2π → ∞.Thus at O( ), for the triadic wave process, we have Lkpq (e k e p + e k e q − e p e q ) (using Equation (33) and def.(38) as shown below) (39) = C J ( constant in τ).
The second equality follows from the following calculation.Applying the product rule of differentiation, we have i . The first term on the r.h.s. is calculated as follows: − Lnmk : 0 (Wick's contraction rule [7]) Lkmn e m e n δ m p δ n q + δ m q δ n p δ k mn (using def.(38), = − Lkpq e p e q , where k = p + q.
Here, the integral is expressed as a sum to make the calculation explicitly clear for the reader.Note that dimensional consistency demands ∑ m,n ↔ Λ 4 dmn where Λ 4 appears here because of the relation δ m,p δ n,q + δ m,q δ n,p = Λ 4 δ m p δ n q + δ m q δ n p in the def.(38).Similarly, 1 Λ 2 c * k (i∂ τ c p )c q = Lkpq e k e q and 1 Λ 2 c * k c p (i∂ τ c q ) = Lkpq e k e p with k = p + q.Thus follows the second equality in Equation ( 40).The r.h.s. of Equaiton ( 40) is independent of time because we are interested in stationary solutions.The solution to the differential Equation ( 40) is , where κ is a constant.
By taking the limit T τ → ∞ for 1, ω s T τ ∼ 1 −T τ e iφ 1 (ω)τ dτ T τ →∞ −−−→ 0 when φ 1 (ω) = 0. Thus in the long time limit, the estimated average J kpq is expressed as follows: Thus, we have a reasonable estimate for J kpq in Equation (37) that is obtained as above at O( ).Recall that the O( ) solutions described in Section 4 are the slow R-RHD waves that enable this closure procedure.We now plug in this estimate J kpq for J kpq in Equation (37) to close the kinetic system.The occurrence of the singularity (φ 1 (ω) = 0) is averted by circumventing the pole by adding a small term i ˜ to the denominator.Then we use the identity: x−x 0 dx − iπ f (x 0 ), where P refers to the principal value.For the special case when f (x) = 1, the preceding identity is concisely written as where it is generally understood that the integral is taken over the domain with ˜ → 0 + .The imaginary part of the l.h.s., thus, yields the delta function term on the r.h.s. of the identity) and substitute the resulting term for J kpq in Equation (37) to obtain the closed form of the kinetic equation: | Lkpq | 2 (e p e q − e k e p − e k e q )δ k,pq δ ω k ,ω p ω q − 2| Lpkq | 2 (e k e q − e p e k − e p e q )δ p,kq δ ω p ,ω k ω q dpdq. ( Recall that δ(φ(ω)) ≡ δ ω k ,ω p ω q because φ(ω) := ω k − ω p − ω q .Equation ( 43), that describes the evolution of the energy spectral density, can be rewritten as a single integral with the interaction operator proportional to the square of Lkpq by using the Zakharov-Kuznetsov conformal transformation [30,40].5.2.2.Dimensional Consistency of the kinetic Equation (43) In Equation (43) Therefore, the derived kinetic Equation ( 43) is dimensionally consistent.

Physical Realizability of Spectrum Prescribed by Equation (43) and Comparison with Quasi-Normal Type Closures
Quasi-Normal family of closures (including modified versions like EDQNM) are well known statistical closures which are widely used in analytical study of turbulence [41][42][43][44].However, as pointed out by Orszag [42], the original Quasi-Normal closure results in an unphysical negative energy spectrum.This issue was shown to be related to the presence of time-history integrals in the evolution equation of the correlation function.This was corrected in modified versions of the closure, such as Quasi-Normal Markovian and Eddy Damped Quasi-Normal Markovian (EDQNM) closures [37,42,45].EDQNM is widely used in turbulence predictions and works well in the absence of waves [46].However, as shown by Bowman et al. [46], Bowman and Krommes [47] and noted in [48,49], in the presence of waves, especially Rossby, inertial and drift waves, EDQNM is unphysical in the sense that the energy spectrum, obtained by performing the closure, becomes negative for dynamically important Fourier modes (see Figures 1 and 2 in [46] and Figure 1 in [47]).This unphysical energy condition, in the presence of waves or mode coupling, is shown to stem from the violation of the realizability criterion θ kpq (t) ≥ 0, ∀t ≥ 0. This is because of the application of the fluctuation-dissipation ansatz that must be invoked in order to capture the two time correlation information within a single parameter, the triad interaction time θ kpq [50].Conceptually, this loss of realizability can be checked for a test case by assuming the total damping coefficient to be constant in time but complex valued (e.g., purely imaginary), whence θ kpq can be shown to be oscillatory about zero.In fact, Bowman [50] explicitly shows that the realizability condition is violated for an actual stochastic problem involving three-wave interactions (see ch. III in [50] for details).As a remedy, a realizable Markovian closure was proposed by Bowman et al. [46], Bowman and Krommes [47], Bowman [50].
In contrast, in this paper, a wave turbulence closure [7,30] is used as explained in Section 5.2.1.In this context, it is essential to check the physical realizability of the energy spectrum permitted by the kinetic system derived above.Following the strategy prescribed by Orszag [42] (cf.p. 313, ch.IV), we envision a scenario where, starting with a positive initial condition for the energy spectrum, we assume e k > 0 ∀k − {k * } and ∀τ ≤ τ * , e k * > 0 ∀τ < τ * , and e k * = 0 at τ = τ * .
Here ∀ should be read as for all.The above assumptions mean that physically we have a situation where the energy in a certain mode k * is null at time τ * .Specifically, note that e k * = 0 and e p , e q > 0 at τ = τ * .The goal is to investigate if the energy becomes negative for this specific wave mode k * at times τ > τ * .By inspection, it is clear that δ p,k * q δ ω p ,ω k * ω q e p e q dpdq > 0 (44) because e p , e q τ=τ * > 0 and the fact that the terms in the square parentheses are positive.The above inequality is obtained by substituting e k * = 0 in the r.h.s. of Equation ( 43) when evaluated at time τ * .
This means that the mode with null energy receives energy from the interacting modes of the triad and the energy of that mode becomes positive, i.e., the nonlinear transfer to the mode k * is positive and energy is dumped into this mode.The above argument can be extended to any wave mode and for all times τ > τ * in a straight forward manner.Thus, the spectrum prescribed by Equation ( 43), and obtained by the wave turbulence closure described in the previous section, is non-negative and physically realizable.This can be attributed to the absence of any time-history integrals in Equation (43).
In fact, the separation of time scales, captured by relation (15), enables us to take the long time limit and thereby time dependent oscillatory terms are eliminated as shown in the short paragraph following Equation (41).It is essential to note that the above inequality does not comport with the stationary state of the spectrum for the obvious reason that the derivative is non-zero.So the spectrum eventually adjusts itself to the stationary state permitted by the kinetic system derived earlier and for which the power law solutions are given in a subsequent section.

Invariants of the Closed Kinetic Equations
It can be shown easily that the total interaction energy given by Equation ( 43) is conserved, i.e., ∂ t e k dk = 0 using the result of Appendix A.3.The proof in the appendix is essentially a statement of conservation of energy in each wave triad.The form of the interaction coefficient ensures the convergence of the collision integral appearing on the r.h.s. of Equation ( 43), thereby it meets the first criteria for the realizability of a Kolmogorov spectrum, i.e., stationary energy spectrum solution to (43) [36].This point is elaborated in Section 5.2.6.

Kolmogorov Solution of the Kinetic Equations
The exact steady state solution of the kinetic Equation ( 43), as power laws, is obtained by assuming locality of the scale-by-scale energy transfer.This is illustrated in Nazarenko [7], Zakharov et al. [30] and earlier papers referenced therein.
Owing to the constrain imposed by the AT limit k Z k ⊥ , the turbulence cascade under investigation is in the k ⊥ direction and k Z = B r (k 0 ) (i.e., k Z is within a relatively small r vicinity of k 0 ) serves as a tracking parameter to ensure that we do not violate the condition: The four possible stationary solutions for the anisotropic spectrum, e k ∼ k −x i Z k −y i ⊥ , ∀k z = 0, are listed as follows [7,30]: (i) x 1 = 1 and y 1 = −1 (Rayleigh-Jeans (RJ) spectrum).(ii) x 2 = 1 and y 2 = 0. (iii) x 3 = (1 + u) and y 3 = (2 + v), where 2u and 2v are respectively the powers of k Z and k ⊥ in | Lkpq | 2 .Clearly, in our case, u = 0, v = 1.Thus, x 3 = 1 and y 3 = 3. (iv) x 4 = 1 and y 4 = 7/2.This solution corresponds to the constant flux in the z-component of the momentum.
Solution (iii), above, corresponds to the only constant energy flux solution, a necessary requirement of Kolmogorov's theory and of primary concern in this paper.The constancy of energy flux can be easily verified by noting that ; whereby on integrating with respect to k ⊥ and demanding constant energy flux in the perpendicular direction, we can extract the aforementioned solution.Here, Π denotes the flux of energy.The direction of the flux can be ascertained by comparing the power law solutions of the RJ spectrum with that of the constant energy flux solution, i.e., since −y 3 < −y 1 , the energy cascade is direct and the energy flux is towards smaller scales ( p. 139 in ref. [7]).Thus, the exact solution for the Kolmogorov-Zakharov-Kuznetsov spectra with constant energy flux is as follows: This result is in agreement with experimental and computational simulations [17,22,23].Additionally, LaCasce [51] has discussed the observation of a k −3 spectrum in the context of barotropic vorticity model using data obtained from several numerical simulations.Moreover, since lim k ⊥ →∞ e k dk converges for the spectrum (45), the Kolmogorov spectrum obtained above has finite capacity [7].The limit k ⊥ → ∞ is further discussed in the next section.

Locality of Interactions and Convergence of Collision Integral
As stated in the previous section, the Kolmogorov solutions obtained above are based on the assumption of locality of interactions.To validate the applicability of the locality assumption, it is sufficient to comment on the convergence of the collision integral on the r.h.s. of Equation ( 43) for the concerned power law solution [7,30,36].However, analyzing the convergence of the collision integral is a delicate task [52] and a rigorous commentary is not within the scope of this paper.However, in this section, we will argue on the convergence of the collision integral given by the r.h.s. of Equation ( 43) based on inspection of the terms in the integrand, the surface over which the integration is performed and numerical integration of the same within the region of validity of the AT limit of R-RHD.
The integral in Equation ( 43) is in six dimensions, viz., (• • •)dp x dp y dp Z dq x dq y dq Z if we further consider the decomposition of the perpendicular modes, k ⊥ = (k x , k y ).We have four delta functions, viz., three over the wave numbers and one over the frequency.Integrating over each delta function reduces the dimension of integration by one.So by substituting q = ±(k − p) in the respective integrands, we can completely eliminate integration over the three q wave numbers.Consequently, re-writing the frequency resonance condition in terms of the substituted expressions for q Z and q ⊥ , we finally have a surface integral over p x and p y dimensions where the surface is defined by p Z (p x , p y ).This is illustrated concisely as follows.
, where + is taken for the first integrand and − for the second.
Here S ≡ p Z (p x , p y ) defines the surface over which the integration is carried out when the collision integral is projected on the (p x , p y ) plane.Interestingly, and conveniently, the surface of integration turns out to be identical for both the first and second integrands when projected on the (p x , p y ) plane.
By inspection, the surface S is singular at k ⊥ = 0 and (k x − p x ) 2 + (k y − p y ) 2 − p 2 x + p 2 y = 0. Clearly, k ⊥ = 0 is outside the validity of AT because we must have k Z k ⊥ , k Z > 0, so we must not be concerned about the pole at k ⊥ = 0.The second condition requires more careful inspection.It defines the equation of a plane in the (p x , p y ) coordinates because it can be expressed as Figure 3 illustrates this singular plane on the surface S for two distinct limiting wavenumber configurations: α = 1 =⇒ p x + p y = k and α 1 =⇒ p y = k/2 (α 1 is just the opposite case).Equation (47) can also be re-written as follows, This implies that for the first integral corresponding to the collision integral in Equation ( 43), we have q ⊥ = 1 2 k ⊥ because of δ k,pq , whence we have p ⊥ • q ⊥ = 0 and thereby Lkpq ≡ 0. Similarly for the second integral we have q ⊥ = − 1 2 k ⊥ because of δ p,kq from which we have k ⊥ • q ⊥ = 0 and thereby Lpkq ≡ 0. This means that the integrand is identically zero on the singular plane on S and hence the integration along this singularity may be safely circumvented without blow up.Similarly, averting p Z = 0 demands k Z = 0 (always true) and k ⊥ = 1 2 p ⊥ .The latter condition is ensured because Lkpq , Lpkq ≡ 0 when k ⊥ = 1 2 p ⊥ based on very similar arguments as mentioned above.There are two other regions in Figure 3 that stand out for inspection, viz., (p x , p y ) = (0, 0) and k x = p x , k y = p y for the value of the integrand is high at these points due to the large values of Lkpq and Lpkq .These correspond to p ⊥ = 0 and q ⊥ = 0 respectively and for reasons discussed above are outside the region of validity of AT (R-RHD).Hence these points do not contribute to the integral.In general, the anisotropic limit 1 ensures that the region of integration is typically farther away from the singular plane defined by (47) where the value of the integrand is diminishingly small as can be seen from the contour plots in Figure 3.So we know at the very least that we are safe in terms of blow up of the integrand.However, we must still investigate whether or not the accumulation of these small values of the integrand increase in an unbounded fashion over a reasonably large radial wavenumber in the k ⊥ plane.In order to check this we resort to numerical evaluation of the collision integral (The singular points are carefully excluded for they are found to be inconsequential as discussed above) for a range of perpendicular wave numbers.The general form of the collision integral in the kinetic Equation ( 43) upon projection on the (p x , p y ) plane is The analytic expression for the Jacobian term + 1 is too messy to be written explicitly here but can be inspected symbolically using the computer program attached in the Appendix A.5.Most importantly, this term shares the location of the poles with that of the surface S = p Z (p x , p y ).Hence no additional caution is required for this term.Since we are interested in analyzing the convergence of this integral, we have ignored the constant multiplier on the r.h.s. of Equation (43).We numerically integrate the function (the integrand in expression ( 49)) over the surface S defined by Equation ( 46) by two different methods, viz., adaptive quadrature method and adaptive Monte Carlo method.The numerical integration is performed over S by excluding the singular points described in the previous paragraph as the singularities were either outside the region of validity of AT (R-RHD), i.e., outside the region where k Z /k ⊥ 1 is true or occurred at points where the integrand is identically zero.The computer program used to generate the plots in this section is given in Appendix A.5.By analyzing Figure 4, it is evident that by staying within the region of validity of AT (R-RHD), as has been elaborated above in this section, the collision integral is convergent.Moreover, it is clear from the results of the numerical integration that the convergence holds for reasonably large values of k ⊥ thereby permitting an energy cascade to develop.Figure 4 shows that the collision integral on the k ⊥ plane plateaus away from the boundary of integration (at the boundary, the numerical integration is less accurate) and is very nearly zero as is expected for a stationary solution.
In this section, we have provided direct evidence of the convergence of the collision integral I(k ⊥ , k Z ) by numerical computation to validate the applicability of the locality assumption for seeking the Kolmogorov constant energy flux solution.For a rigorous analysis of convergence of collision integrals in a three wave system, one may refer to the work by Connaughton [52].Additionally, it is important to note that analysis of stability and evolutionary non-locality of concerned spectra subject to extra perturbative forcing in this anisotropic limit is beyond the scope of the current paper and the interested reader is referred to the work of Balk and Nazarenko [36] and Nazarenko [7] for details.

Non-Zero Helicity Dynamics: Interplay of Energy and Helicity
In this section, we deduce a general set of coupled equations for the two invariants of the flow.This is done by formally extending the symmetrical system of the previous section.

Coupled Equations for Energy and Helicity
Recall that the assumption, However, on relaxing such an assumption, a coupled set of equations for energy and helicity may be arrived at by using Equation (28) in Equation (43) (cf.e k in Equation ( 43) is actually e + k ≡ e − k ).Thus the closed form coupled energy-helicity equation becomes, The individual evolution equation for e k (and h k ) follows by adding (and subtracting) the two set of equations expressed concisely by Equation ( 50) and is given as follows: and Equations ( 51) and ( 52) clearly reveal the coupled dynamics of energy and helicity in a rapidly rotating fluid flow.The possibility of recovering the total energy and helicity spectrum from the simpler zero helicity case relies on extending the functional on the right hand side of Equation ( 43) to the non-zero helicity case.
In essence, we have taken a special case of the kinetic equations which has the functional form as the right hand side of Equation ( 43)) where h k = 0, and extended it to the more general case where h k = 0.In this general case, the domain of f (k) is still the positive real line because the inequality |h k | ≤ k ⊥ e k implies that h k and e k are not independent variables and e k ± h k k ⊥ ≥ 0 is always true.Hence, we have extended the applicability of the kinetic equation from the simpler symmetric case (where c * k = c − k because h k = 0) to the more general case with non-trivial helicity.This way we have circumvented tedious algebraic computations involving multiple correlation functions (to account for the departure in mirror symmetry in the fully helical case) by first, reducing the Hamiltonian system by invoking reflection symmetry and then formally extending the symmetrical system to the more general helical case.The applicability and implication of this procedure is explained in more details in [38].

Generalized Solution of Energy and Helicity Spectra
Note that Equation ( 28) are consistent with the definitions, e k = e + k + e − k and h k = k ⊥ (e + k − e − k ).Clearly, the power law solutions do not change for the non-zero helicity case because the homogeneity of the interaction coefficient and the linear dispersion relation remains the same.Thus, ⊥ that is consistent with earlier findings based on numerical simulations [21].The cylindrically symmetric solutions are: ⊥ , and The power law solution obtained here pertains to the perpendicular cascade within the slow manifold, the reader is referred to Bellet et al. [53] for power law solutions of the spectrum in the axial direction.It is important to note that the solutions given by Equation ( 53) are different from the ones discussed in Pouquet and Mininni [24] as the latter assume isotropy in their arguments to show that the sum of the powers of the wave numbers for E k and H k is equal to 4. The solution set of Galtier [8] belongs to the regime where the effect of fast inertial waves is dominant.In the recent work of Galtier [9], it has been shown that within the theory developed in the earlier work of Galtier [8], the power law solutions can be generalized to the empirical form presented in [24].In contrast, the analysis presented here captures the dynamics in the slow manifold where the flow is highly anisotropic and the effect of the fast inertial waves is sub-dominant, and accounts for a distinct separation of temporal scales, t and τ as is necessary [16] and explained in the introductory section here.In such a regime in the slow manifold that is devoid of fast inertial wave modulation and where the dynamic is well captured by the R-RHD equations, the stationary power law solutions, that are entirely dependent on the functional form of the interaction term L and the dispersion relation ω k , are given by Equation (53) and is in agreement with results from numerical simulations (e.g., the work of Mininni and Pouquet [21], Teitelbaum and Mininni [22]).It may be interesting to point out here that an identical k −2 ⊥ spectrum for E ± ⊥ has been reported in the case of anisotropic MHD turbulence of Alfven waves in the presence of strong external magnetic field (where dominant modes are k z k ⊥ ) [7,34,54].The analytical results for the anisotropic MHD case were developed based on the reduced MHD model referenced in [7].Thus it may not be surprising to draw an analogy between the anisotropic MHD wave turbulence and the anisotropic hydrodynamic case discussed in the current paper where the strong magnetic field plays the role of rapid rotation, and the polarized Alfven waves are congruent to the helical inertial waves discussed here (please see the Table 1 in Appendix A.4).
The important point to note here is that rotating turbulence encompasses several different and distinct dynamical regimes, each with distinct set of solutions.Hence, the importance of using reduced equations for distinct asymptotic limits as has been explained by Nazarenko and Scheckochihin [27] (see Appendix A of [27]).Within such a distinct limit of rapid rotation, the R-RHD equations have been derived [14,15,26,27] and a multiple scales perturbation method has been applied in this paper for the asymptotically reduced equations.In the following sections, we attempt to stitch together the important results of weak and strong turbulence of different dynamical regimes in order to clarify the picture of turbulent cascade that has evolved based on recent research literature.

Hierarchy of Slow Manifolds in Anisotropic Wave Turbulence Diverges from the Critical Balance Route towards Isotropy
The discussion in this section is motivated by the turbulence cascade picture presented by Nazarenko and Scheckochihin [27].We modify the turbulent cascade schematic based on the results presented here.
The critical balance with polarization alignment argument presented by Nazarenko and Scheckochihin [27] is an attempt to explain the k −2 ⊥ energy spectra observed in several numerical simulations of rotating turbulence [17,21,22,55].However, as is evident from Figure 5 and the corresponding sketch in [27], the critical balance with polarization alignment leads to a departure from anisotropy and is a path to the recovery of isotropic scales that are prevalent above the Zeeman wavenumber [23].It must be emphasized that anisotropy is dominant in rotating turbulent flows.It is in this light, we believe that the energy spectra k −2 ⊥ is obtained for anisotropic wave turbulence in the regime of R-RHD defined by AT (ω f T τ ω s T τ 1, ω s T f 1 and k Z k ⊥ ).In this anisotropic regime, the slow inertial wave frequency is much smaller than the fast inertial wave frequency.This solution prevails as the flow traverses a hierarchy of slow manifold regions with successively decreasing k Z k ⊥ and is the anisotropic wave turbulence solution for the energy spectra within the slow manifold.In summary, there seems to be a bifurcation of the energy spectral solution at the critical balance wavenumber (see Figure 5), two distinct spectra evolve, each with a k −2 ⊥ energy spectra: one leading towards the isotropic scale via critical balance and hence falls within the realm of strong turbulence (Strong wave turbulence is fairly complicated (see secs.3.2 and 14.6 in [7]) because the wave amplitudes can be O(1) or higher [56] and the perturbative expansion may not be readily possible.Its description is incomplete as stated by Nazarenko [7]), and the other towards highly anisotropic horizontal scales with k Z k ⊥ sustained by small amplitude resonating wave modes.Here, k i is the isotropic wavenumber, k ⊥c is the classical critical balance wavenumber, k 0 is the injection wavenumber corresponding to an initial wave field.Three distinct regimes are shown: (i) WT (Galtier) corresponding to the wave-turbulence regime with ω f T τ 1, ω f T s 1, (ii) CB w/ pol.(i.e., critical balance with polarization alignment) as explained in [27] leading towards isotropy, and (iii) AT (R-RHD) corresponding to the anisotropic wave turbulence dynamics of the R-RHD equations with ω f T τ ω s T τ 1, ω s T f 1 and k Z k ⊥ .As we move along the horizontal axis from left to right, the flow traverses a hierarchy of slow manifolds with successively rescaled (decreasing) k Z /k ⊥ wave number ratio.Also shown are possible explanation of 2D-3D coupling by non-resonant interactions.The AT (R-RHD) does not explain inverse cascade phenomena as the energy flux obtained is direct.Here k 0 is the wavenumber corresponding to the initial condition.

Comparison with Weak Turbulence Theory of Galtier
In this section, we contrast the theory presented here with earlier works [8,9].The main distinctive features are listed below.
1.The governing equations on which the wave turbulence theory [8,9] is developed are the Navier-Stokes equations (i.e., Equation ( 1)).In these equations, the Rossby number Ro appears explicitly and is used as an order parameter in the perturbation analysis and the fast inertial wave dynamic is dominant.In contrast, the governing equations for the analysis presented here are the R-RHD (i.e., Equations ( 5) and ( 6)).The asymptotic limit of infinitesimally small Ro is already accounted for in the multiple scales analysis to derive the R-RHD and hence do not explicitly appear in the R-RHD equations.This means that the effect of fast inertial waves is sub-dominant here and the R-RHD equations are hence suitably applicable for the slow manifold dynamics.2. The dynamical regime where the theory [8,9] is valid is ω f T τ 1, ω f T s 1.This point has been elaborated in great detail in the work of Nazarenko and Scheckochihin [27] (see Section 2 in [27]).The dynamical regime of the R-RHD is ω f T τ ω s T τ 1, ω s T f 1 and k z k ⊥ .It is the limit in which k z is so small that the turbulent dynamic is populated by slow inertial waves with dispersion relation 1 (also see Figure 5 above).
In other words, within the slow manifold, k z is so small that k Z k ⊥ Ro 1.The perturbative analysis presented here is applied to the R-RHD where the smallness (weakness) of the wave amplitude is measured by 1 and the slow dispersive three wave system undergoes weak non-linear exchanges at the asymptotic order 2 as explained earlier in Sections 4 and 5. 3.In wave turbulence theory of Galtier [8] and Galtier [9], the kinetic equations for energy and helicity are derived by using multiple correlation functions to capture the energetics as well as the absence of symmetry due to helicity.This makes the calculations tediously lengthy.In contrast, the derivation for the helicity kinetic equation is presented here as a natural extension of the symmetrical non-helical system and bypasses the use of calculations using multiple correlation functions.This simpler approach follows a more general philosophy of Hamiltonian reduction exploiting symmetries in the system [25] and their natural extension to understanding asymmetrical phenomena.4. Finally, the spectral power laws obtained for the energy cascade are distinctly different as explained concisely in Figure 5.
Despite the fundamental differences in the region of validity of the two theories, they present a more detailed recipe to better understand turbulent energetics and cascades for rotating flows.A simple schematic towards this goal is presented in Figure 5 above to highlight the key findings in this field.

Comparison with Weak Turbulence Theory of Newell: Cumulant Hierarchy vs. Wave Amplitude Hierarchy
The fundamental difference between the multiple scales perturbation technique employed in this paper to derive the kinetic equations and the weak wave turbulence theory reviewed by Newell et.al [39] is explained in the following paragraphs.
Weak turbulence calculations, as reviewed by Newell et.al. [39], are computed by taking the closure in the limit T → ∞.The third and higher order cumulants survive at the advective timescale (τ NL ) in the presence of fast inertial waves.This implies prevalence of non-Gaussian statistics.This is the reason for writing a hierarchal system for the cumulants in the perturbation approach, cf.Newell et.al. [39], rather than for the Fourier amplitudes.
The dynamics explained here elapse at the much slower advective time scale (τ NL ) compared to the faster wave timescale T → ∞.At this slower time scale, it is assumed that the third (and higher odd) order cumulants are sub-dominant thereby corroborating the suitability of Gaussian statistics for the wave field.This aspect is detailed by Nazarenko [7] (see p. 64 in Section 5.6).This justifies the applicability of Wick's theorem in the derivation of the kinetic equations.Consequently, this is also why we employed a hierarchal expansion for the wave amplitude in Equation (18).

Coupling between Wave and 2D Modes Through Non-Resonant Interactions
Neither the theory presented in [8] nor the model presented in this manuscript address the issue about the possibility of wave modes pumping the 2D modes.Note that Lkpq = 0 when k Z = 0 because of the fact that p ⊥ = q ⊥ .Thus, within the framework of purely resonating wave triads, it is not possible to establish the coupling between wave and purely 2D modes.This should not be surprising because the theory developed is that of dispersive waves that are in resonance.However, a rather phenomenological explanation may be insightful for explaining the coupling between wave and 2D modes, as is presented below.
Suppose that in Equation ( 22), the triadic resonance condition is modified such that φ(ω) = s k ω k − s p ω p − s q ω q = δ ω , δ ω = 0.This condition represents non-resonant (or near-resonant) interaction of the three wave modes.For small δ ω τ, Taylor expansion implies, i.e., Thus computing the fast time integral after Taylor-expanding the exponential term reveals that the higher order slow secular terms (cf.presence of τ in the higher order terms) are inherently embedded in the full system that is not restricted to resonating triadic wave interactions only.These terms account for the coupling with the 2D modes.This is evident by re-writing Equation ( 22) with the higher order terms as follows, The slow secular behavior of the non-resonant term (after integration) can be efficiently moderated by tuning the damping parameter δ ω .The kinetic equation, that is constructed from the above amplitude equation as explained before, can then be decomposed into two parts with contributions from resonating and non-resonating modes considered separately, as follows: where T(k, τ) denotes nonlinear transfer of energy to mode k.Thus, retaining only non-resonant interactions, it is possible to establish the aforementioned coupling in the slow manifold, k Z = 0. Note that in the case of non-resonant interactions, due to the absence of the frequency delta function, a stationary Kolmogorov (constant flux) solution cannot be obtained.The reader is also referred to the works of Janssen [57] and Annenkov and Shrira [58] for a detailed theory of quasi-resonant interactions in four wave dispersive systems.

Conclusions
In conclusion, it is important to emphasize an important point, that of the asymptotic dynamical regime to which the fluid system belongs.In the context of this paper, we have restricted our analysis to the highly anisotropic regime of rapid rotation (i.e., infinitesimally small Rossby number) within the slow manifold where k z is infinitesimally small.It has been shown that the application of a multiple scales perturbation method within this regime yields a k −3 ⊥ law for the anisotropic energy spectrum, e k (cf.equivalently a k −2 ⊥ spectrum for the cylindrically symmetric spectrum, E k ⊥ ) that is in agreement with results from experimental and numerical simulations [17,[21][22][23] as has been stated earlier.An asymptotically reduced system spans a hierarchy of slow manifold regimes and thereby captures the gradual transfer of energy towards the 2D modes.Interestingly, a similar power law solution can also be obtained by applying a critical balance phenomenology (where fast inertial wave time scale balances the nonlinear advection time scale) to the system of rotating turbulence as has been shown in Nazarenko and Scheckochihin [27].This is the realm of strong turbulence where the nonlinear interactions are strong, meaning ω f T τ NL ∼ O(1).However, the anisotropic limit of rapidly rotating turbulence is farther away from modes where critical balance holds, this has been shown through numerical simulations in Di Leoni et al. [59].In addition to the discussion in Section 6.3 above, the reader is referred to [27,38] for a detailed discussion on a wave turbulence and critical balance schematic of the energy cascade process.It is important to note that in the analysis presented in this paper, any physical artifact induced by boundary condition is not considered.Interested readers are referred to the work of [60] that describes discrete boundary effects on wave turbulence formalism.
In summary, we have constructed a multiple scales kinetic model for an asymptotically reduced set of equations that is valid in the limit of rapid rotation and in the anisotropic slow manifold regime.Stationary solutions of invariant quantities have been obtained that are consistent with experimental and simulation data reported in recent work.A coupled set of equations has been derived explaining the nature of the inter-dynamics of the two global invariants of the system, viz., energy and helicity.This has been done by extending the symmetrical non-helical system to the more general helical case where the reflection symmetry is broken.This procedure is novel in the sense that it bypasses construction of multiple correlation functions to account for the departure in mirror symmetry in the helical case.This analytical study will serve a useful reference point for theoretical understanding of atmospheric phenomena of planets that require a better knowledge of anisotropic wave dynamics.Moreover, detailed understanding of anisotropic multi-scale wave-eddy interactions, along the lines of the work presented here, will be useful for engineering new parametrization of mesoscale eddy flux in closure schemes for ocean circulation models [61].
Appendix A.2. Derivation of the Amplitude Equation: Supplementary Calculations Recall from Section 2.1 that the wave field is proportional to terms that evolve at advective time scale τ and the exponential term that elapses at the inertial time scale t (cf.since the R-RHD limit is ωτ = τ t 1, we chose τ = t, 1, also t fast inertial time scale that has been filtered out by the R-RHD), Ψ advective scale τ e iΦ(k,s k ω k t) This separation of scales between τ and t (note τ and t are now independent variables) allows us to average (integrate) out terms varying at the inertial timescale t ∼ O( 1 ), unless the terms are multiplied by terms that are functions of ω ∼ O( ), and leaves us with terms that are dependent on the advective scale τ alone.This enables us to write an evolution equation for the small amplitude of the form ∂ τ c s k k (τ) = r.h.s of Equation (22).This is standard procedure in multi-scale perturbation techniques and the interested reader is referred to texts by Hinch [31], Bender and Orzag [32], Kevorkian and Cole [33].To elucidate that the averaging is done over a large time limit, consider, e.g., = 1 1000 , consequently τ = t = 1 1000 t; this means that for τ to elapse 1 unit, t must elapse 1000 units, i.e., within the scope of τ, t is already very large.The average of a function, say f (t) defined over the domain D is given as We apply the solvability condition 1 2 (h 2k ) = 0 to Equation ( 21) to obtain the amplitude equation.The solvability condition involves an inner product (denoted by the operator • above) that includes a projection onto the helical basis (that we denote by angle brackets here, i.e., 1  2 h −s k k , • in wavenumber space) and a time averaging (that we denote by normal integration symbol) as shown below We operate each term of Equation ( 21) with the operator (64) defined above.
1st term: 1 2 (h 2k dt = 1 T T 0dt = 0.The 0 integrand follows from application of the solvability criterion which is basically a Fredholm alternative in the context of the adjoint problem explained in Section 4. This makes the left hand side of Equation ( 21) null upon application of the inner product.k (τ) can be factored out of the integration over t because t and τ are independent variables as has been explained above.Also note that the integrand is equal to one because we have used the fact that 3rd term: For the third term we interchange the order of integration, i.e., we swap the operations 1 T T (•)dt and 1  2 h s k k , • .Note that all terms except the exponential term are functions of τ (and not t) and hence can be factored out of the integration over t.So we now concentrate only on the averaging of the exponential term, comprising of the ω terms, as follows.Recall, for sake of brevity, we use the following notation φ(ω) := (s k ω k − s p ω p − s q ω q ).Now, in the limit of large t ∼ O( 1 ) as → 0, φ(ω) ∼ O( ) becomes proportionally small thereby φ(ω)t = φ 1 (ω) τ = φ 1 (ω)τ ∼ O(1) in the leading order approximation.Here φ 1 (ω) ∼ O (1).This means we have e p e q δ k,pq δ ω k ,ω p ω q dk = 0. To show this, we approximate the delta function with a limiting exponential function as: δ ω k ,ω p ω q ≈ lim σ→0 e − ω k −ωp−ωq e p e q δ k,pq δ ω k ,ω p ω q dk The above equation is zero on account of the integral e k e p δ k,pq δ ω k ,ω p ω q dk = 0 and 1 k 2

Figure 1 .
Figure 1.Anisotropic energy spectra of rotating helical turbulence simulations published in [21,22].(a) e k ∼ k −3 ⊥ and (b)E k ⊥ ∼ k ⊥ e k ∼ k −2⊥ comport with the spectra obtained analytically in the current article as is explained in later sections.The inertial range in the above plots span roughly over an order of magnitude in the log scale, i.e., roughly 70 to 90 wave numbers.The line showing the −5/3 slope is just for reference as is clearly mentioned in[21].

Figure 2 .
Figure 2. Helical wave basis: ( k ⊥ , , k ) forms a right-handed coordinate system with k , k ⊥ = k ,  = 0 where  = k ×k ⊥ k 2 The total derivative is given by d dt = d dt + d dτ dτ dt and since τ = t, we have dτ dt =

2 1 ,
we have the average of the exponential term approach zero because 1 2T τ T τ (e p e k−p − e k e p − e k e k−p ) − 2| Lpk(p−k) | 2 (e k e p−k − e p e k − e p e p−k ) dS = p Z =0,p ⊥ 0 | Lkp(k−p) | 2 (e p e k−p − e k e p − e k e k−p ) −2| Lpk(p−k) | 2 (e k e p−k − e p e k − e p e p−k )

Figure 3 .
Figure 3. Contour plot showing the surface of integration S : p Z (p x , p y ) of the collision integral (without the constant π) in Equation (43) projected on the (p x , p y ) plane.Colors represent the value of the integrand.Figures on the left are in logarithmic scale.(a) and (b) correspond to α = 1, k = 20 while (c) and (d) correspond to α = 0.1, k = 100 with k Z = 1 fixed such that k Z /k ⊥ 1. (e) is top view of (a).The singular regions are either outside the region of validity of AT (R-RHD) or correspond to identically null value of the integrand as explained in this section.

Figure 4 .
Figure 4. Numerical integration of the collision integral in Equation (43) normalized by the constant π.The plots show I(k ⊥ , k Z ) vs. k ⊥ = (k x , k y ) for fixed k Z = 1.The integration is computed at the nodes k x , k y ∈ [10, 150] in intervals of 20.(a) Method: Gauss Kronrod quadrature; (b) method: adaptive Monte Carlo.

Figure 5 .
Figure5.A sketch of cascade paths for rotating turbulence shows the different flow regimes and the corresponding energy spectra.Here, k i is the isotropic wavenumber, k ⊥c is the classical critical balance wavenumber, k 0 is the injection wavenumber corresponding to an initial wave field.Three distinct regimes are shown: (i) WT(Galtier)  corresponding to the wave-turbulence regime with ω f T τ 1, ω f T s 1, (ii) CB w/ pol.(i.e., critical balance with polarization alignment) as explained in[27] leading towards isotropy, and (iii) AT (R-RHD) corresponding to the anisotropic wave turbulence dynamics of the R-RHD equations with ω f T τ ω s T τ 1, ω s T f 1 and k Z k ⊥ .As we move along the horizontal axis from left to right, the flow traverses a hierarchy of slow manifolds with successively rescaled (decreasing) k Z /k ⊥ wave number ratio.Also shown are possible explanation of 2D-3D coupling by non-resonant interactions.The AT (R-RHD) does not explain inverse cascade phenomena as the energy flux obtained is direct.Here k 0 is the wavenumber corresponding to the initial condition.

k dt = ∂ τ c s k k 1 T
T 1dt = ∂ τ c s k k (τ).Note that the term ∂ τ c s k

k
= 1 as has been explained earlier in Section 4.

1 First we show that 1 k 2 ⊥
(ω)τ dt = e iφ 1 (ω)τ 1 T T dt = e iφ 1 (ω)τ = e iφ(ω)t .(65)Thismay look misleading but note that φ(ω)t remains O(1) and relatively constant even for large t and hence the exponential term can be factored out of the integral and 1T T dt = 1.All other variables besides the exponential term that make up the third term of Equation (21) are functions of τ and upon being operated by 1 2 h s k k , • result in the right hand side of Equation (22).Appendix A.3.Energy Conservation in Wave Triads