1. Introduction
In the oceans, convective transport is driven, in part, by the presence of large-scale gradients of temperature and of salinity. As water is heated, its density decreases, resulting in an upward buoyancy force. Conversely, when the salt content of a fluid parcel increases, its density increases and yields a downward buoyancy force. The fact that these phenomena take place simultaneously and that salt diffuses substantially slower than heat creates a phenomenologically rich physical system referred to as doubly diffusive convection (Schmitt Reference Schmitt1994; Radko Reference Radko2013). Due to its importance in sub-planetary-scale oceanography, but also as a model problem for pattern formation in fluids, doubly diffusive convection is often studied in a planar fluid layer. However, it is also of relevance in a number of astrophysical problems (Garaud Reference Garaud2018). For example, doubly diffusive instabilities are important flow mechanisms in the Earth’s core (Trümper et al. Reference Trümper, Breuer and Hansen2012; Monville et al. Reference Monville, Vidal, Cébron and Schaeffer2019) and doubly diffusive layering has been shown to inhibit vertical transport in the subsurface oceans of Jupiter’s and Saturn’s icy moons (Wong et al. Reference Wong, Hansen, Wiesehöfer and McKinnon2022).
The initial discovery of steady spatially localised convection consisting of rolls surrounded by large areas of still, conductive fluid (Ghorayeb & Mojtabi Reference Ghorayeb and Mojtabi1997), subsequently coined convectons (Blanchflower Reference Blanchflower1999), sparked interest that resulted in the exhaustive characterisation of their bifurcation diagram in natural doubly diffusive convection, i.e. in vertically extended domains where convection is driven by horizontal gradients of temperature and salinity (Bergeon & Knobloch Reference Bergeon and Knobloch2008; Beaume et al. Reference Beaume, Bergeon and Knobloch2013a ). Convectons were also studied in binary fluid convection (the doubly diffusive convection that takes place in a horizontal fluid layer driven by vertical gradients) in the presence (Mercader et al. Reference Mercader, Batiste, Alonso and Knobloch2009, Reference Mercader, Batiste, Alonso and Knobloch2011) and in the absence (Beaume, Bergeon & Knobloch Reference Beaume, Bergeon and Knobloch2011) of the Soret effect. The bifurcation diagram in the latter case shows two branches of convectons which we recomputed in figure 1. The solutions along these branches possess different parity, with either an odd or an even number of rolls, respectively corresponding to an even or an odd number of local extrema of the temperature perturbation profile. The solution branches neatly organise themselves to undergo a series of saddle-node bifurcations between two well-defined values of the thermal Rayleigh number. This behaviour is called snaking and is accompanied by a spatial growth mechanism, whereby a convecton nucleates one roll on either of its sides every time the branch is continued upward past two saddle nodes. This can be observed by comparing successive left saddle-node states in figure 1. An important consequence of the snaking structure is that, for a given range of parameter values, a large number of spatially localised states of different sizes coexist. For an exhaustive description of the formation of spatially localised structures, we refer the reader to the recent reviews by Knobloch (Reference Knobloch2015) and by Bramburger, Hill & Lloyd (Reference Bramburger, Hill and Lloyd2025) and to early work on homoclinic snaking in the Swift–Hohenberg equation (SHE) (Burke & Knobloch Reference Burke and Knobloch2006).
(Left) Bifurcation diagram of the branches of spatially localised convection in planar doubly diffusive convection. The solution branches are represented via their domain-integrated kinetic energy
$E$
as a function of the thermal Rayleigh number
${\textit{Ra}}_T$
. (Right) Representation of the solutions at the saddle nodes with corresponding labels on the bifurcation diagram. The flow is represented using the temperature deviation/perturbation from the conductive state with red (blue) indicating positive (negative) values. Parameter values are: domain aspect ratio,
$20.1643$
; solutal Rayleigh number,
${\textit{Ra}}_S = 500$
; Prandtl number,
$Pr = 1$
; inverse Lewis number,
$\tau = 1/15$
. These results were originally found in Beaume et al. (Reference Beaume, Bergeon and Knobloch2011).
Following Mercader et al. (Reference Mercader, Batiste, Alonso and Knobloch2011) and Beaume et al. (Reference Beaume, Bergeon and Knobloch2011, Reference Beaume, Bergeon and Knobloch2013a , Reference Beaume, Knobloch and Bergeonb ), a good understanding of the formation of spatially localised structures in doubly diffusive convection is now available for simple planar geometries. For this reason, recent effort has veered away from the idealised set-ups initially used. The effect of spatial symmetry breaking by a change in the solutal boundary condition in natural doubly diffusive convection revealed that convectons became asymmetric and started travelling (Lo Jacono, Bergeon & Knobloch Reference Lo Jacono, Bergeon and Knobloch2017). Spatial symmetry breaking was also explored in binary fluid convection through the inclination of the horizontal fluid layer which resulted in dramatic changes to the solution space, including the discovery of a new type of spatially localised solution (Mercader et al. Reference Mercader, Batiste, Alonso and Knobloch2019; Alonso, Batiste & Mercader Reference Alonso, Batiste and Mercader2022). Other work investigating the role of temporal symmetry has focused on travelling localised states and their interaction (Watanabe et al. Reference Watanabe, Iima and Nishiura2012, Reference Watanabe, Iima and Nishiura2016). Beyond understanding the role of symmetries, recent investigations into natural doubly diffusive convection have also sought to understand the importance of subcriticality as well as that of the balance between thermal and solutal effects. In particular, Tumelty, Beaume & Rucklidge (Reference Tumelty, Beaume and Rucklidge2023) revealed the persistence of spatially localised states into the supercritical regime, albeit with an altered bifurcation scenario, while Tumelty, Beaume & Rucklidge (Reference Tumelty, Beaume and Rucklidge2025) explored the formation of spatially localised states away from the balanced case where thermal and solutal effects are equally important, identifying the emergence of convectons via cusps in the thermally dominated regime and a smooth connection to domain-filling convection in the solutally dominated one. Understanding how robust convectons are to such departures from the idealised set-ups in which they were initially investigated is fundamental to transfer knowledge to less ideal systems.
It is now natural to turn to other geometries in order to understand spatial localisation more generally. Although the formation of spatially localised states in polar geometries has been investigated for the SHE (McCalla & Sandstede Reference McCalla and Sandstede2010, Reference McCalla and Sandstede2013; Verschueren, Knobloch & Uecker Reference Verschueren, Knobloch and Uecker2021), to the authors’ knowledge, it has yet to be studied in a spherical geometry. The absence of such studies is, however, not withstanding the fact that examples of spherical physical systems susceptible to harbouring such states do exist. For example, the wrinkling of elastic rings has recently been investigated to show the emergence of spatially localised folds (Foster et al. Reference Foster, Verschueren, Knobloch and Gordillo2022) and fingering patterns (Foster & Knobloch Reference Foster and Knobloch2023), and is therefore suggestive of the existence of a spatial localisation mechanism in the spherical extension of the problem. The wrinkling of spherical membranes indeed shows a hysteretic transition between the unpatterned surface and hexagons (Stoop et al. Reference Stoop, Lagrange, Terwagne, Reis and Dunkel2015), much like the scenario observed in the Rosensweig instability that led to the discovery of snaking in magnetically forced ferrofluids (Lloyd et al. Reference Lloyd, Gollwitzer, Rehberg and Richter2015). In an astrophysical context, the geodynamo is believed to originate via a subcritical bifurcation (Rincon Reference Rincon2019) and, at its laminar to turbulent boundary in phase space, also admits unstable travelling wave solutions (Skene, Marcotte & Tobias Reference Skene, Marcotte and Tobias2024). Given the role played by localised flow states in transition to turbulence in shear flows (Duguet, Schlatter & Henningson Reference Duguet, Schlatter and Henningson2009; Schneider, Gibson & Burke Reference Schneider, Gibson and Burke2010; Gibson & Schneider Reference Gibson and Schneider2016; Pershin, Beaume & Tobias Reference Pershin, Beaume and Tobias2019), it is possible that similar states could help to understand the transition to magnetohydrodynamic turbulence in spherical geometries.
To understand the influence of the geometry on spatial localisation in fluids, this paper asks: Do spatially localised solutions of the Navier–Stokes equations exist for convection in a spherical shell? And, if so, what do these solutions look like and how do they form? In order to enable useful links with established literature, we consider doubly diffusive convection and use the same set-up as in Beaume et al. (Reference Beaume, Bergeon and Knobloch2011) with the only change being the geometry: in previous work, a two-dimensional planar fluid layer with spatially periodic boundary conditions was used, while we are posing the equations and boundary conditions on an axisymmetric spherical shell.
In § 2 the problem formulation is presented. We then introduce, in § 3, the types of solution we found in order to provide their spherical representation and establish terminology. We go back to the basics in § 4, where we solve the linear stability problem for the conduction state, followed, in § 5, by the presentation of our main results. Finally, this paper terminates with a discussion in § 6.
2. Formulation
We consider convection driven by competing buoyancy forces generated by temperature and solute concentration gradients between spheres of radii
$R_1$
and
$R_2\gt R_1$
. The forcing gradients are imposed via a constant temperature and solute concentration at the sphere walls. The inner sphere is set to be hotter and richer in solute such that the temperature and solute concentration differences are
$\Delta T = T_1 - T_2 \gt 0$
and
$\Delta S = S_1 - S_2 \gt 0$
, where
$T_1$
and
$S_1$
(respectively
$T_2$
and
$S_2$
) are taken at the wall of the sphere of radius
$R_1$
(respectively
$R_2$
). We assume that the fluid contained within the resulting spherical shell is well modelled under the Boussinesq approximation: its density is taken to be
$\rho _{\!f}$
everywhere except in the gravitational force term, where it varies linearly with the fluid temperature and solute concentration:
\begin{equation} \rho (T,S) = \rho _{\!f} [1 - \beta _T(T - T_2) + \beta _S (S - S_2) ]. \end{equation}
Here, the expansion coefficients
$\beta _T$
and
$\beta _S$
are positive constants and
$\rho _{\!f}$
is the density at the reference temperature
$T_2$
and solute concentration
$S_2$
. The gravity field is spherically symmetric and follows from the assumption that the inner core (
$r\lt R_1$
) density is much greater than that of the spherical shell fluid:
\begin{equation} \boldsymbol{g} = - g_0 \frac {R_1^2}{r^2} \boldsymbol{\hat {r}}, \end{equation}
where
$g_0$
is the acceleration due to gravity on the inner sphere surface. Choosing the length scale
$d = R_2 - R_1$
, the thermal diffusion time scale
$d^2/\kappa$
and the temperature and solutal scales
$\Delta T$
and
$\Delta S$
, the governing equations may be written in the following non-dimensional form:
\begin{align} \frac {1}{Pr} \left ( \frac {{\textrm D} \boldsymbol{u}}{{\textrm D} t} + \boldsymbol{\nabla }P \right ) &= g(r) \hat {\boldsymbol{r}} ({\textit{Ra}}_T T - {\textit{Ra}}_S \, S) + {\nabla} ^2 \boldsymbol{u}, \end{align}
\begin{align} \frac {\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }T &= {\nabla} ^2 T, \end{align}
\begin{align} \frac {\partial S}{\partial t} + {\boldsymbol u} \boldsymbol{\cdot }\boldsymbol{\nabla }S &= \tau {\nabla} ^2 S, \end{align}
\begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} &= 0, \end{align}
where
$r$
is now non-dimensional,
$\boldsymbol u$
is the velocity field,
$P$
is the pressure,
$T$
is the temperature,
$S$
is the solute concentration and
$t$
is the time. In writing these equations, we have introduced a number of parameters:
\begin{align} \text{Prandtl number: } & \textit{Pr} = \displaystyle \frac {\nu }{\kappa }, \end{align}
\begin{align} \text{Lewis number: } & \tau = \displaystyle \frac {\kappa _S}{\kappa }, \end{align}
\begin{align} \text{(Thermal) Rayleigh number: } & {\textit{Ra}}_T = \displaystyle \frac {g_0 \beta _T \Delta T d^3}{\nu \kappa } \gt 0, \end{align}
\begin{align} \text{Solutal Rayleigh number: } & {\textit{Ra}}_S = \displaystyle \frac {g_0 \beta _S \Delta S d^3}{\nu \kappa } \gt 0, \end{align}
\begin{align} \text{Aspect ratio: } & \varGamma = \displaystyle \frac {R_1 \pi }{d}, \end{align}
\begin{align} \text{Gravitational factor: } & g(r) = \displaystyle \frac {R_1^2}{r^2 d^2}, \end{align}
where
$\nu$
is the kinematic viscosity,
$\kappa$
the thermal diffusivity and
$\kappa _S$
the solutal diffusivity. The thermal and solutal Rayleigh numbers measure the strength of the buoyancy force resulting from the imposed temperature and solutal gradients at the sphere boundaries, the Prandtl and inverse Lewis numbers determine the dominant diffusive processes while
$\varGamma$
and
$g(r)$
are geometric parameters. We specify no-slip velocity boundary conditions for
$\boldsymbol{u}$
and assume perfectly conducting spherical shell walls:
\begin{align} \boldsymbol{u}(r=r_1,\theta ,\varphi ) = 0, \quad \boldsymbol{u}(r=r_2,\theta ,\varphi ) = 0, \end{align}
\begin{align} T(r=r_1,\theta ,\varphi ) = 1, \quad T(r=r_2,\theta ,\varphi ) = 0, \end{align}
\begin{align} S(r=r_1,\theta ,\varphi ) = 1, \quad S(r=r_2,\theta ,\varphi ) = 0, \end{align}
where
$r_1 = R_1/(R_2 - R_1)$
(
$r_2 = R_2/(R_2 - R_1)$
) is the location of the inner (outer) sphere wall,
$\theta$
is the polar coordinate and
$\varphi$
is the azimuthal coordinate. We look for longitudinally invariant solutions that are devoid of zonal flow and therefore impose
$\partial /\partial \varphi = 0$
and
$u_\varphi = 0$
. This restriction allows the use of the following pole boundary conditions:
\begin{equation} \frac {\partial S}{\partial \theta } = \frac {\partial T}{\partial \theta } = \frac {\partial u_r}{\partial \theta } = u_{\theta } = 0, \quad \text{for} \quad \theta = 0,\pi , \end{equation}
where
$\theta = 0$
(respectively
$\pi$
) can be thought of as the north (respectively south) pole of our spherical geometry.
The problem solved by Beaume et al. (Reference Beaume, Bergeon and Knobloch2011) in a planar geometry, and whose results are reproduced in figure 1, is different from the one we pose on the spherical shell in three important ways: (i) the system no longer has the top-down reflection symmetry of the plane layer, (ii) the geometry is curved as opposed to flat and (iii) the
$\theta$
boundary conditions are of a nature different from that of the periodic boundary conditions used in the planar problem.
System (2.6)–(2.16) admits the base state solution
\begin{equation} \boldsymbol{u} = 0, \quad T_0 = S_0 = - A_T/r + B_T, \end{equation}
where
$A_T = r_1 r_2/(r_1 - r_2)$
and
$B_T = r_1/(r_1 - r_2)$
. To understand the symmetries of the system, we introduce the convective variables
$\varTheta$
and
$\varSigma$
such that
$T = T_0 + \varTheta$
and
$S = S_0 + \varSigma$
. As discussed by Beltrame & Chossat (Reference Beltrame and Chossat2015), the full three-dimensional system is equivariant under the antipodal symmetry
$\boldsymbol{r} \to -\boldsymbol{r}$
. Using longitudinal invariance, this symmetry reads
\begin{equation} R_{\textit{AP}}: (r,\theta ) \to (r, \pi - \theta ), \quad (u_r, u_{\theta }, \varTheta , \varSigma ) \to (u_r, -u_{\theta }, \varTheta , \varSigma ) \end{equation}
and can be understood as a simple reflection about the equator. Spatial reversibility, here associated with
$R_{\textit{AP}}$
, is an important feature of systems that admit stationary spatially localised solutions (Burke, Houghton & Knobloch Reference Burke, Houghton and Knobloch2009). For example, if a solution is conductive at the north pole and displays convection rolls at the equator, then the reverse connection can be found from the equator to the south pole, yielding a spatially localised convection structure at the equator.
In practice, the fact the flow is two-dimensional and incompressible allows the introduction of a streamfunction
$\psi$
:
\begin{equation} \boldsymbol{u} = \boldsymbol{\nabla }\times \left (0,0,\frac {\psi (r,\theta ,t)}{r} \right ) = \bigg ( \frac {1}{r^2 \sin \theta } \frac {\partial (\psi \sin \theta )}{\partial \theta }, - \frac {1}{r}\frac {\partial \psi }{\partial r}, 0 \bigg ), \end{equation}
which is chosen following Marcus & Tuckerman (Reference Marcus and Tuckerman1987) so that the resulting primitive variable system projects onto a sine/cosine basis. Using this definition of the streamfunction, the primitive variable system (2.6) can then be written as
\begin{align} \frac {\partial D^2 \psi }{\partial t} + \mathcal{J}( \psi , \psi ) = \textit{Pr g}(r) & \bigg ( {\textit{Ra}}_T \frac {\partial \varTheta }{\partial \theta } - {\textit{Ra}}_S \frac {\partial \varSigma }{\partial \theta } \bigg ) + \textit{Pr D}^2 D^2 \psi , \end{align}
\begin{align} \frac {\partial \varTheta }{\partial t} + J( \psi , T_0 + \varTheta ) & = {\nabla} ^2 \varTheta , \end{align}
\begin{align} \frac {\partial \varSigma }{\partial t} +J( \psi , S_0 + \varSigma ) & = \tau {\nabla} ^2 \varSigma , \end{align}
where
\begin{align} J(\psi ,f) & = \frac {1}{r^2 \sin \theta } \frac {\partial (\psi \sin \theta ) }{\partial \theta } \frac {\partial f }{\partial r} - \frac {1}{r^2} \frac {\partial \psi }{\partial r} \frac {\partial f }{\partial \theta }, \end{align}
\begin{align} \mathcal{J}( \psi , \psi ) & = \frac {\partial }{\partial r} \bigg ( ( r^2 u_r) \frac {D^2 \psi }{r^2} \bigg ) + \frac { \partial }{\partial \theta } \bigg ( (r u_{\theta }) \frac {D^2 \psi }{r^2} \bigg ) \end{align}
and
\begin{align} {\nabla} ^2 & = \frac {1}{r^2} \frac {\partial }{\partial r} \bigg ( r^2 \frac {\partial }{\partial r} \bigg ) + \frac {1}{r^2 \sin \theta } \frac {\partial }{\partial \theta } \bigg (\sin \theta \frac {\partial }{\partial \theta } \bigg ), \end{align}
\begin{align} D^2 & = \frac {\partial ^2 }{\partial r^2} + \frac {1}{r^2 \sin \theta } \frac {\partial }{\partial \theta }\bigg (\sin \theta \frac {\partial }{\partial \theta } \bigg ) - \frac {1}{r^2 \sin ^2 \theta }. \end{align}
In this formulation, the spherical shell boundary conditions become
\begin{equation} \psi = \frac {\partial \psi }{\partial r} = \varTheta = \varSigma = 0, \quad \hbox{for} \; \; r = r_1,r_2, \end{equation}
while the pole conditions are
\begin{equation} \psi = \frac {\partial \varTheta }{\partial \theta } = \frac {\partial \varSigma }{\partial \theta }=0, \quad \text{for} \quad \theta = 0,\pi . \end{equation}
Following Beaume et al. (Reference Beaume, Bergeon and Knobloch2011), who studied the formation of spatially localised convection states in doubly diffusive convection over a periodic planar geometry, we set
$\tau = 1/15$
and
$Pr = 1$
. We use the thermal Rayleigh number as the control parameter against which the bifurcation diagrams of steady states will be represented. However, to enable us to fully explain the localised pattern formation at play in our system, we have had to vary also the solutal Rayleigh number.
3. Solutions
We now outline the main types of localised solutions discussed in this paper. We concentrate exclusively on the simplest types of localised solutions constituting straightforward homoclinic or heteroclinic connections in space. More complex connections exist but are beyond the scope of our work.
The first family of spatially localised solutions consists of solutions that respect the symmetry
$R_{\textit{AP}}$
, i.e. their reflection about the equator returns the same solution. Examples of such solutions are shown in terms of the solutal field
$S$
in figure 2. Since they display identical variations in salinity (and in temperature) on either side of the equator, they can be understood as even-parity states with respect to the equator. Consequently, no convection roll can straddle the equator and, instead, a pair of counter-rotating rolls is present and separated by the equatorial plane on which a purely vertical (i.e. aligned with gravity) flow is observed. We found two types of
$R_{\textit{AP}}$
-symmetric localised states. Solutions represented in figure 2(a,b) are characterised by the presence of convection around the equator while the fluid near the poles is mostly still. We refer to these states as equatorial-convectons. The convecton of figure 2(a) displays an upward flow at the equator and is referred to as
$C+$
. Conversely, the convecton of figure 2(b) is referred to as
$C-$
due to the downward flow it produces at the equator. We have also identified localised solutions where the convection is centred around both poles and the equatorial region displays mostly quiescent fluid. Examples of such states are shown in figure 2(c,d). Owing to the imposed longitudinal invariance, the poles create a constraint on the flow and they bear similarity with solutions found in planar geometry with endwalls, where they are known as anticonvectons (Mercader et al. Reference Mercader, Batiste, Alonso and Knobloch2011; Tumelty et al. Reference Tumelty, Beaume and Rucklidge2025). We adopt this terminology here. The anticonvectons in figure 2(c) connect with the
$C+$
convectons (figure 2
a), as shown in § 5.1, so we refer to them as
$A+$
anticonvectons. A similar connection is found between
$C-$
convectons (figure 2
b) and the anticonvectons of figure 2(d), which are referred to as
$A-$
anticonvectons.
Equatorially symmetric spatially localised states of convection represented by the salinity (red for saltier) in a section of the spherical shell passing by the poles. Four different types of solutions are shown: (a)
$C+$
equatorial-convectons, (b)
$C-$
equatorial-convectons, (c)
$A+$
anticonvectons and (d)
$A-$
anticonvectons. These states are computed at
${\textit{Ra}}_S = 400, \varGamma = 8.9224$
for different
${\textit{Ra}}_T$
and are discussed further in § 5.1.
In addition to symmetric convectons, we are also interested in spatially localised solutions that break
$R_{\textit{AP}}$
. The simplest such solutions are shown in figure 3 and consist of states that display convection near one pole only. We call these pole-convectons. We found two families of such states, which we differentiate by the direction of the flow at the pole. We refer to those with an upward flow at the pole as
$+$
pole-convectons (figure 3
a), while those with a downward flow at the pole are referred to as
$-$
pole-convectons (figure 3
b). Similar states consisting of convection near the south pole and conduction near the north pole also exist but will not be reported here: they can be obtained by applying the symmetry
$R_{\textit{AP}}$
onto the solutions presented and are, consequently, dynamically equivalent. Pole-convectons share a similar spatial structure with localised wall states found in planar binary fluid convection (Mercader et al. Reference Mercader, Batiste, Alonso and Knobloch2011).
Symmetry-breaking spatially localised states of convection represented by the salinity (red for saltier) in a section of the spherical shell passing by the poles. Two different types of solutions are shown: (a)
$+$
pole-convecton and (b)
$-$
pole-convecton. These states are computed at
${\textit{Ra}}_S = 150, \varGamma = 10.029$
for different values of
${\textit{Ra}}_T$
and are discussed further in § 5.2.
Before discussing how these six types of solutions arise, we briefly comment on the consequences of imposing longitudinal invariance. As made apparent by figures 2 and 3, an axis of symmetry passing through the sphere’s north and south poles now exists. The existence of this axisymmetry, in addition to making computations easier and simplifying phase space, implies that convection rolls near the poles take a doughnut- or torus-like form and convect fluid away from the spherical shell boundaries over a smaller surface area than a convection roll of a similar cross-section located near the equator. To simplify visualisation, we present our solutions on a plane where the south pole can be thought to be at the left end and the north pole at the right end. While this geometrical feature is easy to visualise in figures 2 and 3, it is not that clear on the flattened representations that follow. The reader should therefore keep in mind the physical differences between pole (outer regions in the flattened figures) and equatorial rolls.
4. Linear theory
To solve the linear problem, we chose a different spherical streamfunction formulation:
\begin{equation} \boldsymbol{u} = \boldsymbol{\nabla }\times \left (0,0,\frac {\varPsi (r,\theta ,t)}{r \sin \theta }\right ) = \bigg ( \frac {1}{r^2 \sin \theta }\frac {\partial \varPsi }{\partial \theta }, \frac {-1}{r \sin \theta }\frac {\partial \varPsi }{\partial r}, 0 \bigg ), \end{equation}
such that the linear system admits separable solutions in
$r$
and
$\theta$
(Beltrame & Chossat Reference Beltrame and Chossat2015). Upon using this streamfunction in system (2.6) and linearising the resulting equations about the conduction state (2.17), one obtains the perturbation system:
\begin{equation} \bigg ( \frac {\partial }{\partial t} \mathcal{M} - \mathcal{L} \bigg ) \boldsymbol{\tilde {X}} = 0, \quad \boldsymbol{\tilde {X}} = \begin{pmatrix} \tilde {\varPsi } & \tilde {\varTheta } & \tilde {\varSigma } \end{pmatrix}^T, \end{equation}
where
\begin{equation} \mathcal{M} = \begin{pmatrix} \displaystyle \frac {1}{Pr} \, \hat {D}^2 & 0 & 0\\ 0 & \mathbb{I} & 0 \\ 0 & 0 & \mathbb{I} \end{pmatrix}, \quad \mathcal{L} = \begin{pmatrix} \hat {D}^2 \hat {D}^2 & {\textit{Ra}}_T \, g(r) \sin \theta \displaystyle \frac {\partial }{\partial \theta } & - {\textit{Ra}}_S \, g(r) \sin \theta \displaystyle \frac {\partial }{\partial \theta },\\ -\displaystyle \frac {T^{\prime}_0}{r^2 \sin \theta } \frac {\partial }{\partial \theta } & {\nabla} ^2 & 0 \\[6pt] -\displaystyle \frac {S^{\prime}_0}{r^2 \sin \theta } \frac {\partial }{\partial \theta } & 0 & \tau {\nabla} ^2 \end{pmatrix}\!, \end{equation}
and
$\hat {D}^2 = D^2 + 1/r^2 \sin ^2 \theta $
. In writing (4.3), we have introduced the radial derivatives of the base state:
$T^{\prime}_0 = {\textrm d} T_0 / {\textrm d} r$
and
$S^{\prime}_0 = {\textrm d} S_0 / {\textrm d} r$
. Owing to the choice of streamfunction (4.1), system (4.2) admits separable eigenfunctions of the form
\begin{align} \tilde {\varPsi } &= {\textrm e}^{\lambda t} \, \varPsi _{\ell }(r) \, \sin \theta \frac {\partial P_{\ell }(\cos \theta )}{\partial \theta }, \end{align}
\begin{align} \tilde {\varTheta } &= {\textrm e}^{\lambda t} \, \varTheta _{\ell }(r) \, P_{\ell }(\cos \theta ), \end{align}
\begin{align} \tilde {\varSigma } &= {\textrm e}^{\lambda t} \, \varSigma _{\ell }(r) \, P_{\ell }(\cos \theta ). \end{align}
In writing these eigenfunctions, we have introduced the temporal growth rate
$\lambda \in \mathbb{C}$
as well as
$P_{\ell }$
, the Legendre polynomial of degree
$\ell$
. Substituting for
$\varPsi _{\ell }(r), \varTheta _{\ell }(r), \varSigma _{\ell }(r)$
and using the fact that
$T^{\prime}_0(r) = S^{\prime}_0(r)$
allows system (4.2) to be reduced to an eighth-order ordinary differential equation:
\begin{equation} \left( \lambda - \tau {\nabla} ^2_{\ell }\right) \left(\lambda - {\nabla} ^2_{\ell }\right) g(r)^{-1} \left ( \frac {\lambda }{Pr} D_{\ell }^2 - D_{\ell }^2 D_{\ell }^2 \right ) \varPsi _{\ell } = \left ( \lambda - {\nabla} ^2_{\ell } \right ) \tilde {Ra} \frac {T^{\prime}_0}{r^2} \ell (\ell +1) \varPsi _{\ell }, \end{equation}
where
$\tilde {Ra} = {\textit{Ra}}_T - {\textit{Ra}}_S$
and
\begin{equation} D^2_{\ell } = \frac {\partial ^2 }{\partial r^2} - \frac {\ell (\ell + 1) }{r^2}, \quad {\nabla} ^2_{\ell } = \frac {1}{r^2} \frac {\partial }{\partial r} \bigg ( r^2 \frac {\partial }{\partial r} \bigg ) - \frac {\ell (\ell +1)}{r^2}. \end{equation}
Neutral curves demarcating the onset of convection via steady state (
$\lambda = 0)$
and oscillatory state (
$\lambda = \pm {\textrm i} \omega$
) bifurcations, where
$\ell$
denotes the Legendre polynomial associated with the latitudinal mode. These curves were computed for
${\textit{Ra}}_S=400$
. According to (4.7), the critical value of
${\textit{Ra}}_T$
varies trivially with
${\textit{Ra}}_S$
as the critical value of
${\textit{Ra}}_T - {\textit{Ra}}_S$
remains constant.
Letting
$\lambda = s \pm {\textrm i} \omega$
and setting
$s=0$
, one can solve for the critical
$\tilde {Ra}$
for a given aspect ratio
$\varGamma$
and Legendre polynomial degree
$\ell$
. Figure 4 shows the resulting neutral stability curves as the aspect ratio
$\varGamma$
varies for fixed
${\textit{Ra}}_S$
. The advantage of the above method is that the Legendre polynomials in
$\theta$
are decoupled, allowing one to track straightforwardly the neutral stability curve of the various latitudinal modes. As
${\textit{Ra}}_T$
increases from zero, the conduction state first loses stability to time-dependent convection via a Hopf bifurcation, indicated by the neutral stability curves for which
$\lambda = \pm {\textrm i} \omega$
. Bifurcations to steady states only come in at larger values of the thermal Rayleigh number.
Although we present results between these two types of bifurcations separately, the dynamics they generate may impact one another. When posed on a flat plane, two-dimensional doubly diffusive convection possesses the symmetries of the group
$O(2) \times Z_2$
which leads to simultaneous emergence of standing and travelling waves (Crawford & Knobloch Reference Crawford and Knobloch1991). For the values of the parameters we chose, Beaume et al. (Reference Beaume, Bergeon and Knobloch2011) found that the travelling waves bifurcate subcritically but soon turn around a saddle node and extend to larger values of the thermal Rayleigh number until they connect with a branch of spatially periodic steady convection, to which they transfer their stability. On the one hand, the equations posed in a spherical shell are not symmetric with respect to horizontal translations or vertical reflections. They simply possess the symmetry
$Z_2 = \{ I, R_{\textit{AP}}\}$
and Hopf bifurcations will only create standing waves. On the other hand, the larger the value of
$\varGamma$
, the more the present spherical system approaches the planar system of Beaume et al. (Reference Beaume, Bergeon and Knobloch2011). In the planar problem, steady-state bifurcations are subcritical and give rise to convectons, which is the observation that forms the basis of our investigation here.
We decided on the value of the aspect ratio
$\varGamma$
in the following way. First of all, we note that equatorial-convectons and anticonvectons respect
$R_{\textit{AP}}$
and are generated by even Legendre modes. We study them in the subspace generated by the even latitudinal modes. We set
$\ell = 10$
to enable a sufficiently large number of convection rolls to develop and fixed
$\varGamma = 8.9224$
. This value was chosen so that the next even-parity steady-state bifurcations are the farthest away: for
${\textit{Ra}}_S = 400$
, the
$\ell = 10$
bifurcation is located at
${\textit{Ra}}_{T,10} \approx 8351.53$
while the bifurcations to
$\ell = 8$
and
$\ell = 12$
modes only occur at
${\textit{Ra}}_{T,8} \approx {\textit{Ra}}_{T,12} \approx 8479.92$
. Pole-convectons break the symmetry and their search requires the full solution space. We set
$\ell = 11$
in order to enable the formation of structures that break
$R_{\textit{AP}}$
. We chose
$\varGamma = 10.029$
, the value of the aspect ratio for which the next steady-state bifurcations (towards
$\ell = 10$
and
$\ell = 12$
) are the furthest away. The steady-state bifurcation to
$\ell = 11$
modes is found at
${\textit{Ra}}_{T,11} \approx 8275.70$
for
${\textit{Ra}}_S = 400$
. The
$\ell = 10$
and
$\ell = 11$
eigenvectors that are responsible for these steady-state bifurcations are shown in figure 5. Due to its parity, the
$\ell = 10$
mode generates either an upward or a downward flow at both poles as well as a flow in the opposite direction at the equator. Conversely, the
$\ell = 4n$
modes, where
$n = 1, 2, \ldots$
, trigger convective structures that display flows in the same direction at both poles and over the equator. The even modes preserve
$R_{\textit{AP}}$
and, thus, generate what may be thought of as transcritical bifurcations. Like all odd modes, the
$\ell = 11$
mode is characterised by an equator-straddling roll and by opposite vertical flows at the poles. They break the
$R_{\textit{AP}}$
symmetry and are thus responsible for the occurrence of pitchfork bifurcations of the conduction state.
(Left) Even-parity mode
$\ell = 10, \varGamma = 8.9224$
and (right) odd-parity mode
$\ell = 11, \varGamma = 10.029$
at onset represented by the perturbations (top panels) in the in-plane velocity, (middle panels) in the temperature and (bottom panels) in the salinity fields. The velocity is shown by arrows indicating the magnitude and direction of the flow while positive (negative) values of the temperature and salinity perturbations are shown in red (blue). The even-parity mode respects the
$R_{\textit{AP}}$
symmetry while the odd-parity one breaks it. Both modes are spatially modulated such that their maximum occurs at the poles and minimum at the equator.
In addition to the symmetry group
$Z_2$
, system (2.6) also possesses a hidden symmetry due to the operator in (4.2) being self-adjoint. As discussed by Beltrame & Chossat (Reference Beltrame and Chossat2015), this hidden symmetry, here due to
$g(r) \sim r^{-2},\ T^{\prime}_0(r) \sim r^{-2},\ S^{\prime}_0(r) \sim r^{-2}$
, implies that the bifurcation to even-parity solutions is a codimension-2 bifurcation corresponding to the collision between the transcritical bifurcation from the conduction state and the saddle-node bifurcation of the supercritically emerging branch. We refer to this transcritical-saddle-node bifurcation as a TSN bifurcation. The unfolding of this bifurcation can be carried out by varying the functional form of
$g(r)$
relative to
$T^{\prime}_0(r)$
and
$S^{\prime}_0(r)$
and is out of the scope of this paper.
Due to the spherical geometry, the eigenfunctions in the linear stability analysis consist of a spatially modulated array of convection rolls, where the polar rolls are the strongest (see figure 5). This spatial modulation is not a consequence of the pole boundary conditions (2.16), which can be thought of as Neumann boundary conditions, but rather because of the presence of sinusoidal terms in the governing equations. Physically, this modulation occurs because the contact plane between rolls near the poles has a smaller surface area than near the equator and, thus, rolls have to be stronger to transfer the same amount of heat or solute. This is not a consequence of axisymmetry but rather of the fact that convection rolls on a sphere are necessarily curved.
Further insight can be obtained by drawing analogies between the work of Lloyd & Sandstede (Reference Lloyd and Sandstede2009), who studied the SHE on an axisymmetric disc, and the spherical shell system studied here in the limit of large
$\varGamma$
. The similarities between these systems arise from the fact that we look for solutions that are axisymmetric and, thus, that the spatial dynamics near the pole resembles that of the axisymmetric disc near its origin. The eigenfunctions of the linear stability analysis around the trivial state of the axisymmetric disc problem are Bessel functions which decay radially while oscillating. In the present spherical shell problem, the counterpart eigenfunctions are Legendre polynomials whose modulus is maximal at the poles. In the limit of large
$\varGamma$
, these Legendre polynomials are closely related to Bessel functions, further emphasising the connection between the spatial modulation observed here and that studied in the axisymmetric disc. We therefore anticipate that, in the limit of large
$\varGamma$
, anticonvectons (see figure 2) and pole-convectons (see figure 3) will have a similar bifurcation scenario to that of the spot-like localised states bifurcating from the zero state in the SHE (Lloyd & Sandstede Reference Lloyd and Sandstede2009).
5. Results
Following Marcus & Tuckerman (Reference Marcus and Tuckerman1987), we solve system (2.20) using a pseudo-spectral method that employs either a sine basis or cosine basis in
$\theta$
together with a Chebyshev basis in
$r$
. We perform differentiation in collocation space for the radial direction and dealiase nonlinear terms using the 3/2-rule for latitudinal modes such that the largest resolved wavenumber is
$N_{\theta }$
. Our numerical scheme relies on the use of transforms to send physical fields into spectral space. These transforms have
$O(N_r^2 N_{\theta }\log (N_{\theta }))$
numerical complexity, where
$N_{\theta }$
(respectively
$N_r$
) refers to the number of sine/cosine modes (respectively Chebyshev collocation points) used. We perform time stepping, steady-state computation and numerical continuation using matrix-free methods following Tuckerman (Reference Tuckerman1989), Edwards et al. (Reference Edwards, Tuckerman, Friesner and Sorensen1994), Beaume (Reference Beaume2017) and Tuckerman, Langham & Willis (Reference Tuckerman, Langham and Willis2019). Our code is validated against the Dedalus code (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020) and is available on github.
We found spatially localised states in doubly diffusive convection in a spherical shell in two different ways. The first method relies on the initialisation of a Newton iteration using the eigenvectors obtained in § 4 to compute branches emerging from the primary bifurcation. We anticipate that these low-amplitude solutions will develop into localised solutions or into domain-filling states from which bifurcations to convectons can be found. This method does not allow us to find states that are disconnected from the conduction state. In order to obtain them, we use time stepping in the following way. We first obtain time-dependent behaviour for values of the thermal Rayleigh number near the primary Hopf bifurcation, then progressively drive the thermal Rayleigh number towards the regime where steady states exist. Most of our simulations approach an upper branch steady state which we then numerically continued. Throughout the remainder of this paper, we use an extension of the notation introduced by Lo Jacono, Bergeon & Knobloch (Reference Lo Jacono, Bergeon and Knobloch2010) to refer to different branches of spatially localised solutions:
$L_{\ell }^{X,\, \pm }$
, where the subscript
$\ell$
indicates the Legendre polynomial degree to which the branch is associated,
$X = A$
for anticonvectons,
$X = C$
for convectons and
$\pm$
differentiates the branches of solutions from the same family, as described in § 3.
Bifurcation from the conduction state giving rise to the branches of anticonvectons
$L_{10}^{A+}$
and
$L_{10}^{A-}$
. The bifurcation diagram shows the volume-averaged kinetic energy
$\mathcal{E}= ({1}/{2V}) \int _{-\pi }^{\pi } \int _{r_1}^{r_2} |\boldsymbol{u}|^2 \, r^2 \sin \theta \, {\textrm d}r \, {\textrm d} \theta$
as a function of the thermal Rayleigh number
${\textit{Ra}}_T$
for
${\textit{Ra}}_S = 400$
and
$\varGamma = 8.9224$
. As explained in the text, the presence of a hidden symmetry implies that this bifurcation is a codimension-2 bifurcation between a transcritical and a saddle-node bifurcation.
5.1. Symmetric localised states
Figure 6 shows the different branches of solutions which emerge at the steady-state bifurcation generated by the
$\ell =10$
mode at
${\textit{Ra}}_T = 8351.53$
for
${\textit{Ra}}_S = 400$
and
$\varGamma = 8.9224$
. As this bifurcation does not break the
$R_{\textit{AP}}$
symmetry, one would expect it to produce two branches of solutions, one supercritical and one subcritical. As explained in § 4, this bifurcation is, in fact, a TSN bifurcation where the saddle node of the supercritically emerging branch occurs at the same location as the transcritical bifurcation from the conduction state (Beltrame & Chossat Reference Beltrame and Chossat2015). The resulting TSN bifurcation produces two subcritical branches bearing symmetric solutions:
$L_{10}^{A+}$
anticonvectons and
$L_{10}^{A-}$
anticonvectons. Figure 7 shows the path of these branches as they are followed either until their termination point or until they reach large-amplitude states. Starting with
$L_{10}^{A+}$
, as the branch is continued to lower
${\textit{Ra}}_T$
, the solution becomes spatially localised (see figure 2
c). However, it does not produce any snaking like that of figure 1: the branch turns around to larger
${\textit{Ra}}_T$
and smaller energy at
${\textit{Ra}}_T \approx 3510$
before turning around again near the
$\ell = 8$
bifurcation point, where we renamed it
$L_8^{A+}$
to reflect the change in the structure of its solutions. The initial stages of the bifurcation scenario of the
$L_{10}^{A-}$
branch are similar and produce anticonvectons with an upflow at the poles (see figure 2
d). Unlike the
$L_{10}^{A+}$
branch, however, the
$L_{10}^{A-}$
branch goes on to larger energies while producing a behaviour akin to snaking. For the sake of simplicity, we hereafter use the word snaking loosely to refer to the behaviour of solution branches that undergo series of saddle-node bifurcations or oscillations in parameter space. By using time stepping at
${\textit{Ra}}_T = 4000$
, we were able to compute a stable large-amplitude equatorial-convecton in its domain-filling form. We numerically continued it to trace the
$L_{10}^{C+}$
branch in figure 7(a). This branch produces snaking but does not connect to the conduction state. Instead, it turns around a saddle node at low amplitude before undergoing complex behaviour (not shown). We were also able to compute another type of equatorial-convectons:
$L_{10}^{C-}$
(see figure 7
b). These convectons produce snaking branches but, instead of extending to large amplitude like
$L_{10}^{C+}$
, they turn back to low amplitude where they connect with the conduction state at the
$\ell = 12$
bifurcation. The behaviour of these localised state branches is not trivial and hints at the presence of imperfect bifurcations, especially in the case of
$L_{10}^{A+}$
and
$L_{10}^{C-}$
.
Bifurcation diagram of the
$L_{10}^{A+}$
anticonvecton and
$L_{10}^{C+}$
equatorial-convecton branches (a) and of the
$L_{10}^{A-}$
anticonvecton and
$L_{10}^{C-}$
equatorial-convecton branches (b). The equatorial-convectons (anticonvectons) are indicated by a solid (dotted) line. Parameters are
${\textit{Ra}}_S = 400$
and
$\varGamma = 8.9224$
.
The bifurcation scenario presented here differs from the one observed in a planar geometry with periodic boundary conditions, where the primary bifurcation is triggered by eigenmodes with discrete translational invariance and where localised states only emerge from a secondary bifurcation that modulates the amplitude of the spatially periodic solution (Beaume et al. Reference Beaume, Bergeon and Knobloch2011). Owing to the spherical geometry, the eigenfunctions are spatially modulated and favour the formation of convection rolls near the poles. Due to this spatial modulation, the aforementioned secondary bifurcation cannot be achieved and the branch of spatially localised states with convection rolls near the poles,
$L_{10}^{A \pm }$
, consequently bifurcates directly from the conduction state (as was observed for the SHE on the axisymmetric disc by Lloyd & Sandstede (Reference Lloyd and Sandstede2009) and Verschueren et al. (Reference Verschueren, Knobloch and Uecker2021)). The branches of spatially localised states with convection rolls near the equator,
$L_{10}^{C \pm }$
, on the other hand, are disconnected from the branch bifurcating from the
$\ell =10$
mode. We attribute this disconnection to the fact that the spatial modulation of equatorial convectons is opposite to that of the eigenfunctions (i.e. the Legendre polynomials shown in § 4).
As the spherical problem yields four distinct branches of symmetric localised states which would be dynamically equivalent if the system were translational-invariant in
$\theta$
, strong similarities are observed with the bifurcation structure of a planar system in the absence of this symmetry (Houghton & Knobloch Reference Houghton and Knobloch2009; Mercader et al. Reference Mercader, Batiste, Alonso and Knobloch2009). To shed further light on the bifurcation structure, we treat the solutal Rayleigh number,
${\textit{Ra}}_S$
, as a homotopy parameter, i.e. we track how the bifurcation diagram changes with
${\textit{Ra}}_S$
. We describe our findings below, together with a description of the solutions belonging to the various localised state branches computed.
5.1.1. Minus branch
Figure 8 shows how the bifurcation diagram of the
$L_{10}^{A-}$
and
$L_{10}^{C-}$
branches changes with the value of the solutal Rayleigh number. An imperfect bifurcation is present at
${\textit{Ra}}_S \approx 404$
that splits the equatorial-convectons from the anticonvectons near their leftmost saddle node, the point where these branches would otherwise connect spatially periodic states in suitable domains (Houghton & Knobloch Reference Houghton and Knobloch2009). For larger
${\textit{Ra}}_S$
, the
$L_{10}^{C-}$
equatorial-convectons are favoured over the
$L_{10}^{A-}$
anticonvectons whose extent decreases rapidly while the equatorial-convecton branch remains the only one present at large energies. When
${\textit{Ra}}_S$
is decreased below the imperfect bifurcation, the opposite behaviour is observed, although the process is much slower: the equatorial-convectons at
${\textit{Ra}}_S = 250$
still reach larger energies than the anticonvectons at
${\textit{Ra}}_S = 420$
. Changing the value of the solutal Rayleigh number also affects the minimum value of the thermal Rayleigh number at which these localised states can be observed: the leftmost saddle node of the dominant structure is found at
${\textit{Ra}}_T \approx 2896.72$
for
${\textit{Ra}}_S = 250$
while it is located at
${\textit{Ra}}_T \approx 3110$
for
${\textit{Ra}}_S = 420$
.
Bifurcation diagrams of the
$L_{10}^{A-}$
anticonvecton (dotted line) and the
$L_{10}^{C-}$
equatorial-convecton (solid line) branches shown by the kinetic energy
$\mathcal{E}$
as a function of the thermal Rayleigh number
${\textit{Ra}}_T$
for different values of
${\textit{Ra}}_S$
.
By looking at both sides of the imperfect bifurcation identified in figure 8, we were also able to compute the full extent of the
$L_{10}^{A-}$
anticonvectons and of the
$L_{10}^{C-}$
equatorial-convecton spatial development. The
$L_{10}^{C-}$
equatorial-convectons are represented together with their bifurcation diagram at
${\textit{Ra}}_S = 450$
in figure 9. Starting from low energy, the branch first extends to lower thermal Rayleigh numbers where the solution takes the form of an equatorial-convecton displaying downflow at the equator and constituted of four convection rolls, leading to downward flow at both ends of the localised structure. This structure is shown in solution (a) in figure 9, where five adjacent regions are visible where the temperature departure from the conduction state
$\varTheta$
alternates sign. Following this branch to larger
$\mathcal{E}$
, we observe that the
$L_{10}^{C-}$
equatorial-convectons grow in a similar way to snaking localised states with two important exceptions: (i) the branch oscillation in parameter space is of smaller amplitude when the outermost rolls display upflow on the outside of the convective structure and (ii) the snaking branch is slanted to the left and the saddle nodes are, therefore, not aligned. Nonetheless, the solution grows in space by successive nucleation of convection rolls after every oscillation of the branch. In the case of
$L_{10}^{C-}$
at
${\textit{Ra}}_S = 450$
, the slant is such that some saddle nodes are missing thus preventing a one-to-one comparison with standard snaking (see figure 1). To complete the picture of the nucleation process, we manually added points (c) and (d) on the
$L_{10}^{C-}$
branch of figure 9 which would roughly correspond to the location of saddle nodes if the branch were not slanted. The equatorial-convecton becomes domain-filling at point (g) before adjusting its wavelength to the size of the domain and leaving the slanted snaking region to go to large amplitude and large thermal Rayleigh number.
(Middle) Bifurcation diagram showing the
$L_{10}^{C-}$
and the
$L_{10}^{C+}$
equatorial-convecton branches via the kinetic energy
$\mathcal{E}$
as a function of the thermal Rayleigh number
${\textit{Ra}}_T$
. Solutions at the saddle nodes (red circles) and at additional points (blue squares) along
$L_{10}^{C-}$
are represented on the left-hand panel and along
$L_{10}^{C+}$
on the right-hand panel. These representations show the temperature departure from the conduction state profile
$\varTheta$
. The solutal Rayleigh number value is
${\textit{Ra}}_S=450$
.
The
$L_{10}^{A-}$
anticonvecton branch is represented at
${\textit{Ra}}_S = 350$
in figure 10. At its lowest-energy saddle node, the
$L_{10}^{A-}$
anticonvectons take the form of identical localised convection states located at both poles comprised of one roll driving an upward flow at the pole. The branch proceeds in a similar fashion to
$L_{10}^{C-}$
at
${\textit{Ra}}_S = 450$
, describing a similar type of slanted imperfect snaking. The
$L_{10}^{A-}$
anticonvectons also grow similarly: one convection roll is added around each polar structure per branch oscillation until a domain-filling,
$10$
-roll state is reached by saddle node 8. The branch then turns around to larger thermal Rayleigh numbers while convection strengthens. This spatially developed convection state connects with that produced by the
$L_{10}^{C-}$
branch at the imperfect bifurcation at
${\textit{Ra}}_S \approx 404$
.
(Left) Bifurcation diagram showing the
$L_{10}^{A-}$
anticonvecton branch via the kinetic energy
$\mathcal{E}$
as a function of the thermal Rayleigh number
${\textit{Ra}}_T$
. (Right) Solutions at the saddle nodes (red circles in the left-hand panel) and at additional points (blue squares in the left-hand panel) along the branch represented by the temperature departure from the conduction state profile
$\varTheta$
. The solutal Rayleigh number value is
${\textit{Ra}}_S=350$
.
Bifurcation diagram of the
$L_{10}^{A+}$
anticonvecton (dotted line) and the
$L_{10}^{C+}$
equatorial-convecton (solid line) branches shown by the kinetic energy
$\mathcal{E}$
as a function of the thermal Rayleigh number
${\textit{Ra}}_T$
for different values of
${\textit{Ra}}_S$
. The branch depicted in blue is shown for completeness as it is involved in an imperfect bifurcation with
$L_{10}^{A+}$
for
${\textit{Ra}}_S \approxeq 198$
. At large amplitude, this branch carries domain-filling solutions, sometimes with a defect. When it is continued to lower amplitude, these solutions display larger defects and may become spatially localised.
5.1.2. Plus branch
The scenario associated with the
$L_{10}^{A+}$
and
$L_{10}^{C+}$
branches can also be elucidated by using the solutal Rayleigh number as a homotopy parameter. Figure 11 shows how these branches behave for several values of
${\textit{Ra}}_S$
. The connection scenario between
$L_{10}^{A+}$
and
$L_{10}^{C+}$
is not as direct as that between
$L_{10}^{A-}$
and
$L_{10}^{C-}$
. It involves two imperfect bifurcations and a third branch of solutions. At
${\textit{Ra}}_S = 400$
, the
$L_{10}^{C+}$
equatorial-convectons extend to large energy while the
$L_{10}^{A+}$
anticonvectons turn around at a saddle node and remain low-energy states. A first imperfect bifurcation takes place at
${\textit{Ra}}_S \approx 303$
where the upper part of the
$L_{10}^{C+}$
branch disconnects into a branch of domain-filling states, leaving the spatially localised states on a branch that returns back to low energy. As
${\textit{Ra}}_S$
is reduced, the branch of equatorial-convectons recedes while the opposite happens to the branch of anticonvectons. A second imperfect bifurcation takes place at
${\textit{Ra}}_S \approx 198$
, where the anticonvecton branch connects with a branch of domain-filling states with
$12$
rolls and allows the
$L_{10}^{A+}$
anticonvectons to extend to high energy at lower values of the solutal Rayleigh number. The
$L_{10}^{C+}$
equatorial-convectons are best described at larger values of the solutal Rayleigh number. They are represented at
${\textit{Ra}}_S = 450$
in figure 9 and behave in a similar way to the
$L_{10}^{C-}$
equatorial-convectons with the exception that they drive an upward rather than a downward flow at the equator. Conversely, the
$L_{10}^{A+}$
anticonvectons are best described at lower values of the solutal Rayleigh number and are shown in figure 12 for
${\textit{Ra}}_S = 175$
. Although the
$L_{10}^{A+}$
anticonvectons grow in a similar fashion to their
$L_{10}^{A-}$
counterparts, the oscillations of their branch as it is followed towards larger energy are not as pronounced as those of
$L_{10}^{A-}$
. Additionally, we find that the convective regions are not as strongly localised, with small oscillations of the temperature field remaining in the quiescent regions of the localised states of figure 12. We attribute both of these differences to the lower value of the solutal Rayleigh number, which implies a larger disparity between the solutal and thermal buoyancy effects and, therefore, a weaker subcriticality.
(Left) Bifurcation diagram showing the
$L_{10}^{A+}$
anticonvecton branch via the kinetic energy
$\mathcal{E}$
as a function of the thermal Rayleigh number
${\textit{Ra}}_T$
. (Right) Solutions at the saddle nodes (red circles in the left-hand panel) and at additional points (blue squares in the left-hand panel) along the branch represented by the temperature departure from the conduction profile
$\varTheta$
. The solutal Rayleigh number value is
${\textit{Ra}}_S = 175$
.
An important feature of the bifurcation scenario presented here and which is different in planar geometries, is the fact that full branches of equatorial-convectons are found at larger values of the solutal Rayleigh number while full branches of anticonvectons are found at lower values of the same parameter. This selection mechanism may be related to the balance of contributions to the buoyancy force, usually characterised by the magnitude of the buoyancy ratio:
$|N| = {\textit{Ra}}_S/{\textit{Ra}}_T$
(for
${\textit{Ra}}_S \leqslant {\textit{Ra}}_T$
). The further away from
$1$
the value of this number, the less competition there is between thermal and solutal effects, and the weaker the subcriticality is expected to be. Anticonvectons are found for relatively low buoyancy ratio magnitudes: for example the full
$L_{10}^{A-}$
branch was computed for
$|N| \approx 0.10$
in figure 10 while the full
$L_{10}^{A+}$
branch was computed for
$|N| \approx 0.06$
in figure 12. These solutions are characterised by localised convection near the poles. This is where the geometry favours the emergence of convection, as shown by the eigenmodes of the linearised problem (see figure 5). Conversely, full branches of equatorial-convectons of figure 9 are computed for larger buoyancy ratio magnitudes,
$|N| \approx 0.14$
and are not found for
$|N| \lt 0.10$
. We attribute this to the fact that, unlike for the anticonvectons, their existence is not facilitated by geometrical properties and, thus, relies critically on the balance between thermal and solutal effects to sustain convection.
The snaking observed for both equatorial-convectons and anticonvectons is also characterised by the fact that, when the outermost rolls display downward flow (i.e. a negative temperature perturbation) on the outside of the convective structure, the branch covers a wider range of thermal Rayleigh number values than when these rolls display upward flow (i.e. a positive temperature perturbation). We attribute this behaviour to the fact that the spherical shell system, unlike the planar one, does not possess any top-down symmetry. Because of this, the nucleation of a new roll is a different process depending on its sense of rotation. A roll yielding downward flow at the limit of the convective structure is associated with the growth of a negative temperature perturbation. In our figures, these perturbations take the form of blue trapezia with a smaller upper side. Conversely, nucleation of the opposite rolls creates a positive temperature perturbation which looks like red trapezia with a longer upper side.
Bifurcation diagram showing the kinetic energy
$\mathcal{E}$
as a function of the thermal Rayleigh number
${\textit{Ra}}_T$
for the
$L^+_{11}$
(dashed) and
$L^-_{11}$
(solid) pole-convecton branches. Parameters are
${\textit{Ra}}_S=150$
and
$\varGamma = 10.029$
.
5.2. Symmetry-breaking localised states
In addition to the symmetric equatorial-convectons and anticonvectons, system (2.6) admits spatially localised solutions that break the equatorial symmetry
$R_{\textit{AP}}$
. These states may originate from pitchfork bifurcations from the conduction state which are generated by marginal eigenmodes of the linear stability analysis with odd parity. Unlike
$R_{\textit{AP}}$
-symmetric solutions, which we sought on the symmetric (i.e. reduced) subspace, symmetry-breaking localised states necessitate the whole search space to be computed.
Bifurcation diagram of the snaking region (middle) of the
$L^-_{11}$
pole-convectons (solid line, right) and of the
$L^+_{11}$
pole-convectons (dashed line, left) showing the solution kinetic energy
$\mathcal{E}$
as a function of the thermal Rayleigh number
${\textit{Ra}}_T$
. Solutions are taken at successive saddle nodes or roll nucleation points and represented in the left- and right-hand panels via the temperature departure from conduction
$\varTheta$
. Parameter values are
${\textit{Ra}}_S = 150$
and
$\varGamma = 10.029$
.
The complete bifurcation diagram of the
$R_{\textit{AP}}$
-breaking localised states produced by the pitchfork bifurcation at
$\ell = 11$
is shown in figure 13 for
${\textit{Ra}}_S = 150$
. The
$\ell = 11$
pitchfork at
${\textit{Ra}}_T \approx 4526$
produces subcritical branches of domain-filling counter-rotating rolls of convection that reach an imperfect bifurcation point at
${\textit{Ra}}_T \approx 4519$
. Following the branch whose energy grows the slowest past the imperfect bifurcation, spatial modulation develops to favour the pole where the structure displays downflow. The resulting solutions are referred to as
$L_{11}^-$
pole-convectons and the corresponding branch extends subcritically to produce slanted snaking between
${\textit{Ra}}_T = 2750$
and
${\textit{Ra}}_T = 3000$
, described further in figure 14. This time, the snaking produced by the
$L_{11}^-$
branch appears more regular. With the exception of the bottom saddle nodes, which may not always align with the others (Chapman & Kozyreff Reference Chapman and Kozyreff2009), the saddle nodes appear to fall within a snaking region that is slightly slanted to the left. As the branch is followed upward along the snaking, typical spatial growth is observed: a counter-rotating roll of convection is nucleated on the outside of the convective structure during each snaking period until the domain is full. Just like for the
$R_{\textit{AP}}$
-symmetric equatorial-convectons and anticonvectons, we observe that the nucleation of a roll driving a downward flow at the edge of the convective structure leads to a wider-amplitude oscillation of the branch in parameter space than the nucleation of a roll associated with an upward flow at the edge of the convective structure. At what would be its termination point, where the pattern fills the domain, the branch undergoes another imperfect bifurcation induced by the domain geometry at
${\textit{Ra}}_T \approx 2674$
: it turns to higher-energy states, passing through saddle node (m) and taking on the role of the domain-filling state branch. The other side of this imperfect bifurcation contains localised states of the opposite phase, i.e. pole-convectons with an upflow at the pole, which we refer to as
$L_{11}^+$
. Following
$L_{11}^+$
on the side of this imperfect bifurcation where its solution energy increases, one obtains a complementary snaking to the one described by
$L_{11}^-$
, with a similar saddle-node alignment and spatial growth mechanism. The branch can be followed all the way down the snaking and towards larger thermal Rayleigh numbers to reveal that it is also involved in the imperfect bifurcation at
${\textit{Ra}}_T \approx 4519$
, where it abruptly turns around towards larger-amplitude states. The
$L_{11}^+$
branch behaviour past that point is complex and beyond the scope of this paper. Following
$L_{11}^+$
in the other direction past the imperfect bifurcation at
${\textit{Ra}}_T \approx 2674$
reveals that the solution undergoes spatial modulation to preserve convection at both poles but with opposite flow direction. As a result, the
$L_{11}^{+}$
pole-convecton turns into symmetry-breaking anticonvectons, as shown in figure 15, and the branch produces a similar bifurcation diagram to that of the symmetric anticonvectons. After this spatially localised interlude, the branch travels to low-amplitude states where it reconnects to the
$\ell = 11$
bifurcation.
(Left) Bifurcation diagram showing the
$L^+_{11}$
branch of figure 14 (dashed line) using the kinetic energy
$\mathcal{E}$
as a function of the thermal Rayleigh number
${\textit{Ra}}_T$
. (Right) Solutions are taken at the points labelled in the left-hand panel and represented by the temperature departure from conduction
$\varTheta$
. Parameter values are
${\textit{Ra}}_S = 150$
and
$\varGamma = 10.029$
. Solutions on the other part of the branch are shown in figure 14.
We have regretfully not been able to compute the bifurcation diagram of symmetry-breaking equatorial convectons (the symmetry-breaking analogue to the anticonvectons shown in figure 15). We are convinced that these states exist despite our best attempts failing to locate them. We have tried using combinations of time stepping, homotopical continuation and aggressive manipulations of the solution space, among other techniques. The difficulties that prevented us from computing these solutions are: (i) their asymmetric nature, (ii) the fact that their solution branch is not connected to the conduction state, (iii) that the structure of their convection rolls is not as simple as in the planar system and (iv) that the pole-convectons that share parity with the symmetry-breaking equatorial convectons do not possess an equatorially centred amplitude envelope in their domain-filling form (the amplitude dip of states 12 and (m) of figure 14 is located closer to the left of the domain/south pole). We do not, however, expect these solutions to produce substantially different behaviour from that of the other localised solutions presented in this paper.
6. Discussion
Motivated by the recent progress made on spatially localised pattern formation in planar fluid systems and, in particular, in doubly diffusive convection, we focused here on an alternative geometry relevant to astrophysical applications: the spherical shell. In order to connect our findings in the spherical shell with established literature on spatial localisation, we chose parameter values used in studies of doubly diffusive convection in planar geometries. As the present work represents the first analysis of spatially localised states in a spherical shell, we restricted ourselves to axisymmetric solutions.
The most distinguishing feature of doubly diffusive convection in an axisymmetric spherical shell is the fact that convection rolls are not straight as in two-dimensional planar studies but are curved in order to wrap around the inner sphere. The solutions bifurcating from the conduction state are therefore different from those of the planar problem, where they consist of homogeneous arrays of identical convection rolls typically referred to as spatially periodic. In the spherical shell, they consist of spatially modulated solutions which become spatially localised as their solution branches are continued away from the bifurcation point. Such solutions only arise via secondary bifurcations in unbounded planar doubly diffusive convection (Beaume et al. Reference Beaume, Bergeon and Knobloch2011). The absence of a primary bifurcation to spatially periodic states is a consequence of the fact that the axisymmetric spherical shell system is not invariant under spatial translations in the latitudinal direction. The fact that spatially periodic states are not created by the primary bifurcation has already been observed in systems with non-periodic boundary conditions and, therefore, in the absence of translation invariance, like in the SHE (Houghton & Knobloch Reference Houghton and Knobloch2009), in binary fluid convection (Mercader et al. Reference Mercader, Batiste, Alonso and Knobloch2011) and in natural doubly diffusive convection (Beaume et al. Reference Beaume, Knobloch and Bergeon2013b ).
By performing numerical continuation, we unravelled the existence of a number of spatially localised state families, which we differentiate by their symmetry and by the location of the localised pattern, i.e. either/both pole(s) or the equator. In particular, we found: equatorial-convectons, where convection is centred around the equator; pole-convectons, where it is centred around one pole; and anticonvectons, where each pole is the centre of a convection zone. All these solutions have analogues in the planar doubly diffusive convection where they exist within similar values of the parameters (Beaume et al. Reference Beaume, Bergeon and Knobloch2011; Mercader et al. Reference Mercader, Batiste, Alonso and Knobloch2011). This is not the case in the spherical shell, where spatially localised convection is preferred at the poles for relatively small values of the solutal Rayleigh number and at the equator for larger values of the same number. This suggests that the balance between the thermal and solutal components of the buoyancy force plays a selection role in the type of localised pattern that develops in the system.
The best theoretical framework for the analysis of doubly diffusive pole convectons is the axisymmetric SHE posed on a disc (Lloyd & Sandstede Reference Lloyd and Sandstede2009; McCalla & Sandstede Reference McCalla and Sandstede2010). The disc on which the SHE is posed corresponds to a hemishell centred on either pole in such a way that pole-convectons (and anticonvectons) correspond to spot patterns and equatorial-convectons correspond to target patterns in the SHE. The bifurcation diagram of axisymmetric states computed by Verschueren et al. (Reference Verschueren, Knobloch and Uecker2021) for the SHE on a disc reveals a number of similarities to our bifurcation diagrams. In particular, spot patterns bifurcate directly from the trivial solution while target patterns are disconnected from it. We observed the same in doubly diffusive convection, where pole-convectons and anticonvectons are found to connect with the conduction state while equatorial-convectons are not always connected with it. The bifurcation diagram of the localised states in the axisymmetric SHE displays snaking that is not straight but tends to have a more pronounced leftward slant when the edge of the localised pattern approaches the centre of the domain, where the roll curvature becomes important. While our domain size is modest by comparison with that used in standard SHE studies, these results suggests that the leftward-slanted snaking that we observe is a direct consequence of the geometry.
We anticipate that the bifurcation scenario of the fully three-dimensional system will support an even larger variety of spatially localised states. In particular, new types of convectons may arise that break the longitudinal invariance. Predictions about their formation and bifurcation diagram can be made in light of the non-axisymmetric solutions found in the SHE by Verschueren et al. (Reference Verschueren, Knobloch and Uecker2021). However, the expected complexity in the fluid system is best approached with knowledge of the axisymmetric bifurcation scenario. This is what we have aimed to achieve here.
Besides providing here the first study of doubly diffusive convectons in a spherical shell, we believe that the scenario presented and the links made with the axisymmetric SHE possess a wide applicability to spatially localised pattern formation in spherical geometries, in fluid dynamics and beyond. Further work should include the computation of the symmetry-breaking equatorial-convectons that we have not been able to compute here. We will also endeavour to extend our findings to larger values of the domain aspect ratio
$\varGamma$
, which will allow a better connection to the planar geometry literature and to Swift–Hohenberg theory.
Acknowledgements
The authors wish to thank J. Mestel for his careful reading and recommendations for improving the initial draft of the manuscript. The authors would also like to thank R. Hollerbach for the involvement in the early stages of the work. P.M.M. wishes to acknowledge the support of the Engineering and Physical Sciences Research Council (grant number EP/V033883/1) as part of the [D
$^{*}$
]stratify project.
Declaration of interests
The authors report no conflict of interest.
Data availability
The numerical code, supplementary material detailing its implementation and data used to generate the figures are available on Zenodo at https://doi.org/10.5281/zenodo.15462897.
Open access
