1. Introduction
The effect of an off-centre gravity location on the flow structures in turbulent thermal convection has not been studied previously. Here, we use spherical Rayleigh–Bénard (RB) convection as an idealized model system to study the effect of an off-centre gravity location on thermal convection in a spherical system. The system consists of a fluid layer enclosed between two spherical shells, heated from the inner sphere and cooled from the outer one. Convection in spherical shells differs from the classical RB convection (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; Lohse & Xia Reference Lohse and Xia2010; Chillà & Schumacher Reference Chillà and Schumacher2012) owing to the curvature of the boundaries, non-uniform gravity, and the geometrical asymmetry between the boundary layer at the inner and outer spheres (Busse Reference Busse1970; Spiegel Reference Spiegel1971; Jarvis, Glatzmaierand & Vangelov Reference Jarvis, Glatzmaierand and Vangelov1995; Tilgner Reference Tilgner1996; Shahnas et al. Reference Shahnas, Lowman, Jarvis and Bunge2008; Deschamps, Tackley & Nakagawa Reference Deschamps, Tackley and Nakagawa2010; O'Farrell, Lowman & Bunge Reference O'Farrell, Lowman and Bunge2013; Gastine, Wicht & Aurnou Reference Gastine, Wicht and Aurnou2015).
Most studies on spherical RB convection focus on the case in which the gravity centre coincides with the geometric one. Busse (Reference Busse1975) used perturbation analysis to show the qualitative difference between convective flow patterns of odd and even spherical harmonic order just above the onset of convection. Subsequently, in an infinite Prandtl number (
$Pr$) fluid with application to mantle convection, Bercovici, Schubert & Glatzmaier (Reference Bercovici, Schubert and Glatzmaier1989) showed that these convective patterns persist for Rayleigh numbers (
$Ra$) up to 100 times larger than the critical value. For 
$Ra>10^5$, the axisymmetric convective patterns break down when the flow starts to show time-dependent behaviour (Iwase & Honda Reference Iwase and Honda1997; Deschamps et al. Reference Deschamps, Tackley and Nakagawa2010; Futterer et al. Reference Futterer, Krebs, Plesa, Zaussinger, Hollerbach, Breuer and Egbers2013; Yanagisawa, Kameyama & Ogawa Reference Yanagisawa, Kameyama and Ogawa2016). In spherical RB convection, the influence of curvature is important. For constant 
$Ra$, the length scale of the preferentially convective patterns decreases with increasing radius ratio (Deschamps et al. Reference Deschamps, Tackley and Nakagawa2010). Yanagisawa & Yamagishi (Reference Yanagisawa and Yamagishi2005) find that the convective rolls become smaller with increasing 
$Ra$ when the radius ratio is kept constant. Gastine et al. (Reference Gastine, Wicht and Aurnou2015) show that for 
$Pr=1$, sheet-like thermal plumes are formed near the inner and outer spheres, and these plumes undergo morphological changes into mushroom-like plumes when they eject from the boundary layers to the bulk. These plume dynamics are similar to the situation for classic RB convection (Puthenveettil & Arakeri Reference Puthenveettil and Arakeri2005; Shishkina & Wagner Reference Shishkina and Wagner2008; Zhou & Xia Reference Zhou and Xia2010; Chillà & Schumacher Reference Chillà and Schumacher2012). However, due to the asymmetries between the hot and cold surfaces in spherical shells, the mushroom-like plumes emitted from the outer sphere are thicker than those emitted from the inner sphere (Gastine et al. Reference Gastine, Wicht and Aurnou2015). Nevertheless, studies have shown that there are similarities between spherical and classical planar RB convection. For example, Gastine et al. (Reference Gastine, Wicht and Aurnou2015) showed that the effective heat transport (characterized by the Nusselt number 
$Nu$) scaling exponent 
$\alpha$ in 
$Nu \sim Ra^\alpha$ in spherical RB convection has a similar 
$Ra$ number dependence as in planar RB convection, where for 
$Ra\lesssim 10^{11}$, typically 
$0.28 \lesssim \alpha \lesssim 0.31$ (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009).
However, the misalignment between the gravity centre and the geometric one in the spherical RB convection set-up has been overlooked to date. In this study, we use three-dimensional direct numerical simulations to investigate the effect of the gravity centre offset on the flow structures and heat transfer in spherical RB convection. The paper is organized as follows. In § 2, we explain the numerical method and discuss the considered parameter regime. In §§ 3.1–3.3, we discuss the large-scale flow pattern induced by the off-centre gravity. The effect of the large-scale structure on heat transfer is shown in § 3.4. In § 3.5, we study the effect of the gravity profile and 
$Ra$ on flow structures and heat transfer. The paper ends with a summary of our findings and conclusions in § 4.
2. Numerical method and parameters
2.1. Set-up of convection in spherical shells
The spherical RB geometry is illustrated schematically in figure 1. Fluid fills a spherical shell between the inner sphere of radius 
$r_i$ and the outer sphere of radius 
$r_o$. The radius ratio between the inner and outer spheres is given by
\begin{equation} \eta=r_i/r_o. \end{equation}
In this study, we consider a fixed aspect ratio 
$\eta =0.3$. The surface temperatures of the inner and outer spheres are kept constant at 
$T_i$ and 
$T_o$, respectively, with 
$T_i > T_o$. No-slip boundary conditions are imposed on both boundaries. Point 
$G$ indicates the gravity centre, which is offset from the geometric centre 
$O$. The gravity shift is defined as
\begin{equation} \varepsilon=|\boldsymbol{OG}|/r_i, \end{equation}
where 
$|\boldsymbol {OG}|$ indicates the displacement of the gravity compared to the geometrical centre; see figure 1. In the absence of rotation (and magnetic field), for a sphere, there are no north and south directions, therefore any diametrically opposite symmetry is equivalent to another. Nevertheless, since the spherical coordinates are naturally clustered around the polar axis, it is beneficial to consider the jet aligned along this axis, referred to as north–south, from a computation point of view. For convenience, we always shift the gravity centre 
$G$ to the south.
Figure 1. Schematic of spherical RB convection in between two concentric spheres with radius ratio 
$\eta =r_i/r_o$ between the inner and outer sphere, and gap size 
$d=r_o-r_i$. A no-slip boundary condition with constant temperature is used on the inner (hot) and outer (cold) sphere. In spherical coordinates the longitudinal, co-latitudinal and radial directions are represented by 
$\hat {\theta }$, 
$\hat {\varphi }$ and 
$\hat {r}$, respectively. The gravity centre 
$G$ is offset from the geometric centre 
$O$, and 
$P$ is an arbitrary fluid point in the spherical shell.
The dynamics of spherical RB convection is controlled by the Rayleigh and Prandtl numbers
\begin{equation} Ra=\frac{\beta g_o\,\Delta T\,d^3}{\kappa \nu}, \quad Pr=\frac{\nu}{\kappa}, \end{equation}
where 
$\beta$ is the thermal expansion coefficient, 
$g_o$ is the surface-averaged gravity at the outer sphere, 
$\nu$ is the kinematic viscosity, and 
$\kappa$ is the thermal diffusivity. We normalize the fields and distances by the length scale 
$d=r_o-r_i$, the temperature difference 
$\Delta T$ between the inner and outer spheres, and the free-fall velocity 
$U=\sqrt {\beta g_o\,\Delta T\,d}$.
2.2. Underlying dynamical equation and numerical method
We solve the Navier–Stokes equations under the Boussinesq approximation in spherical coordinates, which in dimensionless form read
\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}={-}\boldsymbol{\nabla} p+\sqrt{\frac{P r}{R a}}\,\nabla^{2} \boldsymbol{u}+ \boldsymbol{C}_{g} g T, \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}
\begin{gather} \frac{\partial T}{\partial t}+\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} T=\frac{1}{\sqrt{{Ra\,Pr}}}\,\nabla^{2} T, \end{gather}
where 
$\boldsymbol {u}$, 
$p$, 
$T$ and 
$\boldsymbol {C}_{g}g$ denote the fluid velocity, pressure, temperature and gravitational acceleration acting in the three directions. The coefficients 
$C_{g,i}$ (where 
$i=1,2,3$ correspond to the 
$\hat {\theta }$, 
$\hat {r}$ and 
$\hat {\varphi }$ directions, respectively) denote the decomposition of 
$g$, which is detailed in Appendix A. Since we always shift the gravity centre 
$G$ to the south, the three components of the buoyancy coefficient can be written as
\begin{equation} \boldsymbol{C}_{g}g=(0,r+r_i \varepsilon \cos{\varphi},-\varepsilon{r_i}\sin\varphi)\,r_o^{{-}n_g}\,L_{GP}^{n_g-1}, \end{equation}
in which 
$L_{GP}=|\boldsymbol {GP}|$ indicates the distance between a fluid point and the gravity centre; see figure 1. We consider the gravity profile
\begin{equation} g=(L_{GP}/r_o)^{n_g}, \end{equation}
and we use 
$n_g=-2$, 
$-1$, 
$0$ and 
$1$. The value of 
$n_g$ can be inferred by physical considerations on the mass ratio between the nucleus (
$r< r_i$) and the spherical shell (
$r_i< r< r_o$). The gravity is constant (
$n_g=0$) when the mass of the spherical shell is negligible. This assumption is typically used to model the Earth's mantle (Bercovici et al. Reference Bercovici, Schubert and Glatzmaier1989). If the mass is condensed centrally, then one obtains 
$n_g=-2$. Gastine et al. (Reference Gastine, Wicht and Aurnou2015) showed that only for this gravity distribution function is there an analytical relation between the viscous dissipation rate (
$\epsilon _U$) and 
$Nu$, i.e. 
$\epsilon _U = ({3}/({1+\eta +\eta ^2})) ({1}/{Pr})(Nu-1)$. When the density is constant within the fluid layer, and no mass is contained in the inner sphere, gravity is directly proportional to the radial coordinate 
$r$ (
$n_g=1$) (Tilgner Reference Tilgner1996).
The governing equations (2.4) and (2.5) are discretized by a staggered central second-order finite-difference scheme in spherical coordinates. The numerical scheme is based on the method by Verzicco & Orlandi (Reference Verzicco and Orlandi1996), which has been extended recently to spherical coordinates by Santelli, Orlandi & Verzicco (Reference Santelli, Orlandi and Verzicco2021). The advantage of this scheme is that it allows arbitrary non-uniform grids in the radial and co-latitudinal directions. Here, it is worthwhile to mention that singularities at the poles are prevented by introducing a new set of quantities 
$(u_\theta,u_r r^2, u_\varphi \sin \varphi )$, which results in trivial boundary conditions 
$u_\varphi \sin \varphi =0$ at the north pole (
$\varphi =0$) and south pole (
$\varphi ={\rm \pi}$) in the co-latitudinal direction. The code has been validated carefully by Santelli et al. (Reference Santelli, Orlandi and Verzicco2021), and we refer to that work for details on the method. For additional validation results relevant to the flows considered here, we refer the reader to Appendix B.
2.3. Choice of parameters
To study the effect of the gravity centre location on the flow dynamics in spherical RB convection, we considered different gravity distributions (
$n_g \in \{-2,-1,0,1\}$, see (2.7)) and gravity centre offset (
$\varepsilon$ from 
$0$ to 
$0.8$, see (2.2)) for 
$Ra=10^7$ and 
$10^8$. All simulations in this study are for 
$Pr=1$. To ensure that the flow is resolved fully, we place a sufficient number of computational grid points in the bulk and the boundary layers (Stevens, Verzicco & Lohse Reference Stevens, Verzicco and Lohse2010). As we will use spectral analysis to study the flow structures, we use a uniform grid in the longitudinal and co-latitudinal directions. In the radial direction, the grid cells are clustered towards the inner and outer spheres to ensure that the boundary layers are resolved adequately (Shishkina et al. Reference Shishkina, Stevens, Grossmann and Lohse2010).
We calculate the Nusselt number 
$Nu$ from the normalized averaged temperature gradients at the inner and outer spheres as follows:
\begin{equation} Nu_i = \left. -\eta\,\frac{{\rm d} \overline{\langle T \rangle_{\theta,\varphi}}}{{\rm d}r} \right| _{r_i},\quad Nu_o= \left. -\frac{1}{\eta}\,\frac{{\rm d} \overline{\langle T \rangle_{\theta,\varphi}}}{{\rm d}r} \right| _{r_o}. \end{equation}
Here, 
$\langle \cdots \rangle _{\theta,\varphi }$ represents the average over a spherical surface with fixed 
$r$, and 
$\overline {\cdots }$ indicates a time average. Also, 
$Nu_i$ and 
$Nu_o$ as functions of the co-latitudinal direction are defined as
\begin{equation} Nu_i(\varphi) = \left. -\eta\,\frac{{\rm d} \overline{\langle T \rangle_\theta}}{{\rm d}r} \right| _{r_i},\quad Nu_o(\varphi)= \left. -\frac{1}{\eta}\,\frac{{\rm d} \overline{\langle T \rangle_\theta}}{{\rm d}r} \right| _{r_o}, \end{equation}
where 
$\langle \cdots \rangle _{\theta }$ represents the average over the azimuthal direction. The difference in 
$Nu$ values obtained at the inner and outer spheres is always less than 
$0.2\,\%$. Also, we verify that 
$Nu$ calculated at the spheres is the same as the value obtained from the volume-averaged Nusselt number 
$Nu_v$, which is calculated from the surface-averaged Nusselt number 
$Nu_h$ averaged over concentric spherical surfaces with fixed 
$r$, with
\begin{equation} Nu_h(r)=\frac{\sqrt{Ra\,Pr}\,\langle u_rT \rangle_{\theta,\varphi}- \partial_r \langle T \rangle_{\theta,\varphi}}{- \partial_r T_c}, \end{equation}
where 
$T_c(r)=\eta /[(1-\eta )^2r]-\eta /(1-\eta )$ is the purely conductive temperature profile for constant temperature boundary conditions. Subsequently, we perform radial integration and time averaging to obtain 
$Nu_v$ as
\begin{equation} Nu_v= \overline{\dfrac{1}{d}\int_{r_i}^{r_o} Nu_h(r) \,{\rm d}r.} \end{equation}The volume-averaged root-mean-square (r.m.s.) Reynolds number is given by
\begin{equation} Re_{rms} = \sqrt{Ra/Pr}\,\overline{\sqrt{\langle u_i u_i \rangle}},\quad i=1,2,3. \end{equation}
Here, 
$\langle {\cdots }\rangle$ denotes the average over the whole volume of the spherical shell. The details of the simulations considered in this study are summarized in table 1.
Table 1. Details of the simulations. The columns from left to right indicate: 
$Ra$, the gravity profile exponent 
$n_g$ (see (2.7)); the number of grid points in the longitudinal, radial and co-latitudinal directions, 
$N_{\theta } \times N_r \times N_\phi$; the shift of the gravity centre with respect to the geometrical centre 
$\varepsilon$ (see (2.1)); the average heat transfer 
$Nu$ across the inner 
$Nu_i$ and outer 
$Nu_o$ spheres (see (2.8a,b)); 
$Nu_v$ (see (2.11)); and 
$Re_{rms}$ (see (2.12)).
3. Results
Before we discuss the turbulent flows in §§ 3.2–3.5, we first focus on the large-scale circulation formation induced by off-centre gravity in a steady axisymmetric flow obtained at 
$Ra=10^3$. As indicated in (2.6), the radial and co-latitudinal buoyancy coefficients depend on the gravity centre offset 
$\varepsilon$ and the gravity distribution 
$n_g$. To disentangle these effects, for turbulent flows we first examine the effect of the gravity centre offset for constant gravity (
$n_g=0$) at 
$Ra=10^8$ in §§ 3.2–3.4. The effect of 
$n_g$ and 
$Ra$ is studied in § 3.5.
3.1. Dynamics of steady flow
When the gravity centre offsets with the geometric one in spherical RB convection, the crossing between temperature gradient (
$\boldsymbol {\nabla } T$) and gravitational acceleration (
$\boldsymbol {g}$) results in the baroclinic torque, which plays an important role in the generation of vorticity by density variations. Figure 2(a) shows gravitational potential isolines to intersect the isotherms in an axisymmetric flow. The curl of the buoyancy term (
$\boldsymbol {\nabla } \times \boldsymbol {g} T$) denotes the Boussinesq baroclinic torque under the Boussinesq approximation (Schladow, Patterson & Street Reference Schladow, Patterson and Street1989; Heifetz, Maor & Guha Reference Heifetz, Maor and Guha2020). As shown in figure 2(b), the baroclinic torque is strongest in the vicinity of the inner sphere where the angle between 
$\boldsymbol {\nabla } T$ and 
$\boldsymbol {g}$ is the largest. The baroclinic torque drives a northward-directed flow along the inner sphere. This results in a large-scale circulation pattern as shown in figure 2(c). In steady flows, large-scale circulations induced by baroclinic torque are also reported for side-heated cavity (Jiang, Sun & Calzavarini Reference Jiang, Sun and Calzavarini2019), stratified flow in a semicylindrical cavity (Kumar & Pothérat Reference Kumar and Pothérat2020), and natural convection in spherical annuli (Scurtu, Futterer & Egbers Reference Scurtu, Futterer and Egbers2008).
Figure 2. Baroclinic imbalance induces a regular large-scale circulation due to the offset gravity centre towards the south (
$\varepsilon =0.2$, see (2.2)) for a steady flow with 
$Ra=10^3$. (a) The misalignment of the gravitational potentials (blue dashed circles) and isotherms (red circles) in a meridional cut. The direction is indicated by the blue arrows. (b) The corresponding Boussinesq baroclinic torque 
$\boldsymbol {\nabla } \times \boldsymbol {g} T$. Two dashed arrows correspond to the direction of the baroclinic torque. (c) The vector field (coloured by the velocity magnitude) and streamlines (grey curves) of the velocity in the meridian plane.
For different 
$n_g$, figure 3(a) shows the radial (
$C_{g,r}g \hat {r}$) and co-latitudinal (
$C_{g,\varphi }g \hat {\varphi }$) buoyancy coefficients for 
$\varepsilon =0.8$ on the mid-plane of the spherical shell. It shows that 
$g$ increases monotonically from the north to the south pole when 
$n_g<0$. For 
$n_g>0$, the opposite trend is observed, while for 
$n_g=0$, the gravity is constant as function of the co-latitude location. We keep 
$Ra=10^3$ and vary 
$n_g$ from 
$-2$ to 
$1$. The baroclinic torque (
$\boldsymbol {\nabla } \times \boldsymbol {g} T$) for 
$\varepsilon =0.8$ is shown in figure 3(b): 
$\boldsymbol {\nabla } \times \boldsymbol {g} T$ concentrates in the south polar region when 
$n_g=-2$, whereas the strongest 
$\boldsymbol {\nabla } \times \boldsymbol {g}T$ gradually moves to the equator with increasing 
$n_g$. The corresponding large-scale circulation is shown in figure 3(c); the velocity field strengthens close to the inner sphere and its highest value moves from the south to the equator, which is the same trend as for the baroclinic torque. In addition, the magnitude of the velocity field decreases with increasing 
$n_g$ due to the decreasing mean buoyancy power (
$\overline {\langle g u_r T \rangle }$, see Gastine et al. Reference Gastine, Wicht and Aurnou2015) averaged over the spherical shell.
Figure 3. (a) The radial (
$C_{g,r}g \hat {r}$, red arrow) and co-latitudinal (
$C_{g,\varphi }g \hat {\varphi }$, blue arrow) decomposition of buoyancy 
$g\,\boldsymbol {GP}$ (black arrow) in (2.6) for 
$\varepsilon =0.8$. Different panels indicate the buoyancy components on the mid-plane of the spherical shell for 
$n_g=-2, -1, 0, 1$; see (2.7). The lengths of the arrows are scaled by the co-latitudinal averaged gravity in the mid-plane. (b,c) Baroclinic imbalance induces regular large-scale circulation due to the offset gravity centre 
$\varepsilon =0.8$ for a steady flow with 
$Ra=10^3$ and different gravity profiles. (b) The Boussinesq baroclinic torque 
$\boldsymbol {\nabla } \times \boldsymbol {g} T$ in a meridional cut. (c) The vector field (coloured by the velocity magnitude) and streamlines (grey curves) of the velocity in the meridian planes.
3.2. Generation of convective jet and large-scale circulation
Figure 4 shows the temperature field on the meridional cut in figures 4(a–d), the equatorial cut in figures 4(e–h, and the surfaces corresponding to the locations of the outer and inner thermal boundary layer heights in figures 4(i–p). Figures 4(b–d) show that for 
$\varepsilon \ge 0.05$, the jet ejects hot plumes from the inner sphere towards the outer sphere. When 
$\varepsilon$ is larger, e.g. 
$\varepsilon \ge 0.4$, thermal plumes preferably ascend northwards along the inner sphere, as shown in figures 4(o,p). This happens because the plumes that form along the southern part of the inner sphere move upwards, following the large-scale flow pattern induced by the baroclinic torque; see figure 2(c). Consequently, the plumes combine at the northern part of the inner sphere to form the northbound plume towards the outer sphere (figures 4c,d). This creates a convective jet impinging on the north pole of the outer sphere, after which thermal plumes descend southwards along the outer sphere, as shown in figures 4(k,l). As a result, a limited number of plumes are ejected radially from the hot sphere to the cold sphere in the equatorial region, as shown in figures 4(g,h). In addition, we visualize the advection of temperature fluctuation by baroclinic flow 
$\boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {\nabla } T'$ as shown in figures 4(q–t). We find that 
$\boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {\nabla } T'$ is intense where 
$T$ is high, which indicates that baroclinic flow works on the temperature perturbation to promote the generation of thermal plumes on the northern side of the sphere. This confirms that the plumes combine at the northern part of the inner sphere to form the northbound plume towards the outer sphere.
Figure 4. (a–p) Snapshots of the temperature fields for 
$Ra=10^8$, 
$n_g=0$ and 
$\varepsilon =0, 0.05, 0.4, 0.8$. (a–d) Meridional cut and temperature isosurface for 
$T=0.2$. (e–h) Temperature field in the equatorial plane, and (i–l,m–p) temperature field at the outer (inner) thermal boundary layer heights along the inner and outer spheres, displayed in hammer projections. (q–t) Advection of temperature fluctuation by the baroclinic flow 
$\boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {\nabla } T'$ at the outer thermal boundary layer height.
We compare the unsteady flow fields for 
$\varepsilon =0.01$ and 
$0.05$, and different 
$Ra$, in figures 5(a) and 5(b), respectively. The 
$Ra=10^3$ cases are shown as reference in the leftmost panel. For the very small 
$\varepsilon =0.01$ in figure 5(a), the upward convective jet is the strongest for 
$Ra=10^7$. However, the convective jet and meridional circulation are weakened significantly for large 
$Ra=10^8$. For 
$\varepsilon =0.05$, figure 5(b) shows that the convective jet is the strongest for 
$Ra=10^6$. Again, the convective jet weakens with increasing 
$Ra$. In figure 5(c), we compare different 
$\varepsilon$ for 
$Ra=10^8$. As expected, the convective jet is enhanced significantly with increasing 
$\varepsilon$. However, the shape of the large-scale circulation is almost unchanged for 
$\varepsilon \le 0.4$. When 
$\varepsilon$ increases to 0.8, the large-scale circulation becomes a crescent shaped circulation, which results in its centre moving northwards.
Figure 5. The vector field (coloured by the velocity magnitude) and streamlines (grey curves) of the longitudinal and time-averaged velocity in the meridian plane for: (a) 
$\varepsilon =0.01$, and (b) 
$\varepsilon =0.05$, with increasing 
$Ra=10^3$ (steady flow), 
$Ra=10^6$, 
$Ra=10^7$ and 
$Ra=10^8$; (c) 
$Ra=10^8$ and 
$n_g=0$ with increasing 
$\varepsilon =0.1, 0.2, 0.4, 0.8$.
3.3. Modal analysis to identify large-scale structures
When the gravity centre is moved towards the south pole (increasing 
$\varepsilon$), the total turbulent kinetic energy (
$\mathrm {TKE}$), which is defined as
\begin{equation} \mathrm{TKE} = \frac{1}{2}\,\frac{Ra}{Pr}\,\overline{\langle u_i u_i \rangle},\quad i=1,2,3, \end{equation}
first increases slightly (
$\varepsilon \leq 0.4$) and then decreases considerably (
$\varepsilon = 0.8$) (see figure 6a). To study the effect of the off-centre gravity on the distribution of 
$\mathrm {TKE}$, we perform a modal analysis of the flow. In spherical coordinates, a surface spherical harmonic function is the expansion of plane waves, where the order 
$m$ and degree 
$\ell$ represent the wavenumbers along the longitudinal and co-latitudinal directions, respectively. Here, we define the power spectrum on a spherical surface as
\begin{equation} C_{\ell, m}=a^2_{\ell, m}+b^2_{\ell, m}, \end{equation}
where 
$a_{\ell, m}$ and 
$b_{\ell, m}$ are the spherical harmonics expansion coefficients, which are detailed in Appendix C. We define the power spectrum in the longitudinal direction as
\begin{equation} C_{m}=\frac{1}{n-m+1} \sum_{\ell=m}^{n} C_{\ell, m},\quad n=\min(N_\varphi-1,(N_\theta+1)/2), \end{equation}
where 
$N_\varphi$ and 
$N_\theta$ represent the numbers of grid points in the co-latitudinal and longitudinal directions, respectively. Radial integration and temporal averaging are used to obtain 
$\overline {\langle \cdots \rangle }_r$. The large-scale circulation is represented by 
$(m=0,\ell >0)$ modes as demonstrated in Appendix C.
Figure 6. (a) Volume-averaged 
$\mathrm {TKE}(\varepsilon )/\mathrm {TKE}(\varepsilon =0)$ as a function of 
$\varepsilon$ for 
$Ra=10^8$ and 
$n_g=0$. The error bars indicate the corresponding standard deviations. (b) Plot of 
$\langle C_{m=0} \rangle$ of 
$\mathrm {TKE}$ for 
$\ell > 0$, normalized by 
$\langle C_{m=0} \rangle$ for 
$\varepsilon =0$. The inset shows 
$\langle C_{\ell,m=0} \rangle$, normalized by 
$\langle C_{\ell,m=0} \rangle$ for 
$\varepsilon =0$, as a function of the spherical harmonic order 
$\ell$, for different 
$\varepsilon$.
Figure 6(a) shows that the average 
$\mathrm {TKE}$ is almost unchanged for 
$\varepsilon \leq 0.4$, but the distribution among the different modes is changed. As we have seen in figure 5, a large-scale circulation is formed when 
$\varepsilon \geq 0.05$ for 
$Ra=10^8$. Figure 6(b) shows that the axisymmetric (
$m=0$) mode, which indicates the large-scale mode, increases strongly with 
$\varepsilon$. For 
$\varepsilon = 0.8$, both the 
$m=0$ mode and the total 
$\mathrm {TKE}$ decrease sharply, while it is still enhanced compared with 
$\varepsilon =0$. Interestingly, only a small gravity centre offset (
$\varepsilon \geq 0.05$) is required to form the convective jet and the corresponding meridional circulation for 
$Ra=10^8$. The inset of figure 6(b) shows the distribution of the 
$m=0$ mode over the different co-latitudinal structures (
$\overline {\langle C_{\ell, m=0} \rangle _r}$), which shows that the large-scale circulation is dominated by the co-latitudinal modes with 
$\ell \leq 3$. With increasing 
$\varepsilon$, 
$\overline {\langle C_{\ell, m=0} \rangle _r}$ is enhanced for low (
$\ell \leq 3$) and high (
$\ell \geq 70$) wavenumbers, compared with 
$\varepsilon =0$.
3.4. Effect of large-scale flow on heat transfer distribution
As we have shown above, the gravity centre offset influences strongly the flow structures. This raises the question of how the local temperature profiles and heat fluxes are affected. Figure 7(a) shows the mean temperature in the meridian plane. For 
$\varepsilon =0$, the temperature is uniform in the co-latitude direction. However, with increasing 
$\varepsilon$, the effect of the north-pole convective jet, which brings hot fluid from the inner to the outer sphere, becomes more pronounced. The effect of the jet is more visible in figure 7(b), which shows the temperature profiles at selected angular positions in the co-latitudinal direction. Most of the temperature drop occurs in the thermal boundary layer regardless of 
$\varepsilon$ and 
$\varphi$. However, for 
$\varepsilon =0.4$ and 
$0.8$, the positive temperature gradient (
${\rm d}\overline {\langle T \rangle _\theta }/{\rm d}r$), observed at 
$\varphi ={\rm \pi} /2$ and shown in the inset, indicates some thermal stratification in the bulk region.
Figure 7. (a) Mean temperature 
$\overline {\langle T \rangle _\theta }(r,\varphi )$ averaged over time and longitudinal direction shown in the meridian plane for 
$\varepsilon =0.05$, 0.4 and 0.8. (b) Mean temperature profiles (
$\overline {\langle T \rangle _\theta }(r)$) at three co-latitudes for 
$\varphi =0$, 
${\rm \pi} /2$ and 
${\rm \pi}$, respectively. The inset for 
$\varphi ={\rm \pi} /2$ shows the temperature gradient 
${\rm d}\overline {\langle T \rangle _\theta }(r) / {\rm d}r$.
Next, we examine how the local heat fluxes depend on 
$\varepsilon$. Figure 8 shows the Nusselt number along the inner 
$Nu_i(\varphi )$ and outer 
$Nu_o(\varphi )$ spheres. It is observed that the local heat fluxes on the outer sphere decrease from 
$\varphi =0$ (north pole) to 
$\varphi ={\rm \pi}$ (south pole), whereas the opposite trend is observed on the inner sphere. This happens because the hot convective jet impinges on the outer boundary layer, resulting in a thinner thermal boundary layer and larger heat flux than without the jet. On the inner sphere, the formation of the convective jet reduces temperature gradients, and therefore the heat transport, at the north pole of the inner sphere, as shown in figure 7(a). Surprisingly, the formation of the large-scale circulation and corresponding non-uniform heat flux distribution has a limited influence on the overall heat flux, which globally is not sensitive to the large-scale dynamics in the flow (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009) (see inset in figure 8b), although some reduction in the heat transport is observed for 
$\varepsilon =0.8$ when the flow becomes stratified more thermally.
Figure 8. Nusselt number 
$Nu$ on (a) the inner sphere and (b) the outer sphere, as functions of the co-latitudinal direction for 
$Ra=10^8$ and 
$n_g=0$. The inset of (b) shows 
$Nu (\varepsilon ) / Nu (\varepsilon =0)$ as a function of 
$\varepsilon$. The error bars represent standard deviations. We note that the error bars for standard errors are smaller than the size of the symbols, which are representative for the uncertainty in the presented global Nusselt values.
3.5. Effect of 
$Ra$ and gravity profile
So far, we have focused on the effect of the gravity centre offset on the flow structures. In this subsection, instead, we show that convective jet and large-scale flow structures are formed for a wide range of 
$n_g$ and 
$Ra$, and that the control parameters modify significantly the intensity of the large-wavelength structures. However, as we have seen above, the global heat flux depends only weakly on 
$\varepsilon$ in the considered parameter regime.
Figures 9(a–d) show visualizations of the convective jet for 
$Ra=10^7$, 
$\varepsilon =0.4$, and different gravity profiles. The figures show that the convective jet forms with 
$u_r>0$ in the north pole region when the gravity centre shifts to the south. As expected, the jet becomes more turbulent when 
$n_g$ is decreased, as this increases the mean buoyancy power over the spherical shell. As shown previously in figure 3(a), the co-latitudinal buoyancy component 
$C_{g,\varphi }g \hat {\varphi }$ points to the north pole in all circumstances. However, the radial buoyancy coefficient 
$C_{g,r}g \hat {r}$ is strongly dependent on the co-latitude. It is well accepted that the stronger radial velocity component tends to induce more intense thermal plumes. However, the convective jet is still observed in the north pole region rather than the south pole region when 
$n_g=-2$, as shown in figure 9(d). This confirms that the large-scale circulation is a result of the baroclinic imbalance due to the misalignment of the gravity isopotential and isothermal lines, which is indicated by steady flows shown in figures 3(b,c).
Figure 9. (a–d) Instantaneous isosurfaces of 
$T=0.2$ for different gravity distributions from 
$n_g=1$ to 
$n_g=-2$ (see (2.7)) for 
$Ra=10^7$ and 
$\varepsilon =0.4$. Panel (e) shows 
$Ra=10^8$, 
$n_g=0$ and 
$\varepsilon =0.4$, for comparison. The colour indicates 
$u_r$ at the isosurface. (Two movies for 
$Ra=10^7$, 
$n_g=0$ and 
$\varepsilon =0.8$ are provided; see supplementary movies 1 and 2.)
The corresponding 
$Re_{rms}(r)$, presented in figure 10(a), shows that the turbulent intensity is enhanced greatly with decreasing 
$n_g$. This reveals that the flow becomes more turbulent when 
$n_g$ decreases, due to the increased buoyancy power over the volume. Correspondingly, figure 10(b) shows that the 
$\mathrm {TKE}$ spectrum of 
$\overline {\langle C_{m} \rangle _r}$ increases for all wavenumbers. We note that 
$\overline {\langle C_{m=0} \rangle _r}$, which dominates 
$\mathrm {TKE}$ and indicates the strength of the large-scale circulation, increases with decreasing 
$n_g$. With increasing 
$Ra$ from 
$10^7$ to 
$10^8$, the normalized 
$Re_{rms}(r)$ is enhanced as shown in figure 10(c). As shown in figure 10(d), the enhancement happens for all wavenumbers, especially for the high wavenumber structures. Thus for a constant 
$\varepsilon =0.4$, the large-scale circulation is more energetic for higher 
$Ra$ and smaller 
$n_g$, due to the increased mean buoyancy power over the volume and the relatively stronger baroclinic torque as indicated in figures 3(b,c).
Figure 10. Plots of (a,c) 
$Re_{rms}(r)$ as function of the gap width, and (b,d) the 
$\mathrm {TKE}$ spectrum of 
$\overline {\langle C_{m} \rangle _r}$ for 
$\varepsilon =0.4$. (a,b) Here, 
$Ra=10^7$ with different 
$n_g$; (c,d) here, 
$n_g=0$ with different 
$Ra$. Figures are normalized by 
$Ra=10^7$ with 
$n_g=0$. The values of buoyancy power for 
$Ra=10^7$ and 
$n_g= 1, 0, -1, -2$, and 
$Ra=10^8$ and 
$n_g= 0$, are 3.3, 5.55, 12.3, 42.8 and 9.46, respectively.
The generation of large-scale structures can induce highly non-uniform co-latitudinal heat flux distribution. We showed previously that the global heat transfer is relatively insensitive to 
$\varepsilon$ for 
$n_g=0$ and 
$Ra=10^8$. Figure 11 shows that a similar trend is observed for different 
$n_g$ and 
$Ra$. In particular, for 
$\varepsilon \leq 0.8$, the difference between 
$Nu(\varepsilon >0)$ and 
$Nu(\varepsilon =0)$ is less than 
$8\,\%$, which indicates that the effect of the large-scale flow structures on the overall heat transport is relatively limited.
Figure 11. Plots of 
$Nu (\varepsilon ) / Nu (\varepsilon =0)$ as a function of 
$\varepsilon$ for different 
$Ra$ and gravity profiles 
$n_g$; see (2.7).
4. Summary and conclusions
We demonstrated using direct numerical simulations that even a small shift of the gravity centre location can introduce pronounced changes in the flow organization and local heat fluxes in spherical RB convection. In the Rayleigh number (
$Ra$) space explored in the present study, we observe the presence of the baroclinic flow, which induces a large-scale circulation and jet for 
$\varepsilon \ge 0.05$. However, for very small 
$\varepsilon =0.01$ and large 
$Ra=10^8$, the baroclinic flow is weakened significantly due to the dominance of the strong convective flow. In turbulent flows with 
$\varepsilon \geq 0.05$, a strong convective jet on the northern side of the inner sphere is formed, which leads to hot/cold fluid going towards the north/south along the inner/outer sphere. The formation of the large-scale circulations has been confirmed by modal analysis. Due to the large-scale flow organization, the heat transfer at the inner and outer sphere is highly non-uniform. In particular, around the north pole, the heat flux is enhanced on the outer sphere due to the impingement of the convective jet, while the formation of the convective jet at the north pole of the inner sphere reduces temperature gradients, and therefore the local heat transport. However, surprisingly, the global heat flux is relatively insensitive to the gravity centre shift. These findings do not seem to depend much on 
$Ra$ (over our explored range) or the particular gravity profile.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2022.360.
Acknowledgements
The authors thank the three anonymous referees for constructive comments that improved the paper. We also acknowledge the Irene at Très Grand Centre de Calcul du CEA (TGCC) under PRACE project 2019215098, and the national e-infrastructure of SURFsara, a subsidiary of SURF cooperation, the collaborative ICT organization for Dutch education and research. We acknowledge PRACE for awarding us access to MareNostrum at Barcelona Supercomputing Center (BSC), Spain (projects 2020225335 and 2020235589).
Funding
G.W. and R.J.A.M.S. acknowledge financial support from ERC (the European Research Council) Starting grant no. 804283 UltimateRB. This work was sponsored by NWO Science for the use of supercomputer facilities.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Buoyancy components of off-centre gravity
The two-step transformation of a vector 
$\boldsymbol {GP}$ in Cartesian coordinates 
$xyz$ to spherical coordinates 
$\varphi \theta r$ (
$x''y''z''$) is shown in figure 12. A random fixed point 
$G(x_G,y_G,z_G)$ is representing the gravity centre location. The other random point, 
$P(r\sin {\varphi } \cos {\theta },r\sin {\varphi } \sin {\theta }, r\cos {\varphi })$, is a moving point representing a fluid parcel. Initially, both 
$P$ and 
$G$ are in Cartesian coordinates 
$xyz$. Using two steps, we will transfer 
$\boldsymbol {GP}$ to spherical coordinates 
$\varphi \theta r$ (
$x''y''z''$).
(i) Rotating
$xyz$ to $x'y'z'$ about $z$ by an angle $\theta$ using the right-hand side rule. The corresponding rotation matrix is
(A 1)
\begin{equation}
R_{z}(\theta)=\left[\begin{array}{@{}ccc@{}} {\cos \theta} &
{-\sin \theta} & {0}\\ {\sin \theta} & {\cos \theta} & {0}
\\ {0} & {0} & {1} \end{array}\right].
\end{equation}
(ii) Rotating
$x'y'z'$ to $x''y''z''$ about $y'$ by an angle $\varphi$ using the right-hand side rule. The corresponding rotation matrix is
(A 2)
\begin{equation}
R_{y}(\varphi)=\left[\begin{array}{@{}ccc@{}} {\cos \varphi} &
{0} & {\sin \varphi} \\ {0} & {1} & {0} \\ {-\sin \varphi}
& {0} & {\cos \varphi} \end{array}\right].
\end{equation}
After the above two steps we obtain that 
$x''y''z''$ overlaps with 
$\varphi \theta r$. The transformation matrix 
$\mathcal {M}_{\theta \varphi }$ from 
$\boldsymbol {GP}_{xyz}$ to 
$\boldsymbol {GP}_{x''y''z''}$ is defined by
\begin{align} \mathcal{M}_{\theta
\varphi}=R_{z}(\theta) \times R_{y}(\varphi) &=
\left[\begin{array}{@{}ccc@{}} {\cos \theta} & {-\sin \theta} &
{0} \\ {\sin \theta} & {\cos \theta} & {0} \\ {0} & {0} &
{1} \end{array}\right] \left[\begin{array}{@{}ccc@{}} {\cos
\varphi} & {0} & {\sin \varphi} \\ {0} & {1} & {0} \\
{-\sin \varphi} & {0} & {\cos \varphi} \end{array}\right]
\nonumber\\ &= \left[\begin{array}{@{}ccc@{}} {\cos
\varphi} \cos \theta & {-\sin
\theta} & {\sin \varphi} \cos
\theta \\ {\cos \varphi} \sin
\theta & {\cos \theta} & {\sin \varphi} \sin \theta \\ {-\sin \varphi} & {0} & {\cos
\varphi} \end{array}\right].
\end{align}
In the Cartesian coordinate frame, 
$\boldsymbol {GP}_{xyz}=(r\sin {\varphi } \cos {\theta }-x_G,\, r\sin {\varphi } \sin {\theta }-y_G,\, r\cos {\varphi }-z_G)$. Using the above transformations, we find that 
$\boldsymbol {GP}$ in spherical coordinates 
$x''y''z''$ is
$$\begin{gather} \boldsymbol{GP}_{x''y''z''}=\boldsymbol{GP}_{xyz} \times \mathcal{M}_{\theta \varphi}= (- x_G \cos\varphi \cos\theta-y_G \cos\varphi \sin\theta + z_G \sin\varphi,\nonumber\\ x_G \sin\theta - y_G \cos\theta, r - x_G \sin\varphi \cos\theta -y_G\sin\varphi \sin\theta -z_G\cos\varphi). \end{gather}$$
In the dimensionless equations for the velocity 
$\boldsymbol {u}$, the buoyancy term is given by 
$gT\,\boldsymbol {GP}/L_{GP}$, where
\begin{equation} L_{GP}=|\boldsymbol{GP}|=\sqrt{(r \sin\varphi \cos\theta - x_G)^2 + (r\sin \varphi \sin \theta-y_G)^2+(r \cos\varphi-z_G)^2}. \end{equation}
The three components in the 
$\hat {\varphi }$, 
$\hat {\theta }$ and 
$\hat {r}$ directions can be decomposed as
\begin{align} \left. \begin{aligned}
C_{g,\theta} & =(x_G \sin\theta - y_G \cos\theta)/L_{GP},
\\ C_{g,r} & =(r - x_G \sin\varphi \cos\theta
-y_G\sin\varphi \sin\theta -z_G\cos\varphi)/L_{GP}, \\
C_{g,\varphi} & =(- x_G \cos\varphi \cos\theta-y_G
\cos\varphi \sin\theta + z_G \sin\varphi)/L_{GP}.
\end{aligned} \right\} \end{align}
Figure 12. Schematic figure of the two-step transformation of a vector 
$\boldsymbol {GP}$ in Cartesian coordinates 
$xyz$ to a spherical coordinates 
$\varphi \theta r$ (
$x''y''z''$) frame. Circular arrows about the 
$z$- and 
$y'$-axes show rotation axes.
Appendix B. Code validations
We validated extensively our code against the results presented by Gastine et al. (Reference Gastine, Wicht and Aurnou2015), who performed a systematic parameter study for spherical RB convection with 
$Pr=1$ covering radius ratios 
$0.2 \le \eta \le 0.95$ for different gravity profiles in table 2. The table shows that 
$Nu$ and the volume-averaged Reynolds number 
$Re_{rms}$ agree within 
$2\,\%$ with the results from Gastine et al. (Reference Gastine, Wicht and Aurnou2015) for all cases. To ensure that the off-centre gravity is implemented correctly, we verified that the flow dynamics for the same 
$L_{\boldsymbol {OG}}$ (see figure 1) are identical for five different gravity centre locations. The last five rows of table 2 confirm that both 
$Nu$ and 
$Re_{rms}$ obtained from these simulations are consistent.
Table 2. Parameters of the spherical RB simulations with 
$Pr=1$, which are compared to the results (
$Nu$(ref) and 
$Re_{rms}$(ref)) of Gastine et al. (Reference Gastine, Wicht and Aurnou2015). Here, 
$G(\theta,r,\varphi )$ is the gravity centre location, and 
$N_\theta,N_r,N_\varphi$ indicate the numbers of grid points in the longitudinal, radial and co-latitudinal directions, respectively. The last five simulations for 
$Ra=3 \times 10^4$ are used to verify that the coefficients 
$C_{g,i}$ (see (A 6)) are implemented correctly, by putting the gravity centre 
$G$ in five different locations with the same 
$L_{\boldsymbol {OG}}$ (see figure 1).
Appendix C. Spectral analysis in spherical coordinates
This appendix describes the spectral analysis procedure used in this study. For further reading on this topic, we refer to pp. 122–145 of MacRobert (Reference MacRobert1947). In planar configurations, Fourier transformations can be applied when performing spectral analysis. However, for a spherical configuration, the basis wave functions are found by solving the Laplace equation. Analogous to the complex exponential in planar configuration, the spherical harmonics 
${Y}_{\ell }^{m}({\theta }, {\varphi })$ read
\begin{equation} {Y}_{\ell}^{m}({\theta}, {\varphi})={P}_{\ell}^{m}({\varphi})\,{\rm e}^{{\rm i} m {\theta}}, \end{equation}
where 
$\ell$ and 
$m$ (
$-\ell \le m \le \ell$) represent the wavenumbers along a meridian and the equatorial plane, respectively. The polar angle 
$\varphi$ ranges from 
$0$ to 
${\rm \pi}$, and 
$\theta$ is the azimuthal angle, in the range 
$0 \leq \theta \leq 2{\rm \pi}$. Also, 
$P_{\ell }^{m} (\varphi )$ are the associated Legendre functions. For computational purposes, the normalized associated Legendre functions in (C2) are more attractive since 
$P_{\ell }^{m}(\varphi )$ can overflow the computer (Swarztrauber Reference Swarztrauber1979):
\begin{equation} \bar{P}_{\ell}^{m}(\varphi)=\left[\frac{2 \ell+1}{2}\,\frac{(\ell-m) !}{(\ell+m) !}\right]^{1 / 2} P_{\ell}^{m}(\varphi). \end{equation}
The use of the spherical harmonics for approximating functions on a sphere is motivated by the fact that they form a complete system of orthogonal functions. Therefore, any function 
$f(\theta, \varphi )$ that is continuous and has continuous derivatives up to second-order may be expanded in an absolutely and uniformly convergent series
\begin{equation} f(\theta, \varphi)=\sum_{\ell=0}^{\infty} \sum_{m={-}\ell}^{\ell} \bar{P}_{\ell}^{m}(\varphi)\left(a_{\ell, m} \cos (m \theta)+b_{\ell, m} \sin (m \theta)\right), \end{equation}where the expansion coefficients are given by
\begin{gather} a_{\ell, m}=\frac{1}{\rm \pi} \int_{0}^{2 {\rm \pi}} \int_{0}^{\rm \pi} f(\theta, \varphi)\, \bar{P}_{\ell}^{m}(\varphi) \cos (m \theta) \sin \varphi \,{\rm d} \theta \,{\rm d} \varphi, \end{gather}
\begin{gather} b_{\ell, m}=\frac{1}{\rm \pi} \int_{0}^{2 {\rm \pi}} \int_{0}^{\rm \pi} f(\theta, \varphi)\, \bar{P}_{\ell}^{m}(\varphi) \sin (m \theta) \sin \varphi \,{\rm d} \theta \,{\rm d} \varphi. \end{gather}
The spherical harmonics 
${Y}_{\ell }^{m}({\theta }, {\varphi })$ are orthonormal, so
\begin{equation} \frac{1}{2{\rm \pi}} \int_{\varphi=0}^{\rm \pi} \int_{\theta=0}^{2 {\rm \pi}} Y_{\ell}^{m} Y_{\ell^{\prime}}^{m^{\prime} *} \,{\rm d} \varOmega=\delta_{\ell \ell^{\prime}}\,\delta_{m m^{\prime}}, \end{equation}
where 
$\delta _{i,j}$ is the Kronecker delta, and 
$\textrm {d} \varOmega =\sin \varphi \,\textrm {d} \theta \,\textrm {d} \varphi$ is the differential solid angle in spherical coordinates. Coefficients (C4) and (C5) are obtained by solving (C3) on a uniform spherical grid by using the SPHEREPACK library (Adams & Swarztrauber Reference Adams and Swarztrauber1999).
The power spectrum on a spherical surface 
$C_{\ell, m}$ and in the longitudinal direction 
$C_{m}$ is defined in (3.2) and (3.3), respectively. To demonstrate the formation of 
$(m=0, \ell >0)$ modes corresponding to the meridional large-scale circulation, we perform a simulation for 
$n_g=-2$ and 
$\varepsilon =0.4$ with 
$Ra=10^3$ (much lower than the critical Rayleigh number when 
$\varepsilon =0$). The advantage of using this steady flow is a regular large-scale circulation generated without longitudinal dependent structures (
$m \neq 0$), which can be seen from figures 13(a,b). The spectrum of the kinetic energy, i.e. 
$\overline {\langle C_{\ell,m=0} \rangle }_r$, is shown in figure 13(c). Except for the mode of 
$(m=0,\ell =0)$ representing the mean flow, apparently the large-scale circulation is represented by 
$(m=0,\ell >0)$ modes.
Figure 13. Spectral analysis of a regular large-scale circulation in a steady flow field for 
$Ra=10^3$, 
$n_g=-2$ and 
$\varepsilon =0.4$. (a) The streamlines (curves coloured by the velocity magnitude) and vortical structures (isosurfaces of the 
$Q$-criterion, 
$Q=0.006$, coloured black). (b) The equatorial cut, meridional cut and radial surface located at 
$r=r_i+0.05d$ of the velocity magnitude. (c) Plot of 
$\overline {\langle C_{\ell,m=0} \rangle }_r$ of the velocity field as a function of the spherical harmonic order 
$\ell +1$ in logarithmic scale. Also shown is the spherical harmonic basis of 
$\ell =0, 1, 2, 3, 4$ for 
$m=0$.
Open access
