1. Introduction
Since the seminal work of Chandrasekhar (Reference Chandrasekhar1961), who systematically analysed the foundational theory for convective stability, the problem of thermal convective instability in rotating spheres or spherical shells has continued to attract considerable interest, owing to its broad relevance in geophysical and planetary physics as well as its mathematical richness (Roberts Reference Roberts1968; Busse Reference Busse1970; Zhang Reference Zhang1994; Jones, Soward & Mussa Reference Jones, Soward and Mussa2000; Dormy et al. Reference Dormy, Soward, Jones, Jault and Cardin2004; Kaplan et al. Reference Kaplan, Schaeffer, Vidal and Cardin2017). Subsequent developments have extended this framework, delivering robust numerical solutions and progressively complete asymptotic descriptions (Zhang & Liao Reference Zhang and Liao2004; Zhang, Lam & Kong Reference Zhang, Lam and Kong2017). During the early stages of planetary evolution, the interiors of undifferentiated protoplanets and small bodies are thought to consist of unconsolidated and porous matrices, so a porous medium convection approach (Nield & Bejan Reference Nield and Bejan2017) offers a more realistic representation than pure-fluid models. However, studies of porous convection in spherical geometries (Young, Zhang & Schubert Reference Young, Zhang and Schubert2003; Zhang, Liao & Schubert Reference Zhang, Liao and Schubert2005) are scarce, and the majority of existing work has focused on plane models in Cartesian coordinates, leaving porous convection in rapidly rotating spheres largely unexplored so far.
The stability of rotating systems is conventionally characterised by the Taylor number
$\textit{Ta}$
(or the Ekman number
$Ek$
), which quantifies the balance between rotational and viscous effects. It is well established that rotation suppresses convective motion in pure fluids, thereby elevating the critical Rayleigh number
$Ra_c$
required for the onset of convection (Chandrasekhar Reference Chandrasekhar1961). Research shows that rotation also stabilises convection in horizontal porous layers (Malashetty, Swamy & Sidram Reference Malashetty, Swamy and Sidram2010), elevating
$Ra_c$
and increasing the preferred wavenumber, thereby encouraging the development of a greater number of convective cells (Jou & Liaw Reference Jou and Liaw1987; Malashetty et al. Reference Malashetty, Swamy and Sidram2010). In unbounded horizontal porous domains, linear stability theory predicts that the critical Rayleigh number and wavenumber scale with
$\textit{Ta}$
and
$\textit{Ta}^{0.5}$
, respectively (Riahi Reference Riahi1994; Vadasz Reference Vadasz1998). However, the onset criteria become sensitive to the aspect ratio in confined configurations such as rectangular boxes (Vadasz Reference Vadasz1995; Govender Reference Govender2003). Furthermore, both Coriolis and centrifugal forces can jointly govern the stability in rotating porous systems (Vadasz Reference Vadasz1996; Vadasz & Govender Reference Vadasz and Govender2001). Despite these advances in plane models, understanding how rotation influences the onset of convection in spherical geometries remains an open question of both fundamental and practical importance.
This study aims to fill this theoretical gap by analysing the linear stability of thermal convection in a rotating porous sphere. Using Darcy’s law as the governing framework (Nield & Bejan Reference Nield and Bejan2017), we systematically examine the combined influence of spherical geometry and Coriolis forces on convective onset, and contrast the results with established predictions for rotating fluid convection (Sánchez, Garcia & Net Reference Sánchez, Garcia and Net2016a ). Our analysis provides a rigorous theoretical basis for future studies of finite-amplitude convection and its application to early planetary thermal evolution.
2. Mathematical formulation
We consider a sphere of radius
$R$
filled uniformly with an internally heated, fluid-saturated porous medium that rotates at a constant angular velocity
$\varOmega$
. Gravity is taken to vary linearly with radius, while centrifugal effects are neglected in comparison with gravitational effects (Sánchez et al. Reference Sánchez, Garcia and Net2016a
,
Reference Sánchez, Garcia and Netb
). The fluid density is denoted by
$\rho$
, the internal heating rate by
$H$
, the kinematic viscosity by
$\nu$
, the porosity by
$\varphi$
, the solid-matrix permeability by
$K$
, the thermal expansion coefficient by
$\alpha$
, the surface gravitational acceleration by
$g_0$
, the overall thermal conductivity by
$k$
and the thermal diffusion coefficient by
$\kappa$
. The governing equations are non-dimensionalised using
$R$
,
$R^2/\kappa$
and
$\Delta T_c = \rho \textit{HR}^2/(6k)$
as the characteristic scales of length, time and temperature, respectively (Yu & Wang Reference Yu and Wang2025). Under the Boussinesq approximation (Turcotte & Schubert Reference Turcotte and Schubert2014) and Darcy’s law (Nield & Bejan Reference Nield and Bejan2017), and with the Coriolis force included, the dimensionless momentum equation takes the form
\begin{equation} \frac {\textit{Da}}{\varphi \textit{Pr}}\frac {\partial \boldsymbol{u}}{\partial t} +\boldsymbol{u} = -\boldsymbol{\nabla }\!p + Ra \varTheta \boldsymbol{r} - \tau \boldsymbol{e}_z \times \boldsymbol{u}, \end{equation}
where
$\boldsymbol{u}$
is the seepage velocity,
$p$
denotes the pressure and
$\varTheta$
represents the temperature deviation from the conductive reference profile
$\varTheta _c = 1 - r^2$
. The key dimensionless parameters are defined as
\begin{align} \textit{Pr}=\frac {\nu }{\kappa },\quad {\textit{Da}}=\frac {K}{R^2},\quad {Ra} = \frac {g_0 \alpha \Delta T_c K R}{\nu \kappa }=\frac {\rho H g_0 \alpha K R^3}{6\nu \kappa k}, \quad \tau =\sqrt {\textit{Ta}} = \frac {2 \varOmega K}{\varphi \nu } .\end{align}
Since the Darcy number
$\textit{Da}$
at planetary scales is typically of order
$10^{-10}$
or smaller, the inertial term in the momentum equation is neglected in this study for simplification:
\begin{equation} \boldsymbol{u} = -\boldsymbol{\nabla }\!p + Ra \varTheta \boldsymbol{r} - \tau \boldsymbol{e}_z \times \boldsymbol{u}. \end{equation}
To eliminate the pressure gradient, we take the curl of (2.3), which leads to an equation for the evolution of vorticity:
\begin{equation} \tau \frac {\partial \boldsymbol{u}}{\partial z}= \boldsymbol{\nabla }\times (\boldsymbol{u} - Ra\varTheta \boldsymbol{r}). \end{equation}
The mass conservation and energy equations are, respectively,
\begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0, \end{equation}
\begin{equation} \frac {\partial \varTheta }{\partial t} = {\nabla} ^2 \varTheta + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }(\varTheta + \varTheta _c), \end{equation}
with isothermal and impermeable boundary conditions
$u_r |_{r=1} = \varTheta |_{r=1} = 0$
.
A vector potential decomposition
$\boldsymbol{u}=\boldsymbol{\nabla }\times (\boldsymbol{\nabla }\times v \boldsymbol{r}) + \boldsymbol{\nabla }\times w \boldsymbol{r}$
(Schubert, Turcotte & Olson Reference Schubert, Turcotte and Olson2001) is introduced to satisfy (2.5) automatically. Substituting the expression into (2.6) and neglecting nonlinear terms gives the linearised energy equation
\begin{equation} \frac {\partial \varTheta }{\partial t} = {\nabla} ^2 \varTheta + 2 L^2 v. \end{equation}
Then we apply the radial
$\boldsymbol{r} \boldsymbol{\cdot}$
and the radial curl
$\boldsymbol{r} \boldsymbol{\cdot }\boldsymbol{\nabla }\times$
operators to (2.4), obtaining a form similar to that of Sánchez et al. (Reference Sánchez, Garcia and Net2016a
,Reference Sánchez, Garcia and Net
b
):
\begin{equation} L^2 w + \tau \left (\frac {\partial w}{\partial \phi } - Qv\right ) = 0, \end{equation}
\begin{equation} L^2{\nabla} ^2 v + \tau \left (\frac {\partial }{\partial \phi }{\nabla} ^2 v + Qw\right ) + RaL^2 \varTheta = 0, \end{equation}
where the operators are defined as
$L^2 = -( {1}/{\sin \theta }) ({\partial }/{\partial \theta }) (\sin \theta ( {\partial }/{\partial \theta }) ) -{} ( {1}/{\sin ^2 \theta } )({\partial ^2}/{\partial \phi ^2})$
and
$Q = r \cos \theta {\nabla} ^2 - (L^2 + r ({\partial }/{\partial r} )) (\cos \theta ( {\partial }/{\partial r}) - ( {\sin \theta }/{r}){} ({\partial }/{\partial \theta }) )$
, respectively. Consequently, the linear stability analysis decouples in the azimuthal wavenumber
$m$
but remains coupled in the spherical harmonic degree
$l$
. At the onset of convection in a rotating reference frame, the disturbance is a neutral oscillation with a purely imaginary growth rate
${\rm e}^{{\rm i} \omega t}$
, where
$\omega$
is a real number representing the drift frequency and its sign indicates the direction of propagation. We expand the variable
$\varTheta _m$
in normalised complex spherical harmonics, allowing the problem to be separated azimuthally:
\begin{equation} \varTheta _m = \sum _{l=m}^{\infty } a_{m,l} (r) Y_l^m (\theta , \phi ) {\rm e}^{{\rm i} \omega t} .\end{equation}
Similarly,
$v_m$
and
$w_m$
can be expanded in the same manner by
$b_{m,l} (r)$
and
$c_{m,l} (r)$
, respectively. Substituting them into (2.7)–(2.9) yields
\begin{equation} {\rm i} \omega a_{m,l} = D_l a_{m,l} + 2l(l+1)b_{m,l}, \end{equation}
\begin{equation} -\left [1 + \frac {{\rm i} m \tau }{l(l+1)} \right ] c_{m,l}=\tau (X_l^m D_{l-1}^- b_{m,l-1}+ Z_l^m D_{l+1}^+ b_{m,l+1}), \end{equation}
\begin{equation} \left [1 + \frac {{\rm i} m \tau }{l(l+1)}\right ] D_l b_{m,l} + Ra a_{m,l}=\tau (X_l^m D_{l-1}^- c_{m,l-1} + Z_l^m D_{l+1}^+ c_{m,l+1}), \end{equation}
where
$X_l^m$
and
$Z_l^m$
depend on
$l$
and
$m$
, and the operators are
$D_l=( {\partial ^2}/{\partial r^2} )+ ({2}/{r}){}( {\partial }/{\partial r}) - {l(l+1)}/{r^2}$
,
$D_l^-=({\partial }/{\partial r})-( {l}/{r})$
and
$D_l^+=( {\partial }/{\partial r})+( {l+1})/{r}$
. Substituting (2.11)–(2.12) into (2.13) leads to the following recurrence relation:
\begin{equation} A_l^m(\tau ,\omega )a_{m,l-2} + B_l^m(\tau ,\omega )a_{m,l} + C_l^m(\tau ,\omega )a_{m,l+2} = Ra a_{m,l}, \end{equation}
where
$A_l^m$
,
$B_l^m$
and
$C_l^m$
all depend on
$l,\ m,\ \tau$
and
$\omega$
. The radial component of the temperature field
$a_{m,l}(r)$
must satisfy regularity conditions at the centre, requiring it to be the product of
$r^l$
and an even function. It is therefore expanded in Worland polynomials (Livermore, Jones & Worland Reference Livermore, Jones and Worland2007; Sternberg, Marti & Jackson Reference Sternberg, Marti and Jackson2025):
\begin{equation} a_{m,l} = \begin{bmatrix} W_0^l & \cdots & W_n^l \end{bmatrix} \boldsymbol{F}_{m,l}, \quad \boldsymbol{F}_{m,l} = \begin{bmatrix} f_{0,m,l} \\ \vdots \\ f_{n,m,l} \end{bmatrix}\!. \end{equation}
For a given drift frequency
$\omega$
, finding the critical
$Ra$
reduces to the matrix eigenvalue problem
\begin{equation} \boldsymbol{M}_m \boldsymbol{P}_m = Ra \boldsymbol{P}_m, \end{equation}
where the coefficient vector
$\boldsymbol{P}_m$
and the operator matrix
$\boldsymbol{M}_m$
are defined as
\begin{equation} \boldsymbol{P}_m = \begin{bmatrix} \boldsymbol{F}_{m,l} \\ \vdots \\ \boldsymbol{F}_{m,L} \end{bmatrix}\!,\quad \boldsymbol{M}_m = \begin{bmatrix} B_l^m & C_l^m & & & \\[5pt] A_{l+2}^m & B_{l+2}^m & C_{l+2}^m & & \\[5pt] & \ddots & \ddots & \ddots & \\[5pt] & & A_{L-2}^m & B_{L-2}^m & C_{L-2}^m \\[5pt] & & & A_L^m & B_L^m \end{bmatrix}\!. \end{equation}
The resulting eigenvalue problem is solved numerically by the spectral method (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988; Yu & Wang Reference Yu and Wang2025). Boundary conditions are enforced by modifying two rows of the corresponding submatrices in
$\boldsymbol{M}_m$
. Spectral truncation is applied with a maximum spherical harmonic degree
$l_{\textit{max}} = 200$
, a maximum azimuthal order
$m_{\textit{max}} = 200$
and radial truncation number
$n$
satisfying
$ l+2n\leq 202$
. These truncation values were chosen to ensure complete resolution of all dynamically relevant scales, with the adequacy of resolution confirmed through convergence tests.
The eigenvalues of
$\boldsymbol{M}_m$
in (2.17) are generally complex. A shooting method is employed to select the appropriate
$\omega$
such that the relevant eigenvalue becomes real. The onset of physically feasible convection corresponds to the smallest one, which defines the critical Rayleigh number
$Ra_c$
, with the associated azimuthal wavenumber
$m_c$
and drift frequency
$\omega _c$
. The numerical code was originally developed and validated against benchmark results for rotating fluid convection (Zhang & Liao Reference Zhang and Liao2004; Sánchez et al. Reference Sánchez, Garcia and Net2016a
). Following a similar approach, we employ the Arnoldi inverse shift iteration algorithm from ARPACK (Lehoucq & Sorensen Reference Lehoucq and Sorensen1996) to compute a subset of eigenvalues efficiently. The code has been specifically adapted for the rotating porous convection problem studied here. All programs are implemented in Fortran 90 and are available from the authors upon request.
3. Results
In the absence of rotation, the degeneracy of the convective instability yields a critical Rayleigh number
$Ra_c$
that depends solely on the spherical harmonic degree
$l$
, independent of the azimuthal order
$m$
(Schubert et al. Reference Schubert, Turcotte and Olson2001; Young et al. Reference Young, Zhang and Schubert2003). As shown in table 1, the preferred mode in the sphere is
$l = 2$
with a critical value of 91.950. The result agrees well with that of Young et al. (Reference Young, Zhang and Schubert2003) while noting that the Rayleigh number definition used herein is half of theirs. The governing equations and boundary conditions are spherically symmetric, allowing free orientation of the convective structure’s symmetry axis. For comparison with subsequent discussions, figure 1(a–f) depicts the convective structure for the
$Y^2_2$
mode. This structure is equatorially symmetric, forming two identical upwelling (red in figure 1
a–c) and downwelling (blue in figure 1
a–c) regions. The distribution of azimuthal velocity component
$u_\phi$
and the kinetic energy
$E_k$
are consistent (figure 1
d–f). The configuration corresponds to the limit of slow rotation, and thus the degeneracy issue arising from the uncertainty in
$m$
need not be discussed (Schubert et al. Reference Schubert, Turcotte and Olson2001; Young et al. Reference Young, Zhang and Schubert2003).
Critical Rayleigh number
$Ra_{c}$
versus spherical harmonic degree
$l$
in the non-rotating case.
Convective structure of four solutions. (a–f) The
$Y_2^2$
mode:
$\tau = 0$
,
$Ra_c = 91.95$
. (g–l) Equatorially symmetric mode:
$\tau = 8$
,
$m = 3$
,
$Ra_c = 3.774 \times 10^2$
,
$\omega _c = 7.864$
. (m–r) Equatorially antisymmetric mode:
$\tau = 8$
,
$m = 1$
,
$Ra_c = 8.961$
$\times$
$10^2$
,
$\omega _c = 15.443$
. (s–x) Equatorially symmetric mode:
$\tau = 2 \times 10^4$
,
$m = 2$
,
$Ra_c = 3.704 \times 10^8$
,
$\omega _c = 2.977 \times 10^4$
. (a,d,g, j,m,p,s,v) Equatorial cross-sections. (b,e,h,k,n,q,t,w) Certain meridional cross-sections. (c,f,i,l,o,r,u,x) Surface boundaries. Arrows in (a–c,g–i,m–o) show the velocity field projection onto the respective planes. Contours in (d–f, j–l,p–r) represent kinetic energy
$E_k$
. Background colours in (a,b,g,h,m,n) indicate temperature perturbation, in (c,i,o,u) heat flux perturbation and in (d–f, j–l,p–r,v–x) the azimuthal velocity component
$u_\phi$
. In the strongly rotating regime, both the velocity field and kinetic energy contours are highly concentrated near the polar axis. To prevent visual clutter, velocity arrows and kinetic energy contours have been omitted from (s–x), retaining only the background colour fields.
Under the influence of the Coriolis force, rotation significantly alters the stability and flow patterns of porous convection. The unstable modes decouple in azimuthal wavenumber
$m$
but couple in degree
$l$
, and bifurcate into equatorially symmetric (
$l-m$
is even) and antisymmetric (
$l-m$
is odd) families. As rotation strengthens (with increasing
$\tau$
),
$Ra_c$
required for the onset of convection increases accordingly (figure 2
a). We only discovered non-oscillatory solutions (
$\omega = 0$
) in the axisymmetric mode (
$m=0$
). For modes with azimuthal wavenumber
$m \gt 0$
, similar to pure-fluid convection, the system exhibits a drifting instability where the flow pattern maintains its form and propagates (figure 2
b). As rotation sets in, the drift direction initially opposes the rotation direction (
$\omega \lt 0$
). When the rotation rate increases further,
$\omega$
changes from negative to positive and eventually grows continuously with
$\tau$
.
Critical parameters for convection at low rotation rates: (a)
$Ra_{c}$
versus
$\tau$
; (b)
$\omega _c$
versus
$\tau$
. Different colours represent different
$m$
; solid lines correspond to equatorially symmetric modes, dashed lines to antisymmetric modes.
Critical parameters for the preferred mode under strong rotation: (a)
$Ra_{c}$
versus
$\tau$
; (b)
$\omega _{c}$
versus
$\tau$
.
We focus on the minimum
$Ra$
, which corresponds to the physically realisable mode. The preferred equatorially symmetric convective structure, shown in figure 1(g–l) for
$\tau = 8$
and
$m = 3$
, features an equatorial region partitioned into six cells exhibiting
$120^{\circ }$
rotational symmetry. This structure closely resembles that observed in rotating convection of pure fluids (Marti et al. Reference Marti2014). And for the same
$\tau$
, the preferred antisymmetric mode (figure 1
m–r) corresponds to
$m = 1$
and higher
$Ra_c$
. Figure 2(a) demonstrates that the equatorially symmetric mode is preferred across the entire parameter range explored in this study, with no discernible trend suggesting a potential preference for antisymmetric modes even at higher values of
$\tau$
. The
$m = 2$
mode is optimal over the entire parameter range explored with
$m = 3$
being marginally preferred only within a very narrow interval around
$1.5 \lesssim \tau \lesssim 8.1$
.
As indicated by figure 3 under strong rotation, the critical parameters scale as
$Ra_c \propto \tau ^2$
and
$\omega _c \propto \tau ^1$
. In the subsequent section, we provide a concise analysis of the underlying reasons for the validity of these two simple power laws. Figure 1(g–x) displays the convective structure in the strongly rotating regime: perturbations of temperature and kinetic energy are strongly confined to axially elongated columnar regions aligned with the rotation axis. This configuration concentrates nearly all heat and mass transport close to the polar axis – a manifestation of intense polar focusing – in striking contrast to the small-scale structures that accumulate near the equator in traditional rotating pure-fluid convection.
4. Discussion
To better compare thermal convection in porous media with that in pure fluids within a rotating spherical geometry, we first recall the linear stability results for the classic basal heated-plane-layer model. The Rayleigh number
$Ra$
, the Taylor number
$\textit{Ta}$
and the horizontal wavenumber
$a$
are related by (Chandrasekhar Reference Chandrasekhar1961; Riahi Reference Riahi1994)
\begin{align} \text{Pure fluid:} \quad Ra = \frac {(\pi ^2 + a^2)^3 + \pi ^2 Ta}{a^2}. \end{align}
\begin{align} \text{Porous-layer medium:} \quad Ra = \frac {(\pi ^2 + a^2)^2 + \pi ^2 Ta (\pi ^2 + a^2)}{a^2}. \end{align}
Although the spherical curvature and the internally heated configuration employed in this work introduce geometric and physical differences, preventing a strict correspondence between the plane results and the spherical problem, this analogy can still reveal the key behavioural distinctions between the two systems under the influence of rotation.
In the non-rotating case (
$\textit{Ta} = 0$
), the wavenumbers
$a = \pi / \lambda$
that minimise
$Ra$
correspond to aspect ratios
$\lambda = \sqrt {2}$
(pure fluid) (Lord Rayleigh 1916) and
$\lambda = 1$
(porous medium) (Horton & Rogers Reference Horton and Rogers1945; Lapwood Reference Lapwood1948). This indicates that porous-medium convection tends to favour higher wavenumbers compared with pure-fluid convection (Turcotte & Schubert Reference Turcotte and Schubert2014). In spherical geometry, the wavenumber
$a$
no longer varies continuously but is linked to the spherical harmonic degree
$l$
, satisfying an approximate relation with the effective aspect ratio
$\lambda _l$
as
$a = \pi / \lambda _l \sim 2l$
. Substituting the relation into (4.1)–(4.2) explains why the minimum
$Ra$
occurs at
$l=1$
for the pure fluid (Yu & Wang Reference Yu and Wang2019) but at
$l=2$
for the porous medium.
Under rapid rotation (
$\textit{Ta} \gg 1$
), the flow is largely confined to slender columnar regions near the polar axis. The local dynamics is approximately two-dimensional, thereby recovering the scaling behaviour observed in plane porous-layer systems. If the wavenumber is fixed, a reasonable extrapolation based on (4.2) indicates that
$Ra_c \propto Ta$
, a consequence of the constant critical azimuthal wavenumber
$m_c = 2$
in spherical geometry. On the other hand, the sole origin of drifting instability is thermal waves rather than inertial waves because the inertial term is neglected and the momentum equation contains no time derivative. This alignment of the time scales – among the drift of the physical field, the rotation of the reference frame and the energy transport – explains why
$\omega _{c}\propto \tau$
holds. However, the time scales of energy and momentum are not aligned in pure fluid (unless
$\textit{Pr} = 1$
), and the optimal wavenumber can be freely varied to minimise
$Ra$
, leading to more complex asymptotic scaling laws (Sánchez et al. Reference Sánchez, Garcia and Net2016a
).
A crucial question arises: why does the Darcy model favour a fixed
$m_c=2$
, while in pure fluids
$m_c$
tends to infinity as
$\textit{Ta}$
increases? For the pure fluid, rapidly rotating flow is governed by the Taylor–Proudman theorem: variations in vorticity along the rotation axis are strongly suppressed, the flow becomes quasi-two-dimensional, and it organises into Taylor columns aligned with the rotation axis. This implies that under low-viscosity conditions, the presence of inertial waves maintains
$\partial \boldsymbol{u} / \partial z \approx 0$
(Chandrasekhar Reference Chandrasekhar1961). Perturbations are confined within sheet-like structures and they must conform to the spherical geometry, resulting in a horizontal length scale for convection that is much smaller than the vertical scale. The scale of the critical mode is determined by a balance between inertial waves and thermal forcing and exhibits a specific power-law dependence on
$\textit{Ta}$
. In the strong-rotation limit, the critical azimuthal wavenumber
$m_c$
diverges, and the flow exhibits highly localised, slender columnar structures (Busse Reference Busse1970).
In this model, the Darcy drag, which is proportional to velocity, directly provides dissipation and balances the Coriolis force and thermal buoyancy. The reduction in the order of the governing equations substantially enhances damping and disrupts the quasi-geostrophic relation, thereby eliminating the need to form small-scale structures for mediating energy cascades or matching viscous boundary layers. Consequently, the Taylor–Proudman theorem does not apply to the Darcy model (
$\partial \boldsymbol{u} / \partial z$
in (2.4) is no longer a small quantity). The system favours the large-scale mode of lowest energy. This mode couples optimally with the spherical geometry, forming a global columnar structure that spans the polar regions, rather than localised, fine-scale rolls.
From an operator-theoretic perspective, the fundamental bifurcation in mode selection between the two systems under rotation originates from the intrinsic mathematical structure of their governing equations. The absence of the higher-order viscous term renders the rotating Darcy problem a lower-order eigenvalue problem, whose eigenvalue spectrum is discrete and smooth, and the system invariably selects the globally energy-optimal fundamental mode. In contrast, the momentum equation for a pure fluid contains the viscous term, which constitutes a singular perturbation problem in the small-
$Ek$
limit. Its eigenvalue spectrum becomes continuous in wavenumber space. This allows the system to select modes with arbitrarily large wavenumbers in order to optimise heat transport through boundary layers.
5. Conclusion
This study conducts the first linear stability analysis of thermal convection in a rotating porous sphere. Rotation exerts a notable suppressing effect on the onset of convection, which is manifested as a power-law relationship between the critical Rayleigh number
$ Ra_c$
and the Taylor number
$ Ta$
. A new triple balance among Darcy drag, Coriolis force and thermal buoyancy emerges under strong rotation, causing the classical Proudman–Taylor theorem to break down in porous media. This balance mechanism suppresses the development of high-wavenumber disturbances and instead promotes the formation of large-scale columnar structures aligned with the rotation axis, which is quite different from the inertial-wave mechanism that drives small-scale convective cells concentrated near the equatorial region in pure fluids. The system develops an Ekman-like boundary layer in the vicinity of the polar axis rather than near the spherical boundary, while flow in other regions is almost entirely suppressed. This new dynamical balance fundamentally weakens the growth advantage of high-wavenumber modes, thereby stabilising the large-scale convective structure. Consequently, over almost the extensive parameter range explored, the porous system consistently selects the azimuthal wavenumber
$m_c = 2$
. This contrasts sharply with the behaviour in ordinary rotating fluids, where the preferred wavenumber increases indefinitely with the rotation rate.
In the small Darcy number approximation considered in this study, the inertial term in the momentum equation is neglected. Our results therefore apply to the pure Darcy regime where Prandtl number effects are not dominant, corresponding to typical physical conditions such as rapidly rotating planetesimals or undifferentiated asteroids with high permeability. Future work may build upon this foundation by incorporating modest
$\textit{Da}$
and
$\textit{Pr}$
, as well as viscous effects, aiming to understand the transition from porous-medium convection to classical fluid convection and to enable more realistic planetary-scale simulations.
Acknowledgements
The insightful and constructive comments by two anonymous reviewers are highly appreciated. This work was supported by the Funding Project of Director of China Institute of Atomic Energy (grant no. PA010271225776).
Declaration of interests
The authors report no conflict of interest.
Author contributions
The authors jointly designed the study and analysed the results.
