4.1. The $O(1)$ problem for the vertical $z$ structure

The lowest-order problem is very simply built on the assumption that thermal diffusion in the radial direction is negligible, with (3.4a) approximated by ${\partial ^2_z}\theta _0=0$. Integration subject to ${\partial _z}\theta _0=0$ at $z=0$ and $1$ determines

(4.2)\begin{equation} \theta_0 = f(R,T) . \end{equation}

Then neglecting the right-hand side of (3.13b), we see that

(4.3a,b)\begin{equation} R^{{-}1}\psi_0 = {Ra}_c\,g(R,T)\,P(z) , \quad \mbox{where } g=f^\prime \end{equation}

(notation (4.1a)), provided that $P(z)$ solves

(4.4a)\begin{equation} {\mathcal{L}}(P) \equiv {\ddddot P} + 4E^{{-}2}P ={-}1, \end{equation}

(notation (4.1b)) (cf. Dowling Reference Dowling1988, (25)). The boundary conditions (2.9c,d) require

(4.4b)\begin{equation} P(0) = P(1) = {\dot P}(0) = {\ddot P}(1) = 0 . \end{equation}

We summarise the solution in Appendix A. It lacks the simplicity of Dowling's equations (26) and (27), applicable to the case of stress-free boundaries.

On neglecting the right-hand side of (3.13a), we obtain

(4.5a)\begin{equation} R\varpi_0 \equiv 2\,g(R,T)\,W(z) , \end{equation}

on use of (4.2) and (4.3), provided that

(4.5b,c)\begin{equation} {\dot W} ={-} {Ra}_c\,P(z),\quad \mbox{ giving } W ={-} {Ra}_c{\langle{P}\rangle_{{0}}^{{z}}}, \end{equation}

after integration subject to $W(0)=0$ (see (4.8a)), implied by $\varpi =0$ at $z=0$.

On further use of (4.2) and (4.3), the lowest-order approximation of (3.1b) is

(4.6a)\begin{equation} R^{{-}1}\varphi_0 = R^{{-}1}\psi_0 + {\partial_{{{R}}}}\theta_0 ={-} g(R,T)\,{\ddot Q}(z) , \end{equation}

where

(4.6b,c)\begin{equation} {\ddot Q} ={-} {Ra}_c\,P(z)-1,\quad \mbox{giving } {\dot Q} = W(z)-z, \end{equation}

after integration, and without loss of generality, the boundary condition choice ${\dot Q}(0)=0$. Hence, on neglect of the left-hand side of (3.8a), integration of its remaining right-hand side with respect to $R$ implies that $R{\mathcal {G}}_0$ is a constant. Then the boundary conditions (2.9a,b) and (2.10c) establish that ${\mathcal {G}}_0=0$. In turn, substitution of (4.6a) into (3.8b), recalling that ${\partial _z}\theta _0=0$ implies ${\langle {\psi _0\,{\partial _z}\theta _0}\rangle }=0$, yields sequentially

(4.7a–c)\begin{equation} R^{{-}1}{\langle{\varphi_0}\rangle} = {\langle{{\mathcal{G}}_0}\rangle} , \quad {\langle{\ddot Q}\rangle} = 0 , \quad {Ra}_c{\langle{P}\rangle} ={-}1, \end{equation}

on use of (4.6a,b). Performing the integral in (4.7b) gives ${\unicode{x27E6} {{\dot Q}}\unicode{x27E7} }=0$, which, having chosen ${\dot Q}(0)=0$, yields ${\dot Q}(1)=0$. So finally, (4.6c) implies that $W(1)=1$, and in summary,

(4.8a–c)\begin{equation} W(0) = 0 , \quad W(1) = 1 , \quad {\dot Q}(0) = {\dot Q}(1) = 0 . \end{equation}

Our $P,W,Q$ notation is adopted to follow the development in (3.8), (3.10) of Chapman & Proctor (Reference Chapman and Proctor1980).

Finally, we note the useful $z$-average identities

(4.9a,b)\begin{equation} {\langle{{\dot W} \bullet}\rangle} ={-} {\langle{W \dot{\bullet}}\rangle}+{\bullet}(1) ,\quad {\langle{{\ddot Q} \bullet}\rangle} ={-} {\langle{{\dot Q} \dot{\bullet}}\rangle}, \end{equation}

which follow from integration by parts and use of the boundary values (4.8). At this early stage, the emergence of ${\dot Q}$ in (4.6c), as a derivative, appears contrived because its integral $Q(z)$ is determined only up to an arbitrary constant of integration. Nevertheless, the way our solution method unfolds, $Q$ itself appears only within the $z$-average ${\langle {{\ddot Q}Q}\rangle }$, which on integration by parts takes the unique value $-{\langle {{\dot Q}^2}\rangle }$ (see (4.9b). Other useful related results are

(4.10)\begin{equation} \left.\begin{aligned} - {\langle{{\ddot Q}W}\rangle} & = {\langle{{\dot Q}{\dot W}}\rangle}\\ & ={-} {Ra}_c{\langle{P{\dot Q}}\rangle} \end{aligned}\right\} = {\langle{{\dot Q}({\ddot Q}+1)}\rangle} = {\langle{\dot Q}\rangle} . \end{equation}

4.2. The $O(\epsilon ^2)$ problem

Just as for the $O(1)$ problem, we begin our $O(\epsilon ^2)$ study with the heat conduction equation (3.4a), whose right-hand side ${\mathcal {N}}_\theta$ (3.5a) is determined at leading order by two terms, $R^{-1}\,{J}(\psi _0 , \theta _0 )= -{Ra}_c\, g^2{\dot P}$ and ${\partial _{{{R}}}}\varphi _0=-(Rg)^\prime {\ddot Q}$. The ensuing ${\mathcal {N}}_\theta$ may be integrated with respect to $z$ so that the corresponding integral of (3.4a) gives

(4.11a)\begin{equation} {\partial_z}\theta_2 ={-}{Ra}_c\,Pg^2+{\dot Q}\,{\partial_{{{{R}}}}^+} g , \end{equation}

which, in view of (4.4b) and (4.8c), satisfies the boundary conditions ${\partial _z}\theta _2=0$ at $z=0, 1$. Multiplication of (4.11a) by $R^{-1}\psi _0 ={Ra}_c\,Pg$ (see (4.3a)) provides the useful result

(4.11b)\begin{equation} - R^{{-}1}\psi_0\,{\partial_z}\theta_2 = {Ra}_c [{Ra}_c\,P^2g^3-P{\dot Q}g\,{\partial_{{{{R}}}}^+} g]. \end{equation}

Moreover, a further integration of (4.11a), which notes $-{Ra}_c{\langle {P}\rangle _{{0}}^{{z}}}=W$ from (4.5c), yields

(4.12a)\begin{equation} \theta_2 = W g^2 + Q\,{\partial_{{{{R}}}}^+} g + f_2(R,T) , \end{equation}

where $f_2(R,T)$, like $f(R,T)$ introduced in (4.2), is at this stage an unknown function, whose value (not needed by us) is fixed only by closure at a higher order. Indeed, since $Q(z)$ is defined only up to an arbitrary constant $\hat{Q}$, the corresponding contribution $\hat{Q}\,{ {\partial _{{{{R}}}}^+} g}$ may be absorbed by $f_2$. The radial derivative of (4.12a) determines

(4.12b,c)\begin{equation} {\partial_{{{R}}}}\theta_2 = 2W g g^\prime + Q \varDelta g + g_2 ,\quad g_2 = f^\prime_2 , \end{equation}

where we have recalled that ${\partial _{{{R}}}} {\partial _{{{{R}}}}^+} = \varDelta$ (see (3.6a) and (4.1c)).

Our next objective is to solve the inhomogeneous equation (3.13b) for $\psi _2$. The leading-order terms on its right-hand side are determined from

(4.13a)\begin{align} {\mathcal{N}}_\varpi &= 2\sigma^{{-}1}\,{Ra}_c(P{\dot W}-{\dot P}W) g\, {\partial_{{{{R}}}}^+} g - 2W\varDelta g ,\end{align}

(4.13b)\begin{align} {\mathcal{N}}_\psi &= \sigma^{{-}1}[{Ra}_c^2(P{\dddot {P}}g\, {\partial_{{{{R}}}}^+} g- {\dot P}{\ddot P} g\,{\partial_{{{{R}}}}^-} g) + 8E^{{-}2} W{\dot W} R^{{-}1}g^2] \nonumber\\ &\quad - 2\,{Ra}_c\,{\ddot P}\varDelta g - \mu^2 g \end{align}

(notation (4.1c,d)). Together with the additional contribution $-{Ra}_c\,R\,{\partial _{{{R}}}}\theta _2$ (use (4.12b)) from its left-hand side, (3.13b) determines

(4.14)\begin{equation} {\partial^4_z}\psi_2+4E^{{-}2}\psi_2 = R[{\mathcal{N}}_\psi - {Ra}_c\,{\partial_{{{R}}}}\theta_2 - 2E^{{-}2} {\langle{{\mathcal{N}}_\varpi}\rangle_{{z}}^{{1}}}], \end{equation}

with ${\mathcal {N}}_\psi$, ${\mathcal {N}}_\varpi$ given by (4.13). The equation must be solved subject to $\psi _2=\partial \psi _2/\partial z =0$ at $z=0$, and $\psi _2=\partial ^2 \psi _2/\partial z^2 =0$ at $z=1$. The solution may be expressed in the form

(4.15a)\begin{align} R^{{-}1}\psi_2 &= P ({Ra}_c\,g_2+ \mu^2g)+P_{{{D}}} \varDelta g +P_{{{{W}}}} g g^{\prime}\nonumber\\ & \quad +\sigma^{{-}1}[P_{{{PP}}} g\,{\partial_{{{{R}}}}^-} g+ (P^+_{{{{PP}}}} +P^+_{{{{WW}}}}) g\,{\partial_{{{{R}}}}^+} g +P_{{{{WW}}}} R^{{-}1}g^2] . \end{align}

Here, the various $P_\bullet (z)$ functions solve

(4.15b,c)\begin{gather} {\mathcal{L}}(P_{{{D}}}) = (4/E^2) {\langle{W}\rangle_{{z}}^{{1}}}- {Ra}_c(2{\ddot P} + Q) ,\quad {\mathcal{L}}(P_{{{W}}}) ={-} 2\,{Ra}_c\,W , \end{gather}

(4.15d,e)\begin{gather}{\mathcal{L}}(P_{{{PP}}}) ={-} {Ra}_c^2\,{\dot P}{\ddot P} ,\quad {\mathcal{L}}(P^+_{{{PP}}}) = {Ra}_c^2\,P{\dddot {P}} , \end{gather}

(4.15f,g)\begin{gather}{\mathcal{L}}(P^+_{{{WW}}}) ={-} (4\,{Ra}_c/E^2) {\langle{P{\dot W}-{\dot P}W}\rangle_{{z}}^{{1}}} ,\quad {\mathcal{L}}(P_{{{WW}}}) = (8/E^2) W{\dot W}, \end{gather}

subject to the boundary conditions $P_\bullet (0)={\dot P}_\bullet (0)=P_\bullet (1)={\ddot P}_\bullet (1)=0$ of (4.4b). So, on multiplying each of (4.15b–g) by $P(z)$, taking the $z$-average, integrating by parts and noting the property ${\mathcal {L}}(P)=-1$ from (4.4a), we obtain the important result

(4.16)\begin{equation} {\langle{P_\bullet}\rangle}={-}{\langle{P_\bullet {\mathcal{L}} (P)}\rangle}={-}{\langle{P {\mathcal{L}} (P_\bullet)}\rangle} \end{equation}

(an extension of the technique employed in Appendix A of Chapman & Proctor Reference Chapman and Proctor1980).

Armed with the result (4.15a), we may now use (3.13a) to obtain

(4.17)\begin{equation} R\,{\partial_z}\varpi_2 ={-} 2R^{{-}1}\psi_2 - {\langle{{\mathcal{N}}_\varpi}\rangle_{{z}}^{{1}}} , \end{equation}

which, upon integration subject to $\varpi _2=0$ at $z=0$, determines $\varpi _2$. However, that result is not needed to close our problem, as we now demonstrate.

4.3. Closure

The amplitude equation for $f$ follows from (3.8a,b), which at lowest order yields

(4.18a)\begin{gather} {\partial_{{{T}}} } {\langle{\theta_0}\rangle} = {\partial_{{{{R}}}}^+} {\mathcal{G}}_2 \equiv R^{{-}1}(R{\mathcal{G}}_2)^\prime , \end{gather}

(4.18b)\begin{gather}R{\mathcal{G}}_2 = {\langle{\phi_2}\rangle} - {\langle{\psi_0\,{\partial_z}\theta_2}\rangle} = {\langle{\psi_2}\rangle} + R\,{\partial_{{{R}}}}{\langle{\theta_2}\rangle} - {\langle{\psi_0\,{\partial_z}\theta_2}\rangle} . \end{gather}

The terms on the right-hand side of (4.18b) are determined respectively by the mean values of (4.15a), (4.12b) and (4.11b). Collecting them together and noting that the two terms involving $f_2^\prime$ cancel, because ${Ra}_c {\langle {P}\rangle }=-1$ in (4.7c) implies $(1+{Ra}_c {\langle {P}\rangle })f_2^\prime =0$, we are left with

(4.19)\begin{align} {\mathcal{G}}_2 &={-} \mu^2\,{Ra}_c^{{-}1}\,g - {\mathcal{F}}_{{{D}}} \varDelta g-{\mathcal{F}}^{{{W}}}_{{{Q}}} g g^\prime-{\mathcal{F}}_{{{Q}}}\, g\,{\partial_{{{{R}}}}^+} g +{\mathcal{F}}_\theta g^3\nonumber\\ & \quad -\sigma^{{-}1}[{\mathcal{F}}_{{{{PP}}}}\, g {\partial_{{{{R}}}}^-} g + ({\mathcal{F}}^+_{{{{PP}}}}+{\mathcal{F}}^+_{{{{WW}}}}) g {\partial_{{{{R}}}}^+} g +{\mathcal{F}}_{{{{WW}}}} R^{{-}1}g^2]. \end{align}

Here, the coefficients of the terms independent of $\sigma$ are

(4.20a,b)\begin{gather} {Ra}_c^{{-}1}={-} {\langle{P}\rangle} ,\quad {\mathcal{F}}_{{{D}}} ={-} {\langle{P_{{{D}}}}\rangle} - {\langle{Q}\rangle} , \end{gather}

(4.20c,d)\begin{gather}\left.\begin{array}{c@{}} {\mathcal{F}}^{{{{W}}}}_{{{{Q}}}} ={-} {\langle{P_{{{{W}}}}}\rangle}-2{\langle{W}\rangle} ,\\ {\mathcal{F}}_{{{{Q}}}} = {Ra}_c {\langle{P{\dot Q}}\rangle} ={-} {\langle{\dot Q}\rangle}, \end{array}\right\} \quad {\mathcal{F}}_\theta = \left\{\begin{array}{@{}l}{Ra}_c^2{\langle{P^2}\rangle} ,\\ {Ra}_c{\langle{{\dot P}W}\rangle} ,\end{array}\right. \end{gather}

where the reductions in (4.20c,d) have respectively involved (4.10) and (4.5b). The remaining coefficients of the terms proportional to $\sigma ^{-1}$ are

(4.20e,f)\begin{gather} {\mathcal{F}}_{{{{PP}}}} ={-} {\langle{P_{{{PP}}}}\rangle} , \quad {\mathcal{F}}^+_{{{{PP}}}} ={-} {\langle{P^+_{{{{PP}}}}}\rangle} , \end{gather}

(4.20g,h)\begin{gather}{\mathcal{F}}^+_{{{{WW}}}} ={-} {\langle{P^+_{{{{WW}}}}}\rangle} , \quad {\mathcal{F}}_{{{{WW}}}} ={-} {\langle{P_{{{{WW}}}}}\rangle} . \end{gather}

Each of the six ${\langle { P_\bullet }\rangle }=-{\langle {P {\mathcal {L}} (P_\bullet )}\rangle }$ in (4.20b,c,e–h) are evaluated following various integrations by parts and repeated use of (4.5b,c), (4.6b,c) and (4.10), giving

(4.21a)\begin{align} {\mathcal{F}}_{{{D}}} &= 2\,{Ra}_c{\langle{{\dot P}^2}\rangle} - {\langle{{\dot Q}^2}\rangle} - {\mathcal{R}}_c^{{-}1}\!{\langle{W^2}\rangle}\nonumber\\ &= 2\,{Ra}_c{\langle{{\dot P}^2}\rangle} - (1+{\mathcal{R}}_c^{{-}1}){\langle{W^2}\rangle} + 2{\langle{zW}\rangle}-\tfrac13 , \end{align}

(4.21b)\begin{align} \tfrac12{\mathcal{F}}^{{{{W}}}}_{{{{Q}}}} &= {\mathcal{F}}_{{{{Q}}}} ={-}{\langle{\dot Q}\rangle} ={-}{\langle{W}\rangle} + \tfrac12 , \end{align}

(4.21c)\begin{align} \tfrac12 {\mathcal{F}}^+_{{{{PP}}}} &= {\mathcal{F}}_{{{{PP}}}} = \tfrac12\,{Ra}_c^2{\langle{{\dot P}^3}\rangle} , \end{align}

(4.21d)\begin{align} \tfrac23 {\mathcal{F}}^+_{{{{WW}}}} &= {\mathcal{F}}_{{{{WW}}}} ={-} (4/E^2){\langle{{\dot P}W^2}\rangle} ={-} (8/E^2)\,{Ra}_c{\langle{P^2W}\rangle} , \end{align}

where in (4.21a) we have introduced the alternative measure

(4.22)\begin{equation} {\mathcal{R}}_c = E^2\,{Ra}_c/ 4, \end{equation}

of ${Ra}_c$. Aided by the identities (4.21b–d), we may reduce (4.19) to

(4.23a)\begin{equation} {\mathcal{G}}_2 ={-} \mu^2\,{Ra}_c^{{-}1}\,g - {\mathcal{F}}_{{{D}}} \varDelta g-{\mathcal{F}}_{{{{Q}}}\sigma} g(3g^\prime+R^{{-}1}g)-\sigma^{{-}1}2{\mathcal{F}}_{{{{WW}}}} R^{{-}1}g^2+{\mathcal{F}}_\theta g^3 , \end{equation}

where

(4.23b)\begin{equation} {\mathcal{F}}_{{{{Q}}}\sigma} = {\mathcal{F}}_{{{Q}}}+\sigma^{{-}1} ({\mathcal{F}}_{{{{PP}}}}+\tfrac12{\mathcal{F}}_{{{{WW}}}}). \end{equation}

In the non-rotating case $E^{-1}=0$, the $E=\infty$ coefficient ${\mathcal {F}}_{{{{WW}}}}(E)$ in (4.23) vanishes. The other coefficients are linked to $A$–$D$ introduced in (1.2b), and their values follow from Table 1, Case C of Chapman & Proctor (Reference Chapman and Proctor1980), which after appropriate scaling (different units) yields

(4.24a)\begin{gather} {Ra}_c(\infty) = 16/A = 320 , \end{gather}

(4.24b)\begin{gather}{\mathcal{F}}_{{{D}}}(\infty) = B/4 = 58/693 , \end{gather}

(4.24c)\begin{gather}{\mathcal{F}}_\theta(\infty) = C = 760/567 , \end{gather}

(4.24d)\begin{gather}{\mathcal{F}}_{{{{Q}}}\sigma}(\infty) = D/6 = 1/18 + (5/126)\sigma^{{-}1}, \end{gather}

composed of

(4.24e,f)\begin{equation} {\mathcal{F}}_{{{Q}}}(\infty) = 1/18 , \quad {\mathcal{F}}_{{{{PP}}}}(\infty)+\tfrac12{\mathcal{F}}_{{{{WW}}}}(\infty) = 5/126. \end{equation}

For the finite $E$ rotating case, the linear problem (4.4) is addressed in Appendix A by considering the Ekman layer style equations (A3) for velocities ${\mathcal {U}}(z)$, ${\mathcal {V}}(z)$ (see (A2)), which relate to ${\dot P}$, $W$ (see (A1b,c)). Since all the ${\mathcal {F}}_\bullet (E)$ coefficients (4.20d) and (4.21a–d) needed to define ${\mathcal {G}}_2$ in (4.23a) depend on ${\dot P}$ and $W$, we are able to determine their values in terms of ${\mathcal {U}}$ and ${\mathcal {V}}$ in Appendix B. Remarkably, the needed $z$-averages (B1), as well as ${Ra}_c^{-1}=4E^{-2}{\mathcal {R}}_c^{-1}$ (see (4.20a) and (4.22)) determined by (A6), may be expressed entirely in terms of the end point values (A5), (A10) at $z=0$ and $1$ of the linear solution. The derivation of the integral results (B2)–(B5) is relegated to Appendix C.

The explicit formulae assembled in Appendices A–C show that

(4.25a)\begin{equation} {\mathcal{F}}_{{{D}}}\gtrless 0\quad \mbox{when } E \gtrless E_c \approx 0.2274 , \end{equation}

in agreement with the value $(E^{-1}=)\ \varOmega _1=4.3966$ given on p. 1347 of Cox (Reference Cox1998). The positivity

(4.25b)\begin{equation} {\mathcal{F}}_\theta > 0 \end{equation}

is guaranteed by (4.20d). We also have

(4.25c)\begin{equation} {\mathcal{F}}_{{{{WW}}}} < 0 , \end{equation}

increasing monotonically through negative values to zero, as $E\to \infty$. Moreover,

(4.25d)\begin{equation} {\mathcal{F}}_{{{Q}}} > 0, \end{equation}

increasing from zero at $E=0$ monotonically to $1/18$ (see (4.24e)), as $E\to \infty$, while

(4.25e)\begin{equation} {\mathcal{F}}_{{{{PP}}}}+\tfrac12{\mathcal{F}}_{{{{WW}}}}\gtrless 0\quad \mbox{ when } E \gtrless {\bar E} \approx 0.4650 , \end{equation}

specifically increasing from $-0.5$ at $E=0$ monotonically to $5/126$ (see (4.24f)), as $E\to \infty$. All the behaviours (4.25) pertain to the plots of ${\mathcal {F}}_\bullet$ versus $E$ in figure 1. Each plot is restricted to the range $E> E_c$, where ${\mathcal {F}}_{{{D}}} > 0$ (see figure 1b), necessary for the application of our long radial length scale asymptotic assumption. The values of ${Ra}_c(E)$, ${\mathcal {F}}_{{{D}}}(E)$ and ${\mathcal {F}}_\theta (E)$ portrayed in figures 1(a–c) are normalised by their $E\to \infty$ values (4.24a–c).

Figure 1. Plots of (a) ${Ra}_c/{Ra}_c(\infty )$, (b) ${\mathcal {F}}_{{{{D}}}}/{\mathcal {F}}_{{{{D}}}}(\infty )$, (c) ${\mathcal {F}}_\theta /{\mathcal {F}}_\theta (\infty )$, and (d) ${\mathcal {F}}_{{{{Q}}}}$ (solid black), $\tfrac 12 {\mathcal {F}}_{{{{WW}}}}$ (grey) and ${\mathcal {F}}_{{{{PP}}}}+\tfrac 12 {\mathcal {F}}_{{{{WW}}}}$ (dashed black) versus $E$ on the range $E>E_c\approx 0.2274$; ${\mathcal {F}}_{{{D}}}(E_c)=0$ (see (4.25a)).

Since the algebra required to determine the results described in Appendices A–C is so intricate, we undertook a numerical check for some specific values of $E$. That involved the direct numerical solution of (4.4) for $P(z)$, including, of course, ${Ra}_c=-1/{\langle {P}\rangle }$ (see (4.7c)). Whence the values of the other ${\mathcal {F}}_\bullet (E)$ coefficients in (4.20) and (4.21), needed for (4.23), were obtained directly by numerical integration. The results for the selected $E$-values are identified by the dots in figure 1, in perfect agreement with the analytic results.

5.1. Axial and outer boundary layer considerations

In addition to the thermal and kinematic boundary conditions (5.1b,c), the equations of motion (2.8) are subject to stress boundary conditions embedded within (2.9). Relevant to that are the tangential components of velocity

(5.2a,b)\begin{equation} R^{{-}1}\,{\partial_{{{R}}}}\psi_0 = {Ra}_c\,{\partial_{{{{R}}}}^+} g\,P(z) ,\quad R\varpi_0 = 2g\,W(z) , \end{equation}

and the vertical and azimuthal stresses proportional to

(5.3a,b)\begin{equation} {\partial_{{{R}}}} (R^{{-}1}\,{\partial_{{{R}}}}\psi_0) = {Ra}_c\,\varDelta g\,P(z) , \quad R\,{\partial_{{{R}}}}\varpi_0 = 2\, {\partial_{{{{R}}}}^-} g\,W(z) , \end{equation}

all on a cylinder $R= \textrm {const.}$ Their appropriate application almost certainly leads to a viscous layer near the outer boundary $R=1$ of radial extent $1-R=O(\epsilon )$, i.e. in a relatively small roughly square region, not accessible by our asymptotics. Though consideration of this layer is needed to determine the solution in the boundary layer, it ought not to influence the ‘mainstream’ solution elsewhere at leading order, so we consider it no further.

As our solutions of the amplitude equation (5.1) have $g\propto R$ as $R\downarrow 0$, the vertical velocity and angular velocity determined by (5.2) are finite on the axis $R=0$, while in turn the stresses (5.3) vanish there, as required. The outer boundary $R=1$ is more interesting. Consideration of the expression (4.23a) for ${\mathcal {G}}_2$ shows that together, $g(1,T)=0$ and ${\mathcal {G}}_2(1,T)=0$ (see (5.1b,c)) imply

(5.4)\begin{equation} \varDelta g = 0\quad \mbox{ at } R = 1 . \end{equation}

This means that whereas the azimuthal velocity (5.2b) is brought to rest ($g = 0$), as required by (2.9b), the vertical velocity (5.2a) is not (${ {\partial _{{{{R}}}}^+} g}\not =0$), contrary to (2.9b). Interestingly, a similar problem would arise in the case of a stress-free outer boundary. In that case, the vertical stress (5.3a) vanishes ($\varDelta g=0$), while the azimuthal stress (5.3b) does not (${ {\partial _{{{{R}}}}^-} g}\not =0$, essentially $g^\prime \not =0$ again).

Whether the outer boundary is rigid or stress-free, only one (but not both) of the stress boundary conditions can be met, so a boundary layer is required. Interestingly, for the non-rotating problem $E^{-1}=0$, there is no azimuthal flow. So for that case, the problem with a stress free boundary ${\rm \Delta} g=0$ (see (5.3a) and (5.4)) at $R=1$ does not require a boundary layer, whereas the case of a rigid boundary, needing ${ {\partial _{{{{R}}}}^+} g}=0$, does.

We cannot overemphasise our assumption that the $R$ length scale is large compared to $\epsilon$. So whenever solutions of (5.1) vary significantly on that relatively short $\epsilon$ length scale, i.e. the vertical extent, our asymptotic assumption is violated and the solution of (5.1) must be viewed with suspicion. The worth of such solutions can be assessed only by comparison with the DNS of the complete problem, a matter that we address in § 7.

5.2. The thermal energy balance (3.12)

Our understanding of the nature of the convection and flow is aided by consideration of the thermal energy equation (3.12). The fact that $\theta _0$ in (4.2) is independent of $z$, implying ${\partial _z} \theta _0=0$, has important consequences, which include ${\langle {\varphi _0}\rangle }=0$ (see (4.7a)). In turn, the leading-order terms on the right-hand side of (3.12) vanish,

(5.5a,b)\begin{equation} - \epsilon^{{-}2}{{\langle{ {\langle{\theta_0^2}\rangle}}\rangle}} - 2{{\langle{ {\langle{\theta_0\theta_1}\rangle}}\rangle}} = 0 , \quad - {{\langle{ {\langle{R^{{-}1} \varphi_0 ({\partial_{{{R}}}} \theta_0)}\rangle}}\rangle}} = 0, \end{equation}

leaving only $O(\epsilon ^2)$ terms. What remains, involving $\varphi _2=\psi _2+R\,{\partial _{{{R}}}}\theta _2$ (see (3.1b)), is

(5.6)\begin{equation} \tfrac12 \,{\mathrm d}_{{{{T}}} }{{\langle{ {\langle{\theta_0^2}\rangle}}\rangle}} ={-} {{\langle{ {\langle{({\partial_z} \theta_2)^2}\rangle}}\rangle}} - {{\langle{ {\langle{(R^{{-}1}\varphi_0({\partial_{{{R}}}} \theta_2)}\rangle}}\rangle}} - {{\langle{ {\langle{(R^{{-}1}\psi_2+{\partial_{{{R}}}}\theta_2)({\partial_{{{R}}}} \theta_0)}\rangle}}\rangle}} . \end{equation}

Aided by the expressions (4.11a) for ${\partial _z}\theta _2$, (4.12b) for ${\partial _{{{R}}}}\theta _2$ and (4.15a) for $R^{-1}\psi _2$, the right-hand side may be evaluated tediously. A more direct derivation of the result,

(5.7)\begin{equation} \tfrac12 \,{\mathrm d}_{{{{T}}} } {{\langle{ {\langle{f^2}\rangle}}\rangle}} ={-} {{\langle{ {\langle{g{\mathcal{G}}_2}\rangle}}\rangle}} = \mu^2\,{Ra}_c^{{-}1}\!{{\langle{\langle{g^2}\rangle}\rangle}} - {\mathcal{F}}_{{{D}}} {{\langle{ {\langle{({\partial_{{{{R}}}}^+} g)^2}\rangle}\rangle}}} + \sigma^{{-}1}2{\mathcal{F}}_{{{{WW}}}}{{\langle{ {\langle{ R^{{-}1}g^3}\rangle}}\rangle}}-{\mathcal{F}}_\theta {{\langle{ {\langle{g^4}\rangle}}\rangle}} , \end{equation}

follows from evaluating the weighted average ${{\langle { {\langle {f\,{\partial _{{{T}}} }{f}}\rangle }}\rangle }}$ using (5.1a) and integrating by parts. The resulting integral is evaluated using the formula (4.23a) for ${\mathcal {G}}_2$. In it, the term with the coefficient ${\mathcal {F}}_{{{{Q}}}\sigma }$ evaporates because ${{\langle { {\langle {g^2(3g^{\prime }+R^{-1}g)}\rangle }}\rangle }}= {{\langle { {\langle {R^{-1}[Rg^3]^{\prime }}\rangle }}\rangle }} =0$. Evidently, instability is driven by the term $\mu ^2\,{Ra}_c^{-1}\!{{\langle { {\langle {g^2}\rangle }}\rangle }}$ when $\mu ^2>0$, and damped by the term $-{\mathcal {F}}_\theta {{\langle { {\langle {g^4}\rangle }}\rangle }}\ ( <0)$ (see (4.25b)). The diffusive term $-{\mathcal {F}}_{{{D}}}{{\langle { {\langle {({\partial _{{{{R}}}}^+} g)^2}\rangle }\rangle }}}$ damps only when $E>E_c$ (${\mathcal {F}}_{{{D}}}>0$), otherwise when $E< E_c$ (${\mathcal {F}}_{{{D}}}<0$), it drives the instability (see (4.25a)). The sign of ${{\langle { {\langle { R^{-1}g^3}\rangle }}\rangle }}$ in the term $\sigma ^{-1}2{\mathcal {F}}_{{{{WW}}}}{{\langle { {\langle { R^{-1}g^3}\rangle }}\rangle }}$ is important in determining the nature of the convection, as we argue in the next paragraph. Further consequences are highlighted by our weakly nonlinear theory of § 6.2 below.

Typically, the meridional flow consists of a single (horizontally elongated) cell, for which the direction of circulation may be identified by the sign of the $z$-average of the scaled vertical velocity, namely

(5.8)\begin{equation} {\mathcal{W}}(R,T) = {\langle{R^{{-}1}\,{\partial_{{{R}}}} \psi_0}\rangle} = ({\partial_{{{{R}}}}^+} g)\, {Ra}_c{\langle{P}\rangle} ={-} {\partial_{{{{R}}}}^+} g, \end{equation}

(use (4.7c)), evaluated on the axis $R=0$. There, (5.8) determines

(5.9)\begin{equation} {\mathcal{W}}(0,t) ={-} \tfrac12 g^\prime(0,T) \begin{cases} > 0 & \mbox{upwelling} ,\\ <0 & \mbox{downwelling} .\end{cases} \end{equation}

So, for a single cell with $g(R,T)<0\ (>0)$ on $0< R<1$, we have upwelling (downwelling) on the axis. With that scenario ${{\langle { {\langle { R^{-1}g^3}\rangle }}\rangle }}<0\ (>0)$, and since ${\mathcal {F}}_{{{{WW}}}}<0$ from (4.25c), the term $\sigma ^{-1}2{\mathcal {F}}_{{{{WW}}}}{{\langle { {\langle { R^{-1}g^3}\rangle }}\rangle }}>0\ (<0)$ renders the upwelling state to be preferred. This term, however, vanishes in both the infinite Prandtl number limit $\sigma \to 0$ and the non-rotating limit $E\to \infty$ for which ${\mathcal {F}}_{{{{WW}}}}\uparrow 0$.

6.1. The linear problem

The linearised version of the amplitude equations (6.2a,b) are

(6.4)\begin{equation} \partial_{ {\sf{t}} }{{\sf{f}}} ={-} {\partial_{{{{R}}}}^+}(\lambda {\sf{g}} + \varDelta {\sf{g}} ) , \quad \partial_{ {\sf{t}} }{{\sf{g}}} ={-} \varDelta(\lambda {\sf{g}} +\varDelta {\sf{g}}). \end{equation}

Note that ${\sf {f}}$ is determined only up to an arbitrary constant, which we ignore below in order to reduce clutter. The solutions that satisfy the boundary conditions (5.1b,c) are

(6.5a)$$\begin{gather} - j_m^{{-}2}\,{\partial_{{{{R}}}}^+} {\sf{g}} = {\sf{f}} ={-} j_m^{{-}1}\,A_m({\sf{t}})\,{J}_0(\,j_m R) , \end{gather}$$

(6.5b)$$\begin{gather}{\sf{g}} = {\sf{f}}^\prime = A_m({\sf{t}})\,{J}_1(\,j_m R) , \end{gather}$$

(6.5c)$$\begin{gather}\varDelta {\sf{g}} ={-} j_m^2{\sf{g}} , \end{gather}$$

(6.5d)$$\begin{gather}{\sf{G}} ={-} (\lambda {\sf{g}} +\varDelta {\sf{g}} ) = (-\lambda+j_m^2){\sf{g}} , \end{gather}$$

where $j_m$ is the $m$th zero of the Bessel function ${J}_1$, chosen such that ${\sf {g}}(0,{\sf {t}})={\sf {g}}(1,{\sf {t}})=0$ (see (6.5b)), with the consequence ${\sf {G}}(0,{\sf {t}})={\sf {G}}(1,{\sf {t}})=0$ by (6.5d), provided that

(6.6)\begin{equation} {\mathrm d}_{{\sf{t}}}{A_m} = j_m^2 [\lambda - j_m^2 ]A_m \end{equation}

(where ${\mathrm d}_{{\sf {t}}}\equiv {\mathrm d}/{\mathrm d} {\sf {t}}$). The requirement $\lambda = \mu ^2/({Ra}_c\,{\mathcal {F}}_{{{D}}}) > 0$ (see (6.3a)) for instability is met only when ${\mathcal {F}}_{{{D}}}>0$, which requires $E>E_c\approx 0.2274$ (see (4.25a) and figure 1b).

The steady modes ${\mathrm d}_{{\sf {t}}}A_m=0$ correspond to

(6.7)\begin{equation} \lambda = \lambda_m = j_m^2 . \end{equation}

The lowest mode $m=1$ is identified by the first non-zero zero of ${J}_1$, namely

(6.8a,b)\begin{equation} j_1 \approx 3.83171 , \quad \lambda_1 \approx 14.68197 . \end{equation}

Thus, correct to $O(\epsilon ^2)$, the critical Rayleigh number determined by (6.3b) is

(6.9)\begin{equation} {Ra}_c^{\dagger} = {Ra}_c (1 + \epsilon^2 j_1^2 {\mathcal{F}}_{{{D}}}) . \end{equation}

We note that $\psi \propto R{\sf {g}} \propto R\,{J}_1(\,j_1 R)$ is maximised when ${\mathrm d}_R [R\,{J}_1(\,j_1 R)]=0$, or equivalently, ${J}_0(\,j_1 R)=0$. Thus the first zero of ${J}_0$ determines the location

(6.10a)\begin{equation} R_{{{{M}}} 0} \approx 0.6276 \end{equation}

of the maximum, which itself is proportional to $R_{{{{M}}} 0}\,{J}_1(\,j_1 R_{{{{M}}} 0})$, where

(6.10b)\begin{equation} {J}_1(\,j_1 R_{{{{M}}} 0}) \approx 0.5191 . \end{equation}

As a corollary, ${\sf {f}}\propto {J}_1(\,j_0 R)$ reverses sign across $R_{{{{M}}} 0}$.

6.2. Small-amplitude expansion about critical

For our finite amplitude solutions, a useful measure of supercriticality, relative to ${Ra}_c^{\dagger}$ of (6.9), is

(6.11)\begin{equation} \aleph = \dfrac{\lambda-j_1^2}{j_1^2} = \dfrac{{Ra} - {Ra}_c^{\dagger}}{{Ra}_c^{\dagger} - {Ra}_c} . \end{equation}

In the following two subsubsections, we consider a small-amplitude expansion

(6.12a)\begin{gather} j_1^2\aleph = \lambda - j_1^2 = \delta \varLambda_1 + \delta^2 \varLambda_2 +\cdots , \end{gather}

(6.12b)\begin{gather}{}[{\sf{f}}, {\sf{g}}] - \delta[{\sf{f}}_0, {\sf{g}}_0](R,t)= \delta^2 [{\sf{f}}_1, {\sf{g}}_1](R,t) + \delta^3 [{\sf{f}}_2, {\sf{g}}_2](R,t) +\cdots \end{gather}

($\epsilon ^2\ll \delta \ll 1$) for the lowest steady $m=1$ mode (6.5a,b), which solves

(6.12c,d)\begin{equation} \square\,{\sf{g}}_0 = 0 , \quad \mbox{where} \ \square\,\bullet\equiv (\varDelta + j_1^2 )\bullet \end{equation}

(see (6.5c,d)). The objective is to construct the equation governing the slow evolution of the amplitude $A_1(t)$. The positive parameter $\delta \ (\ll 1)$ is chosen at our convenience to aid identification of the terms, which balance at various orders of $\delta \ (>0)$.

6.2.1. The case $\beta /\sigma =O(1)$

For this generic case, we consider only the leading-order terms $\delta \varLambda _1$ and $\delta ^2 [{\sf {f}}_1, {\sf {g}}_1]$ on the right-hand sides of (6.12a,b). Anticipating evolution on the slow time scale $\delta ^{-1}$, we write

(6.13)\begin{equation} A_1({\sf{t}})={\sf{A}}({\sf{T}}_{1}) , \quad {\sf{T}}_{1}=\delta {\sf{t}} , \quad \partial_{ {\sf{t}}}=\delta\,\partial_{ {\sf{T}}_{1}} , \quad {\mathrm d}_{{\sf{t}}}=\delta \,{\mathrm d}_{{\sf{T}}_{1}} . \end{equation}

Then at $O(\delta )$, (6.2a,c) determine

(6.14)\begin{equation} R^{{-}1}[R\,{ {\boldsymbol \square} }\,{\sf{g}}_1 ]^{\prime} ={-} \partial_{ {\sf{T}}_{1}} {\sf{f}}_0 + R^{{-}1}[-\varLambda_1 R{\sf{g}}_0 + \alpha R {\sf{g}}_0(3{\sf{g}}_0^{\prime} +R^{{-}1}{\sf{g}}_0) - (\beta/\sigma){\sf{g}}_0^2]^{\prime} \end{equation}

(notation (6.12d)), where significantly the cubic term $+{\sf {g}}^3$ in (6.2c), being smaller by another factor $O(\delta )$, has been omitted.

We take the radial average of (6.14) weighted by $J_0(\,j_1 R)$ to eliminate the left-hand side and so obtain

(6.15)\begin{equation} j_1^{{-}2} \,{\mathrm d}_{{\sf{T}}_1 }{\sf{A}} = \varLambda_1{\sf{A}} + N_2 (\beta/\sigma){\sf{A}}^2 . \end{equation}

Here, we have used the property

(6.16a)\begin{equation} {{\langle{ {\langle{{\sf{g}}^2_0(3{\sf{g}}_0^{\prime}+R^{{-}1}{\sf{g}}_0)}\rangle}}\rangle}} = {{\langle{ {\langle{R^{{-}1}\,{\mathrm d}_R(R{\sf{g}}_0^3)}\rangle}}\rangle}} = 0 , \end{equation}

as ${\sf {g}}_0(1,t)=0$, to eliminate the term proportional to $\alpha$, and noted that

(6.16b)\begin{equation} {{\langle{ {\langle{{J}^2_0(\,j_1 R)}\rangle}}\rangle}} = {{\langle{ {\langle{{J}^2_1(\,j_1 R)}\rangle}}\rangle}} = \tfrac12\,{J}^2_0(\,j_1) \approx 0.0811 . \end{equation}

In addition, since ${{\langle { {\langle {R^{-1}\,{J}^3_1(\,j_1 R)}\rangle }}\rangle }}\approx 0.0821$, we have

(6.16c)\begin{equation} N_2 = {{{\langle{ {\langle{R^{{-}1}\,{J}^3_1(\,j_1 R)}\rangle}}\rangle}}}/{{{\langle{ {\langle{{J}^2_1(\,j_1 R)}\rangle}}\rangle}}} \approx 1.0124 . \end{equation}

The bifurcation of the steady trivial solutions ${\sf {A}}=0$ of (6.15) at $\varLambda _1=0$ to the neighbouring finite amplitude solutions

(6.17)\begin{equation} {\sf{A}} ={-} \varLambda_1\sigma/(\beta N_2) \gtrless 0 \quad \mbox{for } \varLambda_1 \lessgtr 0 , \end{equation}

since $N_2 (\beta /\sigma )>0$, is transcritical (see Guckenheimer & Holmes Reference Guckenheimer and Holmes1983; Cross & Hohenberg Reference Cross and Hohenberg1993, figure 6). Obviously, the solutions ${\sf {A}}>0$ (upwelling on the axis) for $\varLambda _1<0$ are unstable and will evolve to a large amplitude for which the weakly nonlinear theory developed here no longer applies. We expand on this matter next.

6.2.2. The case $\beta /\sigma =O(\delta )$

To capture the stabilising term $+{\sf {g}}^3$ of (6.2c) omitted in (6.14), we consider the case $\beta /\sigma =O(\delta )$. In practice, this limit restricts our analysis to the case $E\gg 1$ of slow rotation, but nevertheless reveals, in more detail, the nature of possible finite amplitude solutions of (6.2a). Since the magnitude of the term $N_2 (\beta /\sigma ){\sf {A}}^2$ in (6.15) is reduced by a factor $O(\delta )$, we reduce $\aleph$ in (6.12a) by the same amount. Accordingly, we set

(6.18)\begin{equation} \varLambda_1 = 0 , \end{equation}

while lengthening the time scale:

(6.19)\begin{equation} A_1({\sf{t}})={\sf{A}}({\sf{T}}_{2}) , \quad {\sf{T}}_{2}=\delta^2 {\sf{t}} , \quad \partial_{ {\sf{t}}}=\delta^2\, \partial_{ {\sf{T}}_{2}} , \quad {\mathrm d}_{{\sf{t}}}=\delta^2 \,{\mathrm d}_{{\sf{T}}_{2}} . \end{equation}

Then at $O(\delta )$, (6.14) simplifies leaving us with only the $\alpha$ term on the right-hand side. After integration of what remains and application of the end point conditions ${\sf {G}}=0$ at $R=0$ and $1$, where ${\sf {g}}_0=0$ too, we obtain

(6.20a)\begin{equation} { {\boldsymbol \square} }\,{\sf{g}}_1 = \alpha {\sf{g}}_0(3{\sf{g}}_0^{\prime}+R^{{-}1}{\sf{g}}_0), \end{equation}

with solution

(6.20b)\begin{equation} {\sf{g}}_1 = \alpha {\sf{g}}_0 {\sf{f}}_0 ={-} {\sf{A}}^2\alpha j_1^{{-}1}\,{J}_0(\,j_1 R)\,{J}_1(\,j_1 R), \end{equation}

vanishing at $R=0$ and $1$. In this way, correct to the lowest two orders, we have

(6.21a)\begin{gather} {\sf{f}} = \delta{\sf{f}}_0[1 + \tfrac12\alpha\delta{\sf{f}}_0] , \end{gather}

(6.21b)\begin{gather}{\sf{g}} = {\sf{f}}^\prime = \delta{\sf{g}}_0[1 + \alpha\delta{\sf{f}}_0]. \end{gather}

Consideration of the maximum of $R({\sf {g}}_0+\delta {\sf {g}}_1)$ reveals a shift in the linear value $R_{{{{M}}} 0}$ from (6.10a) for the maximum of $\psi$ to $R_{{{M}}}$ given by the solution of

(6.22a)\begin{equation} j_1\,{J}_0(\,j_1 R) = \alpha \delta {\sf{A}}[{J}_0^2(\,j_1 R)-{J}_1^2(\,j_1 R)]. \end{equation}

The Taylor series expansion of ${J}_0(\,j_1 R)$ about $R_{{{{M}}} 0}$, at which ${J}_0(\,j_1 R_{{{{M}}} 0})=0$, reveals the lowest-order result

(6.22b)\begin{equation} R_{{{M}}} - R_{{{{M}}} 0} = \alpha j_1^{{-}2} \delta {\sf{g}}_0(R_{{{{M}}} 0}) = \alpha \delta {\sf{A}} j_1^{{-}2}\,{J}_1(\,j_1 R_{{{{M}}} 0}) , \end{equation}

with ${J}_1(\,j_1 R_{{{{M}}} 0})\approx 0.5191$ (see (6.10b)). The result (6.22b) quantifies the outward (inward) shift of the maximum of $|\psi |$ for solutions that upwell (downwell), ${\sf {A}} >(<)\ 0$, on the axis.

The $O(\delta ^2)$ terms in (6.2a,c) give

(6.23a)\begin{equation} R^{{-}1}[R\,{ {\boldsymbol \square} }\,{\sf{g}}_2 ]^{\prime} ={-} \partial_{ {\sf{T}}_{2}}{\sf{f}}_0 \!+\! R^{{-}1}[-\varLambda_2 R{\sf{g}}_0 \!+\! \alpha(3R({\sf{g}}_0 {\sf{g}}_1)^{\prime}+2{\sf{g}}_0{\sf{g}}_1) - \delta^{{-}1}(\beta/\sigma){\sf{g}}_0^2 - R{\sf{g}}_0^3]^{\prime}; \end{equation}

cf. (6.14). As in § 6.2.1, we take the radial average of (6.23a) weighted by $J_0(\,j_1 R)$ to eliminate the left-hand side. Recalling that ${\sf {g}}_1=\alpha {\sf {g}}_0 {\sf {f}}_0$ (6.20b), evaluation of the term proportional to $\alpha$ is aided by the identity

(6.23b)\begin{equation} {\sf{g}}_0 (3({\sf{g}}_0 {\sf{g}}_1)^{\prime}+2 R^{{-}1}{\sf{g}}_0{\sf{g}}_1) = 2\,{\partial_{{{{R}}}}^+}({\sf{g}}^2_0{\sf{g}}_1) - \alpha{\sf{g}}^4_0 . \end{equation}

In this way, we obtain

(6.24a,b)\begin{equation} j_1^{{-}2} \,{\mathrm d}_{{\sf{T}}_2 }{\sf{A}} = \varLambda_2{\sf{A}} + N_2 (\beta/\sigma)\delta^{{-}1}{\sf{A}}^2 - \varUpsilon N_3{\sf{A}}^3 , \quad \varUpsilon = 1+\alpha^2 , \end{equation}

in which

(6.25)\begin{equation} N_3 = {{{\langle{ {\langle{{J}^4_1(\,j_1 R)}\rangle}}\rangle}}}/{{{\langle{ {\langle{{J}^2_1(\,j_1 R)}\rangle}}\rangle}}} \approx 0.2517 , \end{equation}

since ${{\langle { {\langle {{J}^4_1(\,j_1 R)}\rangle }}\rangle }}\approx 0.02041$.

Equation (6.24a), albeit valid only when $\beta /\sigma =O(\delta )$, reveals the nature of the bifurcation beyond the transcritical regime identified in § 6.2.1 for $\beta /\sigma =O(1)$. As the steady finite amplitude solutions ${\sf {A}}$ satisfy

(6.26)\begin{equation} \varLambda_2 + N_2(\beta/\sigma){\sf{A}} - \varUpsilon N_3{\sf{A}}^2 = 0 , \end{equation}

the transcritical bifurcation at $\varLambda _2=0$, described by (6.17), becomes the tangent to the parabola (6.26). For that, there are two positive ${\sf {A}}$ solutions on $\varLambda _{{min}}<\varLambda _2<0$, which coalesce at $\varLambda _2=\varLambda _{{min}}$ with value ${\sf {A}}={\sf {A}}_{{min}}$, where

(6.27a,b)\begin{equation} \varLambda_{{min}} ={-} \frac12 N_2(\beta/\sigma){\sf{A}}_{{min}} , \quad {\sf{A}}_{{min}} = \frac{N_2(\beta/\sigma)}{4\varUpsilon N_3} \quad ( N_2(\beta/\sigma)>0 ). \end{equation}

Presumably, for $\varLambda _{{min}}<\varLambda _2<0$, only the upper branch ${\sf {A}}>{\sf {A}}_{{min}}$ is stable, whereas for $\varLambda _2>0$, both the positive and negative ${\sf {A}}$ branches are stable. This presumption suggests that large-amplitude solutions of the DNS for the full problem (2.8)–(2.9) exist in the generic case $\beta /\sigma =O(1)$, but are not accessible by the weakly nonlinear theory of § 6.2.1.

6.2.3. The non-rotating case $E^{-1}=0$

Interestingly, the transcritical instability identified by (6.17) degenerates when $\beta /\sigma =0$. That happens in the non-rotating case $E^{-1}=0$, upon which we comment briefly here. It is a special case of that in § 6.2.2 with the quadratic term $N_2 (\beta /\sigma )\delta ^{-1}{\sf {A}}^2$ absent from (6.24a). As a consequence, the bifurcation at $\varLambda _2=0$ is a pure pitchfork. For $\varLambda _2>0$, the two steady-state solutions $\pm |{\sf {A}}|$ are determined by the vanishing of what remains of (6.24a), which, noting $j_1^2\aleph =\delta ^2\varLambda _2$ (see (6.12a)), gives

(6.28)\begin{equation} |\delta {\sf{A}}|^2 = j_1^2\aleph/(\varUpsilon N_3). \end{equation}

We emphasise that the symmetry of the bifurcation (6.28), possessing solutions $\pm |{\sf {A}}|$, is a low order result. Taken to next order, the solution (6.21), in which ${\sf {f}}_0$ and ${\sf {g}}_0$ are proportional to ${\sf {A}}$, is clearly not invariant under the sign change ${\sf {A}}\mapsto -{\sf {A}}$. Moreover, both correction terms $\alpha \delta {\sf {f}}_0$ in (6.21a,b) change sign, as $R$ crosses $R_{{{{M}}} 0}$, at which ${\sf {f}}_0=0$, as noted below (6.10b). The shift of the maximum (6.22b) for each, obtained using (6.24b) and (6.28), is

(6.29)\begin{equation} R_{{{M}}} - R_{{{{M}}} 0} ={\pm} \dfrac{1}{j_1}\sqrt{\dfrac{\alpha^2 \aleph}{(1+\alpha^2)N_3}}\,{J}_1(\,j_1 R_{{{{M}}} 0}) . \end{equation}

7.1. The non-rotating case, $E^{-1}=0$

For the case $E^{-1}=0$, $\sigma =0.1$, ${Ra}=330$, we compare in figure 2 the streamlines obtained from the DNS (figures 2a,c) and the asymptotics (figures 2b,d). Relative to the onset values ${Ra}_c=320$, ${Ra}_c^{\dagger} \approx 323.9321$, the supercriticality (6.11) of the finite amplitude solution is

(7.2a,b)\begin{equation} \aleph \approx 1.54315,\quad \mbox{or equivalently}\quad \lambda \approx 37.33844 . \end{equation}

The remaining ${\sf {G}}$ coefficients in (6.2c) are $\beta =0$ and

(7.2c)\begin{equation} \alpha \approx 0.1659 + 0.11848\sigma^{{-}1} \approx 1.3506 \quad \mbox{for } \sigma=0.1 . \end{equation}

Figure 2. No rotation case ($E^{-1}=0$, $\sigma =0.1$, ${Ra}=330$). A comparison in the $r$–$z$ plane of results obtained from DNS with asymptotics (labelled $A$): (a) $\psi _{{{{DNS}}}}>0$, (b) $\psi _{{{{A}}}}>0$, (c) $\psi _{{{{DNS}}}}<0$, and (d) $\psi _{{{{A}}}}<0$. Colour scale from $-1.3$ (blue) to 0 (green) to $1.3$ (red).

As we stressed in § 6.2.3, the pitchfork bifurcation at ${Ra}_c^{\dagger}$ sheds two solutions: one characterised by upwelling on the axis, near which $\psi >0$ (figures 2a,b); the other characterised by downwelling on the axis, near which $\psi <0$ (figures 2c,d). Though the bifurcation is symmetric with the infinitesimal maximum of $|\psi |$ at $R=R_{{{{M}}} 0}$, on increasing $\lambda$ that maximum shifts outwards for the $\psi >0$ solutions and inwards for the $\psi <0$ solutions. Such behaviour was predicted by (6.22b) for the case of small but finite amplitude motion. However, the value $|R_{{{M}}} - R_{{{{M}}} 0}|\approx 0.2696$ determined from (6.29), though qualitatively plausible, overestimates the shifts visible in figure 2, because $\aleph \approx 1.54315$ from (7.2a) is too large for the small-$\delta$ asymptotics of § 6.2 to provide quantitative accuracy. By contrast, the excellent agreement of the DNS (figures 2a,c) with the numerical solutions (figures 2b,d) of

(7.3)\begin{equation} (R^{{-}1}(R{\sf{g}})^{\prime})^{\prime} ={-} \lambda {\sf{g}} + \alpha {\sf{g}}(3{\sf{g}}^{\prime}+R^{{-}1}{\sf{g}}) + {\sf{g}}^3 , \end{equation}

namely ${\sf {G}}=0$ (see (6.2c)) with $\beta =0$, is most encouraging.

We can make an interesting comparison of the contour plots in figure 2 with those in figure 4(b) of Chapman et al. (Reference Chapman, Childress and Proctor1980) for their internally heated case exhibiting up/down asymmetry like us. We may capture the structure of the steady-state version of their Cartesian asymptotic equation (15) by dropping the curvature terms in (7.3), where, in our § 1 notation, $R$ has become $X$. This leaves ${\sf {g}}^{\prime \prime }=-\lambda {\sf {g}} +3\alpha {\sf {g}}{\sf {g}}^{\prime }+{\sf {g}}^3$, but note the sign reversal in (6.1b). The location $X=X_{{{M}}}$ of their maximum $|\psi |$ occurs at the mid-point $X_{{{{M}}} 0}=0.5$ at onset but, on increasing $\lambda$ shifts towards the downwelling side boundary due to the quadratic nonlinearity $3\alpha {\sf {g}}{\sf {g}}^{\prime }$, exactly as we predict in (6.29) and find (figures 2a–d) for $R=R_{{{M}}}$.

In figure 3, we plot the maximum value of $|\psi |$ on the entire domain, but signed depending on whether the solution describes upwelling ($\psi >0$) or downwelling ($\psi <0$) on the symmetry axis $R=0$. We note that at given ${Ra}$, the amplitude $\max |\psi |$ of the upwelling solution is greater than that for the downwelling solution. This is a finite amplitude effect that the weakly nonlinear calculation ($\delta \ll 1$) in § 6.2.3 could not identify, at any rate to the order taken. Note too that the solution portrayed in figure 2 at ${Ra}_c=330$ is very close to the bifurcation point in figure 3, yet, as already mentioned, is outside the range of validity of our small-$\delta$ weakly nonlinear asymptotics of § 6.2. Bearing those limitations in mind, it is remarkable how well our long length scale (small-$\epsilon$) asymptotic amplitude equation (7.3) works, giving good maximum amplitude up to remarkably large ${Ra}=2000$ and beyond (not portrayed). A partial asymptotic explanation is provided in the second paragraph of § 7.2.3 (below).

Figure 3. Bifurcation diagram for the no-rotation case ($E^{-1}=0$, $\sigma =0.1$), showing $\textrm {sign}(\psi )\max |\psi |$ as a function of ${Ra}$. Black dots indicate DNS; green squares indicate asymptotics.

7.2. The rotating case, $E( >E_c)$ finite

The small-$\delta$ analysis of § 6.2.1 identified a transcritical bifurcation, for which the subcritical branch is presumably unstable. The origin of that instability is encapsulated by the quadratic term $(\beta /\sigma ) R^{-1}{\sf {g}}^2$ in the expression (6.2c) for ${\sf {G}}$. The analysis of § 6.2.2, valid for sufficiently small $\beta /\sigma$, identified possible recovery on a stable steady solution upper branch. Whether or not such a branch exists for finite $\beta /\sigma$ remains a matter of speculation, a consideration that emphasises the importance of the size of $\beta /\sigma$. From a general point of view, complications that limit the validity of the approach are likely, as $E$ decreases towards $E_c$. Furthermore, the importance of the $(\beta /\sigma )$ term must increase with decreasing $\sigma$. In the following subsubsections, we investigate how far we can decrease $E$ and $\sigma$ and yet still obtain useful asymptotic results.

7.2.1. Meridional flow

Inspection of the asymptotic results illustrated in figure 1 shows that the coefficients (4.20), which appear in our expression (4.23) for ${\mathcal {G}}_2$ of our amplitude equation (5.1a), vary measurably, only on decreasing $E$ from $\infty$, at about $E=1$. On further decrease of $E$, the variation becomes more significant. So, as a tentative first step, we consider the case $E=1$, $\sigma =0.3$ (moderately small), ${Ra}=345$. Relative to the critical values ${Ra}_c\approx 325.3612$, ${Ra}_c^{\dagger} \approx 329.0900$, supercriticality is measured by

(7.4a,b)\begin{equation} \aleph \approx 4.26679\quad \mbox{or equivalently},\quad \lambda \approx 77.3269 . \end{equation}

The remaining ${\sf {G}}$ coefficients are

(7.4c–e)\begin{equation} \alpha \approx 0.4868 , \quad \beta \approx 0.1043 , \quad \beta/\sigma \approx 0.34772 . \end{equation}

We illustrate the streamlines for $\psi \gtrless 0$ in figures 4(a–d) following the style of figure 2 and exhibiting many of the same features. To highlight any differences between the asymptotics and DNS, we plot horizontal and vertical cross-sections in figures 4(e,f), respectively. The agreement is almost perfect except for the steep descent curves, $\psi >0$, in figure 4(e) between about $r=8$ and the end $r=10$ (recall that $r=10R$). This may be explained by the outer boundary layer caused by the outer rigid boundary condition ${\partial _r}\psi |_{r=10}=0$ (see (2.9b)), which is not met by the asymptotic solution. A similar, but weaker boundary layer is evident in the more gently sloping ascent curves, $\psi <0$.

Figure 4. Rotating case ($E=1$, $\sigma =0.3$, ${Ra}=345$): (a) $\psi _{{{{DNS}}}}>0$, (b) $\psi _{{{{A}}}}>0$, (c) $\psi _{{{{DNS}}}}<0$, and (d) $\psi _{{{{A}}}}<0$. Colour scales from $0$ (green) to $2.5$ (red) in (a,b), and from $-2.2$ (blue) to $0$ (green) in (c,d). (e) Horizontal cross-sections at $z=0.5$, and (f) vertical cross-sections at $r=5$, of the $\psi$ fields shown in (a–d). DNS ($\psi _{{{{DNS}}}}$) in red; asymptotics ($\psi _{{{{A}}}}$) in black.

Other than the presence of rotation in figure 4, the use of the lower Prandtl number $\sigma =0.1$ in figure 3 is significant as it increases the influence of inertia. We will return to this point in § 7.2.2

We test matters further in figure 5, which addresses the case $E=0.5$ at the same Prandtl number $\sigma =0.3$ but increased Rayleigh number ${Ra}=360$. Relative to the critical values ${Ra}_c\approx 341.4403$, ${Ra}_c^{\dagger} \approx 344.5690$, supercriticality is measured by

(7.5a,b)\begin{equation} \aleph \approx 4.95105,\quad \mbox{or equivalently},\quad \lambda \approx 87.3732 . \end{equation}

The remaining ${\sf {G}}$ coefficients are

(7.5c–e)\begin{equation} \alpha \approx 0.2498 , \quad \beta \approx 0.4453 , \quad \beta/\sigma \approx 1.4843 . \end{equation}

There is not much change in the streamline patterns of figures 5(a–d) from that displayed in figures 4(a–d). Indeed, the cross-sections in figures 5(e,f) compare well with similar right-hand boundary layer discrepancies visible in figure 5(e). However, a more worrying feature of that figure is the small but clearly evident differences outside that layer between $r=0$ and $r=8$, which cannot be explained by boundary layer arguments. Indeed, studies of even more testing cases of smaller $\sigma$ and/or $E$ reveal even greater $\psi$ differences. For them, the key to the failure is linked to the azimuthal motion due to the rotation.

Figure 5. Rotating case ($E=0.5$, $\sigma =0.3$, ${Ra}=360$). Same as figure 4. Colour scales from $0$ (green) to $3$ (red) in (a,b), and from $-2$ (blue) to $0$ (green) in (c,d).

7.2.2. Azimuthal flow

In § 7.2.1 we considered only the meridional flow. The complete solution involves the interaction of the meridional and azimuthal flows through their coupling via the Coriolis force. In this subsubsection, we investigate that interaction by considering the azimuthal velocity $v=r\omega$ (see (2.7a)). We portray the DNS and asymptotic results for $v$ in figures 6(a,b,d,e) and 7(a,b,d,e) for the cases that correspond to figures 4(a–d) and 5(a–d), respectively. However, in style, there is one important new addition in figures 6(c,f) and 7(c,f) that we describe as ‘hybrid’, which for the moment must be ignored together with the extra dashed curves in figures 6(g,h) and 7(g,h). On their omission, the remaining comparisons of the DNS and asymptotics are visibly poor. To avoid possible confusion, we stress that, unlike figures 4(e) and 5(e), where the upper $\psi >0$ curves correspond to the top figures 4(a,b) and 5(a,b), in figures 6(g) and 7(g), the lower $v<0$ curves correspond to the top figures 6(a–c) and 7(a–c). Generalised and expressed succinctly, $v\lessgtr 0$ corresponds to $\psi \gtrless 0$ almost everywhere, with some notable exceptions near the axis.

Figure 6. Rotating case ($E=1$, $\sigma =0.3$, ${Ra}=345$): (a) $v_{{{{DNS}}}}<0$, (b) $v_{{{{A}}}}<0$, (c) $v_{{{{H}}}}<0$, (d) $v_{{{{DNS}}}}>0$, (e) $v_{{{{A}}}}>0$, and (f) $v_{{{{H}}}}>0$. Colour scales from $-0.35$ (blue) to $0_+$ (i.e. a little positive) (green) in (a–c), and from $0$ (green) to $0.6$ (red) in (d–f). (g) Horizontal cross-sections at $z=0.5$, and (h) vertical cross-sections at $r=5$, of the $v$ fields shown in (a–f). DNS ($v_{{{{DNS}}}}$) in red; asymptotics ($v_{{{{A}}}}$) in solid black; hybrid ($v_{{{{H}}}}$) in dashed black.

Figure 7. Rotating case ($E=0.5$, $\sigma =0.3$, ${Ra}=360$). Same as figure 6. Colour scales from $-1$ (blue) to $0_+$ (i.e. some positive) (green) in (a–c), and from $0$ (green) to $0.8$ (red) in (d–f).

The discrepancies visible in the azimuthal flow contour plots in figures 6(a,b,d,e) and 7(a,b,d,e) are brought into sharp focus by comparing the red DNS and black asymptotic curves in figures 6(g) and 7(g), which describe radial cross-sections. Together, they indicate that for the case of up/downwelling on the axis, the asymptotics over/underestimate the (correctly predicted by the DNS) magnitude $|rv|$ of the angular momentum advected away from/towards the rotation axis in the neighbourhood of the upper boundary. This asymptotic failure is a low Prandtl number effect, i.e. the increased role of inertia, exacerbated by curvature effects manifested by the various powers of $r$ in the angular momentum equation (2.8b), which lead to large azimuthal velocity gradients that violate the long radial length scale assumption on which the asymptotics is based. Though the asymptotic trends are not far off the mark near the outer boundary $r=10$, particularly for the $v<0$ (corresponding to $\psi >0$) curves, they are definitely unsatisfactory elsewhere.

The upshot of the above assessment is that the feedback of the azimuthal flow on the meridional flow is relatively weak in the parameter ranges considered. That said, the azimuthal flow clearly influences the meridional flow as evinced by the fact that the critical Rayleigh number is a function of $E$. So we may suppose that though our azimuthal flow predicted by the asymptotics is flawed, it is sufficiently accurate to generate a totally acceptable meridional flow as illustrated in figures 4 and 5.

On the basis that our asymptotically predicted $\psi$ is good, we solved the azimuthal component of the momentum equation, namely (2.8a) for the angular velocity $\omega$, with that $\psi$, subject to the $\omega$ boundary conditions appearing in (2.9). We call the solutions of this linear problem ‘hybrid’ solutions. The hybrid solutions in figure 6 agree very well with the DNS, vindicating the hybrid approach. On the one hand, this indicates that the asymptotics is on the right track, but that its parameter range of validity is limited. For more testing parameter values, the hybrid and DNS $v$ shown in figure 7 continue to compare reasonably well, but discrepancies are beginning to emerge. They can be explained from the evidence in figure 5(e) that the asymptotic $\psi$ results, on which the hybrid $v$-solution builds, are losing a little accuracy. Evidently, on pushing the parameter values much further, the asymptotic $\psi$ will be too poor to enable the construction of useful hybrid results.

An important feature of the asymptotics is that, at lowest order, the vertical $z$ profile is the same for all $r$, though, of course, the profile amplitude changes. With this restriction, if $v$ has only one sign at some $r$, then it cannot exhibit a sign reversal at another $r$. That means that the DNS and hybrid solutions, portrayed in figures 7(a,c), exhibiting a sign reversal across the contour beginning on the axis at $z\approx 0.5$ and terminating on the lower boundary just beyond $r=8$, cannot be described by our asymptotics shown in figure 7(b).

7.2.3. A large ${Ra}$ application

Oruba et al. (Reference Oruba, Davidson and Dormy2017) portray results for the case $E=0.1$, $\sigma =0.5$, with ${Ra}=15\ 000$ in their figure 3(a–d). From our point of view, the parameter values $E=0.1< E_c$ outside the domain of validity for the amplitude modulation equation (5.1), and ${Ra}$ large, are extreme. Nevertheless, it is instructive to make a tentative comparison of the DNS results, recalculated and displayed as $\psi _{{{DNS}}}$ and $v_{{{DNS}}}$ in our figures 8 and 9, with results from a yet more extreme version of our hybrid approach outlined below. Our idea is motivated by the encouraging comparison in figure 3 of our max$|\psi |$ amplitudes for the DNS and asymptotic evaluation of $\psi$ in no-rotation cases at largish ${Ra}$. Their robustness suggests that such $\psi =\psi _{{{{E}}}\infty }$ (for $E^{-1}=0$) might provide a plausible approximation of $\psi _{{{DNS}}}$ at finite rotation (i.e. for $E^{-1}\not =0$) – at any rate from a qualitative point of view.

Figure 8. Rotating case ($E=0.1$, $\sigma =0.5$, ${Ra}=15\ 000$): (a) $\psi _{{{{DNS}}}}$, and (b) $\psi _{{{{E}}}\infty }$. Colour scale from $-100$ to $100$. (c) Horizontal cross-sections at $z=0.5$, and (d) vertical cross-sections at $r=5$, of $\psi /r$ for the $\psi$ fields shown in (a,b): DNS ($\psi _{{{{DNS}}}}/r$) in red; $E^{-1}=0$ ($\psi _{{{{E}}}\infty }/r$) in green.

Figure 9. Rotating case ($E=0.1$, $\sigma =0.5$, ${Ra}=15\ 000$): (a) $v_{{{{DNS}}}}$, and (b) $v_{{{{H}}}}$. Colour scale from $-100$ to $100$. (c) Horizontal cross-sections at $z=0.5$, and (d) vertical cross-sections at $r=5$, of the $v$ fields shown in (a,b): DNS ($v_{{{{DNS}}}}$) in red; hybrid ($v_{{{{H}}}}$) in dashed green.

To assess our hypothesis, we plot asymptotic $\psi _{{{{E}}}\infty }$ results in figure 8 for the parameter values of Oruba et al. (Reference Oruba, Davidson and Dormy2017) but, of course, by definition replace their $E=0.1$ with $E^{-1}=0$. For that case, we recall that ${Ra}_c=320$, ${Ra}_c^{\dagger} \approx 323.9321$, while, for their large ${Ra}=15\ 000$, we have

(7.6a)\begin{equation} \aleph \approx 3732.3350,\quad \mbox{or equivalently},\quad \lambda \approx54 812.7155 , \end{equation}

in place of (7.2a,b). Noting that $\beta =0$ and

(7.6b)\begin{equation} \alpha \approx 0.1659 + 0.11848\sigma^{{-}1} \approx 0.4028 \quad \mbox{for } \sigma=0.5 , \end{equation}

in place of (7.2c), the large-$\lambda$ asymptotic mainstream solution of (7.3) is

(7.6c)\begin{equation} {\sf{g}} \approx \sqrt{\lambda} \approx 234.1212 . \end{equation}

This corresponds to the approximate solution

(7.7a)\begin{equation} 0 = {\mathcal{G}}_2/g \approx{-} \mu^2\,{Ra}_c^{{-}1} + {\mathcal{F}}_\theta g^2 ={-} \mu^2\, {Ra}_c^{{-}1} + {Ra}_c^2{\langle{P^2}\rangle}\,g^2 \end{equation}

of (5.1a), noting (4.23a) and (4.20d). It describes the balance between the buoyant driving and the nonlinear convection of heat as traced via $R^{-1}\,{J}(\psi _0 , \theta _0 )= -{Ra}_c\,g^2{\dot P}$ and (4.11a,b). Moreover, (7.7a) determines the leading-order result

(7.7b)\begin{equation} R^{{-}1}\psi_0 = {Ra}_c\,g\,P(z) \approx \mu\,{Ra}_c^{{-}1/2}\,P(z)/{\langle{P^2}\rangle}^{1/2} \end{equation}

(see (4.3a)) and, on use of $(\epsilon \mu )^2={Ra}-{Ra}_c$ (see (3.3)), equivalently

(7.7c)\begin{equation} \dfrac{\psi_0}{r} \approx \sqrt{\dfrac{{Ra}-{Ra}_c}{{Ra}_c}}\, \dfrac{P(z)}{{\langle{P^2}\rangle}^{1/2}}, \end{equation}

independent of $r$. The result (7.7c) holds everywhere except in the boundary layers, roughly square regions adjacent to the lateral boundaries $r=0$ and $r=10$, where the solution is invalid. Those layers are evident in figure 8(c), which describes the horizontal cross-section $r^{-1}\psi _{{{{E}}}\infty }\ (\propto {\sf {g}})$ (note the factor $r^{-1}$ absent in previous cross-sections). In the mainstream $2\lessapprox r\lessapprox 8$, the agreement of $r^{-1}\psi _{{{{E}}}\infty }$ with $r^{-1}\psi _{{{DNS}}}$ for the rotating case, $E=0.1$, is qualitatively remarkable, in view of the tenuous assumptions made. It suggests that rotation modifies but does not control the meridional flow. From this point of view, figure 8 sheds new light on the no-rotation upper branch results portrayed in figure 3. There, only results up to ${Ra}=2000$ are illustrated, but calculations up to ${Ra}=30\ 000$ ($>15\ 000$, used in figure 8) were also performed. As the percentage errors ceased to change over that considerable extension, those results are not reported here. The same asymptotic–DNS agreement is also evident on the right of figure 8(c), from which the similar sizes of $\max \psi _{{{{E}}}\infty }$ and $\max \psi _{{{{DNS}}}}$ (albeit for $E=0.1$) may be estimated.

We undertook a hybrid calculation (referred to as ‘hybrid-$\{E=\infty \}$’), employing $\psi =\psi _{{{{E}}}\infty }$ derived from the $E^{-1}=0$ asymptotics described above, rather than $\psi _{{{A}}}$, for which the asymptotics is irrelevant in the case of interest, $E=0.1$. The $v_{{{H}}}$ contours for that hybrid calculation are illustrated in figure 9(b). Beyond $r\approx 0.5$, they compare favourably with the DNS illustrated in figure 9(a). A more precise quantitative measure comes from the horizontal cross-section in figure 9(c). The failure of $v_{{{H}}}$, relative to the true $v_{{{DNS}}}$, for $r\lessapprox 5$ is readily traced to the singular behaviour of $r^{-1}\psi _{{{{E}}}\infty } (\propto {\sf {g}})$ for small $r$. The maximum of $v_{{{H}}}$ is located at $r\approx 0.3$, which also measures the width of the $\psi _{{{{E}}}\infty }$ boundary layer visible near $r=0$ in figure 8(c). The presence of strong $v_{{{H}}}$ on $0.3\lessapprox r \lessapprox 5$ stems from the overestimation of the strength of the meridional flow, as measured by $\psi _{{{{E}}}\infty }$ over that domain. Such strong flow advects the angular momentum $rv_{{{H}}}$ and pertains to the proposal in Oruba et al. (Reference Oruba, Davidson and Dormy2017, Reference Oruba, Davidson and Dormy2018) that at large ${Ra}$, angular momentum tends to be constant on streamlines.

At ${Ra}=360$, not far above critical, both $v_{{{DNS}}}$ and $v_{{{H}}}$ portrayed in figures 7(a) and 7(c) are generally negative except for a small region close to the lower boundary but terminating before $r=9$, where $v_{{{DNS}}}$ and $v_{{{H}}}$ are both positive. Interestingly, as ${Ra}$ is increased, that region of positive $v$ expands to largely fill all the space except for a small region near the outer boundary. This feature is illustrated by $v_{{{DNS}}}$ in figure 9(a). Remarkably, in view the almost draconian hybrid hypothesis employed, it is also captured by $v_{{{H}}}$ in figure 9(b). The encouraging agreement vindicates the long horizontal length scale hypothesis for the meridional cell, which is possibly stabilised by the differential rotation caused by angular momentum transport.