1. Introduction
The interface of two immiscible fluids enclosed in a vertically vibrating container upon reaching a certain frequency and acceleration becomes unstable. This unstable phase is known as Faraday instability and results in the formation of nonlinear standing waves at the interface of the two fluids. These waves are termed as Faraday waves and were first described by Faraday (Reference Faraday1831). The onset of Faraday instability in immiscible liquids was theoretically demonstrated by Benjamin & Ursell (Reference Benjamin and Ursell1954) using linear stability analysis of the interface of an ideal fluid governed by the set of equations relevant to the system of Mathieu equations. Linear stability of the immiscible, viscous and finite depth fluid problem for one and two frequency excitation was studied by Kumar & Tuckerman (Reference Kumar and Tuckerman1994) and Besson, Edwards & Tuckerman (Reference Besson, Edwards and Tuckerman1996), respectively, using Floquet analysis. Apart from theoretical investigations, experimental (Gollub & Meyer Reference Gollub and Meyer1983; Douady & Fauve Reference Douady and Fauve1988; Simonelli & Gollub Reference Simonelli and Gollub1989; Douady Reference Douady1990; Müller Reference Müller1993; Edwards & Fauve Reference Edwards and Fauve1994; Binks & van de Water Reference Binks and van de Water1997; Kudrolli, Pier & Gollub Reference Kudrolli, Pier and Gollub1998; Arbell & Fineberg Reference Arbell and Fineberg2002; Westra, Binks & Van De Water Reference Westra, Binks and Van De Water2003; Kityk et al. Reference Kityk, Embs, Mekhonoshin and Wagner2005; Rajchenbach, Clamond & Leroux Reference Rajchenbach, Clamond and Leroux2013; Shao et al. Reference Shao, Bevilacqua, Ciarletta, Saylor and Bostwick2020, Reference Shao, Wilson, Saylor and Bostwick2021) and numerical (Perinet, Juric & Tuckerman Reference Perinet, Juric and Tuckerman2009; Takagi & Matsumoto Reference Takagi and Matsumoto2011; Kahouadji et al. Reference Kahouadji, Périnet, Tuckerman, Shin, Chergui and Juric2015; Takagi & Matsumoto Reference Takagi and Matsumoto2015) studies were also performed by many researchers to understand the Faraday wave patterns at the interface of immiscible liquids.
The majority of the studies on Faraday waves deals with immiscible liquids. However, in the past few decades researchers have shifted their focus to understanding the dynamics of Faraday waves in miscible liquids. Zoueshtiagh, Amiroudine & Narayanan (Reference Zoueshtiagh, Amiroudine and Narayanan2009) performed experiments and two-dimensional numerical simulations to investigate the Faraday instability of diffuse interfaces between pairs of miscible liquids of different densities in a rectangular cell. They observed that, above a certain forcing amplitude, the standing waves that appear on the diffuse interface are highly disorganized and interact with each other leading to the mixing of fluids followed by the disappearance of the waviness on the interface when the two fluids are mixed. They also found the mixing phase to be subharmonic, similar to the subharmonic nature of the Faraday waves in immiscible fluids. Further, this finding was supported by Diwakar et al. (Reference Diwakar, Zoueshtiagh, Amiroudine and Narayanan2015) who used Floquet theory in conjunction with a quasisteady approximation to carry out linear analysis of Faraday instability in miscible fluids. An experimental and numerical framework similar to Zoueshtiagh et al. (Reference Zoueshtiagh, Amiroudine and Narayanan2009) has been used by Amiroudine, Zoueshtiagh & Narayanan (Reference Amiroudine, Zoueshtiagh and Narayanan2012) to report the exponential growth of the mixing layer thickness owing to fingering at the interface, followed by a small growth rate that demonstrates the saturation of mixing. Recently, Gréa & Adou (Reference Gréa and Adou2018) conducted three-dimensional (3-D) numerical simulations and demonstrated that Faraday instability can generate turbulent mixing zones for strong forcing parameter and/or for sufficiently random initial condition at the interface of two miscible fluids. Further, they formulated a system of equations based on the second-order correlation spectra for turbulent quantities and used perturbation analysis to estimate the final size of the mixing zone,
\begin{equation} L_{sat}=\frac{2\mathcal{A}g_0}{\omega^2}(2F+4), \end{equation}
where 
$L_{sat}$ is the final size of mixing zone (defined later in (2.5)) in the saturation state, 
$\mathcal {A}=(\rho _1-\rho _2)/(\rho _1+\rho _2)$ is the Atwood number expressing the density contrast between heavy fluid (
$\rho _1$) and light fluid (
$\rho _2$), 
$g_0$ is the mean acceleration, 
$F$ is the oscillation amplitude and 
$\omega$ is the forcing frequency. Gréa & Adou (Reference Gréa and Adou2018) validated the predicted 
$L_{sat}$ against their simulations for a wide range of parameters within the homogeneous framework and full-inhomogeneous systems of two miscible fluids. They reported that the instability of the diffuse interface begins with a small harmonic phase, but the main instability phase is dominantly subharmonic, where the turbulent mixing zone develops. The irreversible mixing associated with the growth of the mixing zone owing to the harmonic to subharmonic transition was numerically quantified by Briard, Gréa & Gostiaux (Reference Briard, Gréa and Gostiaux2019) using potential energies. The total potential energy (TPE) was split into background potential energy (BPE) and available potential energy (APE). They demonstrated that the BPE, which signifies the measure of irreversible mixing, increases after the transition. The APE, which denotes the fraction of the TPE that can be converted to BPE through irreversible mixing, peaks at saturation and is partially released in the flow as BPE. This increase in the BPE causes the numerically obtained final mixing zone size 
$L$ to exceed the theoretically predicted 
$L_{sat}$ (1.1). To further elucidate the dynamics of the turbulent mixing zone driven by the Faraday instability, Briard, Gostiaux & Gréa (Reference Briard, Gostiaux and Gréa2020) performed experiments with fresh and salty water and supported their finding with direct numerical simulations (DNS) and theoretical predictions. They concluded that when the instability is triggered, a natural wavelength appears at the interface between the two fluids. As the amplitude of this wave increases, well-defined structures break up to produce a turbulent mixing layer. At the saturation of the instability, turbulence is inhibited and a mixing layer of final size consistent with the analytically predicted 
$L_{sat}$ was obtained for a wide range of parameters. Mondal & Kumar (Reference Mondal and Kumar2004), performed a linear stability analysis to investigate the effects of Coriolis force on the Faraday waves in a thin sheet of viscous fluid placed on the vibrating plate. They used Floquet theory to solve the set of equations in the form of Mathieu equations and explained the effect of rotation rate on the unstable and stable regions of the stability diagrams. They found that the subharmonic waves get suppressed with increasing rotation rates resulting in the delay of the onset of the Faraday waves. They also reported the existence of a tricritical point at the onset of the instabilities, where subharmonic, harmonic and superharmonic waves are simultaneously excited.
Many researchers explored the physics associated with the onset of Faraday instabilities and the subsequent turbulent mixing layer in two miscible fluids. However, the effect of the Coriolis force on the onset and saturation of the Faraday instabilities that drive turbulent mixing in miscible fluids is absent from the literature. In stably stratified fluids, the turbulent mixing mechanisms depend on the nature of external forcing and background stratification. The vertical time-varying displacement of the density surfaces in the ocean is mainly associated with tidal oscillations (Garrett & Kunze Reference Garrett and Kunze2007; Sarkar & Scotti Reference Sarkar and Scotti2017). These oscillations generate internal waves that can interact with ocean currents and topography, resulting in instability and turbulent mixing in the deep ocean. In oceanic flows, the Earth's rotation significantly influences the dynamics of mixing processes in the thermocline (Noh & Long Reference Noh and Long1990; Fernando Reference Fernando1991; Fleury et al. Reference Fleury, Mory, Hopfinger and Auchere1991). For example, 
$N/f \approx 4.67$, where 
$N$ is the Brunt Väisälä frequency, in the abyssal Southern Ocean at midlatitude (Nikurashin, Vallis & Adcroft Reference Nikurashin, Vallis and Adcroft2013; Rosenberg et al. Reference Rosenberg, Pouquet, Marino and Mininni2015). The presence of rotation in stratified flows leads to a generation of inertial-gravity waves (Hanazaki Reference Hanazaki2002; Smith & Waleffe Reference Smith and Waleffe2002; Liechtenstein, Godeferd & Cambon Reference Liechtenstein, Godeferd and Cambon2005; Ferrari & Wunsch Reference Ferrari and Wunsch2009; Sarkar & Scotti Reference Sarkar and Scotti2017). When stratification is strong and rotation weak 
$N/f>1$, the motions are primarily controlled by gravity wave dynamics with maximum frequency 
$\approx N$, resulting in the propagation of internal gravity waves. However, when rotation is stronger than stratification 
$N/f<1$, inertial waves with maximum frequency 
$\approx f$ are generated. Both rotation and stratification produce highly anisotropic flow patterns with quasihorizontal velocity fields (Liechtenstein et al. Reference Liechtenstein, Godeferd and Cambon2005; Praud, Sommeria & Fincham Reference Praud, Sommeria and Fincham2006). However, vertically elongated column-like structures are formed in the rotation-dominant (
$N/f<1$) flows, whereas horizontally layered pancake-like structures are formed in stratification-dominant flows (
$N/f>1$). Praud et al. (Reference Praud, Sommeria and Fincham2006) performed experiments to investigate the decaying grid turbulence in a rotating stratified fluid. They found an important effect of rotation that the increase in rotation rate (
$\,f$) inhibits the turbulent kinetic energy decay. Most studies in rotating stratified flows have focused on investigating the impact of the ratio 
$N/f$ on the direct and inverse energy cascades depending on the forcing scales at which the energy is injected into the flow (Smith & Waleffe Reference Smith and Waleffe2002; Deusebio, Vallgren & Lindborg Reference Deusebio, Vallgren and Lindborg2013; Marino et al. Reference Marino, Mininni, Rosenberg and Pouquet2013, Reference Marino, Mininni, Rosenberg and Pouquet2014; Herbert et al. Reference Herbert, Marino, Rosenberg and Pouquet2016; Alexakis & Biferale Reference Alexakis and Biferale2018; Pouquet et al. Reference Pouquet, Rosenberg, Marino and Herbert2018; Li et al. Reference Li, Wan, Wang and Chen2020). Marino et al. (Reference Marino, Mininni, Rosenberg and Pouquet2013) conducted DNS and showed that the interplay between rotation and stratification plays an important role in the inverse energy cascade, where the small scales transfer energy to large scales. They found that the inverse energy cascade is most efficient in the range 
$1/2\leq N/f \leq 2$. These studies have demonstrated how rotation is crucial in stratified turbulent flows. Therefore, it would be interesting to explore the influence of rotation (or ratio 
$N/f$) on stratified turbulent flows when the external forcing is time-periodic, which is fundamentally different from that used in the literature mentioned above. The generation of turbulent kinetic and potential energy depends on the strength of the periodic forcing. This gives us an opportunity to investigate how the ratio 
$N/f$ affects the onset, growth and saturation/non-saturation of the turbulent mixing driven by the Faraday instability in the stably stratified miscible fluids. The Faraday instability, which arises when the interface of two fluids is subjected to time-periodic accelerations has the potential to promote mixing in stably stratified fluids. Therefore, understanding the dynamics of the Faraday instability under the influence of rotation can help to gain some insights into turbulent mixing processes in the deep ocean owing to the time-varying oscillations of density surfaces caused by tidal forcing. The present investigation of rotating miscible fluids with time-periodic accelerations will provide a framework to understand the physics associated with the mixing processes in oceanic flows. We present a linear stability analysis of the Faraday waves in miscible fluids under the influence of the Coriolis force. We aim to investigate the effects of the Coriolis force on the harmonic to subharmonic transition, the onset of the subharmonic instability phase, and the saturation phase of the instability. Our analytical model accounts for the vertical inhomogeneity in the background density profile and the Coriolis effect. We followed the classical approach of Yih (Reference Yih1960) to derive the linear equations of the form of the Mathieu equations and use Floquet theory to solve them. We draw the stability diagrams for the corresponding stable and unstable solutions of the Mathieu equations. We also perform DNS at different forcing amplitudes and Coriolis frequencies to check the validity of our analytical model.
The problem formulation, which includes the governing equations and theoretical framework for linear stability analysis, is provided in § 2. The details of the numerical methodology, and simulation parameters, are given in § 3. The predictions of the linear stability analysis, and the numerical results, are explained in § 4. We conclude § 5.
2. Problem formulation
We present the details of the configuration and the governing equations in § 2.1. Then we perform the linear stability analysis and derive the Mathieu equations governing the stability of a full inhomogeneous inviscid rotating system under vertically periodic forcing in § 2.2. We determine the characteristic frequencies 
$\varOmega _i$ (defined later) of the Mathieu equations of the full inhomogeneous problem for the 
$N>f$ case in § 2.3. Based on the values of 
$\varOmega _i$ for the 
$N>f$ case, we investigate the linear theory for the limiting case of the fully developed mixing zone in § 2.3.1 and the limiting case of the interface in § 2.3.2. Finally, we compute the 
$\varOmega _i$ for the 
$N< f$ and 
$N=f$ cases in §§ 2.4 and 2.5, respectively.
2.1. Governing equations
We consider two miscible fluids of small contrasting densities. The lighter fluid of density 
$\rho _2$ is above the denser fluid of density 
$\rho _1$ in a cubical domain (see figure 1). This domain is subject to periodic vertical oscillations with acceleration 
$g(t)=g_0(1+F\cos {(\omega t))}$, where 
$g_0$ is the mean acceleration, 
$F$ is the oscillation amplitude and 
$\omega$ is the forcing frequency. We are assuming that the density of the mixture is linearly varying with mass concentration, with 
$C(\rho _1)=1$ and 
$C(\rho _2)=0$. The 3-D conservation equations for mass, momentum and concentration field under vertically periodic forcing 
$g(t)$, for unsteady incompressible flow with Boussinesq approximation, are
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial \boldsymbol{U}}{\partial t} + \boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{U} =-\boldsymbol{\nabla} P - f\hat{k} \times \boldsymbol{U} + \nu {\nabla}^2 \boldsymbol{U} - 2\mathcal{A}{C}g(t) \hat{k}, \end{gather}$$
$$\begin{gather}\frac{\partial{C}}{\partial t} + \boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{\nabla} {C} = \kappa {\nabla}^2 {C}. \end{gather}$$
Here 
$\boldsymbol {U}$ represent the velocity vector with components (
$U_1$, 
$U_2$ and 
$U_3$) in the 
$x_1$, 
$x_2$ and 
$x_3$ (vertical) directions, respectively, 
${P}$ is the pressure deviation from the hydrostatic background state, 
$f$ is the Coriolis's frequency, 
$\nu$ (
$\mathrm {m}^2\ \mathrm {s}^{-1}$) is the kinematic viscosity and 
$\kappa$ (
$\mathrm {m}^2 \mathrm {s}^{-1}$) is the diffusion coefficient.
Schematic of the domain containing two miscible liquids subjected to vertically periodic forcing.
We decompose a variable into its mean 
$\langle B \rangle _H$ and fluctuating component 
$b$ following Reynolds decomposition,
\begin{equation} B=\langle B \rangle _H + b, \end{equation}
where 
$\langle B \rangle _H$ is the average in the 
$x_1$ and 
$x_2$ horizontal directions, and 
$\langle b \rangle _H=0$. We consider periodic boundary conditions in both 
$x_1$ and 
$x_2$ directions and assume the quantities are statistically invariant in these directions. Therefore, we can define the horizontal average for any quantity, say 
$B(\boldsymbol {x},t)$, as follows:
\begin{equation} \langle B \rangle_H({x_3},t)=\lim_{\substack{(l_{x_1},l_{x_2}) \to (\infty,\infty) }} \frac{1}{l_{x_1} l_{x_2}} \int_{{-l_{x_1}}/{2}}^{{+l_{x_1}}/{2}} \int_{{-l_{x_2}}/{2}}^{{+l_{x_2}}/{2}} B (\boldsymbol{x},t)\,\mathrm{d}\kern0.7pt x_1\,\mathrm{d}\kern0.7pt x_2. \end{equation}
Substituting Reynolds decomposition (2.2) into the governing equations (2.1a), (2.1b) and (2.1c), with mean velocity field 
$\langle \boldsymbol {U} \rangle _H =0$ and 
$\langle \boldsymbol {u} \rangle _H=\langle c \rangle _H=\langle \,p \rangle _H=0$, we obtain the equations for fluctuating velocity 
$\boldsymbol {u}(\boldsymbol {x},t)$ and concentration 
$c(\boldsymbol {x},t)$ fields as follows:
$$\begin{gather} \frac{\partial u_i}{\partial x_i}=0, \end{gather}$$
$$\begin{gather}\frac{\partial u_i}{\partial t} + u_j \frac{\partial u_i}{\partial x_j} =- \frac{\partial p}{\partial x_i} + f\epsilon_{ij3} u_j\hat{e}_3 - 2\mathcal{A}g(t)c\delta_{i3} + \frac{\partial \langle u_i u_3 \rangle_H}{\partial x_3}+ \nu {\nabla}^2 u_i, \end{gather}$$
$$\begin{gather}\frac{\partial c}{\partial t} + u_j \frac{\partial c}{\partial x_j} =- u_3 \frac{\partial \langle C \rangle_H}{\partial x_3} + \frac{\partial \langle u_3 c \rangle_H}{\partial x_3} + \kappa {\nabla}^2 c. \end{gather}$$ We will compute the mixing zone width 
$L(t)$ from the mean concentration profile (Andrews & Spalding Reference Andrews and Spalding1990) as
\begin{equation} L(t)=6\int_{-\infty}^{+\infty} \langle C \rangle _H (x_3,t) (1-\langle C \rangle _H (x_3,t))\,\mathrm{d}\kern0.7pt x_3. \end{equation}We define the Brunt Väisälä (or stratification) frequency as follows:
\begin{equation} N=\sqrt{-2\mathcal{A} g_0\frac{\partial \langle C \rangle_H}{\partial x_3}} = \sqrt{-2\mathcal{A} g_0 \varGamma}. \end{equation} Here, 
$\varGamma = {\partial \langle C \rangle _H}/{\partial x_3}$ is the vertical gradient of mean concentration, which can be approximated as 
${\partial \langle C \rangle _H}/{\partial x_3}=-1/L$.
2.2. Linear theory for the inhomogeneous inviscid rotating system under parametric oscillations
Using similar formulation, assumptions and analytical techniques as in Briard et al. (Reference Briard, Gostiaux and Gréa2020), we extend their linear stability analysis for the inviscid inhomogeneous system to rotating cases and develop a theoretical model to investigate the effects of initial stratification along with the Coriolis force on the triggering and development of instability. We define the initial stratification by Atwood number (
$\mathcal {A}$) and initial mixing zone size (
$L_0$). We consider a piecewise concentration profile, constant in the bottom and upper pure (unmixed) fluids and linear in the mixing layer, with a vertical gradient of concentration 
$\varGamma =0$, in the unmixed region and 
$\varGamma =\text {constant}$, in the mixed region. We also assume that the fluids are inviscid and the fluctuations in velocity and concentration are small. After neglecting the nonlinear, viscous and diffusion terms from (2.4b) and (2.4c) and then eliminating the pressure term, 
$u_1$ and 
$u_2$ from the resulting equation of fluctuating velocity 
$\boldsymbol {u}(\boldsymbol {x},t)$, we obtain (see details in Appendix A)
$$\begin{gather} \frac{\partial^2}{\partial t^2}(\nabla^2 u_3) =-f^2\frac{\partial^2 u_3}{\partial x_3^2} + 2 \mathcal{A} g(t) \varGamma \nabla^2_{H} u_3, \end{gather}$$
$$\begin{gather}\frac{\partial c}{\partial t} =- u_3\varGamma . \end{gather}$$
Here, 
$\nabla ^2_{H} \equiv \partial ^2/\partial x_1^2 + \partial ^2/\partial x_2^2$ is the horizontal Laplacian operator. To solve (2.7a) we first consider a case constant forcing, 
$g(t)=g_0$ and assume a solution for 
$u_3$ of the form
\begin{equation} u_3(x_1,x_2,x_3,t) = \phi(x_3) \, {\rm e}^{{\rm i} (kx_1+lx_2-\varOmega t)}, \end{equation}
where 
$\phi (x_3)$ is the amplitude, 
$\varOmega$ is the temporal frequency or characteristic frequency, and 
$k$ and 
$l$ are the components of the wavevector in the horizontal directions 
$x_1$ and 
$x_2$, respectively. Substitution of (2.8) in (2.7a) yields (see details in Appendix B)
\begin{equation} \frac{\partial^2\phi}{\partial x_3^2} + {K}^2 \frac{( N^2 -\varOmega^2)}{(\varOmega^2 - f^2 )} \phi =0. \end{equation}
Here, 
$K$ is the horizontal wavenumber which is defined as 
$K=\sqrt {k^2 + l^2}$. Now we focus on solving the case with vertically periodic forcing (see (2.7a)). The horizontal directions 
$x_1$ and 
$x_2$ are assumed to be periodic. Therefore, we apply a Fourier transform on 
$u_3(x_1,x_2,x_3,t)$ with respect to 
$x_1$ and 
$x_2$ such that
\begin{equation} u_3(x_1,x_2,x_3,t) = \hat{u}_3(k,l,x_3,t)\, {\rm e}^{{\rm i} (kx_1+lx_2)}. \end{equation}Substituting the solution (2.10) in (2.7a), we get
\begin{equation} \frac{\partial^2}{\partial t^2} \left( \frac{\partial^2 \hat{u}_3 }{\partial x_3^2} - K^2 \hat{u}_3 \right) =-f^2 \frac{\partial^2 \hat{u}_3}{\partial x_3^2} - 2\mathcal{A} g(t) \varGamma K^2 \hat{u}_3. \end{equation}
We further decompose the amplitude 
$\hat {u}_3(k,l,x_3,t)$ on the basis of previously calculated solutions 
$\phi _i$ of the stationary system (2.9), as follows:
\begin{equation} \hat{u}_3(k,l,x_3,t) =\sum_{i} A_i(k,l,t)\phi_i(k,l,x_3). \end{equation}
Here 
$A_i(k,l,t)$ represents the time-dependent part of the amplitude of vertical velocity field modes evolving as a result of periodic forcing and 
$\phi _i(k,l,x_3)$ corresponds to the spatially (
$x_3$) dependent amplitude of the vertical velocity field modes. Substituting the decomposition (2.12) in (2.11), we get
\begin{equation} \sum_{i} \frac{\partial^2}{\partial t^2} \left( A_i\frac{\partial^2 \phi_i }{\partial x_3^2} - K^2 \phi_i A_i \right) =- \sum_{i}f^2 A_i \frac{\partial^2 \phi_i }{\partial {x_3}^2} - \sum_{i}2\mathcal{A} g(t) \varGamma K^2 A_i \phi_i. \end{equation} Since we are interested in the temporally evolving amplitude (
$A_i$) with periodic forcing, we eliminate 
$\phi _i$ and 
$K$ from (2.13) by substituting (2.9) to obtain (see Appendix C)
\begin{equation} \sum_{i} \frac{\partial^2 A_i}{\partial t^2} = \sum_{i}\left( -\frac{f^2 ( N^2 -\varOmega_i^2 ) + N^2( \varOmega_i^2 - f^2 ) (1+F\cos{\omega t}) }{N^2 - f^2} \right) A_i. \end{equation} Now, we define the amplitude of the concentration field mode 
$a_i$ which is related to the amplitude of vertical velocity mode 
$A_i$ via 
${\partial c}/{\partial t} = - u_3\varGamma$ (from (2.7b)), therefore
\begin{equation} \sum_{i} \frac{\partial a_i}{\partial t} = \sum_{i} - A_i\varGamma. \end{equation} Substituting 
$A_i$ from (2.15) in (2.14) and then integrating the resulting equation with respect to 
$t$ and rearranging the terms, we get
\begin{equation} \frac{\partial^2 a_i}{\partial t^2} + \left( f^2 \left(\frac{N^2 - \varOmega_i^2}{N^2 - f^2}\right) + N^2 \left(\frac{\varOmega_i^2 - f^2}{N^2 - f^2}\right) (1+F \cos(\omega t)) \right) a_i = 0. \end{equation} This system of equations is equivalent to a set of Mathieu equations (Kovacic, Rand & Mohamed Sah Reference Kovacic, Rand and Mohamed Sah2018) of characteristic frequency 
$\varOmega _i$, which governs the stability of our full inhomogeneous inviscid rotating system under vertically periodic forcing. Briard et al. (Reference Briard, Gostiaux and Gréa2020) obtained an equation similar to (2.16) in the absence of rotation (
$\,f=0$) for the full inhomogeneous inviscid system. The discrete values of characteristic frequency 
$\varOmega _i$ can be determined for the cases when 
$N>f$, 
$N< f$ and 
$N=f$, as detailed below.
2.3. Characteristic frequency 
$\varOmega _i$ of the full inhomogeneous problem (2.16) when $N>f$
To solve (2.9) for 
$N>f$, we consider the walls of the domain at 
$x_3 = \pm H/2$. Therefore, the boundary conditions at 
$x_3 = \pm H/2$ are 
$\phi =0$. As we are solving (2.9) for the without periodic forcing case, we assume constant 
$\phi ={u_3}_{top}$ at 
$x_3=L_0/2$ (interface between upper lighter fluid and mixed fluid) and constant 
$\phi ={u_3}_{bot}$ at 
$x_3=-L_0/2$ (interface between mixed fluid and bottom denser fluid). Here 
$L_0$ represents the initial mixing zone size, and 
$H$ represents the domain height, as shown in the figure 1. Equation (2.9) has a finite number of solutions (
$\phi _i$) corresponding to the characteristic frequency 
$\varOmega _i$ for a given horizontal wavenumber (
$K$). The solutions of (2.9) (see Appendix D for the solution steps) in the upper pure fluid region (
$x_3 \geq {L_0}/{2},\varGamma =0$), bottom pure fluid region (
$x_3 \leq {-L_0}/{2},\varGamma =0$) and mixed fluid region (
$\lvert x_3 \rvert \leq {L_0}/{2},\varGamma =-N^2/(2\mathcal {A}g_0)$), respectively, are
$$\begin{gather} \phi_i \left( x_3 \geq \frac{L_0}{2} \right) =-{u_3}_{top}\frac{\displaystyle \sinh{\left( K Y_i \left(x_3-\frac{H}{2} \right) \right)}}{\displaystyle \sinh{\left( K Y_i \left(\frac{H-L_0}{2} \right) \right)} }, \end{gather}$$
$$\begin{gather}\phi_i \left( x_3 \leq- \frac{L_0}{2} \right) =+{u_3}_{bot}\frac{\displaystyle \sinh{\left( K Y_i \left(x_3+\frac{H}{2} \right) \right)}}{\displaystyle \sinh{\left( K Y_i \left(\frac{H-L_0}{2} \right) \right)} }, \end{gather}$$
$$\begin{gather}\phi_i \left( \lvert x_3 \rvert \leq \frac{L_0}{2} \right) = \frac{{u_3}_{top}-{u_3}_{bot}}{\displaystyle 2\sin{\left(K X_i \frac{L_0}{2} \right)} } \sin{(K X_i x_3 )} + \frac{{u_3}_{top}+{u_3}_{bot}}{\displaystyle 2\cos{\left(K X_i\frac{L_0}{2} \right)} } \cos{(K X_i x_3 )}, \end{gather}$$where,
\begin{equation} X_i=\sqrt{\frac{N^2 - \varOmega_i^2}{\varOmega_i^2 -f^2}}, \quad Y_i=\sqrt{\frac{\varOmega_i^2}{\varOmega_i^2 -f^2}}.\end{equation}
Both 
$\phi _i$ and 
$\partial _{x_3}\phi_i$ are continuous at 
${x_3}=\pm L_0/2$. Therefore, we differentiate equation (2.17c) with respect to 
$x_3$ at 
${x_3}=+L_0/2$ and 
${x_3}=-L_0/2$, and equate it to the derivative of (2.17a) at 
${x_3}=+L_0/2$ and derivative of (2.17b) at 
${x_3}=-L_0/2$, respectively. After comparing the coefficients of 
${u_3}_{top}$ and 
${u_3}_{bot}$ we obtain the following equations (see Appendix E):
$$\begin{gather} \tan\left(K \sqrt{\frac{N^2 - \varOmega_i^2}{\varOmega_i^2 -f^2}}\frac{L_0}{2} \right) =- \sqrt{\frac{N^2-\varOmega_i^2}{\varOmega_i^2}} \tanh{\left(K \sqrt{\frac{\varOmega_i^2}{\varOmega_i^2 -f^2}} \left(\frac{H-{L}_0}{2} \right) \right)}, \end{gather}$$
$$\begin{gather}\tan\left(K \sqrt{\frac{N^2 - \varOmega_i^2}{\varOmega_i^2 -f^2}}\frac{L_0}{2} \right) = \sqrt{\frac{\varOmega_i^2}{N^2-\varOmega_i^2}} \frac{1}{\displaystyle \tanh{\left(K \sqrt{\frac{\varOmega_i^2}{\varOmega_i^2 -f^2}} \left(\frac{H-{L}_0}{2} \right) \right)}}. \end{gather}$$ If we assume the vertical height of the domain to be infinite, the terms with 
$\tanh {()}$ can be replaced by 1. We are interested in the solutions of (2.19a) and (2.19b) in terms of 
$\varOmega$ for a given value of 
$N$, 
$f$ and 
$KL_0$. In general, the solutions of (2.19a) and (2.19b) show that for given 
$KL_0$, the possible values of characteristic frequency 
$\varOmega _i$ fall between 
$f$ and 
$N$, where 
$N>f$. To demonstrate this, we assume 
$N=3$ and 
$f=1$ and the solutions show that 
$1<\varOmega _i<3$ (as illustrated in figure 2a). All the 
$\varOmega _i$ become asymptotic for 
$KL_0 \gg 1$ and 
$\varOmega _0 \approx N$ represents the fully developed mixing zone case corresponding to the limiting case of homogeneous stratified turbulence as reported by Gréa & Adou (Reference Gréa and Adou2018) and Briard et al. (Reference Briard, Gostiaux and Gréa2020). However, for 
$KL_0 \ll 1$, from (2.19b), we obtain 
$\varOmega _0^2 \simeq f^2 + N^2 (KL_0/2) \simeq f^2 + \mathcal {A}g_0K$ (see details in Appendix F), which corresponds to the limiting case of the interface. Notably, for 
$f=0$, the frequency of interface reduces to 
$\varOmega _0 \simeq \sqrt {N^2 (KL_0/2)} \simeq \sqrt {\mathcal {A}g_0K}$ and the limiting case of the interface presented by Benjamin & Ursell (Reference Benjamin and Ursell1954) and Briard et al. (Reference Briard, Gostiaux and Gréa2020) is recovered.
Possible values of characteristic frequency 
$\varOmega _i$ of the full inhomogeneous problem (2.16) as a function of horizontal wavenumber 
$K$ and initial mixing zone width 
$L_0$. (a) Here 
$\varOmega _i$, normalized by 
$N$, as a function of 
$KL_0$ for 
$N>f$, where 
$\varOmega _i$ are the solutions of (2.19a) and (2.19b) for rotating 
$\,f=1$ and non-rotating 
$f=0$ cases with 
$N=3$. Solutions of (2.19a) are odd and of (2.19b) are even. For the limiting case of the interface (
$KL_0\ll 1$) plot of 
$\varOmega _0/N$ is shown in the inset. (b) Here 
$\varOmega _i$, normalized by 
$f$, as a function of 
$KL_0$, where 
$\varOmega _i$ are solutions of (2.27a) and (2.27b) with 
$f=3$, 
$N=1$ and 
$KH=100$. Plot of 
$\varOmega _i/f$ for 
$f=3$, 
$N=1$ and 
$KH=1000$ is shown in the inset. Solutions of (2.27a) are odd and of (2.27b) are even.
2.3.1. Linear theory for the limiting case of the fully developed mixing zone (when $N>f$)
We now simplify the full inhomogeneous inviscid problem (general Mathieu equation (2.16)) for the homogeneous case corresponding to the fully developed mixing zone where 
$KL_0\gg 1$, because all the 
$\varOmega _i$ values are constant. First, we rewrite the (2.9) as
\begin{equation} \frac{\partial^2\phi_i}{\partial {x_3}^2} + m^2 \phi_i = 0, \end{equation}where
\begin{equation} m^2=\frac{K^2(N^2-\varOmega_i^2)}{(\varOmega_i^2-f^2)}. \end{equation} According to the approximate Wentzel–Kramers–Brillouin (WKB) solution of (2.20) (see Appendix L for WKB approximation), 
$m$ represents the local wavenumber in the 
$x_3$ (vertical) direction. So, (2.21) represents the dispersion relation for the inertial–gravity waves (Smith & Waleffe Reference Smith and Waleffe2002), and is given by
\begin{equation} \varOmega_i^2 = \frac{K^2N^2+m^2f^2}{K^2+m^2} = \frac{K^2}{\lvert\boldsymbol{K}\rvert^2}N^2+ \frac{m^2}{\lvert\boldsymbol{K}\rvert^2}f^2, \end{equation}
where, 
$\lvert \boldsymbol {K}\rvert =\sqrt {K^2+m^2}$ is the magnitude of wavevector 
$\boldsymbol {K}$. We introduce 
$\tan (\theta _i)=K/m$, where 
$\theta _i$ is the angle between the vertical axis and the wavevector 
$\boldsymbol {K}$, with components 
$(k,l,m)$ in the 
$x_1$, 
$x_2$ and 
$x_3$ (vertical) directions, respectively. After substituting the definition of angle 
$\theta _i$ in (2.21), we obtain the characteristic frequency as follows (see Appendix G):
\begin{equation} \varOmega_i^2=N^2\sin^2{(\theta_i)}+f^2\cos^2{(\theta_i)}. \end{equation}
The above relation shows that 
$\varOmega _i$ depends only on the direction of wavevector 
$\boldsymbol {K}$ and not on the magnitude of 
$\boldsymbol {K}$, which is consistent with the observation that all 
$\varOmega _i$ are constant for 
$KL_0\gg 1$. Now, we substitute the characteristic frequency 
$\varOmega _i$ from (2.23) in (2.16) to get (see Appendix H)
\begin{equation} \frac{\partial^2 a_i}{\partial t^2} + (\,f^2 \cos^2{(\theta_i)} + N^2 \sin^2{(\theta_i)} ( 1 + F \cos{(\omega t)} )) a_i= 0 . \end{equation} We define non-dimensional time as 
$\tau =\omega t$ and replace 
$\theta =\theta _i$ and 
$a=a_i$ such that (2.24) becomes
\begin{equation} \frac{\partial^2 a}{\partial \tau^2} + \left( \frac{f^2 \cos^2{(\theta)}}{\omega^2} + \frac{N^2 \sin^2{(\theta)}}{\omega^2} (1+F\cos{(\tau)} ) \right) a = 0. \end{equation}We solve this Mathieu equation (2.25) using Floquet theory, which provides the solutions to linear differential equations with periodic coefficients. Details of solving (2.25) using Floquet theory are in Appendix K. The resulting 3-D stability diagram for (2.25) is illustrated in figure 3. The stability diagram of the general equation (2.16) for the full inhomogeneous inviscid rotating system is similar to (2.25) except for the eigenvalues.
A 3-D stability diagram obtained by solving the Mathieu equation (2.25) in the parameter space 
$(N^2 \sin ^2(\theta )/\omega ^2,f^2 \cos ^2(\theta )/\omega ^2,F )$. All solutions inside the red and cyan coloured regions are unstable in nature with subharmonic and harmonic responses, respectively. Solutions are stable outside the unstable regions.
2.3.2. Linear theory for the limiting case of the interface (when $N>f$)
In the limiting case 
$KL_0 \ll 1$, where the frequency of the well-defined interface is 
$\varOmega _0^2 \simeq f^2 + N^2 (KL_0/2) \simeq f^2 + \mathcal {A}g_0K$, the full inhomogeneous inviscid problem (see (2.16)) reduces to the interface problem (see derivation in Appendix I), which is written as
\begin{equation} \frac{\partial^2 a}{\partial \tau^2} + \left( \frac{f^2}{\omega^2} + \frac{\mathcal{A}g_0K}{\omega^2} (1+F \cos(\tau) ) \right)a = 0. \end{equation}
The stability diagram for the interface problem is obtained by solving (2.26) using Floquet theory (see Appendix K). We show this stability diagram in the 
$(\mathcal{A}g_0K/{\omega ^2}-F)$ plane at different 
${f^2}/{\omega ^2}$ in a later figure and will discuss it in the results section (§ 4.3).
2.4. Characteristic frequency $\varOmega _i$ of the full inhomogeneous problem (2.16) when $N< f$
To solve (2.9) for 
$N< f$, we use the same boundary conditions and piecewise background concentration profile (for simplicity) as used for the previous case 
$N>f$. The final solutions (
$\phi _i$) (see Appendix J for solution steps) are used to derive an equation for 
$\varOmega _i$ by ensuring that 
$\phi _i$ and 
$\partial _{x_3}\phi_i$ are continuous at 
${x_3}=\pm L_0/2$ and using the same procedure as in the earlier case of 
$N>f$ (Appendix E). The resulting equations are written as
$$\begin{gather} \tan\left(K \sqrt{\frac{\varOmega_i^2 - N^2}{f^2 - \varOmega_i^2}}\frac{L_0}{2} \right) =- \sqrt{\frac{\varOmega_i^2 - N^2}{\varOmega_i^2}} \tan{\left(K \sqrt{\frac{\varOmega_i^2}{f^2 - \varOmega_i^2 }} \left(\frac{H-{L}_0}{2} \right) \right)}, \end{gather}$$
$$\begin{gather}\tan\left(K \sqrt{\frac{\varOmega_i^2 - N^2}{f^2 - \varOmega_i^2 }}\frac{L_0}{2} \right) = \sqrt{\frac{\varOmega_i^2}{\varOmega_i^2 - N^2}} \frac{1}{\displaystyle \tan{\left(K \sqrt{\frac{\varOmega_i^2}{f^2 - \varOmega_i^2 }} \left(\frac{H-{L}_0}{2} \right) \right)}}. \end{gather}$$ We solve (2.27a) and (2.27b) to obtain the characteristic frequency 
$\varOmega _i$. Figure 2(b) illustrates 
$\varOmega _i$, normalized by 
$f$, as a function of 
$KL_0$ for 
$f=3$, 
$N=1$ and 
$KH=100$. The solutions 
$\varOmega _i$ at higher 
$KH$ (
$=1000$) are shown in the inset of the figure 2(b). Note that here 
$KL_0< KH$. For all 
$KL_0$ each 
$\varOmega _i$ is constant with possible values between 
$N$ and 
$f$, i.e. 
$N<{\varOmega _i}< f$ or 
$0.333<{\varOmega _i}/{f}<1$. Moreover, for 
$KH=100$ and 
$1000$ each corresponding 
$\varOmega _i$ are approximately equal. The maximum characteristic frequency 
$\varOmega _0\sim f$ corresponds to the frequency of inertial waves in a stratified rotating fluid, when rotation is stronger and stratification is weak (
$N< f$) (Ferrari & Wunsch Reference Ferrari and Wunsch2009 and Davidson Reference Davidson2013, chapters 3 and 4). Therefore, the full inhomogeneous inviscid problem (2.16) can be simplified for the 
$N< f$ case because all the 
$\varOmega _i$ are constant for all 
$KL_0$ and 
$KH$. We rewrite the expression (2.21) for 
$N< f$ as follows:
\begin{equation} m^2=\frac{K^2(\varOmega_i^2-N^2)}{(\,f^2-\varOmega_i^2)}, \end{equation}
where 
$m$ denotes the local vertical (
$x_3$) wavenumber from WKB approximation (see Appendix L). Notably, the dispersion relation becomes the same as that for 
$N>f$ which is given by (2.22) and (2.23). The observation of the constant values of all 
$\varOmega _i$ for each 
$KL_0$ and 
$KH$ (see figure 2b) agrees with the independence of 
$\varOmega _i$ on the magnitude of wavevector and dependence on the direction (
$\theta _i$) of the wavevector shown by the dispersion relation (2.23). Since the expression of 
$\varOmega _i$ is same for the 
$N< f$ and 
$N>f$ cases, (2.16) for the full inhomogeneous inviscid problem simplifies to (2.25). As a result, the 3-D stability diagram for 
$N< f$ is the same as that for the 
$N>f$ case.
2.5. Characteristic frequency $\varOmega _i$ of the full inhomogeneous problem (2.16) when $N=f$
When the stratification and rotation are of equal strength (
$N=f$), the dispersion relation for the inertial-gravity waves (2.22) gives 
$\varOmega = N = f$. Here, we use the concept of group velocity 
$\boldsymbol {c}_g$, which is defined as a gradient of 
$\varOmega$ in the wavenumber space (
$k,l,m$) and shows the velocity at which energy propagates away from a disturbance in the form of wave packets (Davidson Reference Davidson2013, chapters 3 and 4). Since 
$\varOmega = N = f$, 
$\boldsymbol {c}_g=\partial \varOmega /\partial \boldsymbol {K}_i=0$, indicating that the wave energy does not propagate or there is no dispersal of energy by waves, and there is no anisotropy in the large-scale eddies.
3. Numerical methodology and simulation parameters
The governing equations (2.1a), (2.1b) and (2.1c) are discretized using the finite-difference method on a staggered grid arrangement. We store the velocity fields at the cell faces and the pressure and concentration fields at the cell centres. A second-order central finite-difference scheme is employed to compute all the spatial derivatives. For the advancement in time, we use an explicit third-order Runge–Kutta method except for the diffusion terms, which are solved implicitly using the Crank–Nicolson method (Pal Reference Pal2020; Pal & Chalamalla Reference Pal and Chalamalla2020). The pressure Poisson equation is employed to project a velocity field into a divergence-free space and is solved using a parallel multigrid iterative solver to obtain the dynamic pressure. This numerical solver has been extensively validated and used for several DNS of stratified free-shear and wall-bounded turbulent flows (Brucker & Sarkar Reference Brucker and Sarkar2010; Pal, de Stadler & Sarkar Reference Pal, de Stadler and Sarkar2013; Pal & Sarkar Reference Pal and Sarkar2015; Pal Reference Pal2020; Pal & Chalamalla Reference Pal and Chalamalla2020; Naskar & Pal Reference Naskar and Pal2022a,Reference Naskar and Palb). To control the spurious reflections from the disturbances propagating out of the domain, we use sponge regions near the top and bottom boundaries, where damping functions gradually relax the values of the variables to their corresponding values at the boundary. These damping functions are added on the right-hand side of the momentum equation (2.1b) as explained (Brucker & Sarkar Reference Brucker and Sarkar2010). The sponge region is far away from the mixing zone such that it does not affect the dynamics of the mixing of fluids. Periodic boundary conditions are used in the 
$x_1$ and 
$x_2$ (horizontal) directions for all the variables. In the 
$x_3$ (vertical) direction, the top and bottom boundaries are walls, with Dirichlet and Neumann boundary conditions as follows:
\begin{equation} U_1=U_2=U_3=0, \quad \frac{\partial P}{\partial x_3}=\frac{\partial C}{\partial x_3}=0, \text{at} \ x_3={\pm} H/2. \end{equation}The initial concentration profile (Briard et al. Reference Briard, Gréa and Gostiaux2019, Reference Briard, Gostiaux and Gréa2020) in the present simulations is given as
\begin{equation} C(\boldsymbol{x},t=0 )=\frac{1}{2} \left(1-\tanh{\left(\frac{2x_3}{\delta}\right)}\right), \end{equation}
where the parameter 
$\delta$ is used to change the initial mixing zone size-
$L_0$. All of the cases simulated here have the same 
$\delta =0.035$ and corresponding 
$L_0=0.096$. Broadband fluctuations (Pal & Sarkar Reference Pal and Sarkar2015) are imposed on the initial concentration profile (3.2), with an initial energy spectrum as follows:
\begin{equation} E(\mathcal{K})=(\mathcal{K}/\mathcal{K}_0)^4\, \mathrm{e}^{-2(\mathcal{K}/\mathcal{K}_0)^2}, \end{equation}
where, 
$\mathcal {K}_0=30$. The fluctuations are localized at the centreline and are damped exponentially away from the centreline by multiplying the fluctuating concentration 
$c'$ with a damping function given by
\begin{equation} g(x_3)=\mathcal{B}\,\mathrm{e}^{-\frac{1}{2}({x_3}/{\delta})^2}, \end{equation}
where 
$\mathcal {B}=0.7$ is the maximum intensity of fluctuation at the centreline. The size of the computational domain is 
$l_{x_1}=2{\rm \pi} \ \mathrm {m}$ and 
$l_{x_2}=2{\rm \pi} \ \mathrm {m}$ in the horizontal directions. The height of the domain is 
$l_{x_3}=2H=3.5{\rm \pi} \ \mathrm {m}$ which includes sponge region of thickness 
$H_s=0.64{\rm \pi} \ \mathrm {m}$ at top and bottom boundaries. We have used uniform grids 
$N_{x_1}=N_{x_2}=512$ in the horizontal direction, while in the vertical direction 
$x_3$, non-uniform grids 
$N_{x_3}=512$ are used. The grid is clustered at the vertical centre region of thickness 
$1.53{\rm \pi} \ \mathrm {m}$ with 
$\Delta x_{3_{min}}=\Delta x_1=\Delta x_2$, which is fine enough to resolve all the length scales in the mixing zone. We perform simulations with different parameters listed in table 1 to verify our theoretical solutions. The forcing frequency 
$\omega$ is chosen according to the saturation criterion (1.1) for all simulation cases with Coriolis frequency 
$f=0$, such that 
$L_{sat}\simeq 2.45$ at all forcing amplitudes 
$F$. We use the same 
$\omega$ for rotation cases at the corresponding 
$F$ to compare non-rotation and rotation cases.
Parameters for the simulation: here 
$F$ is the forcing amplitude, 
$\omega$ is the forcing frequency, 
$f$ is the Coriolis frequency and 
$L_{end}$ is the final mixing zone size. The Atwood number 
$\mathcal {A}=0.01$, initial mixing zone width 
$L_0=0.096 \ \mathrm {m}$, gravitational acceleration 
$g_0=10\ \mathrm {m\ s}^{-2}$, kinematic viscosity 
$\nu =1\times 10^{-4}\ \mathrm {m}^2\ {\rm s}^{-1}$ and diffusion coefficient 
$\kappa =1\times 10^{-4}\ \mathrm {m}^2\ {\rm s}^{-1}$ is used for all of the cases (Briard et al. Reference Briard, Gréa and Gostiaux2019).
4. Results
The 3-D stability diagram of our Mathieu equation (2.25) in the homogeneous limit corresponding to the fully developed mixing zone case for 
$N>f$ is shown in figure 3. The stability diagram (figure 3) is also valid for the 
$N< f$ case since the Mathieu equation (2.25) is the same for both 
$N>f$ and 
$N< f$ cases. The stability diagram shows the effect of the Coriolis force on the triggering and saturation/non-saturation of the instability through representative examples (line segments of different colours). The coloured curves correspond to the neutral stability curves, which separate the 
$(N^2 \sin ^2(\theta )/\omega ^2,f^2 \cos ^2(\theta )/\omega ^2 )$-plane into regions of stable and unstable (subharmonic or harmonic) solutions. Inside the red-coloured subharmonic instability tongues, the possible frequencies of the growing solutions are only odd multiple of the 
$\omega /2$, which is half of the external forcing frequency. The solutions inside the cyan-coloured harmonic instability tongues can have frequencies that are integer multiples of 
$\omega$. Other solutions that lie outside the instability tongues are stable.
4.1. Exploration of the stability diagram in the absence of the Coriolis force
Figure 4 shows the stability curves only for initial stratification conditions in the absence of the Coriolis force ( 
$f=0$) in the 
$(N/\omega )^2-F$ parameter plane. We obtained this stability curve similar to that of Gréa & Adou (Reference Gréa and Adou2018) and Briard et al. (Reference Briard, Gostiaux and Gréa2020). The stable regions are separated by the unstable subharmonic (inside the red stability curve) and harmonic (inside the cyan stability curve) regions. Multiple frequencies are excited for a given initial stratification condition 
$N_0$ at a fixed 
$F$ and 
$\omega$, corresponding to an angle 
$\theta$ between 
$0$–
${\rm \pi} /2$, as demonstrated by the horizontal brown line segment delineated by two crosses (
$\times )$ in figure 4. This segment ranges from 
$0$ (left-hand end; 
$\theta =0$) to 
$(N_0/\omega )^2$ (right-hand end; 
$\theta ={\rm \pi} /2$), where 
$(N_0/\omega )^2$ represents the initial stratification condition corresponding to the initial mixing zone width-
$L_0$ (because 
$(N/\omega )^2=2\mathcal {A}g_0/(L \omega ^2) \propto 1/L$, from (2.6)). When the initial condition 
$(N_0/\omega )^2$ is in the first stable region, the entire segment (shown with a diamond in figure 4) lies in the stable region. The corresponding 
$\theta$-modes are stable, and the instability does not trigger. However, if the initial conditions 
$(N_0/\omega )^2$ are in the first subharmonic tongue (segment with circles) the instability is triggered immediately by the unstable subharmonic 
$\theta$-modes, resulting in the rapid growth of 
$L$ and decrease in 
$(N/\omega )^2$, as indicated by the arrows pointing to the left. Therefore, in a non-rotating system 
$f=0$, Gréa & Adou (Reference Gréa and Adou2018) deduced the following sufficient condition for the appearance of the instabilities:
\begin{equation} \left(\frac{N}{\omega} \right)^2>\mathcal{G}(0,F). \end{equation}
Here 
$\mathcal {G}(0, F)=1/(2F+4)$ represents the marginal stability curve, shown by the black dashed curve in figure 4. This curve lies on the left-most boundary of the first subharmonic tongue. Finally, the mixing zone size saturates (
$L_{sat}$) owing to the saturation of Faraday instability when the instability condition, (4.1), is no longer valid. The estimation of the saturated mixing zone size in this situation will be as follows (Gréa & Adou Reference Gréa and Adou2018):
\begin{equation} \left(\frac{N}{\omega}\right)^2=\mathcal{G}(0,F)\ \mathrm{or}\ L_{sat}=\frac{2\mathcal{A}g_0}{\omega^2}\frac{1}{\mathcal{G}(0,F)}. \end{equation}
The black squares indicate the 
$L_{sat}$ on the marginal stability curve. When the initial condition is in the first harmonic tongue (segment with squares) or far away (
$(N_0/\omega )^2 > \mathcal {G}(0, F)$), the mixing zone size grows due to both subharmonic and harmonic instabilities, and 
$(N/\omega )^2$ decreases. The transition from harmonic to subharmonic occurs when no more harmonic 
$\theta$-modes are excited, such that 
$(N/\omega )^2$ will cross the left stability curve of the first harmonic tongue. The corresponding mixing zone size is 
$L_{tr}$ and is demonstrated by the black circle (see figure 4). We will also show this transition using numerical simulations in section § 4.4. Afterwards, the mixing zone size saturates at 
$L_{sat}$ (as shown by the black square on the marginal stability curve) when the instability condition 4.1 is no longer valid. Therefore, the mixing zone size-
$L$ saturates for all cases without rotation.
Stability diagram of the Mathieu equation (2.25) for stratification (
$N$) initial condition in the absence of rotation (Coriolis frequency, 
$f=0$). The red and cyan-coloured stability curves separate the stable and unstable (subharmonic and harmonic) regions. The brown horizontal line segment represent the frequencies excited, corresponding to an angle 
$\theta$ between 
$0$–
${\rm \pi} /2$, for the initial condition 
$(N_0/\omega )^2$ (right-hand end (
$\times$) at 
$\theta ={\rm \pi} /2$) in the different tongue-like zones. Arrows indicate the evolution of 
$(N/\omega )^2$ as the mixing zone size-
$L$ increases.
4.2. Exploration of the stability diagram in the presence of the Coriolis force
Figure 5 shows the top view of the 3-D stability diagram (figure 3) at forcing amplitude 
$F=1$. Various frequencies are excited for given 
$f$ and 
$(N_0/\omega )^2$, corresponding to an angle 
$\theta$ between 0–
${\rm \pi} /2$, as represented by the inclined segments in figure 5. In contrast to the horizontal segment in the ‘without rotation’ cases, the inclined segment in the presence of rotation ranges from 
$[(N/\omega )^2,0]$ to 
$[(N/\omega )^2,(\,f/\omega )^2]$ for a given value of 
$f$. The right-hand end (
$\times$) of the inclined segment (at 
$\theta ={\rm \pi} /2$) corresponds to the initial stratification condition 
$(N_0/\omega )^2$, while the left-hand end (
$\times$) (at 
$\theta =0$) corresponds to the given 
$f$.
Top view of the 3-D stability curve (figure 3) at forcing amplitude 
$F = 1$. The inclined segments represent the frequencies excited, corresponding to an angle 
$\theta$ between 
$0$ to 
${\rm \pi} /2$, for 
$(N/\omega )^2$ (right end (
$\times$) at 
$\theta ={\rm \pi} /2$) and 
$(\,f/\omega )^2$ (left end (
$\times$) at 
$\theta =0$). The initial conditions 
$(N_0/\omega )^2$ for given 
$(\,f/\omega )^2$ are in the first left-most (a) stable region and subharmonic tongue, and (b) harmonic tongue. Arrows indicate the evolution of 
$(N/\omega )^2$ as the mixing zone size-
$L$ increases.
Similar to the ‘without rotation’ case, we define the condition for the occurrence of instabilities in the presence of rotation as follows:
\begin{equation} \left(\frac{N}{\omega} \right)^2>\mathcal{G}(\,f^2/\omega^2, F), \end{equation}
where 
$\mathcal {G}(\,f^2/\omega ^2, F)$ represents the extreme left-hand boundary of the first subharmonic tongue. Beyond 
$\mathcal {G}$, Faraday instability is triggered by the unstable subharmonic or harmonic 
$\theta$-modes, resulting in the growth of 
$L$ and decrease of 
$(N/\omega )^2$. To understand the stability diagram in the presence of rotation, we study the following possible configurations based on the initial condition 
$(N_0/\omega )^2$ and a given 
$(\,f/\omega )^2$.
(i) If the initial condition is in the first left-hand most stable region such that
$(N_0/\omega )^2 < \mathcal {G}(\,f^2/\omega ^2, F)$, then two scenarios, namely $(\,f/\omega )^2<0.25$ and $(\,f/\omega )^2 \geq 0.25$, are possible. For $(\,f/\omega )^2<0.25$, all the excited frequencies lie in the stable region, and the corresponding $\theta$-modes are stable, as depicted by the inclined pink line segment in figure 5(a). Therefore, the instability is not triggered, and the mixing zone size will remain at its initial value of $L(t=0)$. Note that we have defined the instability condition, (4.3), based on the fact that if at least one $\theta$-modes fall in any unstable tongues solution will grow, resulting in instability. Therefore, we define the saturation state of the Faraday instability when the instability condition (4.3) is no longer valid and $(\,f/\omega )^2<0.25$ as follows:
(4.4)For
\begin{gather}
\left(\frac{N}{\omega}\right)^2=\mathcal{G}(\,f^2/\omega^2,
F)\ \mathrm{or}\ L_{sat}=\frac{2\mathcal{A}g_0}{\omega^2}\frac{1}{\mathcal{G}(\,f^2/\omega^2,
F)},\quad (\mathrm{when}\ (\,f/\omega)^2<0.25).
\end{gather}
$(\,f/\omega )^2 \geq 0.25$, some of the excited frequencies are in the first unstable subharmonic tongue (red), and the corresponding $\theta$-modes are unstable as demonstrated by the dark blue line segment (figure 5a). Therefore, instabilities trigger, $L$ grows, and $(N/\omega )^2$ decreases as indicated by the arrows directing to the left. The mixing zone size-$L$ will never reach a saturation state because, at every instant, some frequencies will always be excited in the subharmonic region, and the corresponding $\theta$-modes will be unstable, as depicted by the dashed blue segment in figure 5(a). Therefore, the instability condition 4.3 is always satisfied, and $L$ will continue to grow. However, the growth rate of $L$ might be slow at later stages because the $\theta$-modes that fall in the first unstable tongue decrease as $(N/\omega )^2$ decreases.(ii) The instability condition, (4.3), is always satisfied when
$(N_0/\omega )^2$ is chosen to be in the first subharmonic (red) region, and therefore, instabilities are triggered by the unstable $\theta$-modes. For $(\,f/\omega )^2<0.25$, the most unstable $\theta$-modes are excited in the first subharmonic region, as depicted by the solid orange segment in the figure 5(a). Therefore, $L$ grows and $(N/\omega )^2$ decreases. The instability will achieve a state of saturation when no subharmonic unstable $\theta$-modes are excited, i.e. the instability condition, (4.3), is no longer satisfied. The corresponding mixing zone size, $L_{sat}$, is given by the saturation condition, (4.4). Now $L_{sat}$ may lie anywhere on the left stability curve of the first subharmonic tongue, as shown by the orange ellipse (see figure 5a). For $(\,f/\omega )^2 \geq 0.25$, the solid green segment (see figure 5a) represents the initially excited frequencies. For this case, the mixing zone size-$L$ will never reach a saturation state because some unstable $\theta$-modes are always excited in the subharmonic region, as demonstrated by the dashed green segment.(iii) Now we consider the case when the initial condition
$(N_0/\omega )^2$ is beyond the first subharmonic tongue. For the initial condition in the first harmonic tongue (cyan), the grey segment (see figure 5b) represents the case for $(\,f/\omega )^2<0.25$. Owing to the validity of (4.3), both harmonic and subharmonic instabilities are triggered by the unstable $\theta$ -modes. For some initial time, the most unstable $\theta$-modes ($\theta \simeq {\rm \pi}/2$) are excited in the harmonic region, so the mixing zone size-$L$ grows due to harmonic instability and thus $(N/\omega )^2$ decreases. The growth of $L$ is followed by the harmonic to subharmonic transition, and the corresponding mixing zone size $L_{tr}$ is shown by the solid grey ellipse on the left boundary of the first harmonic tongue (see figure 5b). Afterwards, $L$ grows due to subharmonic instability if $\theta$-modes are excited in the first unstable subharmonic region. We will also demonstrate the harmonic to subharmonic transition, followed by a subharmonic instability, using numerical simulations in section § 4.4. Finally, the mixing zone size saturates at $L_{sat}$ (shown by a dashed grey ellipse) when the saturation condition, (4.4), is satisfied. For $(\,f/\omega )^2 \geq 0.25$, only the harmonic instability is triggered, as shown by the solid black segment in figure 5(b). After some time, during the growth of $L$, the $(N/\omega )^2$ can be in the stable region between the first subharmonic and harmonic tongues, as shown by the dash–dotted black segment (figure 5b). No unstable $\theta$-modes are excited in this stable region, hence $L$ will not grow. Therefore, subharmonic instability will not be triggered in the present case. This suggests that the instability condition 4.3 is not valid if the entire inclined segment lies in any of the stable regions beyond the first subharmonic tongue for a given $(\,f/\omega )^2$ and $(N/\omega )^2$. Now, we recall that the present theoretical model does not include the diffusion term. However, in realistic configurations, $L$ can diffuse and grow owing to the molecular diffusion resulting in a decrease in $(N/\omega )^2$. As a result, the right-hand end of the dash–dotted black segment $(N/\omega )^2$ enters the first subharmonic tongue. The subharmonic instability then leads to the growth of $L$. However, the mixing zone size-$L$ never achieves a state of saturation because, at all times, some frequencies are excited in the subharmonic region as depicted by the dashed black segment in figure 5(b).
Figures 6(a), 6(b) and 6(c) illustrate the stability diagram of the Mathieu equation (2.25) at 
$F=0.75, 2$ and 
$3$, respectively. We can observe that the first subharmonic (red) and harmonic (cyan) tongues become wide with increasing forcing amplitude 
$F$ (see comparison among figure 6a–c). Therefore, the right-hand boundary of the first subharmonic tongue (red) shifts towards the right-hand side, and the left-hand boundary of the first harmonic tongue (cyan) moves towards the left-hand side. As a result, the stable region between these tongues shrinks. As 
$F$ increases more subharmonic frequencies are excited and the possibility that the frequencies fall in the stable region between the first subharmonic and harmonic tongues decreases. Consequently, at higher 
$F$, the instabilities are always triggered for all values of 
$(\,f/\omega )^2$ and 
$(N/\omega )^2$ satisfying the instability condition, (4.3), unlike the lower 
$F$ cases where the entire segment can lie in the large stable region between the first subharmonic and harmonic tongues. Additionally, for all 
$(\,f/\omega )^2$, the excitation of more subharmonic frequencies triggers the instability quickly, resulting in a rapid increase in 
$L$. This signifies that for higher 
$F$, the effect of rotation on triggering the instabilities is negligible. The numerical simulations presented in section § 4.4 will further support these findings.
Top view of the 3-D stability curve (figure 3) at three different values of the forcing amplitude (a) 
$F=0.75$, (b) 
$F=2$ and (c) 
$F=3$. The green segment shown in (a) represents the frequencies excited for the ‘without rotation’ case 
$(\,f/\omega )^2=0$ (left end (
$\times$) at 
$\theta =0$) with the initial condition 
$(N_0/\omega )^2$ (right end (
$\times$) at 
$\theta ={\rm \pi} /2$) in the first harmonic tongue.
To summarize, the mixing zone saturates at 
$L_{sat}$ for 
$(\,f/\omega )^2<0.25$, whereas for 
$(\,f/\omega )^2 \geq 0.25$, the mixing zone never saturates and continues to grow due to the excitation of unstable subharmonic 
$\theta$-modes at all times.
4.3. Effect of Coriolis force on the onset of Faraday instability from a sharp interface
Figure 7(a) demonstrates the stability diagram of the interface problem given by (2.26) in the (
$\mathcal {A}g_0K/\omega ^2,F$)-plane. The marginal stability curves are shown for different values of rotation rates 
$f^2$, normalized by 
$\omega ^2$, ranging from 0 to 0.3. Corresponding to 
$f^2/\omega ^2$, the region bounded by the solid curves represents the subharmonic (
$SH$) instability tongue, whereas the region bounded by the dashed curves represents the harmonic (
$H$) instability tongue. The stable regions are outside the instability tongues, which are separated by the solid and dashed neutral stability curves. The growth rates (
$\sigma \neq 0$) of unstable modes for 
$f^2/\omega ^2=0,0.1$ and 
$0.2$ inside the first subharmonic tongue are depicted in the figure 7(b). These growth rates (
$\sigma$) are determined using Floquet analysis (see Appendix K). Corresponding to 
$f^2/\omega ^2$ the solid lines represent 
$\sigma = 0$, whereas dash–dotted lines represent 
$\sigma \neq 0$ where 
$\sigma$ ranges from 0.02 to 0.14 with an interval of 0.02. Inside the instability tongues, the well-defined interface between two fluids is unstable to any perturbation resulting in the generation of standing waves with wavelength 
$\lambda$ (
$=2{\rm \pi} /K$). For example, at 
$F=0.75$ for 
$f^2/\omega ^2=0$, the subharmonic waves can be excited with wavelength 
$\lambda ^{sh}$ between 
$\lambda ^{sh}_{min}$ and 
$\lambda ^{sh}_{max}$ as indicated by the black circles on the horizontal orange dashed line in figure 7(a). We can observe that the first subharmonic tongue moves upwards at higher rotation rates 
$f^2/\omega ^2 \geq 0.25$. Consequently, the subharmonic Faraday waves are not triggered at lower forcing amplitudes 
$F<2$ and 
$f^2/\omega ^2 \geq 0.25$, signifying that the Coriolis force suppresses the formation of subharmonic waves at relatively higher rotation rates. However, the subharmonic instability can always be excited by increasing the forcing amplitude above a critical value 
$F_c$ which is defined as the minima of the first subharmonic tongue for a given 
$f^2/\omega ^2$. In contrast, harmonic waves can always be triggered because harmonic tongues do not move upwards in the presence of Coriolis force, as shown in figure 7(a).
Panel (a) shows the stability diagram of the interface problem (Mathieu equation) (2.26) in the parameter space 
$(\mathcal {A}g_0K/\omega ^2,F )$ for different rotation rates 
$f^2/\omega ^2$. The tongue-like region bounded by the solid curves of respective 
$f^2/\omega ^2$ correspond to the subharmonic (
$SH$) instability, whereas the region bounded by the dashed curves correspond to the harmonic (
$H$) instability. The stable solutions lies outside the 
$SH$ and 
$H$ instability tongues. Panel (b) is an enlargement of the first left-most subharmonic tongues for 
$f^2/\omega ^2=0, 0.1$ and 
$0.2$ in (a). For the corresponding 
$f^2/\omega ^2=0$, solid curves denote the neutral stability curve for the growth rate 
$\sigma =0$, whereas dash–dotted curves denote the contours of different 
$\sigma$ at intervals of 0.02 ranging from 0.02 to 0.14. Arrow lines indicate the increasing order of 
$\sigma$. Panel (c) demonstrates the shrinkage of the first subharmonic tongue with increasing 
$f^2/\omega ^2$ in the ( 
$f^2/\omega ^2 + \mathcal {A}g_0K/\omega ^2,F$)-plane.
Figure 7(a) demonstrates that the first subharmonic tongue shifts leftwards as 
$f^2/\omega ^2$ increases from 
$0$ to 
$\leq 0.25$. In addition, all inside U-tongues for growth rates ranging from 
$\sigma = 0.02$ to 
$0.14$ shifts upwards as 
$f^2/\omega ^2$ increases from 0 to 0.2 (see figure 7b). For example, the minima of U-tongue for 
$\sigma = 0.02$ moves from 
$F\approx 0.162$ to 
$0.27$ and 
$0.76$ when 
$f^2/\omega ^2$ increases from 
$f^2/\omega ^2=0$ to 0.1 and 0.2, respectively, as indicated by the orange squares in figure 7(b). At constant forcing amplitude 
$F$, the growth rate of the most unstable mode decreases with an increase in rotation rate from 
$f^2/\omega ^2=0$ to 
$<0.25$. For example, at 
$F=0.75$ indicated by the horizontal dashed orange line in figure 7(b), the growth rate of the most unstable mode is 
$\sigma \approx 0.095$, 
$0.06$ and 
$0.02$ for 
$f^2/\omega ^2=0$, 0.1 and 0.2, respectively. These results suggest that at a fixed 
$F$, the Coriolis force suppresses the excitation of the subharmonic waves. As a result, the subharmonic waves will take time to develop and grow, thereby delaying the onset of the subharmonic instability. We conclude that at lower forcing amplitudes, the subharmonic instability suppressing and delaying effect of the Coriolis force increases with an increase in rotation rates and eventually becomes so strong that no subharmonic instability is triggered at higher rotation rates, i.e. when 
$F<2$ and 
$f^2/\omega ^2 \geq 0.25$. However, the subharmonic instability can always be excited by increasing 
$F$ above a critical value 
$F_c$. Here we can recall that if molecular diffusion exists, the subharmonic instability can also be triggered from the diffuse interface when 
$F<2$ and 
$f^2/\omega ^2 \geq 0.25$, but not from the sharp interface, as discussed in the previous section (§ 4.2).
Additionally, the first subharmonic tongue and the inside U-tongues shrink when 
$f^2/\omega ^2$ increases from 0 to 0.2 (see figure 7b). This shrinkage can be clearly observed in the 
$(\,f^2/\omega ^2 + \mathcal {A}g_0K/\omega ^2,F)$-plane shown in figure 7(c). Notably, the shrinkage with increasing 
$f^2/\omega ^2$ is significant at lower 
$F$, and becomes insignificant at higher 
$F$. This is illustrated by the following example: in comparison with 
$f^2/\omega ^2=0$ the shrinkage of subharmonic tongue for 
$f^2/\omega ^2=0.2$ is 
$\approx 80\,\%$ at 
$F=0.75$ and 
$\approx 20\,\%$ at 
$F=3$. It suggests that fewer most unstable modes are available at lower 
$F$ than at higher 
$F$, which can excite and grow the subharmonic waves. According to the previous discussion, at fixed 
$F$, the growth rate of the most unstable mode is less for 
$f^2/\omega ^2=0.2$ than for 
$f^2/\omega ^2=0$, resulting in the suppression and delay of the subharmonic instability by the Coriolis force. This observation is only true at lower 
$F$ owing to the additional fact that shrinkage of the subharmonic tongue results in the availability of fewer unstable modes for 
$f^2/\omega ^2=0.2$ than that for 
$f^2/\omega ^2=0$ which are also responsible for the excitation and growth of the subharmonic waves. In contrast to the lower 
$F$, a large number of most unstable modes are available at higher 
$F$, nearly equivalent to that for 
$f^2/\omega ^2=0$, that can rapidly grow the subharmonic waves. Therefore, we can deduce that at higher forcing amplitudes, the flow is energetic enough to overcome the efficacy of rotation in suppressing and delaying the onset of subharmonic instability.
For each 
$f^2/\omega ^2$, the U-tongues for different 
$\sigma$ move upwards as 
$\sigma$ increases, as indicated by the arrows in figure 7(b). For 
$f^2/\omega ^2=0$ and at lower forcing amplitudes approximately 
$F\leq 1$, the minima of all U-tongues deviate slightly from 
$\mathcal {A}g_0K/\omega ^2=0.25$ as shown by a vertical dashed black line in figure 7(b). As a result, for the ‘without rotation’ case 
$f^2/\omega ^2=0$, at a given 
$F$ the wavenumber 
$K^{sh}$ corresponding to the growth rate of the most unstable mode at the onset of subharmonic instability can be approximated as 
$K^{sh}=\omega ^2/(4\mathcal {A}g_0)$ (Briard et al. Reference Briard, Gostiaux and Gréa2020). In the presence of rotation, the first subharmonic tongue (
$\sigma =0$) and the inside U-tongues (
$\sigma \neq 0$) move leftwards with an increase in 
$f^2/\omega ^2$ until 
$\leq 0.25$. Similar to the case 
$f^2/\omega ^2=0$, for approximately 
$F\leq 1$ the minima of U-tongues depart slightly from 
$\mathcal {A}g_0K/\omega ^2=0.15$ and 
$0.05$ for 
$f^2/\omega ^2=0.1$ and 
$0.2$ cases, respectively. Consequently, for 
$f^2/\omega ^2\leq 0.25$, the approximation of 
$K^{sh}$ corresponding to maximum growth rate can be written as
\begin{equation} K^{sh}=\frac{\omega^2}{\mathcal{A}g_0}\left(\frac{1}{4}- \left(\frac{f}{\omega}\right)^2\right). \end{equation} Substituting 
$\lambda ^{sh}=2{\rm \pi} /K^{sh}$ in (4.5), we get
\begin{equation} \lambda^{sh}=\frac{2{\rm \pi}\mathcal{A}g_0}{\omega^2} \left(\frac{1}{{1}/{4}-\left(\,{f}/{\omega}\right)^2}\right). \end{equation} Therefore, this approximation (4.6) of the wavelength of the initial most amplified mode starting from the well-defined sharp interface signifies that 
$\lambda ^{sh}$ increases (or 
$K^{sh}$ decreases) with an increase in 
$f^2/\omega ^2(\leq 0.25)$. With further increase in 
$f^2/\omega ^2(>0.25)$, 
$\lambda ^{sh}$ starts to decrease (or 
$K^{sh}$ increase) since the first subharmonic tongue begins to shift rightwards (see figure 7a). Similar effects of the Coriolis force on the onset of subharmonic waves in a thin sheet of viscous fluid were reported by Mondal & Kumar (Reference Mondal and Kumar2004), who performed a linear stability analysis by solving the set of Mathieu equations using Floquet theory.
4.4. Numerical results
In the analytical model, the instability of the interface under the influence of rotation is presented in terms of the mixing zone size-
$L(\propto \omega ^2/N^2)$. We present the time evolution of the mixing zone size-
$L$ (see (2.5)) for forcing amplitude 
$F=0.75$ in figure 8(a) obtained from our simulations. The initial stratification condition corresponding to the initial mixing zone size-
$L_0=0.096$ is 
$N_0^2/\omega ^2=4.64$ that satisfies the instability condition (4.3) and ensures the onset of the Faraday instability. For the case without rotation (F075f/
$\omega$0) at 
$\omega t\simeq 58$, 
$L\simeq 0.48$ and begins to oscillate with a small amplitude that rapidly grows (shown in the inset of figure 8a indicated by black arrow) owing to the subharmonic instability resulting in turbulent mixing. According to the stability diagram in figure 6(a) (green line segment), the 
$N^2/\omega ^2\simeq 0.93$ corresponding to 
$L\simeq 0.48$ lies in the stable region between the first subharmonic and harmonic tongues. When 
$L$ increases, 
$N^2/\omega ^2$ decreases and enters the first subharmonic tongue resulting in the rapid growth of 
$L$ and turbulent mixing. After 
$\omega t\simeq 130$, an asymptotic state is achieved, where the mixing zone size oscillates without any further increase, signifying the saturation of instability. The final size of the mixing zone is 
$L_{end}\simeq 3.1$ for F075f/
$\omega$0. In the stability diagram (figure 6a, green line segment), 
$N^2/\omega ^2\simeq 0.144$ corresponding to 
$L_{end}\simeq 3.1$ lies in the first stable region and signifies the saturation state of the instability.
Evolution of the mixing zone size-
$L$ with non-dimensional time 
$\omega t$ for without rotation and rotation cases (a) at 
$F=0.75$, (b) at 
$F=1$, (c) at 
$F=2$ and (d) at 
$F=3$. Horizontal dashed black lines represent 
$L_{sat}$ calculated from the prediction (1.1) for 
$f=0$ cases. Inset labelled by transition region is the enlarged view of small oscillation before the harmonic to subharmonic transition 
$L_{tr}$. Insets indicated by arrows correspond to the onset of subharmonic instability (shown by square box).
The presence of rotation at 
$F = 0.75$ impedes the growth of the instabilities. For the case with 
$f/\omega =0.48$ the subharmonic instability is triggered at 
$\omega t\simeq 340\unicode{x2013}355$, 
$L\simeq 1.08\unicode{x2013}1.12$, 
$N^2/\omega ^2\simeq 0.412\unicode{x2013}0.398$, as illustrated by the case F075f/
$\omega$48 in figure 8(a) (see the inset indicated by a green arrow). When the Coriolis frequency is increased to 
$f/\omega =0.59$, the subharmonic instability begins at 
$\omega t\simeq 635\unicode{x2013}650$, 
$L\simeq 1.48\unicode{x2013}1.52$, 
$N^2/\omega ^2\simeq 0.3\unicode{x2013}0.293$, as depicted by the case F075f/
$\omega$59 in figure 8(a) (see the inset indicated by a red arrow). Therefore, the onset of the subharmonic instability causing turbulent mixing is significantly delayed owing to rotation. The large Coriolis frequency ( 
$f/\omega =0.59$) also suppresses the amplitude of oscillations in 
$L$ in the subharmonic instability phase. For 
$f/\omega =0.48$ ( 
$f^2/\omega ^2=0.23$) the instability saturates (
$L_{end}\simeq 3.25$) and small oscillations in the asymptotic state are damped out by diffusion, whereas for 
$f/\omega =0.59$ ( 
$f^2/\omega ^2=0.35$) the instability never saturates and hence the mixing zone size-
$L$ continues to evolve with oscillations. This result is consistent with our theoretical prediction that the instability saturates for 
$f^2/\omega ^2<0.25$ but never saturates for 
$f^2/\omega ^2\geq 0.25$.
At forcing amplitude 
$F=1$, the time evolution of 
$L$ for ‘without rotation’ (F1f/
$\omega$0) and ‘with rotation’ (F1f/
$\omega$48 and F1f/
$\omega$59) cases are illustrated in figure 8(b). For all the cases, the instability condition, (4.3), is satisfied by the initial stratification condition 
$N_0^2/\omega ^2=4.25$ corresponding to the initial mixing zone size-
$L_0=0.096$ that ensures the onset of the Faraday instability. Before the rapid growth of 
$L$ small oscillations of varying periods (shown in the inset labelled by a transition region of figure 8b) is observed for both the ‘without rotation’ ( 
$f/\omega =0$) and ‘with rotation’ ( 
$f/\omega =0.48$ and 
$f/\omega =0.59$) cases. These small oscillations represent the harmonic phase following which the harmonic to subharmonic transition (Briard et al. Reference Briard, Gréa and Gostiaux2019) occurs. For 
$f/\omega =0.48$ and 
$0.59$, this transition occurs at 
$\omega t\simeq 60\unicode{x2013}70$ and 
$65\unicode{x2013}75$ (see inset of figure 8b labelled by transition region) followed by the stable region until 
$\omega t\simeq 150\unicode{x2013}160$ and 
$290\unicode{x2013}300$, respectively. No harmonic or subharmonic modes are excited in the stable region, and 
$L$ grows due to diffusion. We also observe similar stable regions for F075f/
$\omega$48 and F075f/
$\omega$59 in figure 8(a) (see the Supplementary movies 1 and 2 are available at https://doi.org/10.1017/jfm.2023.744 of the concentration field for 
$F = 0.75$ and 
$1$). The presence of this stable region delays the onset of the subharmonic instability. This verifies our theoretical prediction in § 4.2 that no instability is triggered if the excited frequencies lie in the stable region between the first subharmonic and harmonic tongues for 
$f^2/\omega ^2\geq 0.25$. However, the presence of diffusion causes 
$L$ to grow and eventually subharmonic instability triggers. Based on the numerical results we can also conclude that the same observation is also true for 
$f^2/\omega ^2<0.25$ cases if all the modes are excited in the stable region. Here 
$L$ rapidly grows after the triggering of the subharmonic instability resulting in turbulent mixing. Similar to 
$F = 0.75$, we observe that 
$L$ saturates for F1f/
$\omega$48, but continues to grow for F1f/
$\omega$59.
The effect of rotation on the evolution of 
$L$ at higher forcing amplitudes 
$F=2$ and 
$F=3$ is demonstrated in figures 8(c) and 8(d), respectively. The initial stratification condition is 
$N_0^2/\omega ^2=3.26$ for 
$F=2$ and 
$N_0^2/\omega ^2=2.6$ for 
$F=3$ corresponding to the initial mixing zone size-
$L_0=0.096$ which ensures the onset of Faraday instability. At 
$F=2$, the harmonic to subharmonic transition (
$L_{tr}\simeq 0.5$) occurs at 
$\omega t\simeq 13$ for F2f/
$\omega$0, F2f/
$\omega$48 and F2f/
$\omega$59 as indicated by the black circle in the inset of figure 8(c). Further increasing the forcing amplitude to 
$F=3$, the transition (
$L_{tr}\simeq 0.375$) occurs at 
$\omega t\simeq 6.5$ for the F3f/
$\omega$0, F3f/
$\omega$48 and F3f/
$\omega$59 cases, as shown by the black circle in the inset of figure 8(d). The instability delaying effect of rotation is subdued as the forcing amplitude increases. This observation is consistent with our analytical model that predicts the shrinkage of the stable region between the unstable harmonic and subharmonic tongues with increasing forcing amplitudes (see stability diagrams at 
$F = 2, 3$ in figures 6b and 6c, respectively). For cases with 
$f^2/\omega ^2=0$ and 
$f^2/\omega ^2=0.23$, the instability saturates, and the oscillations in the asymptotic state are damped. However, the instability never saturates for cases with 
$f^2/\omega ^2=0.35$ and 
$L$ continues to evolve with oscillations. Note that at each corresponding forcing amplitude 
$F=0.75,1,2$ and 
$3$, the values of the final mixing zone size 
$L_{end}$ for the corresponding rotation cases 
$f/\omega =0.48$ ( 
$f^2/\omega ^2<0.25$) are very close to the non-rotation cases ( 
$f/\omega =0$). We report 
$L_{end}$ values in table 1. Therefore, for the parameters used in the present simulations, the main difference between the ‘with rotation’ and ‘without rotation’ cases at each 
$F$ comes principally from the growth of the instability.
Fourier transforms of the derivatives of 
$L$ for different time intervals, each starting from 
$t=0$, are performed to capture the harmonic to subharmonic transition in our simulations. For case F1f/
$\omega$0, the normalized Fourier transform of 
${\rm d}L/{\rm d}t$ as a function of the normalized frequency 
$f/(\omega/2{\rm \pi} )$ is depicted in figure 9(a–d). Two peaks of different frequencies are present in 
$L$. The first and the second peaks represent the subharmonic and harmonic regimes, respectively. We can observe that for 
$\omega t \leq 45$, the second peak is higher (see figure 9a), denoting the dominance of the harmonic regime. At 
$\omega t = 45.92$ (figure 9b) the amplitude of both the peaks are similar (see figure 9b) indicating the transition from harmonic to subharmonic regime. Note that the small oscillations in 
$L$ also rapidly grow from 
$\omega t\simeq 45$ (see figure 8b) owing to the onset of the subharmonic instability. When 
$\omega t > 45$ the first peak becomes the dominating peak as illustrated in figures 9(c) and 9(d), signifying the primary role of the subharmonic instabilities in turbulent mixing. The Fourier transform of 
${\rm d}L/{\rm d}t$ for F1f/
$\omega$48 demonstrates the harmonic to subharmonic transition, stable region and the onset of the subharmonic instability. The dominating second peak at 
$\omega t = 65.4$ in figure 9(e) shows small oscillations of 
$L$ (see inset of figure 8b indicated by a green arrow) are in the harmonic region. The amplitude of the first and second peak remains approximately unchanged until 
$\omega t\simeq 150$ (figures 9f and 9g) which represents the stable region. At 
$\omega t\simeq 150\unicode{x2013}160$, the first peak starts increasing, and the second peak starts decreasing. Eventually the first peak becomes the dominant (figure 9h) marking the onset of the subharmonic instability at 
$\omega t\simeq 150\unicode{x2013}160$. An increase in the Coriolis frequency to 
$f/\omega =0.59$ impedes the harmonic to subharmonic transition, expands the stable region, and delays the onset of the subharmonic instability (figures not shown).
Normalized Fourier transform of 
${\dot L}$ as a function of normalized frequency, (a–d) for F1f/
$\omega 0$ case, (e–h) for F1f/
$\omega 48$ case, for different time intervals each starting from 
$\omega t=0$ to 
$\omega t$ mentioned on the top of each figure. The first peak at 
$f /(\omega /2 {\rm \pi})\simeq 1$ corresponds to the subharmonic regime and second peak at 
$f /(\omega /2 {\rm \pi})\simeq 2$ corresponds to the harmonic regime (Briard et al. Reference Briard, Gréa and Gostiaux2019).
Figure 10 depicts the 3-D concentration field in the mixing zone for ‘without rotation’ F1f/
$\omega$0 and ‘with rotation’ F1f/
$\omega$48 cases. The random wave patterns on the interface at 
$\omega t=41.5$ for F1f/
$\omega$0 indicate the inception of the Faraday instability as shown in figure 10(a). The turbulent mixing caused by the subharmonic instability is demonstrated at 
$\omega t=70.6$ for F1f/
$\omega$0 in figure 10(b) followed by the beginning of instability saturation and fully saturated instability in figures 10(c) and 10(d) at 
$\omega t=120$ and 
$\omega t\sim 152$, respectively. For F1f/
$\omega$48, the stable region after the harmonic to subharmonic transition is depicted by figure 10(f) at 
$\omega t=120$. The subsequent figures 10(g) and 10(h) for F1f/
$\omega$48 at 
$\omega t \sim 212$ and 
$329$ show a significant delay in turbulent mixing and saturation of the instability with respect to F1f/
$\omega$0. Movies (1)–(4) of the concentration field at the 
$x_2 = {\rm \pi}$ plane for different forcing amplitudes are included as Supplementary movies, to demonstrate the influence of rotation on the evolution of the mixing zone.
Visualization of 3-D concentration field, (a–d) for F1f/
$\omega$0 case, (e–h) for F1f/
$\omega$48 case, at different time instants. Colour map: blue for 
$C=0$ (lighter fluid) and red for 
$C=1$ (denser fluid). For better visualization of the mixing of fluids, the pure fluids are made transparent.
5. Conclusions
The effect of rotation on the onset and saturation of Faraday instability in two miscible fluids subjected to vertical oscillations is investigated by means of theoretical analysis and DNS. A linear stability analysis is performed in the presence of rotation, yielding a set of Mathieu equations for the full inhomogeneous inviscid rotating problem. These equations are investigated further for the limiting case of the interface and homogeneous case when 
$N>f$. The 
$N< f$ case analysis yields the same Mathieu equations as the homogeneous case for 
$N>f$. The analytical model utilizes Floquet theory to solve a set of Mathieu equations derived from linear stability analysis. On solving the set of Mathieu equations, a 3-D stability diagram is obtained that consists of stable and unstable harmonic and subharmonic regions. In the limiting case of interface, we found that increasing rotation rates suppress and delay the onset of subharmonic instability at lower forcing amplitudes (
$F$). Eventually, no subharmonic instability is triggered at relatively strong rotation rates. However, the subharmonic instability can always be excited by increasing 
$F$ above a critical value 
$F_c$. In contrast, the instability suppressing and delaying effect of rotation diminishes with an increase in 
$F$. The reason for this occurrence is the increase in the area of the first subharmonic unstable tongue resulting in the majority of excited modes being most unstable. We also predicted a wavelength (
$\lambda ^{sh}$) of the subharmonic waves corresponding to the initially most amplified mode starting from the sharp interface, revealing that 
$\lambda ^{sh}$ increases with an increase in 
$f^2/\omega ^2(\leq 0.25)$.
Further, we have performed simulations at different forcing amplitudes and rotation rates to verify the findings of our theoretical model. The harmonic to subharmonic transition followed by the onset of the subharmonic instability, and its saturation, is shown in terms of the time evolution of the mixing zone size 
$L$. For 
$F = 0.75$ and 
$1$, the evolution of the mixing zone size manifests a significantly slow growth under the influence of rotation. This slow growth of 
$L$ represents the stable region between the harmonic to subharmonic transition and before the onset of the subharmonic instabilities. After this stable region, the subharmonic instability triggers, and 
$L$ rapidly grows, resulting in turbulent mixing. At forcing amplitudes of 
$2$ and 
$3$, the vertical vibrations are strong enough to diminish the instability delaying effect of rotation. Therefore, the mixing zone size for the cases with 
$F = 2$ and 
$3$ with rotation grows similar to the non-rotating cases. We also perform the Fourier transform of the rate of change of 
$L$ to capture the harmonic to subharmonic transition and demonstrate the delay in this transition in terms of the peak frequencies.
Finally, we have compared the simulations for 
$f/\omega =0.48$ and 
$f/\omega =0.59$, which correspond to 
$(\,f/\omega )^2 < 0.25$ and 
$(\,f/\omega )^2 \geq 0.25$, respectively, on the 3-D stability diagram, at different forcing amplitudes. We find that for 
$f/\omega =0.48$, the instability saturates, and the mixing zone size asymptotes. However, for 
$f/\omega =0.59$, the instability does not saturate, and the mixing zone size keeps growing. This confirms our theoretical observation that for 
$(\,f/\omega )^2 < 0.25$ the instability saturates, whereas it never saturates for 
$(\,f/\omega )^2 \geq 0.25$ because some 
$\theta$-modes are always excited in the first subharmonic tongue at all times.
The turbulent mixing, attributed to the onset of the subharmonic instability, is associated with exchanges among the TPE, the BPE and the APE. We aim to extend the present investigation to explore the influence of rotation on these energy exchanges through reversible buoyancy flux, irreversible diapycnal mixing and irreversible kinetic energy dissipation to quantify the irreversible mixing efficiency.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2023.744.
Acknowledgements
We would like to thank the anonymous referees whose helpful suggestions improved the quality of this paper.
Funding
We gratefully acknowledge the support of the Science and Engineering Research Board, Government of India, grant no. SERB/ME/2020318. We also want to thank the Office of Research and Development, Indian Institute of Technology Kanpur for the financial support through grant no. IITK/ME/2019194. The support and the resources provided by PARAM Sanganak under the National Supercomputing Mission, Government of India at the Indian Institute of Technology, Kanpur are gratefully acknowledged.
Declaration of interests
The authors report no conflict of interest.
Author contributions
The authors contributed equally to analysing data and reaching conclusions, and in writing the paper.
Appendix A. Steps to derive (2.7a)
We neglect the nonlinear, viscous and diffusion terms from (2.4b) and (2.4c), which yields
$$\begin{gather} \frac{\partial u_i}{\partial t} =- \frac{\partial p}{\partial x_i} + f\epsilon_{ij3} u_j\hat{e}_3 - 2\mathcal{A}g(t)c\delta_{i3} , \end{gather}$$
$$\begin{gather}\frac{\partial c}{\partial t} =- u_3\varGamma . \end{gather}$$Differentiating (2.4a) with respect to time and using (A1a) yields
\begin{equation} \nabla^2_{H} p = f\left(\frac{\partial u_2}{\partial x_1} - \frac{\partial u_1}{\partial x_2} \right) + \frac{\partial^2 u_3}{\partial x_3 \partial t}, \end{equation}
where 
$\nabla ^2_{H} \equiv \partial ^2/\partial x_1^2 + \partial ^2/\partial x_2^2$ is the horizontal Laplacian operator. Now, eliminating the pressure term from (A2) by taking the 
$\nabla ^2_{H}$ of (A1a) for 
${i}=3$ and using (A2) we obtain
\begin{equation} \frac{\partial}{\partial t}(\nabla^2 u_3) =-\frac{\partial}{\partial x_3}\left[\,f \left(\frac{\partial u_2}{\partial x_1}-\frac{\partial u_1}{\partial x_2}\right)\right] - 2 \mathcal{A} g(t) \nabla^2_{H} c. \end{equation}
Next, we eliminate 
$u_1$ and 
$u_2$ from (A3) by differentiating it with respect to ‘
$t$’, and using (A1a) and (A1b) yields (2.7a).
Appendix B. Steps to derive (2.9)
Starting with
\begin{align} &\frac{\partial^2}{\partial t^2} \left( \left(\phi (i^2k^2+i^2l^2)+\frac{\partial^2\phi}{\partial x_3^2}\right)\, {\rm e}^{{\rm i} (kx_1+lx_2-\varOmega t)} \right) =-f^2\frac{\partial^2 \phi}{\partial x_3^2}\,{\rm e}^{{\rm i} (kx_1+lx_2-\varOmega t)} \nonumber\\ &\quad + 2 \mathcal{A} g_0 \varGamma \phi (i^2k^2+i^2l^2)\, {\rm e}^{{\rm i} (kx_1+lx_2-\varOmega t)} , \end{align}
where 
$i^2=-1$, and substituting 
$N^2=-2\mathcal {A}g_0 \varGamma$ (from (2.6)) and horizontal wavenumber 
$K$, which is defined as 
$K=\sqrt {k^2 + l^2}$, in (B1) we get
$$\begin{gather} (-i\varOmega)^2 \left(-\phi K^2 +\frac{\partial^2\phi}{\partial x_3^2}\right) {\rm e}^{{\rm i} (kx_1+lx_2-\varOmega t)} =-f^2\frac{\partial^2 \phi}{\partial x_3^2}\,{\rm e}^{{\rm i} (kx_1+lx_2-\varOmega t)} +N^2 \phi K^2 {\rm e}^{{\rm i} (kx_1+lx_2-\varOmega t)} , \end{gather}$$
$$\begin{gather}-\varOmega^2 \left(-\phi K^2 +\frac{\partial^2\phi}{\partial x_3^2}\right) =-f^2\frac{\partial^2 \phi}{\partial x_3^2} +N^2 \phi K^2 . \end{gather}$$Appendix C. Derivation of (2.14)
Substituting equation (2.9) into (2.13) yields
\begin{align} \sum_{i} \frac{\partial^2}{\partial t^2} \left(-A_i K^2 \left(\frac{ N^2 -\varOmega_i^2 }{\varOmega_i^2 - f^2}\right)\phi_i - K^2 \phi_i A_i \right) &= \sum_{i}f^2 A_i K^2 \left(\frac{ N^2 -\varOmega_i^2 }{\varOmega_i^2 - f^2}\right) \phi_i \nonumber\\ &\quad - \sum_{i}2\mathcal{A} g(t) \varGamma K^2 A_i \phi_i , \end{align}rearranging the terms we get
\begin{align} \sum_{i} - K^2 \left(\frac{ N^2 -\varOmega_i^2 }{\varOmega_i^2 - f^2}+1\right)\phi_i \frac{\partial^2 A_i}{\partial t^2} &= \sum_{i}f^2 A_i K^2 \left(\frac{ N^2 -\varOmega_i^2 }{\varOmega_i^2 - f^2}\right)\phi_i \nonumber\\ &\quad - \sum_{i}2\mathcal{A} g(t) \varGamma K^2 A_i \phi_i . \end{align} Taking the 
$K^2$ common from both sides of the above (C2) which cancels out each other, and after arranging the terms, we get
\begin{equation} \sum_{i} -\left(\frac{ N^2 -f^2 }{\varOmega_i^2 - f^2}\right)\frac{\partial^2 A_i}{\partial t^2}\phi_i = \sum_{i} \left( f^2 \left(\frac{ N^2 -\varOmega_i^2 }{\varOmega_i^2 - f^2}\right) - 2\mathcal{A} g(t) \varGamma \right) A_i \phi_i.\end{equation} Since 
$\phi _i$ forms a basis on which the decomposition is unique, therefore, the coefficients of each 
$\phi _i$ on both sides of the above (C3) would be the same. Hence, we get
\begin{equation} \sum_{i} -\left(\frac{ N^2 -f^2 }{\varOmega_i^2 - f^2}\right)\frac{\partial^2 A_i}{\partial t^2} = \sum_{i} \left( f^2 \left(\frac{ N^2 -\varOmega_i^2 }{\varOmega_i^2 - f^2}\right) - 2\mathcal{A} g(t) \varGamma \right) A_i . \end{equation}Further rearranging the terms of the above (C4) and doing some simplifications, we get
\begin{equation} \sum_{i} \frac{\partial^2 A_i}{\partial t^2} = \sum_{i}\left( -\frac{f^2 \left( N^2 -\varOmega_i^2 \right) - ( \varOmega_i^2 - f^2 ) 2 \mathcal{A}g(t)\varGamma}{N^2 - f^2} \right) A_i. \end{equation}After substituting 
$g(t)=g_0(1+F\cos {\omega t} )$ and 
$\varGamma =-N^2/(2 \mathcal {A}g_0)$ (see (2.6)) in (C5) we get (2.14).
Appendix D. Steps for solving (2.9) when $N>f$
We show the steps to solve (2.9) for the 
$N>f$ case which represents the second-order linear homogeneous equation with constant coefficients. First, we solve this equation in the upper pure fluid regime where vertical gradient of mean concentration 
$\varGamma =0$ and thus 
$N=\sqrt {-2\mathcal {A} g_0 \varGamma }=0$, according to our assumption of piecewise background concentration profile. So (2.9) becomes
\begin{equation} \frac{\partial^2\phi}{\partial x_3^2} - \frac{{K}^2 \varOmega^2 }{(\varOmega^2 - f^2 )} \phi = 0 . \end{equation} Let 
$q^2=K^2\varOmega ^2/(\varOmega ^2 - f^2 )$ and thus the above equation can be written as
\begin{equation} \frac{\partial^2\phi}{\partial x_3^2} - q^2 \phi = 0 . \end{equation} Let 
$\phi ={\rm e}^{rx_3}$ be a solution to this equation, for some as-yet-unknown constant 
$r$. After substituting this assumed solution in (D2) we get 
${\rm e}^{rx_3}(r^2-q^2 )=0$. Since 
${\rm e}^{rx_3}$ is never zero, the above equation is satisfied if and only if 
$r^2-q^2=0$ and roots of this characteristic equation are 
$r=\pm q$, which are real and distinct. Therefore, the general solution of (D2) is
\begin{equation} \phi= {A\,{\rm e}^{qx_3}+B\,{\rm e}^{-qx_3}} . \end{equation} To find out the constants 
$A$ and 
$B$, we consider the boundary condition 
$\phi =0$ at 
$x_3=H/2$ and 
$\phi ={u_3}_{top}$ at 
$x_3=L_0/2$. Substituting the boundary conditions in (D3) and solving for the constants 
$A$ and 
$B$ we get
$$\begin{gather} A = {{-u_3}_{top}\frac{ \exp ({-q H/{2}})}{\exp ({q (H-L_0)/{2}})-\exp ({-q (H-L_0)/{2}})} }, \end{gather}$$
$$\begin{gather}B = {{u_3}_{top}\frac{ \exp ({q H/{2}})}{\exp ({q (H-L_0)/{2}})-\exp ({-q (H-L_0)/{2}})} } . \end{gather}$$ After substituting the constants 
$A$ and 
$B$ from (D4a) and (D4b), respectively, into (D3), we get
\begin{equation} \phi\left(x_3\geq\frac{L_0}{2}\right) =-{u_3}_{top}\frac{\exp ({q (x_3-H/2)})-\exp ({-q (x_3-H/2)})}{\exp ({q (H-L_0)/{2}})-\exp ({-q (H-L_0)/{2}})} . \end{equation}We can write the exponential terms in the form of hyperbolic trigonometric functions as
\begin{equation} \sinh (a)=\frac{\exp (a)-\exp (-a)}{2} . \end{equation}Therefore, we can write
\begin{equation} \phi \left( x_3 \geq \frac{L_0}{2} \right) =-{u_3}_{top}\frac{\displaystyle \sinh{\left( q \left(x_3-\frac{H}{2} \right) \right)}}{\displaystyle \sinh{\left(q \frac{\left( H-L_0 \right)}{2} \right)} } . \end{equation} Now, substituting the value of 
$q$ in the above (D7) we get the final solution as follows:
\begin{equation} \phi \left( x_3 \geq \frac{L_0}{2} \right) =-{u_3}_{top} \frac{ \sinh{\left( K\sqrt{\dfrac{\varOmega^2}{\varOmega^2 - f^2}} \left(x_3-\dfrac{H}{2} \right) \right)}} {\sinh{\left( K\sqrt{\dfrac{\varOmega^2}{\varOmega^2 - f^2}} \left(\dfrac{H-L_0}{2} \right) \right)}}. \end{equation} We follow the above same procedure to find out the solution for (2.9) in the bottom pure fluid region from 
$x_3=-L_0/2$ to 
$x_3=-H/2$. In this region, we consider the vertical gradient of mean concentration 
$\varGamma =0$, thus 
$N=\sqrt {-2\mathcal {A} g_0 \varGamma }=0$ and with boundary condition 
$\phi ={u_3}_{bot}$ at 
$x_3=-L_0/2$ and 
$\phi =0$ at 
$x_3=-H/2$, we get the solution of (2.9) as follows:
\begin{equation} \phi \left( x_3 \leq-\frac{L_0}{2} \right) = {u_3}_{bot} \frac{\displaystyle \sinh{\left( K\sqrt{\frac{\varOmega^2}{\varOmega^2 - f^2}} \left(x_3+\frac{H}{2} \right) \right)}}{\displaystyle \sinh{\left( K\sqrt{\frac{\varOmega^2}{\varOmega^2 - f^2}} \left(\frac{H-L_0}{2} \right) \right)}} . \end{equation} Next, we see the solution of (2.9) in the mixed fluid region between 
$x_3=L_0/2$ to 
$x_3=-L_0/2$ where the vertical gradient of mean concentration is 
$\varGamma =-N^2/2\mathcal {A}g_0$. Let 
$q^2={K}^2 ({( N^2 -\varOmega ^2 )}/{(\varOmega ^2 - f^2 )})$ then (2.9) can be written as
\begin{equation} \frac{\partial^2\phi}{\partial x_3^2} + q^2 \phi = 0 . \end{equation} Again we assume the solution of the above equation of the form 
$\phi ={\rm e}^{rx_3}$ where 
$r$ is an unknown constant. Substituting this solution in (D10) we get 
${\rm e}^{rx_3}(r^2+q^2 )=0$, where 
${\rm e}^{rx_3}\neq 0$, therefore we get the characteristic equation 
$r^2+q^2=0$ which has two complex roots 
$r=\pm iq$. Thus, the general solution of (D10) is
\begin{equation} \phi= {A\,{\rm e}^{{\rm i}qx_3}+B\,{\rm e}^{-{\rm i}qx_3}} . \end{equation} We use Euler's formula 
${\rm e}^{{\rm i}a}=\cos (a)+i\sin (a)$ to write the solution in the form of trigonometric terms and we get
$$\begin{gather} \phi= A\cos (qx_3)+iA\sin (qx_3)+B\cos (-qx_3)+iB\sin (-qx_3), \end{gather}$$
$$\begin{gather}\phi= (A+B)\cos (qx_3)+i(A-B)\sin (qx_3) , \end{gather}$$
$$\begin{gather}\phi= C\cos (qx_3)+iD\sin (qx_3) , \end{gather}$$
where 
$C$ and 
$D$ are the unknown constants. Now we apply the boundary conditions 
$\phi ={u_3}_{top}$ at 
$x_3=+L_0/2$ and 
$\phi ={u_3}_{bot}$ at 
$x_3=-L_0/2$ to the (D12c) and solving for constants 
$C$ and 
$D$, we get
\begin{equation} C=\frac{{u_3}_{top}+{u_3}_{bot}}{\displaystyle 2\cos \left(\frac{qL_0}{2}\right)};\quad D=\frac{{u_3}_{top}-{u_3}_{bot}}{\displaystyle 2i\sin \left(\frac{qL_0}{2}\right)} . \end{equation}Substituting values of these constants in (D12c) we obtain
$$\begin{gather} \phi \left( \lvert x_3 \rvert \leq \frac{L_0}{2} \right)=\frac{{u_3}_{top}+{u_3}_{bot}}{\displaystyle 2\cos \left(\frac{qL_0}{2}\right)}\cos (qx_3)+i\frac{{u_3}_{top}-{u_3}_{bot}}{\displaystyle 2i\sin \left(\frac{qL_0}{2}\right)}\sin (qx_3) , \end{gather}$$
$$\begin{gather}\phi \left( \lvert x_3 \rvert \leq \frac{L_0}{2} \right)=\frac{{u_3}_{top}+{u_3}_{bot}}{\displaystyle 2\cos \left(\frac{qL_0}{2}\right)}\cos (qx_3)+\frac{{u_3}_{top}-{u_3}_{bot}}{\displaystyle 2\sin \left(\frac{qL_0}{2}\right)}\sin (qx_3). \end{gather}$$ Substituting the value of 
$q=K \sqrt {({N^2 -\varOmega ^2})/({\varOmega ^2 - f^2})}$ we finally obtain
\begin{align} \phi \left( \lvert x_3 \rvert \leq \frac{L_0}{2} \right)&= \frac{{u_3}_{top}+{u_3}_{bot}}{\displaystyle 2\cos \left(K \sqrt{\frac{N^2 -\varOmega^2}{\varOmega^2 - f^2}}\frac{L_0}{2}\right)}\cos \left(K \sqrt{\frac{N^2 -\varOmega^2}{\varOmega^2 - f^2}}x_3\right) \nonumber\\ &\quad + \frac{{u_3}_{top}-{u_3}_{bot}}{\displaystyle 2\sin \left(K \sqrt{\frac{N^2 -\varOmega^2}{\varOmega^2 - f^2}}\frac{L_0}{2}\right)}\sin \left(K \sqrt{\frac{N^2 -\varOmega^2}{\varOmega^2 - f^2}}x_3\right) . \end{align}Appendix E. Derivation of the equation for $\varOmega _i$ when $N>f$
We differentiate (2.17a) with respect to 
$x_3$ which yields
$$\begin{gather} \frac{\partial{\phi_i}}{\partial x_3}\left( x_3 \geq \frac{L_0}{2} \right)=-{u_3}_{top}\frac{K Y_i \cosh{\left( K Y_i \left(x_3-\dfrac{H}{2} \right) \right)}}{\sinh{\left( K Y_i \left(\dfrac{H-L_0}{2} \right) \right)} }, \end{gather}$$
$$\begin{gather}\frac{\partial{\phi_i}}{\partial x_3}\left( x_3 \geq \frac{L_0}{2} \right)=-{u_3}_{top}\frac{K Y_i \cosh{\left( K Y_i \left(\dfrac{H}{2}-x_3 \right) \right)}}{\sinh{\left( K Y_i \left(\dfrac{H-L_0}{2} \right) \right)} }; \end{gather}$$
here we have used the trigonometric hyperbolic function identity 
$\cosh {(-\theta _i)}=\cosh {(\theta _i)}$. Further calculating the above derivative of 
$\phi _i$ at 
${x_3}=L_0/2$ we get
\begin{equation} \left(\frac{\partial{\phi_i}}{\partial x_3}\right)_{x_3=L_0/2}=-{u_3}_{top}\frac{K Y_i }{\displaystyle \tanh{\left( K Y_i \left(\frac{H-L_0}{2} \right) \right)} } . \end{equation} Differentiating (2.17c) with respect to 
$x_3$ we get
\begin{align} \frac{\partial{\phi_i}}{\partial x_3} \left( \lvert x_3 \rvert \leq \frac{L_0}{2} \right) &= \frac{{u_3}_{top}-{u_3}_{bot}}{\displaystyle 2\sin{\left(K X_i \frac{L_0}{2} \right)} } K X_i \cos{(K X_i x_3 )} \nonumber\\ &\quad +\frac{{u_3}_{top}+{u_3}_{bot}}{\displaystyle 2\cos{\left(K X_i\frac{L_0}{2} \right)} } (-K X_i ) \sin{(K X_i x_3 )} ; \end{align}
again calculating this derivative at 
$\phi _i$ at 
${x_3}=L_0/2$ we get
\begin{align} \left(\frac{\partial{\phi_i}}{\partial x_3}\right)_{x_3=L_0/2} = ({u_3}_{top}-{u_3}_{bot}) \frac{K X_i}{2\tan{\left(K X_i \dfrac{L_0}{2} \right)}} - ({u_3}_{top}+{u_3}_{bot}) \frac{K X_i}{2} \tan{\left(K X_i \frac{L_0}{2} \right)} . \end{align} We equate (E1c) and (E2b) for continuity of 
$\partial \phi _i/\partial x_3$ at 
$x_3=L_0/2$. First, we compare the coefficients of 
${u_3}_{top}$ and get
\begin{equation} \frac{-K Y_i}{\tanh{\left(K Y_i \dfrac{H-L_0}{2}\right)}} = \frac{K X_i}{2\tan\left(K X_i\dfrac{L_0}{2}\right)}-\frac{K X_i}{2} \tan\left(K X_i\frac{L_0}{2}\right) , \end{equation}
then we compare the coefficients of 
${u_3}_{bot}$ and get
$$\begin{gather} 0 =-\frac{K X_i}{2\tan\left(K X_i\dfrac{L_0}{2}\right)}-\frac{K X_i}{2} \tan\left(K X_i\frac{L_0}{2}\right) , \end{gather}$$
$$\begin{gather}\frac{1}{\tan\left(K X_i\dfrac{L_0}{2}\right)}=-\tan\left(K X_i\frac{L_0}{2}\right) . \end{gather}$$Substituting the above (E3c) in (E3a) we get
$$\begin{gather} \frac{-Y_i}{\tanh{\left(K Y_i \dfrac{H-L_0}{2}\right)}} = \frac{-X_i }{2}\tan\left(K X_i\frac{L_0}{2}\right)-\frac{X_i}{2} \tan\left(K X_i\frac{L_0}{2}\right) , \end{gather}$$
$$\begin{gather}\frac{-Y_i}{\tanh{\left(K Y_i \dfrac{H-L_0}{2}\right)}} =- X_i \tan\left(K X_i\frac{L_0}{2}\right) ; \end{gather}$$rearranging the terms we get
\begin{equation} \tan\left(K X_i\frac{L_0}{2}\right) = \frac{Y_i}{X_i}\frac{1}{\displaystyle \tanh{\left(K Y_i \frac{H-L_0}{2}\right)}} . \end{equation} Substituting the 
$X_i$ and 
$Y_i$ from (2.18a,b) in the above (E3f), we get
\begin{equation} \tan\left(K \sqrt{\frac{N^2 - \varOmega_i^2}{\varOmega_i^2 -f^2}}\frac{L_0}{2} \right) =\sqrt{\frac{\varOmega_i^2}{N^2-\varOmega_i^2}} \frac{1}{ \tanh{\left(K \sqrt{\dfrac{\varOmega_i^2}{\varOmega_i^2 -f^2}} \left(\dfrac{H-{L}_0}{2} \right) \right)}} . \end{equation}Now we substitute
\begin{equation} -\tan\left(K X_i\frac{L_0}{2}\right)= \frac{1}{\displaystyle \tan\left(K X_i\frac{L_0}{2}\right)} , \end{equation}from (E3c) in (E3a) and we get
\begin{equation} \frac{-Y_i}{\tanh{\left(K Y_i \dfrac{H-L_0}{2}\right)}} = \frac{X_i}{2\tan\left(K X_i\dfrac{L_0}{2}\right)}+\frac{X_i}{2\tan\left(K X_i\dfrac{L_0}{2}\right)} ; \end{equation}rearranging the terms we get
\begin{equation} \tan\left(K X_i\frac{L_0}{2}\right) =-\frac{X_i}{Y_i} \tanh{\left(K Y_i \frac{H-L_0}{2}\right)} . \end{equation} Further, substituting the 
$X_i$ and 
$Y_i$ from (2.18a,b) in (E5c) we get
\begin{equation} \tan\left(K \sqrt{\frac{N^2 - \varOmega_i^2}{\varOmega_i^2 -f^2}}\frac{L_0}{2} \right) =-\sqrt{\frac{N^2-\varOmega_i^2}{\varOmega_i^2}} \tanh{\left(K \sqrt{\frac{\varOmega_i^2}{\varOmega_i^2 -f^2}} \left(\frac{H-{L}_0}{2} \right) \right)} . \end{equation}Appendix F. Derivation of $\varOmega _0$ for the limiting case of the interface
For the limiting case of the interface 
$KL_0\ll 1$, (2.19b) reduces to
\begin{equation} K \sqrt{\frac{N^2 - \varOmega_0^2}{\varOmega_0^2 -f^2}}\frac{L_0}{2} =\sqrt{\frac{\varOmega_0^2}{N^2-\varOmega_0^2}}. \end{equation}
Here, we utilized a small-angle approximation 
$\tan (x)\sim x$, and assumed an infinite vertical domain height which implies that in (2.19b) 
$\tanh {()}\sim 1$. Squaring both the sides of an (F1) and rearranging the terms, we get
\begin{equation} \varOmega_0^4 (1-Q^2) + (2N^2 Q^2-f^2) \varOmega_0^2 - N^4 Q^2=0, \end{equation}
where 
$Q = {KL_0}/{2}$. Because 
$Q \ll 1$, above (F2) reduces to
\begin{equation} (\varOmega_0^2)^2 + (2N^2 Q^2-f^2) \varOmega_0^2 - N^4 Q^2=0. \end{equation}The roots of above quadratic equation are
$$\begin{gather} \varOmega_0^2 = \frac{-(2N^2Q^2-f^2) \pm \sqrt{(2N^2Q^2-f^2)^2 + 4 N^4 Q^2}}{2}, \end{gather}$$
$$\begin{gather}\varOmega_0^2 = \frac{f^2 - 2N^2Q^2 \pm \sqrt{f^4 - 4N^2Q^2 f^2 + 4 N^4 Q^4 +4 N^4 Q^2 }}{2}. \end{gather}$$
Since, 
$4 N^4 Q^4 \ll 4 N^4 Q^2$, the above (F4b) reduces to
$$\begin{gather} \varOmega_0^2 = \frac{f^2 - 2N^2Q^2 \pm \sqrt{f^4 - 4N^2Q^2 f^2 + 4 N^4 Q^2 }}{2}, \end{gather}$$
$$\begin{gather}\varOmega_0^2 = \frac{f^2 - 2N^2Q^2 \pm \sqrt{(\,f^2+2N^2Q)^2 - 4N^2Qf^2 - 4N^2Q^2 f^2}}{2}. \end{gather}$$
Here, 
$4N^2Q^2f^2 \ll 4N^2Q f^2$, which yields
\begin{equation} \varOmega_0^2 = \frac{f^2 - 2N^2Q^2}{2} \pm \frac{ \overbrace{\sqrt{(\,f^2+2N^2Q)^2 - 4N^2Qf^2}}^{(\mathrm{I})}}{2}.\end{equation}
Substituting 
$Q=KL_0/2$ in term 
$(\mathrm {I})$ of the above (F6), and rearranging the terms we can rewrite term 
$(\mathrm {I})$ as
$$\begin{gather} (\mathrm{I}) = \frac{f^2}{2}\sqrt{\left(1+\left(\frac{N}{f}\right)^2KL_0\right)^2 - 2\left(\frac{N}{f}\right)^2KL_0}, \end{gather}$$
$$\begin{gather}(\mathrm{I}) = \frac{f^2}{2}\left(1+\left(\frac{N}{f}\right)^2KL_0\right) \sqrt{1 - \frac{\displaystyle 2\left(\frac{N}{f}\right)^2KL_0}{\displaystyle \left(1 +\left(\frac{N}{f}\right)^2KL_0\right)^2}}, \end{gather}$$
$$\begin{gather}(\mathrm{I}) = \frac{f^2+N^2KL_0}{2} \sqrt{1 - \frac{\displaystyle 2\left(\frac{N}{f}\right)^2KL_0}{\displaystyle 1+\left(\frac{N}{f}\right)^4 \left(KL_0\right)^2 + 2\left(\frac{N}{f}\right)^2KL_0}}, \end{gather}$$
$$\begin{gather}(\mathrm{I}) = \frac{f^2+N^2KL_0}{2} \sqrt{1 - \frac{2}{ \underbrace{\frac{1}{\displaystyle \left(\frac{N}{f}\right)^2KL_0}}_{(\mathrm{II})} + \underbrace{\left(\frac{N}{f}\right)^2 KL_0}_{(\mathrm{III})} + \underbrace{2}_{(\mathrm{IV})}}}. \end{gather}$$ As 
$f$ increases from 
$0$ (
$\,f=0$ corresponds to the non-rotation case) to 
$N/f>1$, the factor 
$(N/f)^2$ decreases for a constant 
$N$. Moreover, we can take 
$KL_0$ sufficiently small due to the limiting case 
$KL_0 \ll 1$ so that the product 
$({N}/{f})^2KL_0 \ll 1$. Consequently, term 
$(\mathrm {II}) \gg 1$ and term 
$(\mathrm {III}) \ll 1$ of the above (F7d), and therefore we can neglect terms 
$(\mathrm {III})$ and 
$(\mathrm {IV})$ in comparison with term 
$(\mathrm {II})$, which yields
$$\begin{gather} (\mathrm{I}) = \frac{f^2+N^2KL_0}{2} \sqrt{1 - \frac{2}{\displaystyle \frac{1}{\left(\dfrac{N}{f}\right)^2KL_0}}}, \end{gather}$$
$$\begin{gather}(\mathrm{I}) = \frac{f^2+N^2KL_0}{2} \sqrt{1 - \underbrace{2\left(\frac{N}{f}\right)^2 KL_0}_{(\mathrm{V})}}. \end{gather}$$ Since 
$({N}/{f})^2KL_0 \ll 1$, term 
$(\mathrm {V}) \ll 1$ in (F9) and we can neglect term 
$(\mathrm {V})$ in comparison with 
$1$. Therefore, (F9) reduces to
\begin{equation} (\mathrm{I}) = \frac{f^2+N^2KL_0}{2} = \frac{f^2+2N^2Q}{2} . \end{equation}
Substituting term 
$(\mathrm {I})$ (equation (F10)) in (F6), we get
\begin{equation} \varOmega_0^2 = \frac{f^2 - 2N^2Q^2 \pm (\,f^2+2N^2Q)}{2}. \end{equation}The first root of the above equation (with a ‘minus’ sign) is
\begin{equation} \varOmega_0^2 = \frac{f^2 - 2N^2Q^2 - (\,f^2+2N^2Q)}{2}. \end{equation}
As 
$2N^2Q^2 \ll 2N^2Q$, we get
\begin{equation} \varOmega_0^2 =- N^2Q. \end{equation}
This implies that the root 
$\varOmega _0 = \sqrt {- N^2Q}$ is complex, and therefore we discard this root. The second root of (F11) (with a ‘plus’ sign) is
$$\begin{gather} \varOmega_0^2 = \frac{f^2 - 2N^2Q^2 + (\,f^2+2N^2Q)}{2}, \end{gather}$$
$$\begin{gather}\varOmega_0^2 = f^2 + N^2Q = f^2 + \frac{KL_0}{2}N^2 . \end{gather}$$The above expression (F13b) shows that for a sufficiently small 
$KL_0$ and a fixed 
$N$ and 
$f$, 
$\varOmega _0^2$ is close to 
$f^2$ and starts to increase from 
$f^2$ with an increase in the 
$KL_0$. To elucidate this, we consider 
$N=3$ and 
$f=1$ and show the plot of 
$\varOmega _0/N$ with 
$KL_0$ for the limiting case of the interface (
$KL_0\ll 1$), i.e. (F13b) in the inset of a figure 2(a). Therefore, the expression (F13b) is consistent with the solutions of (2.19a) and (2.19b) as shown in figure 2(a) (for a fixed 
$N=3$ and 
$f=1$) demonstrating that all the 
$\varOmega _i$ are very close to 
$f$ (with 
$\varOmega _i>f$) when 
$KL_0$ is very small. Additionally, for 
$f=0$, (F13b) recovers the limiting case of the interface (i.e. 
$\varOmega _0 = \sqrt {N^2 (KL_0/2)} = \sqrt {\mathcal {A}g_0K}$) reported by Briard et al. (Reference Briard, Gostiaux and Gréa2020) for the ‘without rotation’ case.
Appendix G. Derivation of (2.23)
First, we rearrange the (2.21) as follows:
\begin{equation} \varOmega_i^2-f^2 =\frac{K^2}{m^2}(N^2-\varOmega_i^2). \end{equation} Then we substitute the definition of angle 
$\theta _i$ in (G1) and obtain the characteristic frequency 
$\varOmega _i$ in terms of the angle 
$\theta _i$ as follows:
$$\begin{gather} \varOmega_i^2-f^2 =\tan^2(\theta_i)(N^2-\varOmega_i^2) , \end{gather}$$
$$\begin{gather}\varOmega_i^2(1+\tan^2(\theta_i)) = N^2\tan^2(\theta_i)+f^2 , \end{gather}$$
using the trigonometric identity 
$1+\tan ^2(\theta _i)=\sec ^2(\theta _i)=1/\cos ^2{(\theta _i)}$ in the above (G2b) and rearranging the terms we get
$$\begin{gather} \varOmega_i^2 = \frac{N^2\tan^2(\theta_i)+f^2}{1/\cos^2{(\theta_i)}} , \end{gather}$$
$$\begin{gather}\varOmega_i^2 = N^2\frac{\sin^2(\theta_i)}{\cos^2(\theta_i)}\cos^2(\theta_i) + f^2 \cos^2(\theta_i) . \end{gather}$$Appendix H. Derivation of (2.24)
Rearranging the terms of (2.16), we obtain
\begin{equation} \frac{\partial^2 a_i}{\partial t^2} + \left( \varOmega_i^2 + N^2 \left(\frac{\varOmega_i^2 - f^2}{N^2 - f^2}\right) F \cos(\omega t) \right) a_i = 0. \end{equation} Substituting the characteristic frequency 
$\varOmega _i$ from (2.23) in the above (H1) yields
\begin{align} \frac{\partial^2 a_i}{\partial t^2} + \left( N^2\sin^2{(\theta_i)}+f^2\cos^2{(\theta_i)} + \underbrace{ \frac{N^2( N^2\sin^2{(\theta_i)}+f^2\cos^2{(\theta_i)} - f^2)}{N^2 - f^2}}_{(\mathrm{I})} F \cos{(\omega t)} \right) a_i = 0 . \end{align} Now we simplify the term 
$(\mathrm {I})$ of above (H2) as follows:
$$\begin{gather} (\mathrm{I}) = \frac{N^2( N^2\sin^2{(\theta_i)}+f^2\cos^2{(\theta_i)} - f^2)}{N^2 - f^2} , \end{gather}$$
$$\begin{gather}(\mathrm{I}) = \frac{N^2( N^2\sin^2{(\theta_i)}+f^2\cos^2{(\theta_i)} + (\,f^2\sin^2{(\theta_i)}-f^2\sin^2{(\theta_i)} ) - f^2)}{N^2 - f^2} , \end{gather}$$
$$\begin{gather}(\mathrm{I}) = \frac{N^2( N^2\sin^2{(\theta_i)}+f^2(\cos^2{(\theta_i)} + \sin^2{(\theta_i)})-f^2\sin^2{(\theta_i)}- f^2)}{N^2 - f^2} , \end{gather}$$
using the Pythagorean trigonometric identity 
$\cos ^2{(\theta _i)} + \sin ^2{(\theta _i)} =1$ , we get
$$\begin{gather} (\mathrm{I}) = \frac{N^2( N^2\sin^2{(\theta_i)}+f^2-f^2\sin^2{(\theta_i)}-f^2)}{N^2 - f^2} , \end{gather}$$
$$\begin{gather}(\mathrm{I}) = \frac{N^2( N^2\sin^2{(\theta_i)}-f^2\sin^2{(\theta_i)})}{N^2 - f^2} . \end{gather}$$ Now we substitute the term 
$(\mathrm {I})$ from above (H3e) in (H2), and we get
\begin{align} \frac{\partial^2 a_i}{\partial t^2} + \left( N^2\sin^2{(\theta_i)}+f^2\cos^2{(\theta_i)} + \frac{N^2( N^2\sin^2{(\theta_i)}-f^2\sin^2{(\theta_i)})}{N^2 - f^2} F \cos{(\omega t)} \right) a_i = 0 . \end{align}After doing some simplifications, we get
\begin{align} &\frac{\partial^2 a_i}{\partial t^2} + \left( \frac{f^2 (N^2-f^2)\cos^2{(\theta_i)}}{N^2 - f^2} \right.\nonumber\\ &\quad \left. + \frac{(N^2)^2\sin^2{(\theta_i)}(1+F \cos{(\omega t)} ) - f^2 N^2\sin^2{(\theta_i)}(1+F \cos{(\omega t)} )}{N^2 - f^2} \right) a_i = 0 , \end{align}
\begin{align} &\frac{\partial^2 a_i}{\partial t^2} + \left( f^2 \cos^2{(\theta_i)} + \frac{N^2 (N^2-f^2)\sin^2{(\theta_i)}}{{N^2 - f^2}}(1+F \cos{(\omega t)}) \right)a_i = 0 . \end{align}Appendix I. Derivation of (2.26)
We substitute 
$\varOmega _0^2 \simeq f^2 + N^2 (KL_0/2)$ from the interface limit 
$KL_0 \ll 1$ in the general Mathieu equation (2.16) and replace 
$a_i$ with 
$a$ for clarity. This gives
\begin{align} \frac{\partial^2 a}{\partial t^2} + \left( f^2 \left(\frac{N^2 - f^2 - N^2 (KL_0/2)}{N^2 - f^2}\right) + N^2 \left(\frac{f^2 + N^2 (KL_0/2) - f^2}{N^2 - f^2}\right) (1+F \cos(\omega t) ) \right) a = 0. \end{align} Since 
$KL_0 \ll 1$, 
$N^2 (KL_0/2) \ll N^2$, then (I1) reduces to
$$\begin{gather} \frac{\partial^2 a}{\partial t^2} + \left( f^2 \left(\frac{N^2 - f^2}{N^2 - f^2}\right) + N^2 \left(\frac{N^2 (KL_0/2)}{N^2 - f^2}\right) (1+F \cos(\omega t) ) \right) a = 0. \end{gather}$$
$$\begin{gather}\frac{\partial^2 a}{\partial t^2} + \left( f^2 + \frac{N^4 (KL_0/2)}{N^2 - f^2} (1+F \cos(\omega t) ) \right) a = 0. \end{gather}$$ Substituting 
$N^2=2\mathcal {A}g_0/L_0$ in the above (I2b) and rearranging the terms, we get
\begin{equation} \frac{\partial^2 a}{\partial t^2} + \left( f^2 + \frac{2\mathcal{A}^2g_0^2K }{{\underbrace{2\mathcal{A}g_0 - L_0f^2}}_{(\mathrm{I})}} (1+F \cos(\omega t) ) \right) a = 0. \end{equation} This (I3) can be further simplified by assuming that the interface between two fluids is well-defined (i.e. no wave/interface breaking) and that the interface is very sharp at the onset of instability. Accordingly, we can take a small value of 
$K$ for a well-defined interface and 
$L_0 \ll 1$ for a sharp interface. This implies that 
$KL_0 \ll 1$ which represents the limiting case of the interface. The assumption of the sharp interface is reasonable because, for the non-rotating system, Briard et al. (Reference Briard, Gostiaux and Gréa2020) has recovered the limiting case of the interface (
$KL_0 \ll 1$) reported by Benjamin & Ursell (Reference Benjamin and Ursell1954) for immiscible fluids with sharp interfaces. Therefore, for a moderate value of 
$f$ the product 
$L_0f^2\ll 2\mathcal {A}g_0$ or 
$f^2 \ll N^2$ (since 
$N^2=2\mathcal {A}g_0/L_0$). Here, we used the assumption of small density contrast 
$\mathcal {A} \approx 0.01$, 
$g_0 \approx 10\ \mathrm {ms}^{-2}$ which gives 
$2\mathcal {A}g_0 = 0.2\ \mathrm {ms}^{-2}$. We can verify the approximation 
$L_0f^2\ll 2\mathcal {A}g_0$ by taking 
$L_0=0.003\ \mathrm {m}$ from the experimental data of Cavelier et al. (Reference Cavelier, Gréa, Briard and Gostiaux2022) and 
$f=0.322\ \mathrm {s}^{-1}$ from the present numerical simulations (see table 1 for case F075f/
$\omega$48), which yields 
$L_0f^2\approx 0.0003\ \mathrm {ms}^{-2}$. This implies that 
$L_0f^2\ll 2\mathcal {A}g_0$ 
$(0.0003\ll 0.2)$ and we can neglect 
$L_0f^2$ compared with 
$2\mathcal {A}g_0$. This also holds true for the other 
$f$ values given in table 1. Therefore, using this approximation 
$L_0f^2\ll 2\mathcal {A}g_0$ in term 
$(\mathrm {I})$ of (I3), we obtain
\begin{equation} \frac{\partial^2 a}{\partial t^2} + (\,f^2 + \mathcal{A}g_0K (1+F \cos(\omega t) ) ) a = 0. \end{equation}
Substituting 
$\tau =\omega t$ in (I4) yields the Mathieu-like (2.26).
Appendix J. Steps for solving (2.9) when $N< f$
We rearrange (2.9) as
\begin{equation} \frac{\partial^2\phi}{\partial x_3^2} + {K}^2 \frac{( \varOmega^2 - N^2 )}{(\,f^2 - \varOmega^2) } \phi = 0. \end{equation}
Here, the same solution strategy, boundary conditions and piecewise background concentration profile is used as for 
$N>f$ case in Appendix D. We first solve (J1) for 
$N< f$ in the upper pure fluid regime 
$x_3 \geq {L_0}/{2}$ where 
$\varGamma =0$ and hence 
$N=\sqrt {-2\mathcal {A} g_0 \varGamma }=0$. Therefore, (J1) becomes
\begin{equation} \frac{\partial^2\phi}{\partial x_3^2} + K^2\frac{\varOmega^2}{(\,f^2-\varOmega^2 )} \phi = 0. \end{equation}The general solution of (J2) is
\begin{equation} \phi= C\cos \left(K\sqrt{\frac{\varOmega^2}{f^2 - \varOmega^2}}x_3\right)+iD\sin \left(K\sqrt{\frac{\varOmega^2}{f^2 - \varOmega^2}}x_3\right), \end{equation}
where 
$C$ and 
$D$ are the unknown constants. These constants are obtained by using the boundary condition 
$\phi =0$ at 
$x_3=H/2$ and 
$\phi ={u_3}_{top}$ at 
$x_3=L_0/2$, which yields the following solution of (J3):
\begin{equation} \phi \left( x_3 \geq \frac{L_0}{2} \right) =-{u_3}_{top} \frac{ \sin{\left( K\sqrt{\dfrac{\varOmega^2}{f^2 - \varOmega^2}} \left(x_3-\dfrac{H}{2} \right) \right)}} {\sin{\left( K\sqrt{\dfrac{\varOmega^2}{f^2 - \varOmega^2}} \left(\dfrac{H-L_0}{2} \right) \right)}}. \end{equation} Following the above same procedure, we obtain the solution for (J1) in the bottom pure fluid region 
$x_3=-L_0/2$ to 
$x_3=-H/2$, where 
$\varGamma =0$, thus 
$N=\sqrt {-2\mathcal {A} g_0 \varGamma }=0$ and boundary condition 
$\phi ={u_3}_{bot}$ at 
$x_3=-L_0/2$ and 
$\phi =0$ at 
$x_3=-H/2$. Then the final solution is given as
\begin{equation} \phi \left( x_3 \leq-\frac{L_0}{2} \right) = {u_3}_{bot} \frac{\displaystyle \sin{\left( K\sqrt{\frac{\varOmega^2}{f^2 - \varOmega^2}} \left(x_3+\frac{H}{2} \right) \right)}}{\displaystyle\sin{\left( K\sqrt{\frac{\varOmega^2}{f^2 - \varOmega^2}} \left(\frac{H-L_0}{2} \right) \right)}} . \end{equation} In the mixed fluid region between 
$x_3=L_0/2$ to 
$x_3=-L_0/2$, we consider 
$\varGamma =-N^2/2\mathcal {A}g_0$. Therefore, the general solution of (J1) is given as
\begin{equation} \phi= A\cos \left(K \sqrt{\frac{\varOmega^2 - N^2}{f^2 - \varOmega^2}}x_3\right)+iB\sin \left(K \sqrt{\frac{\varOmega^2 - N^2}{f^2 - \varOmega^2}}x_3\right) , \end{equation}
where 
$A$ and 
$B$ are the unknown constants. We find out these constants by applying the boundary conditions 
$\phi ={u_3}_{top}$ at 
$x_3=+L_0/2$ and 
$\phi ={u_3}_{bot}$ at 
$x_3=-L_0/2$ to (J6). After substituting the expressions for 
$A$ and 
$B$ in (J6), we get
\begin{align} \phi \left( \lvert x_3 \rvert \leq \frac{L_0}{2} \right)&= \frac{{u_3}_{top}+{u_3}_{bot}}{\displaystyle 2\cos \left(K \sqrt{\frac{\varOmega^2 - N^2}{f^2-\varOmega^2}}\frac{L_0}{2}\right)}\cos \left(K \sqrt{\frac{\varOmega^2-N^2}{f^2-\varOmega^2 }}x_3\right) \nonumber\\ &\quad + \frac{{u_3}_{top}-{u_3}_{bot}}{\displaystyle 2\sin \left(K \sqrt{\frac{\varOmega^2-N^2}{f^2-\varOmega^2}}\frac{L_0}{2}\right)}\sin \left(K \sqrt{\frac{\varOmega^2-N^2}{f^2 - \varOmega^2}}x_3\right) . \end{align}Appendix K. Floquet analysis
In this appendix, we briefly discuss the theorems related to Floquet theory, and provide the steps to solve the Mathieu equation, (2.25), by using these theorems. Proofs of the theorems and solution steps are given in Jordan & Smith (Reference Jordan and Smith2007, pp. 308–317).
Theorem K.1 (Floquet's theorem)
Let 
${\mathrm{d}\kern0.7pt \boldsymbol {x}}/{\mathrm {d}t}=\boldsymbol{\mathsf{A}}(t)\boldsymbol {x}$, be a first-order linear differential system, where 
$\boldsymbol{\mathsf{A}}(t)$ is 
$T$-periodic matrix such that 
$\boldsymbol{\mathsf{A}}(t+T)=\boldsymbol{\mathsf{A}}(t), \forall \:t$. This system has at least one non-trivial solution 
$\boldsymbol {x}=\boldsymbol {\chi }(t)$ such that
\begin{equation} \boldsymbol{\chi}(t+T)=\mu \boldsymbol{\chi}(t),\quad \forall \, t, \end{equation}
where 
$\mu$ is a characteristic number or Floquet multiplier.
If 
$\boldsymbol{\mathsf{\phi}} (t)$ is a fundamental solution matrix of system 
$\boldsymbol {\dot x}=\boldsymbol{\mathsf{A}}(t) \boldsymbol {x}$, then 
$\boldsymbol{\mathsf{\phi}} (t+T)$ is also a fundamental matrix, and there exists a non-singular matrix 
$\boldsymbol{\mathsf{C}}$ such that
$$\begin{gather} \boldsymbol{\mathsf{\phi}}(t+T)=\boldsymbol{\mathsf{C}}\boldsymbol{\mathsf{\phi}}(t),\quad \forall \,t, \end{gather}$$
$$\begin{gather}\boldsymbol{\mathsf{C}}=\boldsymbol{\mathsf{\phi}}^{-1}(t)\boldsymbol{\mathsf{\phi}}(t+T), \end{gather}$$
where characteristic numbers, 
$\mu '\text {s}$, are the eigenvalues of 
$\boldsymbol{\mathsf{C}}$. We can obtain the matrix 
$\boldsymbol{\mathsf{C}}=\boldsymbol{\mathsf{\phi}} (T)$ for initial conditions at 
$t=0$ such that 
$\boldsymbol{\mathsf{\phi}} (0)=\boldsymbol{\mathsf{I}}$, where 
$\boldsymbol{\mathsf{I}}$ is the identity matrix.
Theorem K.2 For system 
$\boldsymbol {\dot x}=\boldsymbol{\mathsf{A}}(t) \boldsymbol {x}$, where 
$\boldsymbol{\mathsf{A}}(t)$ is 
$T$-periodic, with characteristic numbers 
$\mu _1,\mu _2,\ldots, \mu _n$, the product of the characteristic numbers is obtained as
\begin{equation} \mu_1\mu_2\dots\mu_n=\exp \left(\int_0^T \mathrm{tr} \{\boldsymbol{\mathsf{A}}(t)\}\right). \end{equation} Now, we define the characteristic exponent, 
$\sigma$ of the system as 
${\rm e}^{\sigma T}=\mu$.
Theorem K.3 Let 
$\boldsymbol{\mathsf{C}}$ have 
$n$ distinct eigenvalues, 
$\mu _i$ and corresponding 
$\sigma _i, i=1,2,\ldots,n$. Then system 
$\boldsymbol {\dot x}=\boldsymbol{\mathsf{A}}(t) \boldsymbol {x}$ has 
$n$ linearly independent solutions of the form
\begin{equation} \boldsymbol{\chi_i}(t)={\rm e}^{\sigma_i t}\boldsymbol{P_i}(t), \end{equation}
here 
$\boldsymbol {P_i}(t)$ are the periodic vector functions with period 
$T$.
Clearly, characteristic exponents 
${\rm e}^{\sigma _i t}$, will determine the behaviour of the solutions (K4).
Now, we apply above theorems to study the nature of solutions of Mathieu (2.25). We can rewrite (2.25) as
\begin{equation} \ddot{a}+( \gamma +\alpha + \beta \cos{\tau} )a=0, \end{equation}where,
\begin{equation} \gamma = \frac{f^2 \cos^2{(\theta)}}{\omega^2};\quad \alpha = \frac{N^2 \sin^2{(\theta)}}{\omega^2};\quad \beta = F\frac{N^2 \sin^2{(\theta)}}{\omega^2}. \end{equation} We begin by defining 
$a=X$ and 
$\dot{a}=Y$, so (K5) can be expressed as a first-order system,
\begin{equation} \left[ \begin{array}{@{}c@{}} \dot{X} \\ \dot{Y} \\ \end{array} \right] = \left[ \begin{array}{@{}cc@{}} 0 & 1 \\ -(\gamma+\alpha+\beta\cos{\tau}) & 0 \\ \end{array} \right] \left[ \begin{array}{@{}c@{}} X \\ Y \\ \end{array} \right].\end{equation}In the notation of Theorem K.1,
\begin{equation} \boldsymbol{\mathsf{A}}(\tau) = \left[ \begin{array}{@{}cc@{}} 0 & 1 \\ -(\gamma+\alpha+\beta\cos{\tau}) & 0 \\ \end{array} \right] .\end{equation} Here, 
$\boldsymbol{\mathsf{A}}(\tau )$ is periodic with period 
$2{\rm \pi}$ and 
$\mathrm {tr}\{\boldsymbol{\mathsf{A}}(\tau )\}=0$. As our system is of size 
$2\times 2$, so we have two characteristic numbers 
$\mu _1$ and 
$\mu _2$. Therefore, from Theorem K.2
\begin{equation} \mu_1\mu_2=\mathrm{e}^0=1 . \end{equation}For the following initial condition:
\begin{equation} \boldsymbol{\mathsf{\phi}}(0) = \left[ \begin{array}{@{}cc@{}} 1 & 0 \\ 0 & 1 \\ \end{array} \right] ,\end{equation}
we obtain the matrix 
$\boldsymbol{\mathsf{C}}=\boldsymbol{\mathsf{\phi}} (T)$ by numerically integrating (with MATLAB) the system (K7) from 
$\tau =0$ to 
$\tau =2{\rm \pi}$. The eigenvalues or characteristic numbers are the solutions of quadratic characteristic equation of 
$\boldsymbol{\mathsf{C}}$
\begin{equation} \mu^2-({sum\ of\ roots})\mu + {(\,product\ of\ roots)}=0, \end{equation}and using the (K9), we get
\begin{equation} \mu^2-\phi(\gamma,\alpha,\beta)\mu + 1=0, \end{equation}
here, 
$\phi (\gamma,\alpha,\beta )$ represents the sum of roots. The solutions 
$\mu _1$ and 
$\mu _2$ of the (K12) are
\begin{equation} \mu_1,\mu_2=\tfrac{1}{2}(\phi\pm\sqrt{(\phi^2-4)}). \end{equation} Thus, different values of 
$\phi$ will determine the behaviour of solutions as follows.
(i) When
$\phi >2$: $\mu _1$, $\mu _2$ are both real and positive with $\mu _1\mu _2=1$ (from (K9)), one of them, say $\mu _1>1$ and $\mu _2<1$. The corresponding characteristic exponents are real and have the form $\sigma _1=\xi >0$, $\sigma _2=-\xi <0$ (because $\sigma _2=\ln ({\mu _2})/T<0$). So, the general solution from Theorem K.3 can be written as
(K14)where\begin{equation} X(\tau)=c_1\,{\rm e}^{\xi \tau}P_1(\tau)+c_2\,{\rm e}^{-\xi t}P_2(\tau), \end{equation}$\xi >0$, thus $\mathrm {e}^{\xi \tau }$ corresponds to the exponential growth in time and solution becomes unstable with harmonic in nature.(ii) When
$\phi <-2$: $\mu _1$, $\mu _2$ are both real and negative with $\mu _1\mu _2=1$. The form of general solution is
(K15)Here,\begin{equation} X(\tau)=c_1\,{\rm e}^{\left(\xi+\frac{1}{2}\mathrm{i}\right) \tau}P_1(\tau)+c_2\,{\rm e}^{\left(-\xi +\frac{1}{2}\mathrm{i}\right) \tau}P_2(\tau). \end{equation}$\mathrm {e}^{\xi \tau }$ leads to the exponential growth of solution and eventually becomes unstable, with subharmonic response.(iii) When
$-2<\phi <+2$: $\mu _1$, $\mu _2$ are complex, we have $\sigma _1=\mathrm {i}\nu$ and $\sigma _2=-\mathrm {i}\nu$, thus we get the solution
(K16)The solutions are bounded at all times and oscillatory but not necessarily periodic, and lies in stable parameter region.\begin{equation} X(\tau)=c_1\,{\rm e}^{(\mathrm{i}\nu) \tau}P_1(\tau)+c_2\,{\rm e}^{(-\mathrm{i}\nu) \tau}P_2(\tau). \end{equation}(iv) When
$\phi =+2$: $\mu _1=\mu _2=1$, one solution is periodic with period $2{\rm \pi}$ and lies on the harmonic tongue of stability curve $\phi =2$ .(v) When
$\phi =-2$: $\mu _1=\mu _2=-1$, one solution is periodic with period $4{\rm \pi}$ and lies on the subharmonic tongue of stability curve $\phi =2$ .
Appendix L. The WKB approximation of (2.20)
We consider the (2.21) and assume that 
$N(x_3)=\sqrt {-2\mathcal {A}g_0(\partial \langle C \rangle _H/\partial x_3)}$ (see (2.6)) is a slowly varying function such that its fractional change over a vertical wavelength is much less than unity, therefore 
$m(x_3)$ is also a slowly varying function. With such slowly varying properties of the medium (here 
$N$), we can approximate the WKB solution (for details see Kundu, Cohen & Dowling Reference Kundu, Cohen and Dowling2015, pp. 743–745) of (2.20) by assuming a solution of the form
\begin{equation} \phi=A(x_3)\,\mathrm{e}^{{\rm i}\psi(x_3)}, \end{equation}
where 
$A$ represents the amplitude and 
$\psi$ the phase. We substitute this solution in (2.20) and equate the real and imaginary parts, and after doing some simplifications (see Kundu et al. Reference Kundu, Cohen and Dowling2015, pp. 743–745), we get
\begin{equation} A = \frac{A_0}{\sqrt{m}}; \quad \psi={\pm} i\int^{x_3} m \,\mathrm{d}\kern0.7pt x_3, \end{equation}
where 
$A_0$ is a constant. Therefore, the WKB solution becomes
\begin{equation} \phi=\frac{A_0}{\sqrt{m}}\,\mathrm{e}^{{\pm} \int^{x_3} m \,\mathrm{d}\kern0.7pt x_3}. \end{equation} Substitution of this WKB solution (L3) in the assumed solution of 
$u_3$ (from (2.8)) yields
\begin{equation} u_3=\frac{A_0}{\sqrt{m}}\,\mathrm{e}^{{\rm i} (kx_1+lx_2 \pm \int^{x_3} m \,\mathrm{d}\kern0.7pt x_3-\varOmega t )}. \end{equation}Now, taking real parts of the above (L4), we obtain the vertical velocity
\begin{equation} u_3=\frac{A_0} {\sqrt{m}}\cos{\left(kx_1+lx_2\pm\int^{x_3} m \,\mathrm{d}\kern0.7pt x_3 -\varOmega t\right)}. \end{equation} Here, the argument of 
$\cos {()}$ represents the ‘phase’ and 
$\partial (\,phase)/\partial {x_3}=m$, therefore 
$m$ is the local vertical wavenumber in vertical direction 
$x_3$ and we can define the dispersion relation as
\begin{equation} m^2=\frac{k^2(N^2-\varOmega^2)}{(\varOmega^2-f^2)}. \end{equation}
Open access
