1. Introduction
Rotating Rayleigh–Bénard convection (RRBC) is a relevant model in many geo- and astrophysical phenomena. In the case of rapid rotation, the bulk of the system is dominated by a balance between Coriolis and pressure gradient forces known as geostrophy. In the ‘upright’ case, meaning gravity and the rotation axis are aligned, stress-free boundary conditions are compatible with geostrophy and no boundary layer occurs at leading order. However, if no-slip boundary conditions are imposed, the dominant balance near the boundary is between viscous and Coriolis forces, in what is known as the Ekman layer. Conservation of mass in this layer gives rise to the physical phenomenon called Ekman pumping (suction). In computation, the Ekman layer is narrow and computationally expensive to resolve. Stress-free boundary conditions do not break down the geostrophy of the interior, and thus are often chosen for computational convenience. However, many studies have been done to explore how the Ekman layer does effect the flow. Chandrasekhar (Reference Chandrasekhar1961) computes critical Rayleigh numbers for convection onset for both cases (stress-free and no-slip) and finds that these values differ for non-zero Ekman numbers (
$E$). Asymptotic studies (Niiler & Bisshopp Reference Niiler and Bisshopp1965; Heard & Veronis Reference Heard and Veronis1971; Zhang & Jones Reference Zhang and Jones1993) show that while non-zero Ekman numbers have a destabilizing effect, the difference is asymptotically small in the limit of rapid rotation, 
$E \to 0$. Studies of nonlinear effects find an increase in heat transport which does not vanish with 
$E$ (Homsy & Hudson Reference Homsy and Hudson1971; Stellmach et al. Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014; Julien et al. Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016; Pachev et al. Reference Pachev, Whitehead, Fantuzzi and Grooms2020). In our recent publication (Tro, Grooms & Julien Reference Tro, Grooms and Julien2024), we discover that the marginal stability curve of the system with Ekman pumping departs significantly from the stress-free curve at small wavenumbers. We find that this effect does not vanish with 
$E\to 0$. In this paper, we describe that departure with analytical, asymptotic solutions.
2. Formulation
RRBC is governed by the incompressible Navier–Stokes equations (iNSEs). Non-dimensionalization by characteristic velocity scale 
$U$, horizontal length scale 
$\ell$, vertical depth scale 
$H$, advective time scale 
$\ell /U$, pressure scale 
$P$ and temperature difference 
$\Delta T$ results in non-dimensional parameters including the Ekman number, 
$E$, and the reduced thermal Rayleigh number, 
$\widetilde {Ra} = Ra\,E^{4/3}$. We also have a Rossby number, 
$Ro$, which is small in the rapidly rotating limit, so we define 
$\varepsilon$ such that
\begin{equation} Ro=E^{1/3}\equiv\varepsilon\ll 1. \end{equation}
We further select 
$U = \nu /\ell$, where 
$\nu$ is the kinematic viscosity, and 
$\ell / H = \varepsilon$. The resulting non-dimensional iNSE system is given by
\begin{gather} \partial_{t}\boldsymbol{u} +\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} +\frac{1}{\varepsilon}\ \hat{\boldsymbol{z}}\times \boldsymbol{u} + \frac{1}{\varepsilon} \boldsymbol{\nabla} p = \nabla^2 \boldsymbol{u} +\frac{\widetilde{Ra}}{\sigma} \theta \hat{\boldsymbol{z}}, \end{gather}
\begin{gather}\partial_{t}\theta +\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}\theta -\varepsilon \hat{\boldsymbol{z}}\boldsymbol{\cdot}\boldsymbol{u} = \frac{1}{\sigma}\nabla^2 \theta, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}
where 
$\boldsymbol {u}$ and 
$p$ are the velocity and pressure fields, and 
$\theta$ is the temperature with the linear background profile removed, that is, 
$\theta = T +z-1$, where 
$T$ is the temperature field. The parameter 
$\sigma =\nu /\kappa$ is the Prandtl number, where 
$\kappa$ is the thermal diffusivity.
We consider a local 
$f$-plane approximation to the spherical shell, and focus on the case where gravity and the rotation axis are aligned. The unit vectors 
$\hat {\boldsymbol {x}}$, 
$\hat {\boldsymbol {y}}$, 
$\hat {\boldsymbol {z}}$, define the coordinate directions and point east, north and radially. The components of 
$\boldsymbol {u}$ are 
$(u,v,w)$ and are respectively in the 
$\hat {\boldsymbol {x}}$, 
$\hat {\boldsymbol {y}}$, 
$\hat {\boldsymbol {z}}$ directions.
The system (2.2) is accompanied by the boundary conditions
\begin{equation} w = 0,\quad \theta = 0,\quad Z = 0,1, \end{equation}
where 
$Z$ is the non-dimensional coordinate in the 
$\hat {\boldsymbol {z}}$ direction, and either
\begin{gather} \text{no-slip:} \quad (u,v) = 0 \quad Z=0,1, \end{gather}
\begin{gather}\text{stress-free:}\quad \hat{\boldsymbol{z}} \boldsymbol{\cdot} \boldsymbol{\nabla} (u,v) = 0\quad Z = 0,1. \end{gather}Using asymptotic matching, the fluid variables of the iNSE may be decomposed into an interior (outer) geostrophic solution, a middle thermal wind layer and boundary (inner) Ekman layer solution. For no-slip boundary conditions, the result for the interior (outer) solution plus the middle thermal wind layer is the composite quasi-geostrophic (CQG) equations given by Julien et al. (Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016),
\begin{gather} \partial_{t}\nabla_\perp^2\psi +J\big[\psi,\nabla_\perp^2\psi\big] -\partial_{Z}w = \nabla_\perp^4 \psi, \end{gather}
\begin{gather}\partial_{t} w +J\big[\psi,w\big] +\partial_{Z} \psi = \nabla_\perp^2 w +\frac{\widetilde{Ra}}{\sigma} \theta', \end{gather}
\begin{gather}\partial_{t} \theta' +J\big[\psi, \theta'\big] +w\big(\partial_{Z}\bar{\varTheta}-1\big) = \frac{1}{\sigma}\nabla^2 \theta', \end{gather}
\begin{gather}w = \delta \frac{\varepsilon^{1/2}}{\sqrt{2}}\nabla_\perp^2 \psi,\quad \theta = 0,\quad\text{on}\ Z = 0,1, \end{gather}
where 
$\delta = 1$ at 
$Z=0$ and 
$\delta = -1$ at 
$Z=1$. Here, 
$\psi$ is the stream function such that 
$\partial _{x}\psi = v$, 
$\partial _{y}\psi = -u$ and 
$\theta = \bar {\varTheta }+\varepsilon \theta '$ has been split into a lateral mean (denoted by the overbar) and fluctuating component (denoted by the prime). The operators 
$\boldsymbol {\nabla } = (\partial _{x},\partial _{y},\varepsilon \partial _{Z})$, 
$\nabla _\perp =(\partial _{x},\partial _{y})$ and 
$J[\psi,{\cdot } ]= \partial _{x}\psi \partial _{y} - \partial _{y}\psi \partial _{x}$. The boundary condition (2.5d) on 
$w$ is known as the pumping boundary condition since it encapsulates the pumping/suction effect of the Ekman boundary layer on the interior solution. Note that in this model, the vertical diffusion of 
$\theta '$ is retained due to the inclusion of the 
${{O}}(\varepsilon ) = {{O}}(E^{1/3})$ thermal wind layer. Without the vertical derivatives, the thermal boundary conditions cannot be satisfied.
We explore the linear stability of both the rescaled iNSE (2.2) and the CQG model (2.5) by substituting the normal mode ansatz
\begin{equation} \boldsymbol{v} = \breve{\boldsymbol{v}}(Z)\exp( s t + {\rm i} \boldsymbol{k}_\perp\boldsymbol{\cdot}\boldsymbol{x}_\perp) \end{equation}
for convective rolls, where we define the wavenumber 
$\boldsymbol {k}_\perp = (k_x,k_y)$. Inputting parameters 
$\varepsilon$, 
$\sigma$, 
$k_\perp = |\boldsymbol {k}_\perp |$ and 
$s=0$ for marginal stability, the result is a differential equation eigenvalue problem for 
$\widetilde {Ra}$. The iNSE problem is then solved numerically (see Tro et al. Reference Tro, Grooms and Julien2024 for the numerical formulation), and we explore asymptotic solutions to the CQG linear stability problem. While our following analysis is general for 
$\sigma = {{O}}(1)$, note that by solving the problem with 
$s=0$, we are neglecting oscillatory convection for 
$\sigma <1$.
Assuming 
$w_0\sim \psi _0\sim \theta _0 \sim k_\perp ^2 \sim \widetilde {Ra}_0\sim 1$, where the subscript denotes the leading order term in an asymptotic expansion, the pumping boundary condition is subdominant at 
${{O}}(\varepsilon ^{1/2})$. This means that the boundary condition at leading order is 
$w_0 = 0$ at 
$Z = 0,1$. This yields the linear stability problem given by Julien & Knobloch (Reference Julien and Knobloch1998), with solution
\begin{gather} w = A\sin(n{\rm \pi} Z)+ {{O}}\big(\varepsilon^{1/2}\big), \end{gather}
\begin{gather}\psi ={-}A\frac{n{\rm \pi}}{k_\perp^4}\cos(n{\rm \pi} Z)+ {{O}}\big(\varepsilon^{1/2}\big), \end{gather}
\begin{gather}\theta = A\frac{\sigma}{k_\perp^2}\sin(n{\rm \pi} Z)+ {{O}}\big(\varepsilon^{1/2}\big), \end{gather}and eigenvalue
\begin{equation} \widetilde{Ra} \sim \frac{k_\perp^6+n^2{\rm \pi}^2}{k_\perp^2}, \end{equation}
which is the same as found by Chandrasekhar (Reference Chandrasekhar1961). The curve (2.8) in the 
$k_\perp$–
$\widetilde {Ra}$ plane is the marginal stability curve, meaning that for parameter values above the curve, there is a positive growth rate 
$s$ and the system is unstable to perturbation. For stress-free boundary conditions in the rapidly rotating limit 
$\varepsilon \to 0$, the marginal stability curve of the iNSE problem will converge to (2.8). However, we discovered that (2.8) does not describe the marginal stability of the no-slip problem at low wavenumbers 
$k_\perp$ (Tro et al. Reference Tro, Grooms and Julien2024). This numerical result is reproduced in figure 1 by the solid blue line, contrasted with the stress-free result shown by the blue dash-dotted line.
Marginal stability of the system in the 
$\widetilde {Ra}$–
$k_\perp$ plane for 
$\varepsilon = 10^{-4}$ and 
$\sigma = 1$. The numerically computed curve for the iNSE system with Ekman pumping is shown in solid blue. The blue dash-dotted curve is the stress-free case, 
$\widetilde {Ra} = k_\perp ^4+n^2{\rm \pi} ^2/k_\perp ^2$. The asymptotic approximations for Cases (i)–(iv) are in red, purple, yellow and green, respectively.
In this paper, we seek an analytical description of the marginal stability curve when Ekman pumping is present, shown numerically by Tro et al. (Reference Tro, Grooms and Julien2024) and plotted in figure 1. The pumping boundary condition remains 
${{O}}(1)$ when 
$w\sim \varepsilon ^{1/2}k_\perp ^2\psi \sim 1$, (
$\psi \sim k_\perp ^{-2}\varepsilon ^{-1/2}$). Under this scaling, we expand all dependent variables and 
$\widetilde {Ra}$ in half-powers of 
$\varepsilon$, and consider only the leading order (
${{O}}(1)$) solution. Four cases occur: (i) 
$k_\perp ^2\sim \varepsilon ^{1/2}$; (ii) 
$\varepsilon ^{1/2}\gtrsim k_\perp ^2 \gtrsim \varepsilon ^2$; (iii) 
$k_\perp ^2 \sim \varepsilon ^2$; and (iv) 
$k_\perp ^2 \lesssim \varepsilon ^2$, each of which is outlined in § 3 below, and each of which describes a section of the marginal stability curve.
3. Results
To ensure pumping remains 
${{O}}(1)$ as we explore 
$k_\perp ^2 \lesssim \varepsilon ^{1/2}$, we begin by rescaling 
$\psi \to k_\perp ^{-2}\varepsilon ^{-1/2}\check {\psi }$, where 
$\check {\psi }\sim 1$. Note that throughout the manuscript, we use the inverted hat (check) to denote the rescaled quantities, including on their asymptotic expansions. Applying this scaling to the linear problem associated with (2.5) (dropping the prime on the fluctuating temperature), yields
\begin{gather} -\partial_{Z} w= k_\perp^2\varepsilon^{{-}1/2} \check{\psi}, \end{gather}
\begin{gather}\partial_{Z}\check{\psi}= k_\perp^2 \varepsilon^{1/2}\frac{\widetilde{Ra}}{\sigma}\theta - k_\perp^4\varepsilon^{1/2} w, \end{gather}
\begin{gather}- w = \frac{1}{\sigma}\big(\varepsilon^2 \partial_{Z}^2 -k_\perp^2\big)\theta, \end{gather}with
\begin{equation} w ={-}\frac{\delta}{\sqrt{2}} \check{\psi},\quad \theta =0,\quad Z = 0,1. \end{equation}
Under this scaling, the horizontal diffusion of 
$w$, the term 
$k_\perp ^4 \varepsilon ^{1/2} w$, will always be less than 
${{O}}(\varepsilon ^{1/2})$ for the 
$k_\perp ^2<1$ cases considered below. We expand the eigenfunctions in powers of 
$\varepsilon ^{1/2}$, but only consider the leading order terms,
\begin{equation} w = w_0 +
{{O}}\big(\varepsilon^{1/2}\big), \quad \psi = \psi_0
+{{O}}\big(\varepsilon^{1/2}\big), \quad \theta = \theta_0 +
{{O}}\big(\varepsilon^{1/2}\big).
\end{equation}
We also expand the eigenvalue 
$\widetilde {Ra}$ in the same way
\begin{equation} \widetilde{Ra} = \widetilde{Ra}_0 + {{O}}\big(\varepsilon^{1/2}\big). \end{equation}3.1. Case (i): 
$k_\perp ^2\sim \varepsilon ^{1/2}$
In this case, vortex stretching is balanced by the horizontal diffusion of 
$\check {\psi }$ in (3.1a). In (3.1c), the horizontal diffusion of 
$\theta$ is dominant over the vertical, and to bring this into balance with 
$w$, we scale 
$\theta \to k_\perp ^{-2}\check {\theta }$. Finally, we rescale 
$\widetilde {Ra}\to \varepsilon ^{-1/2}\widehat {Ra}$ such that buoyancy is in balance with the vertical gradient of 
$\check {\psi }$ in (3.1b). This results in the leading order system
\begin{gather} -\partial_{Z} w_0= k_\perp^2\varepsilon^{{-}1/2} \check{\psi}_0, \end{gather}
\begin{gather}\partial_{Z}\check{\psi}_0 = \frac{\widehat{Ra}_0}{\sigma}\check{\theta}_0, \end{gather}
\begin{gather}- w_0 ={-}\frac{\check{\theta}_0}{\sigma}, \end{gather}with boundary conditions
\begin{equation} w_0 ={-}\frac{\delta}{\sqrt{2}} \check{\psi}_0, \quad Z = 0,1,\ \delta ={\pm}. \end{equation}This has solution
\begin{gather} w_0(Z) = A\left(
\frac{\varepsilon^{1/2}}{\sqrt{2}}\frac{\widetilde{Ra}_0^{1/2}}{k_\perp}
\cos\Big(Z \sqrt{k_\perp^2\widetilde{Ra}_0}\Big)+ \sin\Big(Z
\sqrt{k_\perp^2\widetilde{Ra}_0}\Big) \right) ,
\end{gather}
\begin{gather}\psi_0(Z) = A
\frac{\widetilde{Ra}_0^{1/2}}{k_\perp^3} \left( -\cos\Big(Z
\sqrt{k_\perp^2\widetilde{Ra}_0}\Big)+\frac{\varepsilon^{1/2}}{\sqrt{2}}\frac{\widetilde{Ra}_0^{1/2}}{k_\perp}
\sin\Big(Z \sqrt{k_\perp^2\widetilde{Ra}_0}\Big)\right) ,
\end{gather}
and 
$\theta _0 = \sigma k_\perp ^{-2} w_0$. The eigenvalue 
$\widetilde {Ra}_0$ must satisfy
\begin{equation}
\frac{\sqrt{2\varepsilon}\sqrt{k_\perp^2\widetilde{Ra}_0}}{k_\perp^2}\cos
\Big(\sqrt{k_\perp^2\widetilde{Ra}_0}\Big) +\left(
1-\frac{\varepsilon\widetilde{Ra}_0}{2k_\perp^2}\right)
\sin\Big(\sqrt{k_\perp^2\widetilde{Ra}_0}\Big)=0.
\end{equation}
The numerically computed root to this equation for a range of 
$k_\perp ^2$ is presented in figure 1 by the red dotted line. It agrees well with the underlying blue curve (the numerically computed marginal curve for the iNSE) near 
$k_\perp \approx 0.1 = \varepsilon ^{1/4}$, as expected. If we consider the limit of (3.8) as 
$k_\perp ^2$ becomes smaller than 
$\varepsilon ^{1/2}$, 
$\sin \Big(\sqrt {k_\perp ^2\widetilde {Ra}_0}\Big) \sim \sqrt {k_\perp ^2\widetilde {Ra}_0}$ and 
$\cos\Big(\sqrt{k_\perp ^2\widetilde {Ra}_0}\Big) \sim 1$, we find that 
$\widetilde {Ra}_0 \to 2\sqrt {2/\varepsilon }$, which is consistent with (3.11) below.
Choosing in particular 
$k_\perp ^2 = \varepsilon ^{1/2}$ and the corresponding 
$\widetilde {Ra}$ satisfying (3.8), examples of the vertical profiles (3.7) are plotted in figure 2 by the red dotted line. The underlying blue curves are the eigenfunctions of the numerically computed iNSE marginal stability problem. We see very good agreement between the two curves in the interior of the domain. The analytical asymptotic expressions for 
$w$ and 
$\psi$ derived from the CQG model do not match the numerical solutions from the iNSE in the boundary layer because the CQG model parametrizes rather than resolves the Ekman boundary layer. The temperature boundary condition is not satisfied due to the subdominance of the vertical diffusion of 
$\theta _0$, but there is a boundary layer of size 
$\varepsilon /k_\perp \sim \varepsilon ^{3/4}$ that allows 
$\theta _0 = 0$ to be satisfied.
Eigenfunction profiles for Case (i), 
$k_\perp ^2\sim \varepsilon ^{1/2}$. We have selected as an example 
$\varepsilon = 10^{-3}$ and 
$k_\perp = \varepsilon ^{1/4}$. Numerically computed eigenfunctions to the iNSE problem are in solid blue and the analytical asymptotic approximations (3.7) are in dashed red.
3.2. Case (ii): $\varepsilon ^{1/2}\gtrsim k_\perp ^2 \gtrsim \varepsilon ^2$
In this case, the horizontal diffusion of temperature is still dominant over the vertical diffusion, so we take 
$\theta \to k_\perp ^{-2} \check {\theta }$ and again, to retain the buoyancy term, we take 
$\widetilde {Ra} \to \varepsilon ^{-1/2}\widehat {Ra}$. This yields the system
\begin{gather} \partial_{Z}w_0 = 0, \end{gather}
\begin{gather}\partial_{Z}\check{\psi}_0 = \frac{\widehat{Ra}_0}{\sigma}\check{\theta}_0, \end{gather}
\begin{gather}w_0= \frac{\check{\theta}_0}{\sigma}, \end{gather}
with boundary conditions the same as above in (3.6). The horizontal diffusion of 
$\psi _0$ is now subdominant to vortex stretching, so 
$\partial _{Z}w_0 = 0$, which gives a constant for 
$w_0$. The solution is
\begin{gather} w_0 =A, \end{gather}
\begin{gather}\psi_0=A \frac{\widetilde{Ra}_0}{k_\perp^2}\left( Z-\frac{1}{2}\right) , \end{gather}
\begin{gather}\theta_0= \frac{\sigma A}{k_\perp^2}, \end{gather}with eigenvalue
\begin{equation} \widehat{Ra}_0 =2\sqrt{2}\quad \text{i.e.}\ \widetilde{Ra}_0 = \frac{2\sqrt{2}}{\varepsilon^{1/2}}. \end{equation}
The constant estimate for 
$\widetilde {Ra}$ given in (3.11) is shown in figure 1 in purple. We can see that it aligns well with the flat portion of the numerical iNSE curve. We also compare the eigenfunctions (3.10) with the numerical iNSE in figure 3. Again, we see good agreement in the interior of the domain, but as the asymptotic result is designed to remove the boundary layer, they do not match in the Ekman boundary region. As in Case (i), there is a thermal boundary layer of size 
$\varepsilon /k_\perp$ which allows the temperature boundary condition to be satisfied. In this case, the boundary layer is still asymptotically thin since 
$1 \gtrsim \varepsilon /k_\perp \gtrsim \varepsilon ^{3/4}$.
Eigenfunction profiles for Case (ii), 
$\varepsilon ^{1/2}\gtrsim k_\perp ^2 \gtrsim \varepsilon ^2$, with 
$\varepsilon = 10^{-4}$ and 
$k_\perp \approx 0.0203$. Numerically computed eigenfunctions to the iNSE problem are in solid blue, and the analytical asymptotic approximations (3.10) are in dashed red.
3.3. Case (iii): $k_\perp ^2 \sim \varepsilon ^2$
In this case, the vertical diffusion term of the temperature equation is now in balance with the horizontal, and we again use the scalings 
$\theta \to k_\perp ^{-2}\check {\theta }$ and 
$\widetilde {Ra} \to \varepsilon ^{-1/2}\widehat {Ra}$. This gives us the leading order problem
\begin{gather} \partial_{Z}w_0= 0, \end{gather}
\begin{gather}\partial_{Z}\check{\psi}_0= \frac{\widehat{Ra}_0}{\sigma}\check{\theta}_0, \end{gather}
\begin{gather}- w_0 = \frac{1}{\sigma} \big(\varepsilon^2k_\perp^{{-}2}\partial_Z^2 - 1\big) \check{\theta}_0, \end{gather}with boundary conditions
\begin{equation} w_0 ={-}\frac{\delta}{\sqrt{2}} \check{\psi}_0, \quad \check{\theta}_0=0,\quad Z = 0,1, \delta ={\pm} 1. \end{equation}This has solution
\begin{gather} w_0 = A, \end{gather}
\begin{gather}\psi_0 =\frac{A \widetilde{Ra}_0}{k_\perp^2} \left( Z-\frac{1}{2} -\frac{\varepsilon}{k_\perp}\frac{{\rm e}^{k_\perp Z/\varepsilon}}{1+{\rm e}^{k_\perp{/}\varepsilon}}+\frac{\varepsilon}{k_\perp}\frac{{\rm e}^{{-}k_\perp Z/\varepsilon}}{1+{\rm e}^{{-}k_\perp{/}\varepsilon}} \right) , \end{gather}
\begin{gather}\theta_0 = \frac{\sigma A}{k_\perp^2}\left( 1-\frac{{\rm e}^{k_\perp Z/\varepsilon}}{1+{\rm e}^{k_\perp{/}\varepsilon}}-\frac{{\rm e}^{{-}k_\perp Z/\varepsilon}}{1+{\rm e}^{{-}k_\perp{/}\varepsilon}} \right) , \end{gather}with
\begin{equation} \widetilde{Ra}_0 = \sqrt{\frac{2}{\varepsilon}} \left( \frac{1}{2} +\frac{\varepsilon}{k_\perp}\frac{1}{1+{\rm e}^{k_\perp{/}\varepsilon}}-\frac{\varepsilon}{k_\perp}\frac{1}{1+{\rm e}^{{-}k_\perp{/}\varepsilon}} \right) ^{{-}1}. \end{equation}
Here, the second-order derivative allows for the temperature boundary condition to be satisfied, and the thermal boundary layer width is 
${{O}}(1)$.
In figure 1, the 
$\widetilde {Ra}_0$ curve is plotted as a yellow dotted line. It agrees well with the underlying blue curve around the corner near 
$k_\perp \approx 10^{-4} = \varepsilon$, as well as the region surrounding this value. If we consider (3.15), where 
$k_\perp \gg \varepsilon$, it follows that 
$\widetilde {Ra}_0 \to 2\sqrt {2/\varepsilon }$, which aligns with (3.11). We may also consider 
$k_\perp \ll \varepsilon$, and using a few terms in the series 
${1}/{(1+{\rm e}^{x})}\approx {1}/{2}+{x}/{2}-{x^3}/{48}$, we find that 
$\widetilde {Ra}_0 \sim 24\sqrt {2}\varepsilon ^{3/2}k_\perp ^{-2}$. See below that this agrees with (3.18).
We also plot the eigenfunctions in figure 4, where we again see good agreement with the numerical iNSE linear stability solution.
Eigenfunction profiles for Case (iii), 
$k_\perp ^2 \sim \varepsilon ^2$, with 
$\varepsilon = 10^{-3}$, 
$k_\perp = 9\times 10^{-3}$. Numerically computed eigenfunctions to the iNSE problem are in solid blue, and the analytical asymptotic approximations (3.14) are in dashed red.
3.4. Case (iv): $k_\perp ^2 \lesssim \varepsilon ^2$
In this case, we lose the horizontal diffusion term and the vertical diffusion of temperature is dominant. This implies the scalings 
$\theta \to \varepsilon ^{-2}\check {\theta }$ and 
$\widetilde {Ra} \to \varepsilon ^{3/2}k_\perp ^{-2}\widehat {Ra}$, and yields the leading order system
\begin{gather} \partial_{Z}w_0= 0, \end{gather}
\begin{gather}\partial_{Z}\check{\psi}_0= \frac{\widehat{Ra}_0}{\sigma}\check{\theta}_0, \end{gather}
\begin{gather}- w_0 = \frac{1}{\sigma} \partial_Z^2 \check{\theta}_0, \end{gather}with solution
\begin{gather} w_0 = A, \end{gather}
\begin{gather}\psi_0 =A \frac{\widetilde{Ra}_0}{\varepsilon^2} \left( -\frac{Z^3}{6}+\frac{Z^2}{4}-\frac{1}{24} \right) , \end{gather}
\begin{gather}\theta_0 = \frac{\sigma A}{2\varepsilon^2 }Z(1-Z), \end{gather}and eigenvalue
\begin{equation} \widehat{Ra}_0 = 24\sqrt{2} \quad \text{i.e.} \ \widetilde{Ra}_0 = \frac{24\sqrt{2}\varepsilon^{3/2}}{k_\perp^2}. \end{equation}
This curve is plotted in green in figure 1, and is in good agreement with iNSE in the small 
$k_\perp$ limit. The eigenfunctions shown in figure 5 are also in good agreement with iNSE.
Eigenfunction profiles for Case (iv), 
$k_\perp ^2 \lesssim \varepsilon ^2$, with 
$\varepsilon = 10^{-3}$, 
$k_\perp = 10^{-4}$. Numerically computed eigenfunctions to the iNSE problem are in solid blue, and the analytical asymptotic approximations (3.17) are in dashed red.
4. Discussion and conclusions
Through the four cases considered above, we have demonstrated that the linear stability of the rapidly rotating iNSE system with no-slip boundary conditions can be well represented by leading order asymptotic approximations. Linear stability theory for stress-free boundary conditions has been established for decades. Our previous results of numerical computation show a departure from that theory for no-slip boundary conditions (Tro et al. Reference Tro, Grooms and Julien2024). We have considered analytical asymptotic solutions to the linear stability problem for small wavenumbers 
$k_\perp ^2 \lesssim 1$. The transition of 
$k_\perp ^2$ to values smaller than 
$\varepsilon ^{1/2}$ leads to a loss of the lateral viscous term that normally balances vortex stretching, and this combined with the pumping boundary conditions leads to a constant, non-zero value for 
$w$ in the interior of the domain. The distinction between subsequent cases is whether the vertical and horizontal diffusion of temperature are in balance, or one is dominant over the other. At the smallest scales, for 
$k_\perp ^2 < \varepsilon ^2$, vertical diffusion dominates.
The condition that the pumping boundary conditions become 
${{O}}(1)$ is closely tied to the condition that 
$\partial _{Z}w=0$. At large values of 
$k_\perp$, 
$k_\perp ^4 \psi$, does not become subdominant to 
$\partial _{Z}w$. Scaling the eigenfunctions such that 
$w\sim 1$ in the interior of the domain, by the vortex-stretching balance, 
$1\sim w \sim k_\perp ^4 \psi$. This implies 
$k_\perp ^2\psi \sim k_\perp ^{-2}$, making it impossible for pumping 
$\sim \varepsilon ^{1/2}k_\perp ^2 \psi \sim \varepsilon ^{1/2}/k_\perp ^2$ to become 
${{O}}(1)$ for large 
$k_\perp$.
The CQG model is derived under the assumption that the aspect ratio of horizontal to vertical scales, 
$\ell /H$, is 
${{O}}(\varepsilon )$ in the limit of 
$\varepsilon \to 0$. In Cases (iii) and (iv), when we consider 
$k_\perp \leq \varepsilon$, we are violating that assumption and the validity of the CQG model. To justify this, we show in the Appendix that if we begin from the iNSE system, we can derive the same systems of (3.12) and (3.16). In essence, any terms that could arise when vertical scales are not assumed smaller than horizontal scales remain subdominant due to the smallness of 
$k_\perp$.
From our results, we can estimate the value of 
$\widetilde {Ra}$ above which the departure from stress-free occurs. The stress-free marginal curve is shown in figure 1 as a blue dash-dotted line. We can see for 
$\widetilde {Ra}$ below the purple line, 
$2\sqrt {2/\varepsilon }$, the no-slip and stress free curves match up. However, if 
$\widetilde {Ra}$ is above 
$2\sqrt {2/\varepsilon }$, we see a much wider range of wavenumbers contributing to the instability. From (3.18), we may estimate that wavenumbers
\begin{equation} \frac{24\sqrt{2}\varepsilon^{3/2}}{\widetilde{Ra}} \leq k_\perp^2 \leq \widetilde{Ra}^{1/2}\quad \text{for}\ \widetilde{Ra} > \frac{2\sqrt{2}}{\varepsilon^{1/2}}, \end{equation}produce modes with a positive growth rate. This is significantly different from the range for the stress-free case,
\begin{equation} \frac{{\rm \pi}^2}{\widetilde{Ra}} \leq k_\perp^2 \leq \widetilde{Ra}^{1/2}, \end{equation}
which does not change with 
$\varepsilon$. Whether these modes affect direct numerical simulations of the fully nonlinear problem remains unexplored.
Funding
This work was supported by the National Science Foundation (grant no. DMS-2308337).
Declaration of interests
The authors report no conflict of interest.
Appendix. Asymptotic justification for Cases (iii) and (iv)
Beginning with the rescaled iNSE, (2.2), we consider the scalings of Cases (iii) and (iv). Because 
$\psi$ is the streamfunction such that 
$\partial _{x}\psi = v$, 
$\partial _{y}\psi = -u$ and 
$\psi = p$, the rescaling 
$\psi \to k_\perp ^{-2}\varepsilon ^{-1/2}\check {\psi }$ implies that
\begin{equation} u\to k_\perp^{{-}2}\varepsilon^{{-}1/2}\check{u}, \quad v\to k_\perp^{{-}2}\varepsilon^{{-}1/2}\check{v}, \quad p \to k_\perp^{{-}2}\varepsilon^{{-}1/2}\check{p}. \end{equation}
For Case (iii), we take 
$\theta \to k_\perp ^{-2}\check {\theta }$ and 
$\widetilde {Ra} \to \varepsilon ^{-1/2}\widehat {Ra}$ and for Case (iv), 
$\theta \to \varepsilon ^{-2}\check {\theta }$ and 
$\widetilde {Ra} \to \varepsilon ^{3/2}k_\perp ^{-2} \widehat {Ra}$. Applying this scaling to the linear terms of (2.2a) with the normal mode ansatz (2.6) substituted in, and 
$s=0$, we get
\begin{gather} -\check{v}+k_x\check{p} = \varepsilon ({-}k_\perp^2 +\varepsilon^2 \partial_{Z}^2 ) \check{u}, \end{gather}
\begin{gather}\check{u}+k_y \check{p}= \varepsilon ({-}k_\perp^2 +\varepsilon^2 \partial_{Z}^2) \check{v}, \end{gather}
\begin{gather}0 = k_\perp \varepsilon^{3/2} \left( k_\perp^2 \varepsilon^{1/2}\big({-}k_\perp^2 +\varepsilon^2 \partial_{Z}^2\big) w + \frac{\widehat{Ra} }{ \sigma} \check{\theta} - \partial_{Z}\check{p}\right) , \end{gather}with (2.2b) becoming
\begin{gather} 0 = k_\perp \varepsilon^{3/2}\left( w+\frac{1}{\sigma}\big(\varepsilon^2k_\perp^{{-}2}\partial_{Z}^2 -1\big) \check{\theta}\right) \quad \text{for Case (iii)}, \end{gather}
\begin{gather}0 = k_\perp \varepsilon^{3/2} \left( w + \frac{1}{\sigma}\partial_{Z}^2 \check{\theta}\right) \quad \text{for}\ {\rm Case}\ \text{(iv)}, \end{gather}and (2.2c)
\begin{equation} \boldsymbol{k}_\perp\boldsymbol{\cdot} \check{\boldsymbol{u}}_\perp + k_\perp^2\varepsilon^{3/2} \partial_{Z}w =0, \end{equation}
where 
$\check {\boldsymbol {u}}_\perp = ( \check {u},\check {v})$ and recall 
$\boldsymbol {k}_\perp = ( k_x, k_y)$. We then expand the rescaled variables as follows:
\begin{gather} \check{u} = \check{u}_0 +k_\perp \varepsilon^{3/2}\check{u}_1 + \cdots, \end{gather}
\begin{gather}\check{v} = \check{v}_0 +k_\perp \varepsilon^{3/2} \check{v}_1 +\cdots, \end{gather}
\begin{gather}\check{p} = \check{p}_0 +k_\perp \varepsilon^{3/2} \check{p}_1 +\cdots, \end{gather}
\begin{gather}w = w_0 +\cdots, \end{gather}
\begin{gather}\check{\theta} = \check{\theta}_0 +\cdots. \end{gather}Substituting in these expansions, the leading order system is
\begin{gather} -\check{v}_0 +k_x\check{p}_0 = 0, \end{gather}
\begin{gather}\check{u}_0 +k_y\check{p}_0 = 0, \end{gather}
\begin{gather}k_x \check{u}_0 +k_y \check{v}_0 = 0, \end{gather}
which returns to us the streamfunction 
$\check {\psi }_0 = \check {p}_0$, with 
$k_x\check {\psi }_0 = \check {v}_0$ and 
$k_y\check {\psi }_0 = -\check {u}_0$ as expected.
To determine the next order problem, we must consider the scaling of 
$-k_\perp ^2 +\varepsilon ^2 \partial _{Z}^2$ for Cases (iii) and (iv), where 
$k_\perp ^2 \leq \varepsilon ^2$:
\begin{equation} {-}k_\perp^2 +\varepsilon^2 \partial_{Z}^2 \sim \max\big\{ k_\perp^2,\ \varepsilon^2 \big\} \sim \varepsilon^2. \end{equation}
Note also that 
$\boldsymbol {k}_\perp \sim k_\perp$. Then, the 
${{O}}(k_\perp \varepsilon ^{3/2})$ terms are
\begin{gather} -v_1 + k_x p_1 = 0 , \end{gather}
\begin{gather}u_1 +k_y p_1 =0 , \end{gather}
\begin{gather}k_x u_1 + k_y v_1 ={-}\partial_{Z}w_0, \end{gather}which has the solvability condition that
\begin{equation} \partial_{Z}w_0 = 0, \end{equation}and from (A3),
\begin{gather} 0 = w_0+\frac{1}{\sigma}\big(\varepsilon^2k_\perp^{{-}2}\partial_{Z}^2 -1\big) \check{\theta}_0 \quad \text{for}\ {\rm Case}\ \text{(iii)}, \end{gather}
\begin{gather}0 = w_0 + \frac{1}{\sigma}\partial_{Z}^2 \check{\theta}_0 \quad \text{for}\ {\rm Case}\ \text{(iv)}, \end{gather}and finally from (A4),
\begin{equation} \partial_{Z} \check{\psi}_0 = \frac{\widehat{Ra}_0}{\sigma}\check{\theta}_0. \end{equation}
Open access
