3.1. $n=0$ , $\mathit{Re}_{F}^{-1}={\it\lambda}_{0}\mathit{Ro}_{F}^{-1}$

For $n=0$ the viscous term is of the same order as the rotation term and the operator that needs to be inverted to obtain the first-order term in $\boldsymbol{v}$ is $(\mathbb{L}-{\it\lambda}_{0}{\it\Delta})$ . This operator is positive definite and can always be inverted to obtain the zeroth-order solution. In detail for the TG forcing (1.5) it is obtained that

The next order correction follows in a similar manner without a qualitative change the of the flow behaviour. We note that the stationary solution in this limit follows the scaling $\langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}\rangle _{V}=\mathit{Ro}_{F}^{2}\langle \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{v}\rangle _{V}=3\mathit{Ro}_{F}^{2}(1+27{\it\lambda}_{0}^{2}q^{4})^{-1}$ . That implies that the relation between $(\mathit{Re}_{U},\mathit{Ro}_{U})$ and $(\mathit{Re}_{F},\mathit{Ro}_{F})$ for the laminar flow is given by

Furthermore, $\langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{f}\rangle =9q^{2}\mathit{Ro}_{F}{\it\lambda}_{0}(1+27{\it\lambda}_{0}^{2}q^{4})^{-1}$ and thus

These results (3.8), (3.9) provide the map from the $(\mathit{Re}_{U},\mathit{Ro}_{U})$ space (figure 1 a) to the $(\mathit{Re}_{U},\mathit{Ro}_{U})$ and $(\mathit{Re}_{D},\mathit{Ro}_{D})$ space shown in figure 2 but only for the laminar solutions.

3.2. $n=1$ , $\mathit{Re}_{F}^{-1}={\it\lambda}_{1}\mathit{Ro}_{F}^{1}$

For this value of the ordering index, equation (3.2) becomes

Since $\overline{\boldsymbol{f}}=0$ the zeroth-order solution can be written as $\boldsymbol{v}^{(0)}=\boldsymbol{v}_{\boldsymbol{f}}^{(0)}+\boldsymbol{v}_{2D}^{(0)}$ where

and $\boldsymbol{v}_{2D}^{(0)}$ is an arbitrary 2D3C field. The velocity field in (3.11) is the same as the one in (3.7) with ${\it\lambda}_{0}=0$ . The 2D3C field is determined at next order for which we have

Averaging over the $z$ direction we obtain by detailed calculation that

where $\boldsymbol{w}_{2D}^{(0)}=\boldsymbol{{\rm\nabla}}\times \boldsymbol{v}_{2D}^{(0)}$ and $\boldsymbol{w}_{\boldsymbol{f}}^{(0)}=\boldsymbol{{\rm\nabla}}\times \boldsymbol{v}_{\boldsymbol{f}}^{(0)}$ . We are thus left with

Multiplying by $\boldsymbol{v}_{2D}^{(0)}$ and horizontal averaging we get $\langle (\boldsymbol{{\rm\nabla}}\boldsymbol{v}_{2D}^{(0)})^{2}\rangle _{V}=0$ and thus the undetermined 2D3C field at first order becomes $\boldsymbol{v}_{2D}^{(0)}=0$ . A 2D3C flow can possibly appear at higher order but we will proceed no further than obtaining $\boldsymbol{v}^{(0)}$ .

The relations given in (3.8) simplify (by setting ${\it\lambda}_{0}=0$ ) to

for the Reynolds numbers based on the root-mean-square velocity. Similarly (3.9) simplify to $\mathit{Re}_{D}=3^{2/3}\,\mathit{Re}_{F}^{1/3}\mathit{Ro}_{F}^{2/3}$ and $\mathit{Ro}_{D}=3^{2/3}\,\mathit{Re}_{F}^{-1/3}\mathit{Ro}_{F}^{4/3}$ for the Reynolds numbers based on the energy injection rate.

3.3. $n=2$ , $\mathit{Re}_{F}^{-1}={\it\lambda}_{2}\mathit{Ro}_{F}^{3}$

For this ordering index the expansion follows as before but viscosity is not present at first order and the solvability condition obtained in (3.12) only restricts $\boldsymbol{v}_{2D}$ at being a solution of the Euler equations $\mathbb{P}[\boldsymbol{v}_{2D}^{(0)}\times \boldsymbol{w}_{2D}^{(0)}]=0$ . The first-order correction is then

To proceed at second order and simplify the mathematical complexity that increases rapidly with the order of the expansion we will assume at this point that $\boldsymbol{v}_{2D}=0$ and verify a posteriori the validity of this assumption. To the next order we obtain

Averaging over $z$ the nonlinear term drops to zero. This can be realized by noting that the zeroth-order components ( $\boldsymbol{v}^{(0)},\boldsymbol{w}^{(0)}$ ) only contain wavenumbers $\pm q$ in the $z$ direction while the first-order component ( $\boldsymbol{v}^{(1)},\boldsymbol{w}^{(1)}$ ) has only $\pm 2q$ wavenumbers. Their product will thus have only $\pm q$ and $\pm 3q$ wavenumbers and will thus average to zero. We obtain thus $\overline{{\rm\nabla}^{2}\boldsymbol{v}^{(0)}}=0$ in agreement with our original assumption that $\boldsymbol{v}_{2D}=0$ . Up to this scaling therefore the stationary solution obtained has zero projection to the 2D3C fields.

3.4. $n=3$ , $\mathit{Re}_{F}^{-1}={\it\lambda}_{3}\mathit{Ro}_{F}^{5}$

The assumption that $\boldsymbol{v}_{2D}=0$ ceases to be true for this ordering. If we assume that $\boldsymbol{v}_{2D}^{0}=0$ , then at third order where the dissipation term is present we will obtain

where $\boldsymbol{v}^{(0)}$ , $\boldsymbol{v}^{(1)}$ and $\boldsymbol{v}^{(2)}$ are obtained from (3.11), (3.16) and (3.17), respectively, with ${\it\lambda}_{2}=0$ . Averaging over $z$ the nonlinear term will not become zero and thus will need to be balanced by the viscosity:

Thus, $\boldsymbol{v}_{2D}=0$ is no longer an acceptable solution. Calculation of the nonlinear term (3.19) leads to

with

Thus, even if the initial data have a zero projection to the 2D3C fields the nonlinearity at third order will force a 2D3C component to grow to zeroth-order amplitude.

3.5. $n>4$ , $O(\mathit{Re}_{F}^{-1})<O(\mathit{Ro}_{F}^{5})$

At higher values of $n$ the nonlinearity at the fourth-order term has non-zero projection to the 2D3C fields and viscosity is not present to provide a balance. Thus, no solution can be obtained by the expansion. Physically this would imply that even if initial conditions were prepared carefully to avoid any instabilities the 2D3C component of the flow would grow with time to values where the original scaling $\boldsymbol{u}=\mathit{Ro}_{F}\boldsymbol{v}\ll 1$ would not be valid and the expansion would fail. This also indicates that the quasi-2D condensates that were shown in figure 1 will show up at this order. This gives a lower bound on the unknown exponent ${\it\alpha}\geqslant 1/5$ since at this order there is direct injection of energy in the 2D3C flow.

4.1. Energy stability

The energy stability method provides a sufficient criterion for stability of a stationary flow subject to perturbations of arbitrary amplitude. It is based on the energy evolution equation of the perturbation. Denoting the stationary solution as $\boldsymbol{u}_{f}$ and without taking any asymptotic limit yet we can write the velocity field $\boldsymbol{u}$ as the sum of the stationary solution $\boldsymbol{u}_{f}$ and a fluctuating part $\boldsymbol{u}_{p}$ . Multiplying (1.1) by the fluctuating part $\boldsymbol{u}_{p}$ and space averaging leads to the equation

The stationary solution is then called energy stable if the right-hand side is negative definite for all incompressible fields $\boldsymbol{u}_{p}$ . Note that the Coriolis term is absent in (4.2) and the same equation is obtained if $\mathit{Ro}_{F}=\infty$ . However, there is dependence on $\mathit{Ro}_{F}$ through the functional form of $\boldsymbol{u}_{f}$ . Stability of the stationary solution is then guaranteed if

in which case the energy of the perturbation $\boldsymbol{u}_{p}$ is bounded by

by Gronwall’s inequality. (Detailed derivations of the inequalities used in this section can be found in Doering & Gibbon (Reference Doering and Gibbon1995).) The energy stability analysis then reduces to solving the involved Euler–Lagrange equations and determining ${\it\mu}(\mathit{Re}_{F},\mathit{Ro}_{F})$ . The energy stability boundary is then determined by the ${\it\mu}(\mathit{Re}_{F},\mathit{Ro}_{F})=0$ curve in the parameter space. Here it is not attempted to solve the Euler–Lagrange equations that are beyond the point of this work but estimates are given based on simple inequalities and the previous approximations of steady-state flow.

We begin by noting that Poincaré’s inequality indicates that $\langle |\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{p}|^{2}\rangle _{V}\geqslant \langle |\boldsymbol{u}_{p}|^{2}\rangle _{V}$ (where the fact that the flow is defined in a periodic box of non-dimensional size $2{\rm\pi}$ has been used). Hölder’s inequality also indicates that

where by $\Vert \boldsymbol{{\rm\nabla}}\boldsymbol{u}_{f}\Vert _{\infty }$ we indicate the maximum of the modulus of the strain tensor $\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{f}$ over space. These estimates together with the condition (4.3) indicate that stability of the flow $\boldsymbol{u}_{f}$ is guaranteed if $\Vert \boldsymbol{{\rm\nabla}}\boldsymbol{u}_{f}\Vert _{\infty }\leqslant \mathit{Re}_{F}^{-1}$ .

This criterion for stability can be calculated using the asymptotic solution we obtained in the previous section. The expression (3.7) gives $\boldsymbol{u}_{f}$ up to an order $\mathit{Ro}_{F}^{-2}$ correction. For ${\it\lambda}_{0}=0$ (3.7) reduces to the same zeroth-order solution as in the $n=1$ and $n=2$ case (see (3.11)), thus this choice is valid for all cases for which $\mathit{Re}_{F}^{-1}\leqslant O(\mathit{Ro}_{F}^{3})$ . For higher ordering the zeroth-order 2D3C field $\mathit{Ro}_{F}^{-1}\boldsymbol{v}_{2D}$ should also be included. Substituting thus (3.7) it is obtained

Substituting in $\Vert \boldsymbol{{\rm\nabla}}\boldsymbol{u}_{f}\Vert _{\infty }\leqslant \mathit{Re}_{F}^{-1}$ and using ${\it\lambda}_{0}=\mathit{Ro}_{F}/\mathit{Re}_{F}$ the stability of the flow calculated for the approximate solution (3.7) is guaranteed if

This leads to the following estimates that guarantee stability. The flow is stable if

and for any value of $\mathit{Ro}_{F}$ stability is guaranteed if

These conservatives estimates of the energy stability boundaries are shown in figure 1 with a solid grey (orange online) line. For $\mathit{Ro}_{F}\gg 1$ the estimate for stability is given by $\mathit{Re}_{F}\leqslant \sqrt{3q/2}$ while the numerical simulations indicate a critical value for $\mathit{Re}_{F}=14.0$ which is considerable larger than the value obtained in (4.9). For $\mathit{Ro}_{F}\ll 1$ (4.8) leads to the condition for stability $\mathit{Re}_{F}\leqslant (4q\mathit{Ro}_{F})^{-1}$ while from the numerical simulations stability is shown in the same limit for $q=2$ if $\mathit{Re}_{F}\lesssim 10\mathit{Ro}_{F}^{-1}$ . The analytical estimates for the stability boundaries (4.8) and (4.9) are thus very conservative which is expected for the rough estimates used here. These estimates however do not aim at obtaining the exact values but rather to give an understanding of the shape of the stability boundary that they clearly capture. They also prove that for any value of $\mathit{Re}_{F}$ , no mater how large, at sufficiently large $\mathit{Ro}_{F}$ the flow will re-laminarize.

4.2. Linear stability

The absence of energy stability does not necessarily imply instability. This is determined by linear stability theory that investigates the evolution of infinitesimal perturbations and is examined in this section. In contrast to energy stability, linear stability does not guarantee stability since a linearly stable system can still be unstable to finite-amplitude perturbations. However, absence of linear stability always guarantees instability. The two methods are thus complementary. For $\mathit{Ro}_{F}\rightarrow 0$ the energy stability indicates that stability is determined at $\mathit{Re}_{F}=O(\mathit{Ro}_{F}^{-1})$ . We thus begin the investigation of linear stability for this ordering. The linear equation for an infinitesimal fluctuation $\boldsymbol{v}_{p}$ on the basic stationary state $\boldsymbol{v}_{f}$ (given at this order by (3.11)) reads

Unlike in the energy stability method, in linear stability the effect of rotation is not removed and the smallness of $\mathit{Ro}_{F}$ can be used to find the growth rate of the perturbations in an asymptotic way. Like before we write $\boldsymbol{v}_{p}=\boldsymbol{v}_{p}^{(0)}+\mathit{Ro}_{F}^{2}\boldsymbol{v}_{p}^{(1)}+\cdots$ . At zeroth order we obtain

The general solution of (4.11) is then given by

where $\boldsymbol{v}_{2D}$ is a 2D3C velocity field that depends only on the ( $x,y$ ) coordinates and on the slow dynamical time ${\it\tau}$ . Here $\boldsymbol{v}_{w}^{(0)}$ are inertial waves whose amplitudes also vary on the long timescale. The general expression for the inertial waves in triple periodic boxes is given by

where $s=\pm 1$ , $\text{c.c.}$ stands for complex conjugate and represents the wavenumbers with negative $k_{z}$ . Here $\boldsymbol{h}_{\boldsymbol{k}}^{s}$ are the helical basis in Fourier space described by Cambon & Jacquin (Reference Cambon and Jacquin1989) and Waleffe (Reference Waleffe1992, Reference Waleffe1993) and $\boldsymbol{h}_{\boldsymbol{k}}^{s}$ are eigenfunctions of the curl operator with $\text{i}\boldsymbol{k}\times \boldsymbol{h}_{\boldsymbol{k}}^{s}=s|\boldsymbol{k}|\boldsymbol{h}_{\boldsymbol{k}}^{s}$ and thus the index $s$ determines the sign of the helicity of the mode. According to (3.4) then

On the periodic box used here for any wavevector $\boldsymbol{k}=[k_{x},k_{y},k_{z}]$ with $k_{\bot }=\sqrt{k_{x}^{2}+k_{y}^{2}}\neq 0$ , $\boldsymbol{h}_{\boldsymbol{k}}^{s}$ is defined as

for the $k_{\bot }=0$ case.

On the fast timescale no instability exists since rotation cannot transfer energy to the perturbation field. The slow timescale evolution of the two fields $\boldsymbol{v}_{2D}^{(0)}$ and $\boldsymbol{v}_{w}^{(0)}$ is obtained by a solvability condition at the next order

At this point we need to define the short time average $\langle g\rangle _{{\it\tau}}$ as

with $-1\leqslant {\it\sigma}<1$ (i.e. the time average is over a timescale much smaller than ${\it\tau}$ and much bigger than $t^{\prime }$ ). With this definition the short time average of a fast oscillating function $g(t)=\text{e}^{\text{i}t^{\prime }}=\text{e}^{\text{i}t/\mathit{Ro}_{F}}$ becomes $\langle g\rangle _{{\it\tau}}=O(\mathit{Ro}_{F}^{1+{\it\sigma}})\ll 1$ and the short time average of a slow oscillating function $g(t)=\text{e}^{\text{i}{\it\tau}}=\text{e}^{\text{i}t\mathit{Ro}_{F}}$ becomes the identity $\langle g\rangle _{{\it\tau}}=1$ . Note that the result is independent of the choice of ${\it\sigma}$ .

To obtain then the evolution of $\boldsymbol{v}_{2D}^{(0)}$ we perform a short time average and an average over the $z$ direction. The left-hand side, and the advection term on the right averages to zero. The vorticity advection term also averages to zero because the vertical average will eliminate all Fourier modes of $\boldsymbol{v}_{p}$ with $k_{z}$ different from that of the forcing $k_{z}\neq q$ . The remaining terms oscillate with frequency $\pm k_{z}/|\boldsymbol{k}|\neq 0$ and will be eliminated by the time average. We are thus left with

which is a diffusion equation and thus any 2D3C perturbation will decay. Therefore, at this order the stationary flow is linearly stable to two-dimensional perturbations.

To obtain the evolution of $\boldsymbol{v}_{w}^{(0)}$ we multiply (4.16) with $\boldsymbol{h}_{\boldsymbol{k}}^{-s_{k}}\text{e}^{-\text{i}(\boldsymbol{k}\boldsymbol{x}+{\it\omega}_{\boldsymbol{k}}t^{\prime })}$ space average and short time average over to obtain the evolution equation for complex amplitude $H_{\boldsymbol{k}}^{s_{k}}$ :

Here $V_{\boldsymbol{q}}^{s_{q}}$ is the complex amplitude of the helical modes of the stationary flow (3.11). It is defined as $V_{\boldsymbol{q}}^{s_{q}}=\langle \text{e}^{-\text{i}\boldsymbol{q}\boldsymbol{x}}\boldsymbol{h}_{\boldsymbol{q}}^{-s_{q}}\boldsymbol{v}_{f}\rangle _{V}$ where $\boldsymbol{q}=(\pm q,\pm q,\pm q)$ is one of the eight forcing wavenumbers. The summation is over all wavenumbers $\boldsymbol{p}\in \mathbb{Z}^{3}$ and $\boldsymbol{q}$ such that $\boldsymbol{k}=\boldsymbol{q}+\boldsymbol{p}$ . The coupling tensor $\unicode[STIX]{x1D60A}_{\boldsymbol{k},\boldsymbol{p},\boldsymbol{q}}^{s_{k},s_{q},s_{p}}$ is given by

The short time average leads to different possibilities that we discuss here in some detail. First if $|{\it\omega}_{p}-{\it\omega}_{k}|=O(1)$ , then the averaged term $\langle \text{e}^{\text{i}({\it\omega}_{p}-{\it\omega}_{k})t^{\prime }}\rangle _{{\it\tau}}$ , becomes of smaller order and thus can be neglected. We refer to these terms as non-resonant terms. If $|{\it\omega}_{p}-{\it\omega}_{k}|=0$ , then the averaged term leads to an order-one contribution $\langle \text{e}^{\text{i}({\it\omega}_{p}-{\it\omega}_{k})t^{\prime }}\rangle _{{\it\tau}}=1$ and needs to be taken into account.

The third intermediate option is when the frequency difference is very small but non-zero. For ${\rm\Delta}{\it\omega}=|{\it\omega}_{p}-{\it\omega}_{k}|=O(\mathit{Ro}_{F}^{2})$ , then the average $\langle \text{e}^{\text{i}({\it\omega}_{p}-{\it\omega}_{k})t^{\prime }}\rangle _{{\it\tau}}$ still remains an order-one quantity. Here two conflicting arguments can be put forward for the role of these quasi-resonances. On the one hand since the value of ${\rm\Delta}{\it\omega}$ for any two wavenumbers $\boldsymbol{p},\boldsymbol{q}$ is fixed and independent of $\mathit{Ro}_{F}$ in the limit $\mathit{Ro}_{F}\rightarrow 0$ the quasi-resonances can be neglected. On the other hand, from the infinity of possible pairs $(\boldsymbol{k},\boldsymbol{p})\in \mathbb{Z}^{6}$ one can always find a pair of wavenumbers $\boldsymbol{k},\boldsymbol{p}$ such that ${\rm\Delta}{\it\omega}=O(\mathit{Ro}_{F}^{2})$ for any value of $\mathit{Ro}_{F}$ . Thus, there are always wavenumbers for which quasi-resonances are important. The resolution of course comes from the fact that not all wavenumbers are available. Any finite value of ${\it\lambda}_{1}$ will introduce a cutoff wavenumber $k_{{\it\lambda}}$ such that all wavenumbers $\boldsymbol{k}$ with $|\boldsymbol{k}|>k_{{\it\lambda}}$ will be damped by viscosity. Thus, pairs from all remaining wavenumbers $|\boldsymbol{k}|\leqslant k_{{\it\lambda}}$ will have finite quasi-resonance frequency ${\rm\Delta}{\it\omega}$ and in the limit $\mathit{Ro}_{F}\rightarrow 0$ can be neglected. Note the importance in this argument that ${\it\lambda}_{1}$ is finite. If the viscous term appeared at next order (i.e.  $\mathit{Re}_{F}=O(\mathit{Ro}_{F})^{-3}$ ) this argument would break down and quasi-resonances that appear for large enough $|\boldsymbol{k}|$ would need to be taken into account.

The cutoff wavenumber $k_{{\it\lambda}}$ can be estimated by the balance of the shear of the basic flow $\langle (\boldsymbol{{\rm\nabla}}\boldsymbol{v}_{f})^{2}\rangle _{V}^{1/2}\propto q$ with the viscous damping ${\it\lambda}_{1}k_{{\it\lambda}}^{2}$ that leads to the estimate

Thus, for small values of $\mathit{Ro}_{F}$ quasi-resonances can be neglected if ${\rm\Delta}{\it\omega}\gg O(\mathit{Ro}_{F}^{2})$ for all Fourier wavenumbers in a sphere of radius $k_{{\it\lambda}}$ . In what follows we give estimates for ${\rm\Delta}{\it\omega}$ and the density of the quasi-resonances as well as solutions for the exact resonances to estimate the validity of these assumptions in our simulations.

For exact resonances of (4.19) the following resonance conditions need to be satisfied for the wavenumbers $\boldsymbol{k},\boldsymbol{p},\boldsymbol{q}$ :

where $\boldsymbol{q}$ is one of the eight forcing wavevectors $(\pm q,\pm q,\pm q)$ that are located at the corners of a $2q$ -length cube centred at the origin. The sign $s=s_{p}s_{k}=\pm 1$ in (4.21b ) indicates whether the coupling is between waves travelling in the same direction ( $s=+1$ ) or opposite direction ( $s=-1$ ). Waves travelling in the same direction implies that the two wavenumbers lie in the same cone $\boldsymbol{k},\boldsymbol{p}$ of opening angle ${\it\theta}$ such that $\tan ({\it\theta})=k_{z}/|k|=p_{z}/|p|={\it\omega}$ while for opposite travelling waves the coupling is between wavenumbers $\boldsymbol{k},\boldsymbol{p}$ that are on opposite cones. Figure 3(a) indicates two such couplings. The projection of the three vectors $\boldsymbol{k},\boldsymbol{q},\boldsymbol{p}$ in the plane perpendicular to the rotation axis is shown in figure 3(b). It is obvious that the three projected vectors $\boldsymbol{k}_{\bot },\boldsymbol{q}_{\bot },\boldsymbol{p}_{\bot }$ can satisfy (4.21a ) if the following triangle inequalities hold $|\,|\boldsymbol{p}_{\bot }|-|\boldsymbol{k}_{\bot }|\,|\leqslant |\boldsymbol{q}_{\bot }|\leqslant |\boldsymbol{k}_{\bot }|+|\boldsymbol{p}_{\bot }|$ . The extreme cases being when $\boldsymbol{k}_{\bot },\boldsymbol{p}_{\bot }$ are parallel or antiparallel to $\boldsymbol{q}_{\bot }$ . Using (4.21b ) we then end up with the following bounds for the allowed wavenumbers lying on the same cone:

while for opposite cones the bounds are

Figure 3.

(a) Two possible couplings of the wavevectors $\boldsymbol{k},\boldsymbol{p}$ (red online) with the forcing wavenumber $\boldsymbol{q}$ (blue online). (b) The $(k_{x},k_{y})$ -plane where two circles of radius $|\boldsymbol{p}_{\bot }|$ and $|\boldsymbol{q}_{\bot }|$ have been drawn. In order for the forcing wavevector $\boldsymbol{q}_{\bot }=(q,q,0)$ to be able to form a triangle with two wavenumbers $\boldsymbol{p}_{\bot }$ and $\boldsymbol{k}_{\bot }$ its two ends need to be placed at the edges of the two circles as shown by the (blue online) arrow. Thus, $\boldsymbol{q}_{\bot }$ needs to be longer than the distance $|p_{\bot }-k_{\bot }|$ and shorter than the distance $|p_{\bot }+k_{\bot }|$ .

The triangle inequalities however do not guarantee the existence of resonances because the wavenumbers, $\boldsymbol{k},\boldsymbol{q},\boldsymbol{p}$ are discrete and the (4.22) need to be solved on a discrete lattice. General solutions of such problems however are non-trivial and have been solved but for very few cases (Bustamante & Hayat Reference Bustamante and Hayat2013). However, there are two simple solutions that can be found by simple inspection. The first is simply when $\boldsymbol{p}=m\boldsymbol{q}/q$ for any integer $m\in \mathbb{Z}$ . Then $\boldsymbol{k}$ lies in the same cone as $\boldsymbol{p}$ and is given by $\boldsymbol{k}=(m\pm 1)\boldsymbol{q}/q$ and the resonance condition (4.22b ) is satisfied exactly. This condition however leads to zero nonlinearity as it corresponds to coupling of parallel shear layers that are exact solutions of the Euler equations. The second type of exact resonances is obtained for

true for every $s_{x},s_{y},s_{z}=\pm 1$ and $m\in \mathbb{Z}$ . This solution couples modes in opposite cones and results in non-zero coupling. Note that these modes only exist when $q$ is even since $q/2$ that appears in (4.25) has to be an integer. For $q=2$ , that is examined here, these two solutions were the only exact solutions found by numerically solving (4.22). For larger values of $q$ more exact solutions were found but not a simple analytical formula for them. Figure 4(a) shows in the plane $k_{z},k_{\bot }=\sqrt{k_{x}^{2}+k_{y}^{2}}$ and for $q=2$ the location of the exact resonances of the first type $\boldsymbol{k}=m\boldsymbol{q}/q$ by (red online) squares and of the second type (4.25) by (red online) triangles. The dashed lines indicate the bounds (4.23) and (4.24).

Figure 4.

(a) The location of the exact resonances $\boldsymbol{k}=m\boldsymbol{q}/q$ by (red online) squares, the exact resonances (4.25) by (red online) triangles, quasi-resonances by (blue online) circles for a frequency tolerance ${\rm\Delta}{\it\omega}=10^{-3}$ in the plane $k_{z},k_{\bot }=\sqrt{k_{x}^{2}+k_{y}^{2}}$ . The dashed lines indicate the bounds (4.23) and (4.24). (b) Number of modes $N_{{\rm\Delta}{\it\omega}}(k)$ that satisfy the quasi-resonance conditions up to a frequency ambiguity ${\rm\Delta}{\it\omega}$ on the spherical shell $k<|\boldsymbol{k}|\leqslant k+1$ .

The modes with the smallest value of $|\boldsymbol{k}|$ and $|\boldsymbol{p}|$ that fall in the set (4.25) for $q=2$ are obtained for $m=0$ and $m=-1$ . For example, a simple pair of coupled modes is given by $\boldsymbol{k}=(2,0,1)$ and $\boldsymbol{p}=(0,-2,-1)$ for $\boldsymbol{q}=(2,2,2)$ . Keeping only these modes we can write a minimal model for their linear evolution:

where $\mathscr{M}_{k}=\unicode[STIX]{x1D60A}_{\boldsymbol{k},\boldsymbol{q},\boldsymbol{p}}^{+,+,-}V_{\boldsymbol{q}}^{+}+\unicode[STIX]{x1D60A}_{\boldsymbol{k},\boldsymbol{q},\boldsymbol{p}}^{+,-,-}V_{\boldsymbol{q}}^{-}$ and $\mathscr{M}_{p}=\unicode[STIX]{x1D60A}_{\boldsymbol{p},\boldsymbol{q},\boldsymbol{k}}^{-,+,+}V_{\boldsymbol{q}}^{+}+\unicode[STIX]{x1D60A}_{\boldsymbol{p},\boldsymbol{q},\boldsymbol{k}}^{-,-,+}V_{\boldsymbol{q}}^{-}$ . This type of three-mode model was used by Waleffe (Reference Waleffe1993) to predict the transfer of energy in Fourier space. Note however that in this case one of the three modes is always the non-oscillating forced mode. The model predicts instability when ${\it\lambda}_{1}$ is smaller than a critical value ${\it\lambda}_{1c}=3/20\sqrt{5}\simeq 0.065\dots$ implying instability for $\mathit{Re}_{F}\geqslant 14.9\dots \mathit{Ro}_{F}^{-1}$ . The stability from the numerical was estimated closer to $\mathit{Re}_{F}\geqslant 20\mathit{Ro}_{F}^{-1}$ indicating that ignoring the higher exact resonances underestimated the stability boundary and/or that sufficiently large $\mathit{Ro}_{F}$ (so that quasi-resonances can be neglected) has not be reached yet in our simulations.

As discussed before in the $\mathit{Ro}_{F}\rightarrow 0$ limit only the exact resonances need to be retained. However, since numerical simulations always operate on finite values of $\mathit{Ro}_{F}$ beside the exact resonances we also need to consider quasi-resonances for which (4.22b ) is satisfied up to an ${\rm\Delta}{\it\omega}=O(\mathit{Ro}_{F}^{2})$ accuracy. The resonance conditions (4.22) define a surface embedded in the $\mathbb{R}^{3}$ space of $\boldsymbol{k}$ where they are satisfied. Allowing for a ${\rm\Delta}{\it\omega}$ ambiguity implies that modes within a distance ${\rm\Delta}k$ normal to the resonance surface should also be considered. Here ${\rm\Delta}k$ is estimated by Taylor expansion ${\rm\Delta}{\it\omega}\simeq ({\rm\Delta}\boldsymbol{k})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}_{\boldsymbol{k}}({\it\omega}_{\boldsymbol{k}}-{\it\omega}_{\boldsymbol{k}+\boldsymbol{q}})\propto ({\rm\Delta}k)q/k^{2}$ , where the last relation is valid in the $q\ll |\boldsymbol{k}|$ limit. The volume $V_{k}$ in Fourier space occupied by these allowed wavenumbers inside a sphere of radius $k$ is then proportional the area of the resonance surface times the thickness ${\rm\Delta}k$ , thus $V_{k}\propto k^{2}{\rm\Delta}k$ . The total number of quasi-resonances on a spherical shell of radius $k$ and width ${\it\delta}k=1$ will scale as $N_{{\rm\Delta}{\it\omega}}(k)=n_{k}\partial _{k}V_{k}$ thus

This scaling is valid provided that ${\rm\Delta}k\simeq {\rm\Delta}{\it\omega}k^{2}/q$ remains small. For large values of $k$ however, ${\rm\Delta}k$ becomes so large that all of the modes inside the spherical shell satisfy the resonance conditions (4.22). The scaling will thus transition to $N_{{\rm\Delta}{\it\omega}}(k)=4{\rm\pi}k^{2}$ . Quasi-resonances were sought for numerically and the results are shown in figure 4. Figure 4(a) shows the location of the quasi-resonances by (blue online) circles for a frequency tolerance ${\rm\Delta}{\it\omega}=10^{-3}$ . Figure 4(b) indicates number of modes that satisfy the quasi-resonance conditions up to a frequency ambiguity ${\rm\Delta}{\it\omega}$ on the spherical shell $k<|\boldsymbol{k}|\leqslant k+1$ . The notation ${\rm\Delta}{\it\omega}=\infty$ implies all pairs of modes in this spherical shell are counted independent of resonances. The scalings $k^{2}$ and $k^{3}$ predicted in the previous paragraph are shown as a reference. It is clear that quasi-resonances increase very rapidly and unless the viscous cutoff wavenumber $k_{{\it\lambda}}$ is sufficiently small, very high rotation rates will be required to justify their neglect.

5.1. Reduced equations

The previous section determined the values of the parameters for which a dynamical behaviour will be observed in the fast rotating limit. The properties of this dynamical behaviour however can only be determined when nonlinearities are taken into account. To understand the nonlinear behaviour the asymptotic expansion in the fast rotating limit is examined including the nonlinear terms. Following the same steps as we did in the previous section we begin from (4.1) with $n=1$ and expand $\boldsymbol{v}=\boldsymbol{v}^{(0)}+\mathit{Ro}_{F}^{2}\boldsymbol{v}^{(1)}+\cdots$ . At zeroth order we obtain the same linear equation for $\boldsymbol{v}$ as we did for $\boldsymbol{v}_{p}$ in (4.11). The general solution of which is then given by

As before $\boldsymbol{v}_{f}^{(0)}$ is the zeroth-order stationary solution (3.11), $\boldsymbol{v}_{2D}^{(0)}$ is a velocity field that depends only on the ( $x,y$ ) coordinates and on the slow dynamical time ${\it\tau}$ , $\boldsymbol{v}_{w}^{(0)}$ are inertial waves whose amplitude also varies on the long timescale and their exact form is given by (4.13).

To next order we then have

To obtain then the evolution of $\boldsymbol{v}_{2D}^{(0)}$ we perform a short time average and an average over the $z$ direction. The left-hand side averages to zero and we are left with

The only quadratic terms that remain non-zero after the vertical average are

The first two correspond to the terms that appear in linear theory and we have already shown that lead to zero contribution after the short time average. The third term can be written using the helical mode expansion (4.13) as

The vector $\boldsymbol{D}_{\boldsymbol{k},\boldsymbol{p}}^{s_{k},s_{p}}$ is given by

(where $\boldsymbol{r}=\boldsymbol{k}+\boldsymbol{p}$ and $\mathbb{P}_{\boldsymbol{k}}[\boldsymbol{g}_{\boldsymbol{k}}]=\boldsymbol{k}\times \boldsymbol{k}\times \boldsymbol{g}_{\boldsymbol{k}}/|\boldsymbol{k}^{2}|$ is the projection operator $\mathbb{P}$ in Fourier space). A remarkable result of Waleffe (Reference Waleffe1993) was that for exact resonances ${\it\omega}_{k}+{\it\omega}_{p}=0$ this term (5.5) also averages to zero. This can be realized by noting that $k_{z}+p_{z}=0$ and ${\it\omega}_{k}+{\it\omega}_{p}=0$ implies that $|\boldsymbol{k}|=|\boldsymbol{p}|$ and $s_{p}=s_{k}$ and therefore the vector $\boldsymbol{D}_{\boldsymbol{k},\boldsymbol{p}}^{s_{k},s_{p}}$ in (5.6) becomes zero and does not effect the 2D3C flow. At this order then the only nonlinear term left for the evolution of the 2D3C flow is the coupling with itself, that leads to

which is the unforced 2D Navier–Stokes equation. As result at this order the nonlinear evolution equation for the 2D3C flow completely decouples from the rest of the flow, and since there is no forcing term in (5.7), $\boldsymbol{v}_{2D}^{(0)}$ will decay to zero at long times.

In the simulations however for all dynamical regimes a 2D3C flow was present. In order to explain the presence of a 2D3C flow in the simulation we have to evoke either quasi-resonances and higher-order terms in the expansion or a breakdown of the fast rotating limit. We consider then quasi-resonances that could lead to the generation of a quasi-2D flow and allow for a frequency tolerance ${\rm\Delta}{\it\omega}\ll 1$ over which a violation of the resonance condition is allowed. For fixed $k_{z}=-p_{z}$ the coupling triads reside in the four-dimensional space defined by $(k_{x},k_{y},p_{x},p_{y})$ . In this space a four spherical shell of unit thickness and of radius $r=(k_{\bot }^{2}+p_{\bot }^{2})^{1/2}$ , has $N\propto r^{3}$ triads out of which the quasi-resonance condition confines the allowed triads to reside in a subspace of thickness ${\rm\Delta}k_{\bot }$ that limits their number to $N_{{\rm\Delta}{\it\omega}}\propto r^{2}{\rm\Delta}k$ . The thickness ${\rm\Delta}k$ is obtained from the dispersion relation ${\rm\Delta}{\it\omega}\simeq {\rm\Delta}k_{\bot }\partial _{k_{\bot }}{\it\omega}=-{\rm\Delta}k_{\bot }k_{z}k_{\bot }/k^{3}$ . Thus, ${\rm\Delta}k_{\bot }\propto {\rm\Delta}{\it\omega}k^{3}/k_{z}k_{\bot }\simeq {\rm\Delta}{\it\omega}r^{2}/k_{z}$ for $k_{\bot }\gg k_{z}$ . The number of these allowed modes then will scale as $N_{{\rm\Delta}{\it\omega}}(r)\propto r^{2}{\rm\Delta}k\propto r^{4}/k_{z}{\rm\Delta}{\it\omega}$ (where locality between the modes $k_{\bot }\sim p_{\bot }\sim r$ has been assumed). The number of quasi-resonances therefore grows with $k$ even faster than the quasi-resonances of the linear term in the previous section. Another noteworthy point is that the smallest values of $k_{z}$ have the largest number quasi-resonances ( $k_{z}=1$ for our case).

At figure 5(a) we plot the exact resonances, and the quasi-resonances for $k_{z}=1,4$ and ${\rm\Delta}{\it\omega}=0.01$ on the $k_{\bot },p_{\bot }$ plane. Owing to their large number only 1 in 100 points has been plotted. Figure 5(b) of the same figure shows the density $N_{{\rm\Delta}{\it\omega}}(r)$ for different values of $k_{z}$ . These modes however will contribute to the 2D3C dynamics with coupling coefficients of the order $\boldsymbol{D}_{\boldsymbol{k},\boldsymbol{p}}^{s_{k},s_{p}}=O({\rm\Delta}k)=O({\rm\Delta}{\it\omega})$ . Their contribution to the $\boldsymbol{v}_{2D}^{(0)}$ evolution will be of the order of $N_{{\rm\Delta}{\it\omega}}\boldsymbol{D}_{\boldsymbol{k},\boldsymbol{p}}^{s_{k},s_{p}}=O({\rm\Delta}{\it\omega}^{2})$ , and thus the same order as the next-order term in the expansion.

Figure 5.

(a) The location of resonances for the modes $\boldsymbol{k},\boldsymbol{p}$ . Exact resonances are shown by (red online) squares, and the quasi-resonances by (blue online) circles ( $k_{z}=-p_{z}=4$ ) and (purple online) triangles ( $k_{z}=-p_{z}=1$ ) for a frequency tolerance ${\rm\Delta}{\it\omega}=10^{-2}$ . (b) Number of modes $N_{{\rm\Delta}{\it\omega}}(r)$ that satisfy the quasi-resonance conditions up to a frequency ambiguity ${\rm\Delta}{\it\omega}$ on the spherical shell $r<\sqrt{k_{\bot }^{2}+p_{\bot }^{2}}|\leqslant r+1$ .

In order thus to capture this contribution the expansion will need to be carried at next order introducing a second slow timescale and increasing $\mathit{Re}_{F}$ to the scaling $\mathit{Re}_{F}={\it\lambda}_{2}\mathit{Ro}_{F}^{3}$ . The process that ‘lives’ on this second timescale although slower, on the long time limit that we investigate here could be dominant and explain the 2D3C dynamics observed in the simulations.

To obtain the evolution equation for the amplitude $H_{\boldsymbol{k}}^{s_{k}}$ of the inertial waves we multiply (5.2) with $\boldsymbol{h}_{\boldsymbol{k}}^{-s_{k}}\text{e}^{-\text{i}(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}+{\it\omega}_{k}t)}$ and space–time average. The right-hand side will then drop out and we are going to be left with the evolution equation for the complex amplitude:

The first term couples the inertial waves to the stationary flow. It is the same term that was studied in the linear stability section and is the term that injects energy into inertial waves. The second term corresponds to the term that couples the 2D3C flow to the inertial waves. This coupling term does not exchange energy with the 2D3C part of the flow as was discussed before. It can, however, redistribute the energy among the inertial wave modes modes $\boldsymbol{k},\boldsymbol{p}$ with $k_{z}=p_{z}$ and $k_{\bot }=p_{\bot }$ (see Waleffe Reference Waleffe1993). The third term couples all of the inertial waves with each other. In the fast rotating limit only exact resonances with will survive at this order after the short time average. For the wavenumbers $\boldsymbol{k}_{1},\boldsymbol{k}_{2},\boldsymbol{k}_{3}$ the following resonance conditions hold:

Analysing the statistical properties of the flow through these exact resonances is the subject of weak wave turbulence theory that has been discussed in Galtier (Reference Galtier2003) and leads to the prediction for the energy spectrum $E(k_{\bot })\propto k_{\bot }^{-5/2}$ (with the possibility of less steep spectra when the helicity flux dominates (Galtier Reference Galtier2014)).

5.2. Intermittent bursts

Given the theoretical background discussed in the previous section we can try to interpret the observed results in the simulations. As discussed in § 2 the striking feature of the flow close to the linear stability boundary is bursts of energy followed by a slow decay. The flow during this slow decay phase is dominated by a large 2D3C component. Figure 6(a) shows a series of such burst for a typical run in this range ( $\mathit{Re}_{F}=1000,\mathit{Ro}_{F}=0.125$ ). The dark (blue online) curve shows the evolution of the total energy and the light grey (red online) line shows the evolution of the energy of the 2D3C component of the flow. Figure 6(b) shows the same signal focused on a single burst that has been shifted to $t=0$ . It can be seen that at the birth of the burst the inertial waves (the non-2D3C component of the flow) become unstable and increase to a high amplitude. The inertial waves then drive the increase of the 2D3C component that follows after. After the burst the 2D3C flow dominates and slowly decays exponentially by a rate proportional to the viscosity. As a result for large values of $\mathit{Re}_{F}$ the decay is slow. Thus, in order to get well-converged averages, runs of long durations are needed.

More detail on the energy distribution during a burst can be seen in figure 7 where a greyscale (colour online) image of the energy spectrum on the $k_{\bot },k_{z}$ plane is shown. The four images correspond to the four times indicated by the vertical dashed lines in figure 6(b). The forcing wavenumber on these plots ( $k_{z}=2,k_{\bot }=2\sqrt{2}\simeq 3$ ) is indicated by a square while the first unstable mode predicted by the linear theory ( $k_{z}=1,k_{x}=2,k_{y}=0,\rightarrow k_{\bot }=2$ , see (4.26)) is denoted by the diamond. The evolution of the burst can then be described as follows: during the decay phase most of the energy is in the $k_{z}=0$ plane (2D3C flow) and a small part on the forcing wavenumber. As the amplitude of the 2D3C flow decays the amplitude of the forcing mode increases approaching the solution (3.11). As this solution is approached a point is reached that the $(2,0,1),(0,-2,-1)$ modes becomes unstable and starts to grow. When their amplitude becomes large enough (peak of the burst) the system becomes strongly nonlinear and allows for the violation of the resonant conditions. It thus couples to a large number of modes including modes with $k_{z}=0$ . The energy of the unstable mode is transferred in all of these modes and cascades to the dissipation scales. This is true for all but the 2D3C modes that due to their quasi-two-dimensionality do not cascade their energy to the small scales but in the large scales. They thus form a large-scale condensate at $k_{z}=0$ , $k_{\bot }=1$ that suppresses the instability. Afterwards it decays on the slow diffusive timescale since there is no further injection of energy to sustain it. This process is then repeated.

Figure 6.

(a) A series of bursts for ( $\mathit{Re}_{F}=1000,\mathit{Ro}_{F}=0.125$ ). The dark (blue online) curve shows the evolution of the total energy and the light grey (red online) line shows the evolution of the energy of the 2D3C component of the flow. (b) The same signal focused on a single burst that has been shifted to $t=0$ .

Figure 7.

Greyscale images of the energy spectrum on the $k_{\bot },k_{z}$ plane. The four images correspond to the four times indicated by the vertical dashed lines in figure 6(b).

This behaviour of a fast increase followed by a slow decay is typical in dynamical systems that pass through a hyperbolic point for which the rate of attraction at the attracting manifold is small. The simplest version perhaps of such a model is written as

where ${\it\epsilon}\ll 1$ . This model has closed flow lines given by $C=\ln (yx^{{\it\epsilon}})-(x^{2}+y^{2})/2$ shown in figure 8, with $x=y=0$ being a hyperbolic unstable point and $x=\sqrt{{\it\epsilon}}$ and  $y=1$ a neutrally stable point. Trajectories passing near the unstable point $(0,0)$ slowly approach it along the $y$ -axis, until the $x$ -mode becomes unstable leading a sudden increase of energy followed again by the slow decay. This procedure leads to the bursts shown in figure 8(b).

Figure 8.

(a) Phase-space trajectories of the model (5.10) and (b) time evolution of a particular solution of the same model.

A more realistic model can perhaps be written taking into account the dynamically most relevant modes in the system. Based on the results of the linear study we consider two complex amplitudes $x,y$ that represent the two unstable modes $(2,0,1)$ and $(0,-2,-1)$ , respectively, with frequency ${\it\omega}_{x}={\it\omega}_{y}={\it\omega}$ one complex amplitude $z$ for the forcing mode $(2,2,2)$ with frequency ${\it\omega}_{z}\neq {\it\omega}$ and one real amplitude $u$ representing the 2D3C field not affected by rotation:

The coupling of the modes $x,y,z$ follows the coupling of the $(2,0,1)$ , $(0,-2,-1)$ , $(2,2,2)$ modes from the Navier–Stokes equation while the coupling with 2D3C mode $u$ that is through strong multi-mode interaction is modelled as a higher-order nonlinearity. The nonlinear term conserves the energy $|z|^{2}+|x|^{2}+|y|^{2}+|u|^{2}$ . In the absence of the modes $x=y=u=0$ the $z$ -mode saturates at the value $z=\text{i}\,f_{0}/{\it\omega}_{z}=O(1)$ . However, it is an unstable solution to the two resonant modes $x,y$ when ${\it\lambda}<|\,f_{0}/{\it\omega}_{z}|$ . The different frequencies of the $z$ -mode and $x,y$ modes will make the nonlinear term $zy$ and $z^{\ast }x$ fast oscillating and thus ineffective to saturate the instability at $x,y=O(1)$ values. Saturation comes at order $x,y=O(1/{\it\epsilon})$ , but before this happens at order $1/\sqrt{{\it\epsilon}}$ the mode $u$ becomes unstable that absorbs almost all energy and decays on a timescale $1/{\it\lambda}$ as observed in the simulations.

A comparison of the results of the model with the DNS can be seen by comparing figure 9 that shows the evolution of the energy of the model with figure 6 that shows the results from DNS. Although the full complexity of the DNS is not recovered due to the simplicity of the coupling to the 2D3C modes in the model, the basic features are reproduced. It is mentioned here that the results of the simulations can be reproduced qualitatively by a finite-mode model because a finite box size is considered (i.e.  $q=2=O(1)$ ). In the case of long boxes different two-scale asymptotic expansions could be considered that would lead to reduced partial differential equations as found in many different set-ups in rotating convection (see for example Julien, Knobloch & Werne Reference Julien, Knobloch and Werne1998; Julien et al. Reference Julien, Knobloch, Milliff and Werne2006a ,Reference Julien, Knobloch, Rubio and Vasil b ).

Figure 9.

The evolution of different components of the energy of the model (5.11).

5.3. Two-dimensional condensates and fully turbulent flows

As discussed in § 2 in the parameter space for $\mathit{Re}_{F}>300$ and $1\geqslant \mathit{Ro}_{F}\geqslant \mathit{Re}_{F}^{-{\it\alpha}}$ the flow forms 2D3C condensates. This class of flows describes the fate of rotating turbulence obtained in the limit $\mathit{Re}_{F}\rightarrow \infty$ for any fixed value of $\mathit{Ro}_{F}$ with $\mathit{Ro}_{F}<0.4$ . The value of ${\it\alpha}$ depends on which order an injection of energy to the 2D3C flow first appears. Here ${\it\alpha}<1$ (i.e.  $n=1$ ) for which case the bursts are observed. If second-order terms couple inertial waves with the 2D3C flow and inject to it energy, then ${\it\alpha}=1/3$ (i.e.  $n=2$ ). If not, ${\it\alpha}=1/5$ (i.e.  $n=3$ ) in which case the stationary flow alone injects energy into the 2D3C component.

Flows in this state are characterized by large values of energy that is concentrated in the largest scales with $k_{z}=0$ . The evolution of the total energy and of the 2D3C component of the energy for the run with $\mathit{Ro}_{F}=0.25$ and $\mathit{Re}_{F}=500$ is shown in figure 10(a). After a short time ( $t<100$ ) during which the energy appears to saturate at a value close to unity the 2D3C flow starts to grow linearly until saturation is reached at very large values of energy. The non-2D3C part of the energy that consist of the energy of the inertial waves and the forced modes comprises only a small part of the energy. This behaviour was observed for all runs that showed a quasi-2D condensate behaviour.

The saturation amplitude is much larger than unity. Figure 10(b) shows the saturation level of the $U^{2}=\langle \langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}\rangle _{V}\rangle _{T}$ as a function of the Rossby $\mathit{Ro}_{F}$ for a fixed value of the Reynolds number $\mathit{Re}_{F}=333$ . Here $U^{2}$ varies discontinuously as $\mathit{Ro}_{F}$ is increased, both at the critical value that it transitions from isotropic turbulence to the condensate state at $\mathit{Ro}_{F}=0.39$ and at the critical value that it transitions from the condensate state to the intermittent bursts behaviour at $\mathit{Ro}_{F}=0.22$ . Thus, the transition from the isotropic turbulent state to this condensate state is found to be subcritical.

Figure 10.

(a) The evolution of the total energy and of the 2D3C component of the energy for the run with $\mathit{Ro}_{F}=0.25$ and $\mathit{Re}_{F}=500$ . (b) The saturation level of the $U^{2}$ as a function of the Rossby $\mathit{Ro}_{F}$ for a fixed value of the Reynolds number $\mathit{Re}_{F}=333$ .

The quasi-2D3C behaviour of the flow can also be seen by looking at the energy spectra. In figure 11 we show the 2D energy spectrum $E_{2D}(k_{z},k_{\bot })$ (a) and the 1D energy spectrum $E_{1D}$ (b) compensated by $k^{5/3}$ defined as

of the run with $\mathit{Ro}_{F}=0.20$ and $\mathit{Re}_{F}=2000$ that corresponds to the run with the largest $\mathit{Re}_{F}$ and smallest $\mathit{Ro}_{F}$ that displayed the quasi-2D3C condensate behaviour. The square in 11(a) indicates the location of the forcing. The dashed lines indicate the location of the circles $k_{z}^{2}+k_{\bot }^{2}=10,20,40,80$ . An isotropic spectrum would be constant along these lines. In the large scales the $k_{z}=0$ modes have significantly larger amplitude and energy is concentrated in $k_{\bot }=\sqrt{2}$ . The small scales appear more isotropic with the exception of the $k_{z}\gg k_{\bot }$ modes that appear to be quenched and deviate stronger from isotropy. The dominance of the large-scale modes can be observed in the 1D energy spectrum. A large amount of energy is concentrated at small wavenumbers followed by a steep drop. The small scales on the other hand display a power-law behaviour with an index slightly smaller than the Kolmogorov prediction $-5/3$ , and much bigger than the wave turbulence prediction $-5/2$ . The inset in 11(b) shows the time-averaged energy flux with the dark solid line, while the grey lines show the instantaneous energy flux. The flux is positive (forward) and almost constant for the wavenumbers $\boldsymbol{k}$ in the range $\sqrt{3}q<|\boldsymbol{k}|\lesssim 30$ , while it is slightly negative for the small wavenumbers ( $|\boldsymbol{k}|<\sqrt{3}q$ ), but with large fluctuations of both signs. It is expected that at even higher values of $\mathit{Re}_{F}$ the isotropic Kolmogorov spectrum will be even better established in the small scales.

Figure 11.

The 2D energy spectrum $E_{2D}(k_{z},k_{\bot })$  (a) and the 1D energy spectrum $E_{1D}$  (b). The inset shows the time-averaged energy flux with the dark solid line, while the grey lines show the instantaneous energy flux.

A visualization of the flow in the condensate state is depicted in figure 12 where the $z$ -component of the vorticity is shown. Figure 12(a) shows the computational box viewed from the top, while the same result is shown viewed from the side in figure 12(b) where the opacity has been reduced so that the structures inside the computational box can be seen. The blue color indicates that vorticity is parallel to the rotation (the flow corotates) while the red color indicates that the vorticity is antiparallel to the rotation. Clearly both large and small scales coexist in the flow. A large-scale, 2D corotating columnar structure can be seen. Opposite to this columnar structure a second columnar structure rotating in the opposite direction can be seen. This second structure (best seen in figure 12(b)) is fluctuating strongly having more small scales. The persistence of corotating structures in rotating turbulence has been observed in experiments (Hopfinger et al. Reference Hopfinger, Gagne and Browand1982; Morize et al. Reference Morize, Moisy and Rabaud2005; Gallet et al. Reference Gallet, Campagne, Cortet and Moisy2014) and has been discussed in various works (Hopfinger & van Heijst Reference Hopfinger and van Heijst1993; Bartello et al. Reference Bartello, Metais and Lesieur1994; Gence & Frick Reference Gence and Frick2001; Sreenivasan & Davidson Reference Sreenivasan and Davidson2008; Staplehurst et al. Reference Staplehurst, Davidson and Dalziel2008; Gallet et al. Reference Gallet, Campagne, Cortet and Moisy2014). Perhaps it is not surprising that structures with vorticity anti-aligned to rotation are more responsive to 3D perturbations. Following Bartello et al. (Reference Bartello, Metais and Lesieur1994) the local Rossby number $\mathit{Ro}_{loc}(\boldsymbol{x})$ (defined by the local rotation rate $2{\it\Omega}_{local}=2{\it\bf\Omega}+\boldsymbol{w}(\boldsymbol{x},t)$ ) will be increased and thus are less likely to show a quasi-2D3C behaviour.

Figure 12.

A visualization of the flow where the $z$ -component of the vorticity is shown. (a) The computational box viewed from the top and (b) the same flow shown viewed from the side where the opacity has been reduced so that the structures inside the computational box can be seen.

With this observation one can easily understand how saturation is reached in the large scales. The fast rotation leads to a quasi-two-dimensionalization of the system that leads to an inverse cascade of energy. As energy of the 2D3C part of the flow increases it reaches the largest scale of the system where it forms two counter-rotating vortexes. As the amplitude of these two vortices is increased a point is reached that the vertical vorticity of the counter-rotating vortex will cancel the effect of rotation locally when an order-one value of the local Rossby number is reached. Then the Taylor–Proudman constrain that leads to quasi-2D flows, the inverse cascade and the pile-up of energy to the large scales is broken and energy starts to flow back to the small scales. Such a mechanism of course implies that saturation is reached when the eddy turnover time of the condensate is the same order with the rotation period or more simply $\mathit{Ro}_{U}\sim 1$ . Indeed in figure 13 we plot $\mathit{Ro}_{U}$ from all of the runs that showed a condensate behaviour as a function of $\mathit{Ro}_{F}$ (figure 13 a) and as a function of $\mathit{Re}_{F}$ (figure 13 a). The value of $\mathit{Ro}_{U}$ appears to be independent of both $\mathit{Ro}_{F}$ and $\mathit{Re}_{F}$ , thus the scaling $\mathit{Ro}_{U}\sim 1$ is verified. This result implies that this state is unlikely to be described by an $\mathit{Ro}_{U}\rightarrow 0$ expansion as the eddy turnover time is the same order as the rotation period. It also implies that condensates saturate by reaching a state of marginal inverse cascade (Seshasayanan, Benavides & Alexakis Reference Seshasayanan, Benavides and Alexakis2014).

Figure 13.

Plots of $\mathit{Ro}_{U}$ as a function of $\mathit{Ro}_{F}$ (a) and $\mathit{Re}_{F}$ (b) for all of the runs that lead to a quasi-2D condensate. The $\mathit{Ro}_{U}=1$ behaviour shown here can also be seen in figure 2(a).