## 1. Introduction

Rotating turbulent flows are ubiquitous in geophysical and astrophysical systems such as stellar interiors, planetary cores, oceans and atmospheres. In a large number of numerical simulations and experiments (see the review by Godeferd & Moisy (Reference Godeferd and Moisy2015)), rotating turbulence is observed to develop a strong anisotropy and to spontaneously form vortices that are invariant along the rotation axis. The latter correspond to a first-order balance between the Coriolis force and pressure gradients and are called ‘geostrophic modes’. Yet, the systematic observation of strong geostrophic modes is at odds with various evidence suggesting that rotating turbulence could also be dominated by inertial waves that are sustained by the restoring action of the Coriolis force. Recent numerical (Le Reun *etal.* Reference Le Reun, Favier, Barker and Le Bars2017) and experimental (Le Reun, Favier & Le Bars Reference Le Reun, Favier and Le Bars2019; Brunet, Gallet & Cortet Reference Brunet, Gallet and Cortet2020) studies have shown that injecting energy in waves solely creates a turbulent state comprising of inertial waves only when the forcing amplitude is sufficiently small, i.e. a discrete version of inertial wave turbulence (Galtier Reference Galtier2003). It is only at larger forcing amplitudes that a secondary instability leads to the classical geostrophic-dominated turbulence. Asymptotic theories describing rotating turbulence in the limit of vanishing forcing amplitude and dissipation also suggest that waves could dominate over geostrophic modes in such a regime (Bellet *etal.* Reference Bellet, Godeferd, Scott and Cambon2006; Sagaut & Cambon Reference Sagaut and Cambon2018). Hence, although bi-dimensionalisation in the form of geostrophic eddies has been commonly observed, it may not be the only equilibrium state of rotating turbulence, be it at moderate (Yokoyama & Takaoka Reference Yokoyama and Takaoka2017; Favier, Guervilly & Knobloch Reference Favier, Guervilly and Knobloch2019) as well as small (van Kan & Alexakis Reference van Kan and Alexakis2019) forcing amplitudes. In addition, the nature of the forcing seems fundamental in determining the equilibrium state of rotating turbulence.

These results altogether call for a better understanding of the fundamental processes by which waves give rise to balanced geostrophic modes. The studies of Le Reun *etal.* (Reference Le Reun, Favier and Le Bars2019) and Brunet *etal.* (Reference Brunet, Gallet and Cortet2020) suggest that such a transfer occurs through an instability. Although wave-to-wave interactions are primarily governed by triadic resonance (Vanneste Reference Vanneste and Grimshaw2005; Bordes *etal.* Reference Bordes, Moisy, Dauxois and Cortet2012), they cannot account for wave-to-geostrophic transfers (Greenspan Reference Greenspan1969), at least in the asymptotic limit of vanishing velocity amplitude and dissipation. Several alternative mechanisms, outside the framework of Greenspan's theorem, have been proposed. Four-mode interactions can transfer energy from waves to geostrophic flows, either directly (Newell Reference Newell1969; Smith & Waleffe Reference Smith and Waleffe1999) or through an instability mechanism (Kerswell Reference Kerswell1999; Brunet *etal.* Reference Brunet, Gallet and Cortet2020). The growth rate of such an instability scales as $Ro^2$, with $Ro$ the dimensionless wave amplitude or Rossby number. It develops over longer time scales than triad-type interactions between waves. The other inviscid mechanism that has been proposed to account for wave-geostrophic transfer is quasi-resonant triadic interaction (Newell Reference Newell1969; Smith & Waleffe Reference Smith and Waleffe1999), that is, a triad between waves whose frequencies do not exactly satisfy the resonance condition (Bretherton Reference Bretherton1964). Their presence and their role in the bi-dimensionalisation of rotating turbulence have been assessed by several numerical studies (Smith & Lee Reference Smith and Lee2005; Alexakis Reference Alexakis2015; Clark di Leoni & Mininni Reference Clark di Leoni and Mininni2016). While it has been shown that such triads can directly transfer energy from two pre-existing waves to geostrophic modes, we show that this transfer can arise spontaneously through an instability mechanism. More precisely, we show with direct numerical simulations (DNS) and theoretical analysis that there exists a linear mechanism by which a single inertial wave drives exponential growth of geostrophic modes through near-resonant triadic interaction.

## 2. The stability of a single inertial wave

### 2.1. Governing equations and numerical methods

Let us consider an incompressible fluid rotating at rate $\Omega \boldsymbol {e}_z$. We investigate the stability of a single plane inertial wave with wavevector $\boldsymbol {k}$ and eigenfrequency $\omega_{k}^{s_k}$. Its amplitude is proportional to $\boldsymbol {U}_w$ with

Here, $\boldsymbol {h}_{\boldsymbol {k}}^{s_k}$ is a helical mode, that is, when $\boldsymbol {k}$ is not parallel to the axis of rotation $\boldsymbol {e}_z$ (Cambon & Jacquin Reference Cambon and Jacquin1989; Waleffe Reference Waleffe1992)

where $s_k = \pm 1$ is the sign of the helicity of the plane wave. If $\boldsymbol {k}$ is parallel to $\boldsymbol {e}_z$, $\boldsymbol{h}_{\boldsymbol{k}}^{s_k} = \boldsymbol {e}_x + \textrm {i}\, s_k \boldsymbol {e}_y$. Here, $\boldsymbol {U}_w$ automatically satisfies the incompressibility condition since $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {U}_w \propto \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {h}_{\boldsymbol {k}}^{s_k} = 0$ and $\boldsymbol {U}_w$ satisfies the linearised rotating Euler equation,

provided that $\omega_{k}^{s_k}$ and $\boldsymbol {k}$ are related by the dispersion relation of inertial waves

$\theta$ being the angle between the wavevector $\boldsymbol {k}$ and the rotation axis $\boldsymbol {e}_z$ (ranging from $0$ to ${\rm \pi}$), and $k = \vert \boldsymbol {k} \vert$. We solve for the time evolution of perturbations $\boldsymbol {u}$ to the wave $\boldsymbol {U}_w$ maintained at a constant amplitude via the following set of equations:

*a*)\begin{gather} \partial_t \boldsymbol{u} + Ro \left(\boldsymbol{U}_w \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} +\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{U}_w\right) + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + 2 \boldsymbol{e}_z \times \boldsymbol{u} ={-} \boldsymbol{\nabla} {\rm \pi} + E \nabla^2 \boldsymbol{u}, \end{gather}

*b*)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

where time is scaled by $\Omega ^{-1}$ and length by the domain size $L$. We have introduced the Ekman number $E = \nu /(L^2 \Omega )$, $\nu$ being the kinematic viscosity, and an input Rossby number $Ro$ quantifying the dimensionless amplitude of the plane wave.

Equations (2.5) are solved numerically in a triply periodic cubic box using the code ‘Snoopy’ (Lesur & Longaretti Reference Lesur and Longaretti2005). The dynamics of the perturbation flow $\{ \boldsymbol {u}, {\rm \pi} \}$ is determined with pseudospectral methods, that is, $\{ \boldsymbol {u}, {\rm \pi} \}$ is decomposed into a truncated sum of Fourier modes $\{ \hat {\boldsymbol {u}}_{\boldsymbol {q}},\hat {{\rm \pi} }_{\boldsymbol {q}} \} \textrm {e}^{ \textrm {i}\, \boldsymbol {q} \boldsymbol {\cdot } \boldsymbol {x}}$. A wavevector $\boldsymbol {q}$ writes $2 {\rm \pi} /L ( n_x, n_y, n_z )$ where $n_{x,y}$ are integers varying from $-N$ to $N$ and $n_z$ is an integer ranging from $0$ to $N$ because of Hermitian symmetry. In the following, $N = 96$; higher resolutions have been tested and yield the exact same results. The temporal dynamics of the modes $\hat {\boldsymbol {u}}_{\boldsymbol {q}}$ is solved using a third-order Runge–Kutta method. Note that the size of the box $L$ is artificial and we thus expect our results to depend on the intrinsic Rossby number based on the imposed wavelength $k Ro$ rather than $Ro$ alone.

### 2.2. Numerical results

Keeping the Ekman number to $E = 10^{-6}$, two simulations of the stability of the inertial wave $\boldsymbol {k} = 2 {\rm \pi} [ 4,0,8 ]$ (with $s_k = 1$) are performed at low ($Ro = 2.83 \times 10^{-3}$) and moderate ($Ro = 2.83 \times 10^{-2}$) wave amplitudes. They are both initiated with a random noise comprising wavenumbers ranging between $0$ and $15 {\rm \pi}$, with very small initial amplitude. The use of spectral methods allows the separation of kinetic energy of the perturbation flow $\boldsymbol {u}$ into a two-dimensional component $\mathcal{E}_{\text{2-D}}$, accounting for all modes $\boldsymbol {q}$ with $q_z=0$, and a complementary three-dimensional component $\mathcal {E}_{\text{3-D}}$. In addition, performing a Fourier transform in space and time allows the projection of the kinetic energy in the sub-space of the dispersion relation of inertial waves to draw the spatiotemporal spectrum $\mathcal {E}(\theta , \omega )$ (Yarom & Sharon Reference Yarom and Sharon2014; Le Reun *etal.* Reference Le Reun, Favier, Barker and Le Bars2017). Note that, in these maps, $\theta$ is restricted to $[0, {\rm \pi} /2]$ since the flow is real and only the wavevectors $\boldsymbol {q}$ with $q_z \geq 0$ are simulated due to the Hermitian symmetry. Moreover, the spectra are symmetrised with respect to $\omega =0$ and the maps are shown as a function of $\vert \omega \vert$ to be more compact.

The kinetic energy time series and maps are shown for both simulations in figure 1. At low imposed wave amplitude (figure 1*a*,*b*), three-dimensional perturbations dominate the growth of the instability. The energy map displays two spots aligned on the dispersion relation with negative frequencies $\omega _{1,2}$ such that $\omega_{k}^{s_k} + \omega _1 +\omega _2 = 0$ which is indicative of several waves undergoing triadic resonance with the imposed mode. The growth of two-dimensional modes is delayed and their growth rate is approximately two times larger than the rate of three-dimensional modes. Removing the nonlinear term $\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u}$ (see (2.5)) suppresses the growth of two-dimensional modes which are thus not unstable themselves, at least on the time scale of the growth of waves. In fact, we find that $\mathcal{E}_{\text{2-D}} \simeq 3 \times 10^{-2} (k \mathcal {E}_{\text{3-D}})^2$ (see figure 1*a*) in the growth phase, which suggests that two-dimensional modes’ growth is due to direct forcing by nonlinear interaction of two growing waves involved in triadic resonances, with close frequencies and opposed vertical wavenumbers. This mechanism corresponds to the direct excitation of geostrophic modes by two waves identified by Newell (Reference Newell1969) and Smith & Waleffe (Reference Smith and Waleffe1999).

At larger wave amplitude (figure 1*c*,*d*), this picture changes: two- and three-dimensional modes grow at the same rate from the start of the simulation. The vertical vorticity field of the growing perturbation and its vertical average are displayed in figure 2(*a*,*b*). The frequency of two-dimensional modes ($\theta ={\rm \pi} /2$) is close to $0$ according to figure 1(*d*), they are thus geostrophic, that is, slow as well as invariant along the $z$ axis. As it may be seen in figure2(*b*), the geostrophic unstable flow contains the wavevector $\boldsymbol {p} = 2 {\rm \pi} [0,5,0 ]$ (and $-\boldsymbol {p}$ for hermiticity), which grows along with other three-dimensional structures whose energy location in the $(\theta , \omega )$ space is reminiscent of triadic resonant instability, but with significantly more spreading. Note that the properties of the transition and the instability we report are robust to changes in the aspect ratio of the box, which discards any spurious effect of the discretisation (Smith & Lee Reference Smith and Lee2005).

To further characterise the growth of geostrophic modes, we carry out simulations with the same imposed wave, but we use as initial condition only one wavevector, $\boldsymbol {p} = 2 {\rm \pi} [0,p_{y},0 ]$ with $p_{y} \in \{ 1,3,5 \}$ and $s_{p} = \pm 1$, instead of random noise. The Ekman number is set to $10^{-8}$ to discard effects of viscosity in the mode growth. The growth rate of geostrophic modes is reported in figure 2(*c*) where $Ro$ is systematically varied. At large wave amplitude, the growth rate increases linearly, but it goes to zero at a finite value of $Ro \sim 10^{-2}$. This value is very large compared with the viscous damping $k^2 E \sim 10^{-5}$ which plays no role here, hence suggesting an inviscid mechanism.

On the kinetic energy map $\mathcal {E}(\theta ,\omega )$ (figure 2*d*) two particular spots appear, one at the location of geostrophic modes and the other at the angle of the vectors closing the triads between $\boldsymbol {k}$ and $\pm \boldsymbol {p}$, which is indicative of triad-type interactions. This may seem *a priori* in contradiction with the result of Greenspan (Reference Greenspan1969). Yet, our results are outside Greenspan's theorem framework in two aspects. First, our study is necessarily limited to finite Rossby number. Second, the frequency of the closing mode $-\boldsymbol {k} \pm \boldsymbol {p}$ (figure 2*d*) is outside its eigenfrequency. In fact, the discrepancy between the observed frequencies and eigenfrequencies of the closing mode corresponds to the sum of the eigenfrequencies of the three modes, ${\rm \Delta} \omega _{kpq}$. Moreover, as $Ro$ is decreased, the fastest growing mode $\boldsymbol {p}_g$ in figure 2(*c*) goes from $p_y = 5$ with ${\rm \Delta} \omega _{kpq}=0.22$ to $p_y = 1$ with ${\rm \Delta} \omega _{kpq}=0.01$: the triad is nearly resonant and draws closer to exact resonance. Two mechanisms may be considered to explain the growth of geostrophic modes. One hypothesis is that the growth of the geostrophic modes is a by-product of a classical resonance between, say, $\boldsymbol {k}$ and $-\boldsymbol {q}$ that forces $\boldsymbol {q}$ (by Hermitian symmetry) and hence $\boldsymbol {p}$. However, in this case, the mode closing the wave triad would have an eigenfrequency $-\omega_{k}^{s_k} + \omega_{q}^{s_q} = {\rm \Delta} \omega_{kqp} -2 \omega_{k}^{s_k}$ which is outside the range of inertial waves for small ${\rm \Delta} \omega _{kpq}$ since $\omega_{k}^{s_k} \simeq 1.78$. Instead, our results point towards a near-resonant triadic instability that transfers energy from an inertial wave to a $z$-invariant geostrophic mode.

## 3. Near resonance of geostrophic modes: theoretical approach

### 3.1. The low Rossby number limit

To provide theoretical insight into the inertial-wave destabilisation observed in the numerical simulations, we turn to linear stability analysis. The flow is $\boldsymbol {U} = \boldsymbol {u} + Ro\, \boldsymbol {U}_w$ where $Ro \,\boldsymbol {U}_w$ is the base flow ($Ro$ being the finite Rossby number) and $\boldsymbol {u} \ll Ro\, \boldsymbol {U}_w$ is now an infinitesimal perturbation. We consider a Cartesian domain with periodic boundary conditions and we decompose $\boldsymbol {U}$ into a superposition of plane waves $\boldsymbol {h}_{p}^{s_p} \exp \textrm {i} (\,\boldsymbol{p} \boldsymbol {\cdot } \boldsymbol {x}-\omega _{p}^{s_p} t)$ with time-dependent amplitudes $b_{p}^{s_p} (t)$. The Euler equation governing $\boldsymbol {U}$ is then equivalent to a set of ordinary differential equations governing the amplitudes $b_{p}^{s_p}$ (Smith & Waleffe Reference Smith and Waleffe1999)

with

*a*,

*b*)\begin{align} C_{pkq}^{s_p s_k s_q} \equiv \frac{1}{2} (s_{q} q - s_k k) \, \boldsymbol{h}_{p}^{s_p*} \boldsymbol{\cdot} \left(\boldsymbol{h}_{k}^{s_k *} \times \boldsymbol{h}_{q}^{s_q *} \right) \quad \mbox{and} \quad {\rm \Delta} \omega_{k p q} \equiv \omega_{k}^{s_k} + \omega_q^{s_q} + \omega_p^{s_p}, \end{align}

and where $k,p = \vert \boldsymbol {k} \vert , \vert \boldsymbol {p} \vert$. Here, ${\rm \Delta} \omega _{k p q}$ is the sum of the eigenfrequencies of the three modes involved in the triad. In general, maximum energy transfer between the three modes is ensured when the oscillations due to the detuning in the right-hand side of (3.1) are cancelled, that is, when ${\rm \Delta} \omega _{k p q} \rightarrow 0$. If all three modes $\boldsymbol {k}$, $\boldsymbol {p}$ and $\boldsymbol {q}$ are inertial waves, this leads to the well known mechanism of triadic resonance (Bretherton Reference Bretherton1964; Vanneste Reference Vanneste and Grimshaw2005; Bordes *etal.* Reference Bordes, Moisy, Dauxois and Cortet2012). When the mode $\boldsymbol {k}$ is imposed with an amplitude $Ro$ and helicity $s_k$, $\boldsymbol {p}$ and $\boldsymbol {q}$ grow exponentially with a rate proportional to $k Ro$. This picture is changed when one of the modes, say $\boldsymbol {p}$, is geostrophic (i.e. $\omega _p^{s_p} = 0$ and $p_z = 0$, regardless of $s_p$). The spatial interaction condition $\boldsymbol {k} + \boldsymbol {p} + \boldsymbol {q} = \boldsymbol {0}$ forces $k_z = - q_z$ and the resonance condition becomes

It may be fulfilled only when $s_q = s_k$ and $q - k \rightarrow 0$. At exact resonance, these conditions impose the location of $\boldsymbol {p}$ on a circle centred on $-\boldsymbol {k}_{\perp } = -[k_x,k_y ]$ with radius $\vert \boldsymbol {k}_{\perp }\vert$. Moreover, in the governing equation for $\dot {b}_p^{s_p}$, the slow oscillation terms involve a coupling coefficient $C_{p k q}^{s_p s_k s_k} \propto k - q \rightarrow 0$. At exact resonance, the coupling coefficient vanishes: there is no energy transfer from waves to geostrophic modes, as proved by Greenspan (Reference Greenspan1969).

Nevertheless, wave-to-geostrophic transfer is still possible when the detuning ${\rm \Delta} \omega _{kpq}$ is small but non-zero (Newell Reference Newell1969; Smith & Waleffe Reference Smith and Waleffe1999; Alexakis Reference Alexakis2015), and we investigate instabilities of geostrophic modes driven by this mechanism. Let us assume that the wave $\boldsymbol {k}$ with helicity $s_k$ is imposed with a small constant amplitude $Ro$, and that $\boldsymbol {p}$ is geostrophic. To infer from (3.1) the time evolution of the geostrophic mode amplitude, we proceed to an asymptotic expansion using a two-time method involving a fast time $\tau = t$ and a slow time $T$. The hierarchy between them must be a power of $kRo$, the intrinsic Rossby number based on the imposed wavelength. Because the wave-to-geostrophic transfer coefficient $C_{ p k q }^{s_p s_k s_k}$ vanishes as ${\rm \Delta} \omega _{kpq} \rightarrow 0$, we find via a heuristic analysis, detailed in appendix A, that $T = (kRo)^2 t$, instead of $(k Ro) \,t$ for classical wave triads. In addition, it imposes the condition that the amplitude of $\boldsymbol {p}$ be smaller by a factor $kRo$ compared with the mode closing the triad $\boldsymbol {q} = -\boldsymbol {k} -\boldsymbol {p}$. The amplitudes of the modes interacting with the imposed wave ($\,\boldsymbol {p}$ and $\boldsymbol {q}$ with both helicity signs $s_{p,q}$) are thus expanded as

where the $B_i^j$ are all $O(1)$. The hierarchy between orders is imposed by the need to match the slow time derivative of the leading order with the next order in the multiple scale expansion. Since the slow derivation introduces a factor $(kRo)^2$, there must be a $(kRo)^2$ hierarchy between the first and second orders.

To find the equations governing the leading-order coefficients, we follow the method of Bretherton (Reference Bretherton1964) and inject the ansatz (3.4) into (3.1). First, the leading-order coefficients are found to be independent of $\tau$ and their long-time evolution with $T$ is determined at next order. As noted by Bretherton (Reference Bretherton1964), for fast triads with detuning larger than $O((kRo)^2)$, the imposed wave only drives fast and bounded oscillations of the second-order terms, and the leading-order terms must be zero to avoid secular growth. It is only for slow triads, i.e. when ${\rm \Delta} \omega _{kpq} = O((kRo)^2)$, that the exponential term in (3.1) drives slow oscillations that contribute to secular growth. Such a condition on the detuning is consistent with the numerics: the fastest growing modes at $Ro = 1 \times 10^{-2}$ ($kRo \sim 0.6$) in figure 2(*c*) is $\boldsymbol {p} = 2{\rm \pi} [0,3,0]$ for which ${\rm \Delta} \omega _{kpq} \simeq 0.1 \sim 0.3 (k Ro)^2$. The condition ${\rm \Delta} \omega _{kpq} = O((kRo)^2)$ may be fulfilled only for the modes $\boldsymbol {p}$ with $s_p = \pm 1$ and $\boldsymbol {q}$ with $s_q = s_k$ and, as noted earlier, the wave-to-geostrophic transfer coefficient $C_{ p k q }^{s_p s_k s_k}$ is also $O((kRo)^2)$. Cancelling the secular growth terms at second order gives the following amplitude equations :

*a*,

*b*)\begin{equation} \partial_T B_{q 1}^{s_q} = \sum_{s_p} C_{qkp}^{s_k s_k s_p}\, B_{p1}^{s_p \,*}\, \textrm{e}^{\textrm{i}\, \frac{{\rm \Delta} \omega_{kpq}}{(kRo)^2} T }\quad \mbox{and}\quad \partial_T B_{p 1}^{s_p} = \frac{C_{pkq}^{s_p s_k s_k}}{(kRo)^2} \, B_{q1}^{s_q\,*}\, \textrm{e}^{\textrm{i}\, \frac{{\rm \Delta} \omega_{kpq}}{(kRo)^2} T }, \end{equation}

the rescaled quantities $C_{ p k q }^{s_p s_k s_k}/(kRo)^2$ and ${\rm \Delta} \omega _{kpq}/(kRo)^2$ being $O(1)$. These equations have exponentially growing solutions and the complex growth rate of the instability normalised by the rotation rate is

Note that the product of coupling coefficients $C_{pkq}^{s_p s_k s_k} C_{qkp}^{s_k s_k s_p\, *}$ is real. The expression of the real part of growth rate, $\sigma _k(\,\boldsymbol{p};Ro)$, is consistent with the numerical findings: for a given near-resonant triad, it drops to zero at a finite value of $Ro$ and it is proportional to $Ro$ at larger $Ro$. However, in the small Rossby number limit, because the detuning ${\rm \Delta} \omega _{kpq}$ and the coupling coefficient $C_{ p k q }^{s_p s_k s_k}$ are both $O((kRo)^2)$, the maximum geostrophic growth rate remains $O((kRo)^2)$ at most.

The $(kRo)^2$ scaling governing the maximum growth rate of geostrophic near resonance is found quantitatively by expanding the frequency detuning, the transfer coefficients and then $\sigma _k(\,\boldsymbol{p};Ro)$ in the neighbourhood of exact resonance. The fact that the growth rate is non-zero only close to the exact resonance is illustrated qualitatively in figure 3(*a*). We thus introduce $\boldsymbol {p} = \boldsymbol {p}_0 + (kRo)^2 \boldsymbol {\delta } \boldsymbol {p}$ where $\boldsymbol {\delta }\boldsymbol {p}$ is an $O(1)$ vector in the geostrophic plane. Consider the closing modes $\boldsymbol {q}_0 = -(\boldsymbol {k} +\boldsymbol {p}_0)$ and $\boldsymbol {q} = -(\boldsymbol {k} +\boldsymbol {p})$, then $\boldsymbol {q} = \boldsymbol {q}_0 - (kRo)^2 \boldsymbol {\delta } \boldsymbol {p}$. The imposed mode $\boldsymbol {k}$ is left unperturbed. Since ${\rm \Delta} \omega _{kpq} =C_{p_0kq_0}^{s_p s_k s_k} = 0$ at exact resonance, the leading orders of the frequency detuning and the transfer coefficients are found from (3.2*a*,*b*) and (3.3) to be $O((kRo)^2)$. The perturbation of the wavevectors is thus consistent with the asymptotic expansion. At leading order,

*a*,

*b*)\begin{equation} \frac{{\rm \Delta} \omega_{kpq} }{(kRo)^2} \simeq s_k \omega_{k}^{s_k} \frac{\boldsymbol{q}_0 \boldsymbol{\cdot}\boldsymbol{\delta} \boldsymbol{p} }{k^2} \quad\mbox{and} \quad \frac{C_{pkq}^{s_p s_k s_k}}{(kRo)^2} \simeq \frac{1}{2} s_k \left(\boldsymbol{h}_{p_0}^{s_p\,*} \boldsymbol{\cdot} \left( \boldsymbol{h}_k^{s_k\,*} \times \boldsymbol{h}_{q_0}^{s_k\,*} \right)\right) \frac{\boldsymbol{q}_0 \boldsymbol{\cdot}\boldsymbol{\delta} \boldsymbol{p}}{k}, \end{equation}

where we have used that $q_0 = k$ and $s_q = s_k$. In the growth rate $\sigma _k(\,\boldsymbol{p};Ro)$, the product of the coupling coefficients is

Therefore, at leading order in powers of $Ro$, the growth rate is

where $X \equiv (\boldsymbol {q}_0 \boldsymbol {\cdot }\boldsymbol {\delta } \boldsymbol {p})/k^2$ and $\mathcal {C}_k(\,\boldsymbol{p}_0)$ is the sum in the right-hand side of (3.8). When $\mathcal {C}_k(\,\boldsymbol{p}_0) >0$, the growth rate reaches an optimum at $X = \mathcal {C}_k /(2 {\omega_{k}^{s_k}}^2)$ with value $(kRo)^2 \mathcal {C}_k (\,\boldsymbol{p}_0) /\vert 4 \omega_{k}^{s_k}\vert$, which remains to be maximised over all exactly resonant wavevectors $\boldsymbol {p}_0$. The coefficient $\mathcal {C}_k(\,\boldsymbol{p}_0)$ is shown in figure 3(*b*) for several wavevectors $\boldsymbol {k}$ with different frequencies $\omega_{k}^{s_k}$. When plotted against $p_0/k_{\perp }$, all the curves $\mathcal {C}_k$ collapse on a master curve that reaches a maximum value of $1$ at $p_0 \rightarrow 0$, regardless of the helicity sign $s_k$ of the imposed wave. As $(kRo) \rightarrow 0$, the unstable geostrophic mode becomes large-scale and stems from the interaction between $\boldsymbol {k}$ and $\boldsymbol {q} \simeq -\boldsymbol {k}$. In the low Rossby number limit, the maximal growth rate is then

We confirm that $k Ro$, the Rossby number based on the wavelength, is the relevant parameter to describe the growth rate of the instability. While each geostrophic mode follows the law (3.6), all the growth rate curves lie below a $(k Ro)^2$ upper envelope. More details are given below in § 3.3 where we compare the law (3.10) to the exact computation of the growth rate curves from (3.6).

### 3.2. The moderate to large Rossby number regime

The $Ro^2$ law governing the growth rate in the small Rossby number limit cannot hold as the wave amplitude is increased since the growth rate has an upper bound following an $Ro$ law. This is proved directly by multiplying (2.5) by $\boldsymbol {u}$ and integrating over the fluid domain $V$, thus giving

where $\left \lVert \boldsymbol {\cdot }\right \rVert _n$ denotes the $L_n$-norm. In virtue of Hölder's inequality (Gallet Reference Gallet2015),

This upper bound applies in particular to the low $Ro$ scaling (3.10), which thus holds up to $k Ro \sim 1$ at most. Note that, beyond that point, non-triad type instabilities (shear, centrifugal) may add to near resonance in driving the dynamics of the flow excited by the maintained wave.

### 3.3. Comparison with exact computation

In figure 3, we sample the growth rate curves given by (3.6) of many geostrophic modes in near-resonant interaction with $\boldsymbol {k} = 2 {\rm \pi} [4,0,8]$, as functions of the Rossby number. We notice that each growth rate curve follows an $Ro$ scaling at sufficiently large $Ro$, which is consistent with the asymptotic expansion carried out in § 3.1. However, because for all modes the growth rate vanishes at a finite value of $Ro$ (which decreases to 0 close to exact resonance), the upper envelope is an $Ro^2$ law that perfectly matches the theoretical prediction (3.10). This instability is thus fundamentally different from four-modes interactions for which each growth rate curve follows a $Ro^2$ law (Kerswell Reference Kerswell1999; Brunet *etal.* Reference Brunet, Gallet and Cortet2020). Besides, in agreement with the upper bound (3.12), we observe the $Ro^2$ maximum growth rate law to break down above $k Ro \sim 1$. Beyond that, the near-resonant growth rate derived from (3.6) remains below the upper bound and follows an $Ro$ law. However, additional instabilities (shear, centrifugal, etc.) may also drive the growth of geostrophic modes in this regime.

### 3.4. Finite size and viscous effects

Finite size domain translates into discretisation of the modes. Since the exactly resonant geostrophic modes lie on a finite radius circle, they cannot be approached with arbitrarily low detuning ${\rm \Delta} \omega _{kpq}$ as $Ro \rightarrow 0$. Hence, discretisation implies the existence of a finite value of $Ro$ below which the near-resonant instability vanishes (Bourouiba Reference Bourouiba2008). Regarding viscous effects, at finite but low Ekman number, i.e. when $(k Ro)^2 \gg k^2 E$, the low Rossby number scaling is unaltered since the near-resonant modes $\boldsymbol {p}$ and $\boldsymbol {q}$ have at most similar wavenumbers to $\boldsymbol {k}$. The $O(Ro)$ upper bound on the growth rate is also unaltered by the inclusion of viscosity.

## 4. A refined model: the double near-resonant triad

Although promising, the model of the previous section needs to be refined to fully account for numerical results. It misses by a factor two the growth rate curves of figure2(*c*) and (3.6) predicts a frequency ${\rm \Delta} \omega _{kpq} /2 \neq 0$ for the geostrophic modes, while it is $0$ according to figure 2(*d*). We must include in our model that not only $\boldsymbol {k}$, but also $- \boldsymbol {k}$ (with the same helicity sign $s$) is imposed due to Hermitian symmetry. Consider two triads $\boldsymbol {k} + \boldsymbol {p} + \boldsymbol {q} = \boldsymbol {0}$ and $-\boldsymbol {k} + \boldsymbol {p} + \boldsymbol {q}' = \boldsymbol {0}$ (see figure 4) with helicity signs $s_q = s_{q'} = s_k$ and $s_p = \pm 1$. As shown schematically in figure 4(*a*), the exact resonance circles associated with $\pm \boldsymbol {k}$ coincide at $\boldsymbol {p} = 0$. In a neighbourhood of this point (the dark grey intersection in figure 4*a*), a geostrophic mode $\boldsymbol {p}$ is in near resonance with both $\pm \boldsymbol {k}$ with both detunings ${\rm \Delta} \omega _{kpq}$ and ${\rm \Delta} \omega _{-kpq'}$ small. It is all the more important that $\boldsymbol {p} = \boldsymbol {0}$ corresponds to the largest near-resonant growth rate as $Ro \rightarrow 0$.

In a very similar way to the two-time asymptotic expansion of § 3.1, we can retrieve the amplitude equations governing the slowly varying envelopes $B_{p1}^{s_p}$ and $B_{Q1}^s$, with $\boldsymbol {K}= \pm \boldsymbol {k}$, $\boldsymbol {Q} = -\boldsymbol {K} - \boldsymbol {p}$ (i.e. $\boldsymbol{Q} = \boldsymbol{q}$ or $\boldsymbol{q}'$) and $s_p = \pm 1$. Similarly to (3.5*a*,*b*), the four amplitude equations are

The characteristic polynomial $\mathcal {P}(\sigma )$ of this set of four linear differential equations in terms of the Rossby number $Ro$ and the detunings is

with $S_{K}$ and $P_0$ two real coefficients defined as

*a*,

*b*)\begin{align} S_K = \sum_{s_p} C_{pKQ}^{s_p s_k s_k} C_{QKp}^{s_k s_k s_p\, *} \quad\mbox{and}\quad P_0 = \left(\sum_{s_p} C_{q k p}^{ s_k s_k s_p} C_{p -k q'}^{ s_p s_k s_k} \right) \left(\sum_{s_p} C_{q' -k p}^{ s_k s_k s_p} C_{p k q}^{ s_p s_{{+}k} s_k}\right). \end{align}

In general, the growth rate $\sigma _k(\,\boldsymbol{p};Ro)$ is found by numerical computation of the roots of $\mathcal {P}$. Nevertheless, an estimate of the maximum growth rate can be obtained in the limit $\vert p_x \vert \ll \vert p_y \vert$, which is relevant to our DNS. In this case, symmetries impose ${\rm \Delta} \omega _{kpq} = -{\rm \Delta} \omega _{-kpq'}$, $S_k = S_{-k}$ and $P_0 = S_{k}^2$. The polynomial $\mathcal {P}$ then has a purely real root,

and the frequency of the growing geostrophic mode is $0$. We show in figure 4(*b*) the excellent agreement between this theoretical law and the DNS data of figure 2(*a*). An expansion similar to § 3.1 reveals that the maximum growth rate follows the exact same $(kRo)^2$ law as in the case of the single triad (3.10). This is confirmed by the systematic computation of the growth rate in the double triad case, as shown in figure4(*c*). Note that although this instability is driven by two imposed modes, it remains different from the four-mode interaction mechanism detailed in Brunet *etal.* (Reference Brunet, Gallet and Cortet2020). The latter features intermediate, non-resonant modes that are absent in the double triad mechanism. Moreover, the near-resonant growth rate (4.4) is proportional to $Ro$ when $Ro$ is sufficiently large, whereas it always follows an $Ro^2$ law in the case of four-mode interaction.

We compare in figure 4(*c*) our theoretical predictions with the geostrophic growth rate extracted from DNS initiated with a large-scale noise, as in § 2, down to $Ro = 5 \times 10^{-3}$. We use the same imposed wavevector as previously ($\boldsymbol {k} = 2 {\rm \pi} [4,0,8 ]$) but also 1.5 and 2 times longer wavevectors to confirm that the growth rate is a function of the intrinsic Rossby number $k Ro$. Such a dependence is expected since in the infinite-domain limit, $L \rightarrow \infty$, $L$ becomes irrelevant and $1/k$ is the only remaining length scale, and the associated Rossby number is $k Ro$. As shown in figure 4(*c*), the numerical growth rate coincides with the maximum near-resonant growth rate in the low $Ro$ regime but also in the moderate $Ro$ regime where it is proportional to $kRo$. Note that for the three lowest $k Ro$ points, the noise has been implemented on two-dimensional modes ($\,p_z = 0$) to facilitate the isolation of the geostrophic instability. It delays the growth of two-dimensional modes with non-zero frequency by direct forcing at the lowest values of $Ro$. Despite their rapid growth, the latter are subdominant in the saturation of wave-driven flows, and do not prevent the long-term growth of unstable geostrophic modes under the mechanism examined here.

Lastly, our analysis allows us to understand the transition in the stability of the geostrophic flow observed in the numerical study (see figure 1). Here, $Ro = 2.83 \times 10^{-2}$ ($k Ro \simeq 1.6$) corresponds to the transition zone where the geostrophic growth rate is $O(kRo)$, as the wave-only triadic resonances (see figure 4*c*). When the Rossby number is decreased by one order of magnitude, as shown in figure 4(*c*), the growth rate of geostrophic modes scales as $(kRo)^2$, i.e. a factor $k Ro$ smaller than the growth rate of wave-only triadic resonances. This is why the geostrophic instability is not observed in the simulation at $Ro = 2.83 \times 10^{-3}$.

## 5. Conclusion

By means of numerical simulations and theoretical analysis, we have described a new instability mechanism by which inertial waves excite $z$-invariant geostrophic modes. We have proved that this instability is driven by near-resonant triadic interaction and derived its theoretical growth rate. When normalised by the global rotation rate $\Omega$, the growth rate follows a $(kRo)^2$ law and a $kRo$ law at small and moderate wave amplitude, respectively, $k$ being the imposed wavenumber. It translates into $\Omega ^{-1}$ and $\Omega ^0$ laws, respectively, for the dimensional growth rate. The near-resonant instability completes the picture proposed by Brunet *etal.* (Reference Brunet, Gallet and Cortet2020) where another inviscid geostrophic instability based on two imposed modes and four-mode interaction is detailed. Although of different nature, both instabilities have a $(kRo)^2$ growth rate in the limit of small Rossby numbers which makes them possibly difficult to distinguish. On the one hand, in the linear growth phase, the eigenmode of the four-mode instability consists of waves and geostrophic flow of comparable amplitude, whereas the geostrophic flow is $kRo$ times smaller than the waves in the present mechanism. On the other hand, the near-resonant instability achieves a growth rate of order $kRo$ for $kRo \gtrsim 1$. Our analysis may thus explain the experimental results of Le Reun *etal.* (Reference Le Reun, Favier and Le Bars2019): they found at moderate Rossby number a geostrophic instability with a growth rate proportional to $Ro$, which thus matches the moderate Rossby number law derived here.

Our work paves the way for new studies dealing with the energy transfers from waves to geostrophic modes in rotating fluids, in particular in rotating turbulence. Firstly, it remains to be seen under which conditions the near-resonant instability may excite slow nearly geostrophic modes with small but non-zero frequencies. These modes have been proved to be important for the development of anisotropy in rotating turbulence (Smith & Lee Reference Smith and Lee2005) even in the asymptotic limit of small Rossby numbers (van Kan & Alexakis Reference van Kan and Alexakis2019). With a heuristic analysis similar to the one developed in § 3.1 and appendix A, we can predict again a $(kRo)^2$ growth rate, but a detailed analysis is difficult and remains to be done. Lastly, the instability we describe in the present article could play an important role in the rotating turbulence dynamics that remains to be fully deciphered. Since wave–wave triadic interactions grow at rate $O(kRo)$ while wave–geostrophic interactions have a $O((kRo)^2)$ growth rate, we may speculate that, at sufficiently low forcing amplitudes, rotating turbulence may remain purely three-dimensional as the wave amplitudes would remain below the threshold of geostrophic instabilities. Such states have been observed by Le Reun *etal.* (Reference Le Reun, Favier, Barker and Le Bars2017), Le Reun *etal.* (Reference Le Reun, Favier and Le Bars2019) and Brunet *etal.* (Reference Brunet, Gallet and Cortet2020), but a systematic investigation of their existence in various realisations of rotating turbulence remains to be carried out. In the future, our study could help bridge the gap between finite Rossby number experiments and simulations of turbulence (Godeferd & Moisy Reference Godeferd and Moisy2015) and asymptotic models of rotating turbulence (Galtier Reference Galtier2003; Bellet *etal.* Reference Bellet, Godeferd, Scott and Cambon2006; van Kan & Alexakis Reference van Kan and Alexakis2019).

## Acknowledgements

The authors acknowledge funding by the European Research Council under the European Union's Horizon 2020 research and innovation programme through grant no. 681835-FLUDYCO-ERC-2015-CoG and FLAVE 757239. T.L.R. is supported by the Royal Society through a Newton International Fellowship (grant reference NIF∖R1∖192181). The computations were performed using the HPC resources of GENCI-IDRIS (grant no. A0080407543).

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. Heuristic justification of the $(kRo)^{\textbf{2}}$ slow time scale

To investigate near-resonant instability, we proceed to an asymptotic expansion of the spectral version of the Euler equation using a two-timing method involving a fast time $\tau = t$ and a slow time $T$. The hierarchy between $\tau$ and $T$ may be derived heuristically from an analysis of the amplitude equation in a triad involving a geostrophic mode. Let us assume the amplitudes $b_{p,q}$ to scale as $\varepsilon _{p,q}$ and that the slow time is given by the small detuning, ${\rm \Delta} \omega _{kpq} = O(\beta )$. We further assume that $\beta$, $\varepsilon _{p,q} = O(kRo)$ at least. The two amplitude equations stemming from (3.1) governing $b_{p,q}$ give the following scaling relationships:

*a*,

*b*)\begin{equation} \beta \varepsilon_p = k Ro \beta \varepsilon_q \quad \mbox{and}\quad \beta \varepsilon_q = k Ro \varepsilon_p, \end{equation}

where we have used that the wave-to-geostrophic coupling coefficient is proportional to the detuning, hence $\beta$, as explained in § 3.1 below (3.3). We have also used that the coupling coefficients scale as the wavenumber $k$. For the hierarchy between $\varepsilon _{p,q}$ and $\beta$ to hold, the determinant of the system (A\,1*a*,*b*) with unknowns $\varepsilon _{p,q}$ must vanish. This condition is satisfied when $\beta = (kRo)^2$. Then, equations (A\,1*a*,*b*) yield $\varepsilon _p = kRo \varepsilon _q$. Therefore, we chose the amplitudes $b_p^{s_p}$ and $b_{q}^{s_q}$ to scale as $(kRo)^2$ and $kRo$, respectively, at leading order.