2.1. Quasi-two-dimensionalisation and inertial waves

In this section we discuss the theoretical results underpinning the present study. A fundamental property of rotating flows is the fact that they support inertial wave motions, whose restoring force is the Coriolis force (Greenspan Reference Greenspan1968). Inertial waves have the peculiar anisotropic dispersion relation

(2.1)\begin{equation} \omega^{s_{\boldsymbol {k}}}({\boldsymbol {k}}) = 2s_{\boldsymbol {k}} \Omega k_{\parallel}/k, \end{equation}

where ${s_{\boldsymbol {k}}}={\pm }1$, $\Omega$ is the rotation rate, $k_{\parallel }$ is the component of ${\boldsymbol {k}}$ along the rotation axis and $k=|{\boldsymbol {k}}|$. Similarly, we define $k_{\perp }$ as the magnitude of the component of ${\boldsymbol {k}}$ perpendicular to the rotation axis. In the remainder of this article, parallel and perpendicular will always refer to the rotation axis. Inertial waves in fast-rotating turbulence are important for understanding the direction of the energy cascade, as will be discussed below. The form of (2.1) shows that motions which are invariant along the axis of rotation, i.e. which are 2-D with three components (2D3C), have zero frequency and are thus unaffected by rotation. This allows decomposition of the flow into two components, the 2D3C modes which are not directly affected by rotation, forming the slow manifold, and the remaining 3-D modes which are affected by the rotation, forming the fast manifold (Buzzicotti et al. Reference Buzzicotti, Clark Di Leoni and Biferale2018b). In the limit ${Ro} \to 0$, it can be shown that only resonant interactions remain present (Waleffe Reference Waleffe1993; Embid & Majda Reference Embid and Majda1998; Chen et al. Reference Chen, Chen, Eyink and Holm2005). Resonant interactions are those interactions between wavenumber triads $({\boldsymbol {k}},{\boldsymbol {p}},{\boldsymbol {q}})$ satisfying

(2.2)\begin{gather} {\boldsymbol {k}} + {\boldsymbol {p}} +{\boldsymbol {q}} = 0, \end{gather}

(2.3)\begin{gather} \omega^{s_{\boldsymbol {k}}}( {\boldsymbol {k}}) +\omega^{s_{\boldsymbol {p}}}( {\boldsymbol {p}}) + \omega^{s_{\boldsymbol {q}}}({\boldsymbol {q}}) = 0, \end{gather}

where $\omega ^{s_{\boldsymbol {k}}}({\boldsymbol {k}}), \omega ^{s_{\boldsymbol {p}}}({\boldsymbol {p}})$ and $\omega ^{s_{\boldsymbol {q}}}({\boldsymbol {q}})$ are given by (2.1). When only resonant interactions are present in the system, it can further be shown that any triad including modes from both the fast and slow manifolds leads to zero net energy exchange between the two manifolds. Thus, with only resonant interactions, the slow and fast manifolds evolve independently from each other without exchanging energy, and there is inverse energy transfer in the perpendicular components of the slow manifold. This decoupling may lead to an inverse energy cascade for the quasi-2-D part of the flow. In fact, it can be proven that for finite Reynolds number $Re\equiv U\ell /\nu$ (where $\nu$ is viscosity, $U$ is root-mean-square velocity and $\ell$ is a forcing length scale) and finite $H$, the flow will become exactly 2-D as $Ro\to 0$ (Gallet Reference Gallet2015).

On the other hand, in the limit of large domain height $H$, very small values of $k_{\parallel }$ are possible, such that quasi-resonant triads, for which (2.3) is only satisfied to $O({Ro} )$, can transfer energy between the slow and fast manifolds. Thus the inverse energy transfer in the slow manifold may be suppressed by interaction with quasi-resonant 3-D modes. Asymptotically, for infinite domains and $k_{\parallel } / k_{\perp } \ll 1$, wave turbulence theory predicts a forward energy cascade and an associated anisotropic energy spectrum (Galtier Reference Galtier2003).

There are thus two mechanisms at play in the energy transport: the dynamics of the slow manifold transferring energy to the large scales and the 3-D interactions transferring energy to the small scales. Which of these two processes dominates depends on the two non-dimensional parameters, the ratio $h=H/\ell _{in}$, where $\ell _{in}$ is the forcing length scale, and the Rossby number $Ro=\epsilon _{in}^{1/3} / (\ell _{in}^{2/3} \Omega )$ based on $\ell _{in}$ and the forcing velocity scale $(\epsilon _{in}\ell _{in})^{1/3}$ resulting from the energy injection rate $\epsilon _{in}$. The main criterion is whether or not 2-D modes are isolated from 3-D modes due to fast rotation. The coupling of 2-D and 3-D motions will be strong enough to stop the inverse cascade if the fast modes closest to the slow manifold ($k_{\parallel } \sim H^{-1}$, $k_{\perp } \sim \ell _{in}^{-1}$) are ‘slow’ enough to interact with the 2-D slow manifold. This implies that their wave frequency $\omega = 2\Omega k_{\parallel }/k_{\perp } \sim 2\Omega \ell _{in}/H$ is of the same order as the nonlinear inverse time scale $\tau _{nl}^{-1} \sim \epsilon _{in}^{1/3} \ell _{in}^{-2/3}$. This leads to the following prediction for the critical height $H_c$, where the transition takes place:

(2.4)\begin{equation} h_c=\frac{H_c}{\ell_{in}} \propto \Omega \epsilon_{in}^{-1/3} \ell_{in}^{2/3} = {Ro}^{-1}. \end{equation}

Importantly, the predicted critical rotation rate and height are linearly proportional. The criterion (2.4) suggests that the non-dimensional control parameter of the transition in the limit of large $h$ and small ${Ro}$ is given by

(2.5)\begin{equation} \lambda = \frac{1}{h \times {Ro}} = \frac{\ell_{in}^{5/3} \Omega }{\epsilon_{in}^{1/3} H}. \end{equation}

2.2. Multiscale expansion

In the present paper, we will explore the regime of simultaneously large $h$ and small ${Ro}$. Brute-force simulations at small ${Ro}$ are costly, since very small time steps are required to resolve the fast waves of interest. Rather, we exploit an asymptotic expansion first introduced in Julien, Knobloch & Werne (Reference Julien, Knobloch and Werne1998), which allows one to test the prediction (2.4) and to investigate the properties of the transition to a split cascade. The expansion is based on the constant-density Navier–Stokes equation in a system rotating at a constant rate $\boldsymbol {\Omega } = \Omega \hat {e}_{\parallel }$, which in its dimensional form reads

(2.6)\begin{gather} {\partial}_t {\boldsymbol {u}} + {\boldsymbol {u}}\boldsymbol{\cdot} \boldsymbol\nabla {\boldsymbol {u}} + 2\Omega \hat{e}_{\parallel}\times{\boldsymbol {u}}= - \boldsymbol\nabla p + \nu \nabla^2{\boldsymbol {u}}+ {\boldsymbol {f}}, \end{gather}

(2.7)\begin{gather} \boldsymbol\nabla \boldsymbol{\cdot} {\boldsymbol {u}} =0 , \end{gather}

with time $t$, velocity ${\boldsymbol {u}}$, pressure (divided by constant density) $p$, kinematic viscosity $\nu$ and forcing ${\boldsymbol {f}}$. The domain considered here is the cuboid of dimensions $2\pi L\times 2 \pi L \times 2 \pi H$, depicted in figure 1, with periodic boundary conditions. For any vector ${\boldsymbol {f}}$, we define the parallel and perpendicular components as ${\boldsymbol {f}}_{\parallel } = ({\boldsymbol {f}}\boldsymbol{\cdot} \hat {e}_{\parallel }) \hat {e}_{\parallel } = F_{\parallel } \hat {e}_{\parallel }$ and ${\boldsymbol {f}}_{\perp } = {\boldsymbol {f}}-{\boldsymbol {f}}_{\parallel }$.

Figure 1.

The long, rapidly rotating box domain.

We choose to consider a stochastic forcing injecting energy at a constant mean rate into both perpendicular and parallel motions $\langle\,{\boldsymbol {f}}_{\perp } \boldsymbol{\cdot} {\boldsymbol {u}}_{\perp } \rangle =\langle\, f_{\parallel } u_{\parallel } \rangle = \epsilon _{in}/2 \Rightarrow \langle\, {\boldsymbol {f}}\boldsymbol{\cdot} {\boldsymbol {u}}\rangle = \epsilon _{in}$, where $\langle \boldsymbol{\cdot} \rangle$ denotes an ensemble average over infinitely many realisations. The forcing is chosen to be 2-D (independent of the parallel direction), for simplicity, and filtered in Fourier space to act only on a ring of perpendicular wavenumbers ${\boldsymbol {k}}$ centred on $|{\boldsymbol {k}}| = k_f=1/\ell _{in}$. A similar 2-D forcing has been widely used in previous studies of the transition toward an inverse cascade (e.g. Smith et al. Reference Smith, Chasnov and Waleffe1996; Celani et al. Reference Celani, Musacchio and Vincenzi2010; Deusebio et al. Reference Deusebio, Boffetta, Lindborg and Musacchio2014). In the present case, it makes sense to restrict the forcing to the 2-D modes because for any forcing with finite correlation time, the injection of energy to the $k_{\|}\ne 0$ modes would be suppressed when the limit $\Omega \to \infty$ is taken due to the high wave frequencies. Thus, most of the energy would be injected into the $k_{\|}=0$ modes. We note, however, that in general the transition to an inverse cascade can depend on the choice of forcing. Recent work in thin-layer turbulence by Poujol, van Kan & Alexakis (Reference Poujol, van Kan and Alexakis2020) suggests that a 3-D forcing, which includes non-zero parallel wavenumbers, is less efficient at generating an inverse cascade and delays the onset. A related problem of interest, which we do not address in the present study, is to investigate the transfer of energy to the 2-D manifold in the case when only the 3-D modes are forced. This has recently been studied experimentally and theoretically (Le Reun, Favier & Le Bars Reference Le Reun, Favier and Le Bars2019; Brunet, Gallet & Cortet Reference Brunet, Gallet and Cortet2020, and Le Reun et al. Reference Le Reun, Gallet, Favier and Le Bars2020).

The forcing imposes a time scale $(\ell _{in}^2/\epsilon _{in})^{1/3}$ as well as a length scale $\ell _{in}$, and thus a velocity scale $(\epsilon _{in} \ell _{in})^{1/3}$. However, the typical scale of parallel variations is $H$, rather than $\ell _{in}$. As detailed in appendix A, non-dimensionalisation using these scales reveals that the non-dimensional control parameters are indeed the Rossby number $Ro = (\epsilon _{in} \ell _{in})^{1/3}/(\Omega \ell _{in})$ and the rescaled domain height $h=H/\ell _{in}$, in addition to the Reynolds number ${Re} = (\epsilon _{in}\ell _{in}^4)^{1/3}/\nu$ and the rescaled domain width $\Lambda = L/\ell _{in}$. We consider tall boxes $h = \epsilon ^{-1}\gg 1$, under the influence of fast rotation, $Ro=O(\epsilon )\ll 1$, such that $\lambda =(h Ro)^{-1}= O(1)$ (independent of $\epsilon$), while ${Re} =O(1)$ and $\Lambda =O(1)$. In particular, the expansion assumes that ${Re} \ll Ro^{-1}$, ${Re} \ll H/\ell _{in}$, $L/\ell _{in} \ll Ro^{-1}$ and $L \ll H$. We note that this limiting procedure does not correspond to the weak turbulence limit, for which the limit $h\to \infty$ is taken before the limit ${Ro} \to 0$. The method of multiple scales or a heuristic derivation (see appendix A) can be used to obtain a set of asymptotically reduced equations for the parallel components of velocity $u_{\parallel }$ and vorticity $\omega _{\parallel } = (\boldsymbol \nabla \times {\boldsymbol {u}})\boldsymbol{\cdot} \hat {e}_{\parallel }$ whose dimensionless form reads

(2.8)\begin{gather} \partial_t u_{\parallel} + {\boldsymbol {u}}_{\perp}\boldsymbol{\cdot}\boldsymbol\nabla_{\perp} u_{\parallel} + 2\lambda \partial_{\parallel}\nabla^{-2}_{\perp}\omega_{\parallel} = \frac{1}{{Re}} \nabla^2_{\perp} u_{\parallel}+ f_{\parallel}, \end{gather}

(2.9)\begin{gather} \partial_t \omega_{\parallel} +{\boldsymbol {u}}_{\perp} \boldsymbol{\cdot} \boldsymbol\nabla_{\perp} \omega_{\parallel} - 2\lambda \partial_{\parallel} u_{\parallel} = \frac{1}{{Re}} \nabla_{\perp}^2\omega_{\parallel}+ f_{\omega}, \end{gather}

where $\partial _{\parallel }$ is the partial derivative in the parallel direction, $\boldsymbol\nabla_{\perp } = \boldsymbol \nabla -\hat {e}_{\parallel } \partial _{\parallel }$ and $f_{\omega } = (\boldsymbol \nabla \times {\boldsymbol {f}})\boldsymbol{\cdot} \hat {e}_{\parallel }$. The perpendicular components ${\boldsymbol {u}}_{\perp }$ are divergence-free to leading order, ${{\boldsymbol {\nabla} }_{\perp }{\boldsymbol{\cdot} u}_{\perp} =0}$, which permits us to write them in terms of a stream function, ${\boldsymbol {u}}_{\perp } = \hat {e}_{\parallel } \times \boldsymbol \nabla \psi$, where $\psi$ is such that $\omega _{\parallel } = \nabla_{\perp }^2 \psi$. These non-dimensional asymptotic equations are valid in the domain $2\pi \Lambda \times 2 \pi \Lambda \times 2 \pi$. Importantly, in (2.8) and (2.9), all the information about $H$ and $\Omega$ is contained in the parameter $\lambda$, which is defined by (2.5). This implies that if a transition from a direct to a split energy cascade is captured in these asymptotic equations, the single control parameter of the transition indeed is given by $\lambda$ (in the limit of large $Re$ and $\Lambda$), as predicted by the heuristic arguments in § 2.1. Variants of the asymptotic equations (2.8) and (2.9) have been extensively used in the past, in particular for studying rotating turbulence (Nazarenko & Schekochihin Reference Nazarenko and Schekochihin2011) and rapidly rotating convection (adding energy equation) (Sprague et al. Reference Sprague, Julien, Knobloch and Werne2006; Grooms etal. Reference Grooms, Julien, Weiss and Knobloch2010; Julien et al. Reference Julien, Knobloch, Rubio and Vasil2012a,Reference Julien, Rubio, Grooms and Knoblochb; Rubio et al. Reference Rubio, Julien, Knobloch and Weiss2014), as well as dynamos driven by rapidly rotating convection (adding the energy and magneto-hydrodynamic induction equations) (Calkins et al. Reference Calkins, Julien, Tobias and Aurnou2015).

Equations (2.8) and (2.9) are closely related to well-known models in geophysical fluid dynamics. In particular, since the leading-order perpendicular velocity is in geostrophic balance and advection is purely perpendicular, the model bears a resemblance to the classical quasi-geostrophic approximation valid in thin layers (see e.g. Pedlosky Reference Pedlosky2013). Indeed, (2.8) and (2.9) have been referred to as generalised quasi-geostrophic equations (Julien et al. Reference Julien, Knobloch, Milliff and Werne2006). A great advantage of the reduced equations over the full Navier–Stokes equations is that they can be efficiently integrated numerically, as explained below. Note that in the expansion, as a consequence of the Taylor–Proudman constraint applied in the limit ${Ro} \to 0$, $h\to \infty$, fast variations in the parallel direction are eliminated, i.e. $k_{\parallel } \ll k_{\perp }$ in terms of dimensional wavenumbers. Equations (2.8) and (2.9) retain inertial waves with the dispersion relation, in non-dimensional form:

(2.10)\begin{equation} \omega({\boldsymbol {k}}) = {\pm}2 \lambda \frac{k_{\parallel}}{k_{\perp}}. \end{equation}

The rescaled wavenumbers $k_{\perp } \geq 1/\Lambda$ and $k_{\parallel }\geq 0$ and $\lambda$ are all $O(1)$ (independent of $\epsilon$), and therefore only the inertial waves with order-one frequencies, i.e. those on the parallel scale of the layer depth, are retained in the reduced equations. Fast inertial waves, i.e. those whose parallel scale is comparable to the perpendicular scale and whose frequencies are comparable to the rotation rate $\Omega$, are filtered out. For this reason, the asymptotic reduction gives a significant improvement in efficiency (the filtering of fast inertial waves here is similar to the filtering of fast inertia–gravity waves in the classical quasi-geostrophic approximation). We perform direct numerical simulations of (2.8) and (2.9) to show that, as predicted by the theory outlined above, there is indeed a transition from a direct to an inverse energy cascade in this extreme parameter regime.

2.3. A homochiral triadic instability

In this subsection, we discuss a linear instability mechanism present in the asymptotically reduced governing equations (2.8) and (2.9), which will be helpful for the interpretation of the direct numerical simulation results. For concreteness, we work in the canonical basis and choose $\hat {e}_{\parallel }=\hat {e}_z$. The Fourier-transformed governing equations then read, in the absence of forcing or dissipation,

(2.11)\begin{gather} \partial_t\hat{u}_{\parallel}({\boldsymbol {k}}) - 2 \textrm{i} \lambda\frac{k_{\parallel}}{k_{\perp}^2} \hat{\omega}({\boldsymbol {k}}) =- \sum_{{\boldsymbol {p}}+{\boldsymbol {q}}+{\boldsymbol {k}}=0} [p_xq_y-q_xp_y] \frac{\hat{\omega}_{\parallel}^*({\boldsymbol {p}})}{p_{\perp}^2} \hat{u}_{\parallel}^*({\boldsymbol {q}}), \end{gather}

(2.12)\begin{gather} \partial_t\hat{\omega}_k - 2\textrm{i} \lambda k_{\parallel} \hat{u}_k = -\sum_{{\boldsymbol {p}}+{\boldsymbol {q}}+{\boldsymbol {k}}=0} [p_xq_y-q_xp_y]\frac{\hat{\omega}_{\parallel}^*({\boldsymbol {p}})}{p_{\perp}^2} \hat{\omega}_{\parallel}^*({\boldsymbol {k}}) , \end{gather}

where it was used that $u_{\parallel } \equiv \sum _{\boldsymbol {k}}\hat {u}_{\parallel }({\boldsymbol {k}}) \textrm {e}^{\textrm {i} {\boldsymbol {k}}\boldsymbol{\cdot} {\boldsymbol {x}}}$, similarly $\omega _{\parallel } \equiv \sum _{\boldsymbol {k}}\hat {u}_{\parallel }({\boldsymbol {k}}) \textrm {e}^{\textrm {i} {\boldsymbol {k}}\boldsymbol{\cdot} {\boldsymbol {x}}}$, in addition to $\omega _{\parallel } = \nabla_{\perp }^2 \psi$ and ${\boldsymbol {u}}_{\perp } = \hat {e}_{\parallel } \times \boldsymbol \nabla \psi$. It is helpful to reformulate the dynamics in a helical basis as in Waleffe (Reference Waleffe1993). The linearisation of (2.11) and (2.12) may be written as

(2.13a,b)\begin{equation} \partial_t {\boldsymbol {x}}({\boldsymbol {k}}) = L({\boldsymbol {k}}) {\boldsymbol {x}}({\boldsymbol {k}}),\quad L({\boldsymbol {k}})= \textrm{i} \omega({\boldsymbol {k}}) \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \end{equation}

with ${\boldsymbol {x}}({\boldsymbol {k}})=(\hat {u}_{\parallel }({\boldsymbol {k}}), \omega _{\parallel }({\boldsymbol {k}})/k_{\perp })$ and $\omega ({\boldsymbol {k}}) = 2 \lambda k_{\parallel }/k_{\perp }$, identical to (2.10) up to the sign. The eigenvectors of $L({\boldsymbol {k}})$ corresponding to the eigenvalues $\pm \omega ({\boldsymbol {k}})$ are given by

\[ \left\{ \begin{pmatrix} 1 \\ 1\end{pmatrix}, \begin{pmatrix} 1 \\ -1\end{pmatrix} \right\}. \]

They correspond to inertial waves of positive and negative helicity, respectively, with dispersion given by (2.10). Representing ${\boldsymbol {x}}({\boldsymbol {k}})$ in this eigenbasis leads to the new variables $Z^{s_{\boldsymbol {k}}}_{\boldsymbol {k}} = \hat {u}_{\boldsymbol {k}} + s_{\boldsymbol {k}} ({\hat {\omega }_{\boldsymbol {k}}}/{k_{\perp }})$, where ${s_{\boldsymbol {k}}}={\pm }1$. The full nonlinear system may then be written entirely in terms of $Z_{\boldsymbol {k}}^{s_{\boldsymbol {k}}}$.

We consider a (not necessarily resonant) triad (${\boldsymbol {k}},\ {\boldsymbol {p}},\ {\boldsymbol {q}}$) of rescaled wavenumbers with ${\boldsymbol {k}}=(k_f=1,0,0)$ being the forcing wavenumber, implying $p_x=-k_f-q_x$, $p_y=-q_y$, $p_z=-q_z$. We choose the forcing-scale mode $Z^{+}_{\boldsymbol {k}}=u_0$, $Z^{-}_{\boldsymbol {k}}=0 \Leftrightarrow \hat {{\boldsymbol {u}}}_k =u_0 (0,-i/2,1/2)$, i.e. the positively helical flow ${\boldsymbol {u}}=u_0 (0,\sin (k_f x), \cos (k_f x))$. We take the modes at ${\boldsymbol {p}}$ and ${\boldsymbol {q}}$ to be small-amplitude inertial waves, and perform a linear stability analysis of this configuration for the homochiral case $s_{\boldsymbol {p}}=s_{\boldsymbol {q}}=1=s_{\boldsymbol {k}}$ (the other cases do not give relevant results). We thus determine the growth rate $\sigma ({\boldsymbol {q}})$ (${\boldsymbol {p}}$ is uniquely determined by ${\boldsymbol {q}}$) of the two inertial-wave modes, whose temporal evolution is given by $Z_{\boldsymbol {q}}^+,Z_{\boldsymbol {p}}^+\propto \exp (\sigma t)$.

Figure 2(a) shows that the maximum of $\sigma$ occurs for small wavenumbers at $q_{\perp }\approx k_f/2$, while figure 2(b) indicates that this maximum is located at $q_{\parallel }=0$.

Figure 2.

($a$) Maximum (over $q_{\parallel }$) of the real part of the growth rate as a function of $q_{\perp }$. The maximum is found for $q_{\perp } \approx k_f/2$. ($b$) Maximum (over $q_{\perp }$) of the real part of the growth rate as a function of $\lambda q_{\parallel }$ with $k_f=1$. There is a monotonic decay with $q_{\parallel }$ up to $\lambda q_{\parallel }\approx 0.35$ and vanishing values beyond this point.

Thus the 2-D base flow becomes unstable to smaller perpendicular wavenumbers $q_{\perp }$ and for a range of parallel wavenumbers $|q_{\parallel }| \lesssim 0.35/\lambda$.

Since the rescaled layer height is given by $2\pi$, the minimum parallel wavenumber is $q_{\parallel }^{\min } = 1$. Thus the $q_{\parallel }\ne 0$ modes are unstable only if $\lambda \lesssim 0.35$. For $\lambda$ larger than this value all $q_{\parallel }\ne 0$ wavenumbers are stable. For large values of $\lambda$ therefore the 2-D modes $q_{\parallel }=0$ are expected to decouple from the 3-D modes $q_{\parallel }\ne 0$ and an inverse cascade is expected. Conversely, for small values of $\lambda$ the 3-D modes become unstable and can possibly redirect energy back to small scales. We find the signature of this triadic instability in the direct numerical simulation results which we discuss in § 4.