2.1. Description

A fluid of constant kinematic viscosity $\nu$ is enclosed between two spheres (see below). The inner one is rigidly rotating with a constant angular velocity $\varOmega _c$ while the outer shell supports a prescribed tangential surface stress $\tau _{*}(\theta )$; $R$ and $\eta R$ are the radii of the outer and inner shells, respectively, and $\theta$ is the colatitude. We sketch out this model in figure 1. We consider the driving stress to be sufficiently weak that the rotation period of the system is much shorter than the typical turnover times associated with the flow in the rotating frame of reference. We thus consider the nonlinear terms to be negligible, which we justify a posteriori in appendix A. The dimensionless equation of vorticity and mass conservation equation in the rotating frame with angular velocity $\varOmega _c$, and in the context of an ‘anelastic’ flow, read

(2.1)\begin{equation} \left. \begin{gathered} \frac{\partial}{\partial t} (\boldsymbol{\nabla} \times \rho {\boldsymbol{u}})+ \boldsymbol{\nabla} \times ({\boldsymbol{e}}_z \times \rho {\boldsymbol{u}}) = E \boldsymbol{\nabla} \times {\boldsymbol{f}}_{visc}, \\ \boldsymbol{\nabla} \boldsymbol{\cdot} \rho {\boldsymbol{u}} = 0, \end{gathered} \right\} \end{equation}

where we have used $R$ as the length scale, $(2\varOmega _c)^{-1}$ as the time scale and $\rho _c$, the density at the inner shell, as the density scale. Here, $\boldsymbol{e}_z$ is the unit vector in the z-direction, $\boldsymbol{u}$ is the non-dimensionalised velocity vector, $E$ is the Ekman number defined as

(2.2)\begin{equation} E=\frac{\nu}{2 \varOmega_c R^2}, \end{equation}

and

(2.3)\begin{equation} {\boldsymbol{f}}_{visc}= \rho (\Delta {\boldsymbol{u}} + \tfrac{1}{3}\boldsymbol{\nabla} \textrm{div} {\boldsymbol{u}} ) + 2(\boldsymbol{\nabla} \rho \boldsymbol{\cdot} \boldsymbol{\nabla}){\boldsymbol{u}} + \boldsymbol{\nabla} \rho \times (\boldsymbol{\nabla} \times {\boldsymbol{u}}) - \tfrac{2}{3} \boldsymbol{\nabla} (\rho \textrm{div} {\boldsymbol{u}}), \end{equation}

is the dimensionless viscous force. This expression of the viscous force is meant to represent the case where the envelope is pervaded by some small-scale turbulence with constant diffusive properties as often used in stellar physics (Brandenburg et al. Reference Brandenburg, Jennings, Nordlund, Rieutord, Stein and Tuominen1996; Käpylä, Mantere & Brandenburg Reference Käpylä, Mantere and Brandenburg2012). Centrifugal effects are neglected altogether, and the flows are considered axisymmetric.

Figure 1.

Schematic view of the system: the inner shell of radius $\eta R$ rotates with an angular velocity $\varOmega _c$ while the outer one of radius $R$ rotates differentially at $\varOmega _s(\theta )$ by means of a prescribed tangential surface stress $\tau _{*}(\theta )$. The dashed lines correspond to the edges of the tangent cylinder $\mathcal {C}$, circumscribing the inner sphere.

This system is then completed by boundary conditions. On the outer shell, we impose a specified stress, which represents the torque resulting from the stellar wind. The simplest dimensionless expression of such a stress, which we take equatorially symmetric, is

(2.4)\begin{equation} \tau(\theta)=\sigma_{r \phi}/\rho = -A\sqrt{\frac{3}{4{\rm \pi}}} \sin \theta, \end{equation}

where $[\sigma ]$ is the dimensionless stress tensor and $A$ is a positive constant that sets the amplitude of the prescribed (braking) stress. Although the wind imposes a small radial flow at the surface, we shall ignore it in this first approach, and impose a vanishing normal velocity on the outer shell. Finally, we impose no-slip boundary conditions at the inner boundary, as justified by the important viscosity jump at the stellar envelope/core interface. Boundary conditions thus read

(2.5a)\begin{gather} u_r = u_{\theta} = u_{\phi} = 0 \quad \textrm{at} \ r=\eta, \end{gather}

(2.5b)\begin{gather} \left. \begin{gathered} u_r = r\frac{\partial}{\partial r}\left( \frac{u_{\theta}}{r}\right) = 0,\\ \frac{\sigma_{r\phi}}{\rho}=r \frac{\partial}{\partial r}\left( \frac{ u_{\phi}}{r}\right) = -A\sqrt{\frac{3}{4{\rm \pi}}} \sin \theta \end{gathered} \right\} \quad \textrm{at} \ r=1. \end{gather}

2.2. Numerical method

We use a spectral decomposition and expand the fields in spherical harmonics for the angular part and Chebyshev polynomials for the radial part (see Rieutord & Valdettaro Reference Rieutord and Valdettaro1997). We write

(2.6)\begin{equation} {\boldsymbol{q}}=\rho {\boldsymbol{u}} = \sum^{+\infty}_{l=0} \sum^{+l}_{m=-l} u_{m}^{l}(r) \boldsymbol{R}^{m}_l + v_{m}^{l}(r) \boldsymbol{S}^{m}_l + w_{m}^{l}(r) \boldsymbol{T}^{m}_l , \end{equation}

where

(2.7a–c)\begin{equation} \boldsymbol{R}^{m}_l = Y^{m}_l {\boldsymbol{e}}_r, \quad \boldsymbol{S}^{m}_l = \boldsymbol{\nabla} Y^{m}_l, \quad \boldsymbol{T}^{m}_l = \boldsymbol{\nabla} \times \boldsymbol{R}^{m}_l, \end{equation}

with $Y^{m}_l$ being the usual normalised scalar spherical harmonic function. Because the ${\boldsymbol {q}}$-flow is divergenceless both in the incompressible and anelastic models, $v_{m}^{l}$ can be expressed as a function of $u_m^l$ only (Rieutord Reference Rieutord1987). Projecting the vorticity equation (2.1) on $\boldsymbol {R}^{m}_l$ and $\boldsymbol {T}^{m}_l$ for an axisymmetric flow ($m=0$), we obtain the following system of equations for the radial parts

(2.8)\begin{equation} \left. \begin{gathered} \frac{\partial w^l}{\partial t} - \left[A_{l-1}^{l} r^{l-1}\frac{\partial}{\partial r}\left(\frac{r u^{l-1}}{r^{l-1}}\right) + A_{l+1}^{l} r^{-l-2}\frac{\partial}{\partial r}(r^{l+3} u^{l+1})\right] = E h^l,\\ \frac{\partial}{\partial t}( \Delta_l r u^l) + B_{l-1}^{l} r^{l-1}\frac{\partial}{\partial r}\left(\frac{w^{l-1}}{r^{l-1}}\right) + B_{l+1}^{l} r^{-l-2}\frac{\partial}{\partial r}(r^{l+2} w^{l+1}) =-\frac{E}{r} \left(\,f^l - \frac{\partial r g^l}{\partial r} \right) ,\\ v^l=\frac{1}{r \varLambda} \frac{\partial r^2 u^l}{\partial r} , \end{gathered} \right\} \end{equation}

where $\varLambda =l(l+1)$. Here, $A_{l-1}^{l}$, $A_{l+1}^{l}$, $B_{l-1}^{l}$ and $B_{l+1}^{l}$ are the coupling coefficients defined as

(2.9a,b)\begin{equation} A_{l-1}^{l} = A_l^{l-1}=\frac{1}{l}\sqrt{\frac{1}{(2l-1)(2l+1)}} , \quad B_{l-1}^{l} = B_l^{l-1}=\sqrt{\frac{l(l^2-1)}{(2l-1)(2l+1)}}. \end{equation}

The viscous force is projected on the spherical harmonics basis as well, namely

(2.10)\begin{equation} {\boldsymbol{f}}_{visc} = \sum^{+\infty}_{l=0} \sum^{+l}_{m=-l} f_{m}^{l}(r) \boldsymbol{R}^{m}_l + g_{m}^{l}(r) \boldsymbol{S}^{m}_l + h_{m}^{l}(r) \boldsymbol{T}^{m}_l , \end{equation}

where

(2.11)\begin{equation} \left. \begin{gathered} f^l = 2 \frac{\partial \rho}{\partial r} \frac{\varLambda b v^l - 2 bu^l}{r} + \frac{\varLambda \rho}{r} z^l + \frac{\partial }{\partial r} \left( \frac{4 \rho}{3} d^l \right),\\ g^l=\frac{1}{r}\frac{\partial}{\partial r} (r \rho z^l) + \frac{2}{r}\frac{\partial \rho}{\partial r} (bu^l - bv^l) + \frac{4 \rho}{3 r} d^l,\\ h^l= \rho \Delta_l b w^l + r\frac{\partial \rho}{\partial r} \frac{\partial}{\partial r} \left( \frac{b w^l}{r}\right), \end{gathered} \right\} \end{equation}

with

(2.12a,b)\begin{equation} z^l= \frac{1}{r}\frac{\partial r b v^l}{\partial r} - \frac{b u^l}{r}, \quad d^l = \frac{1}{r^2}\frac{\partial r^2bu^l}{\partial r} - \frac{\varLambda b v^l}{r}, \end{equation}

and where $b(r)=1/\rho (r)$ is the inverse density function.

Finally, the boundary conditions in the rotating frame of reference are projected on the spherical harmonics and read

(2.13)\begin{equation} \left. \begin{array}{ll} u^l = v^l = w^l = 0, & \textrm{at} \ r=\eta, \\ u^l = r\dfrac{\partial }{\partial r}\left( \dfrac{bv^l}{r}\right) = 0 , \quad r \dfrac{\partial }{\partial r}\left( \dfrac{ bw^l}{r}\right) = - A\delta_{l1}, & \textrm{at}\ r=1, \end{array} \right\} \end{equation}

where $\delta _{ij}$ is the Kronecker symbol.

3.1. The transient phase

As to assess the relevance of the steady-state approximation for geophysical and astrophysical applications, it is important to determine the time scales the flow characteristics are governed by. To do so, we measure the time at which the evolution of the total angular momentum in the rotating frame may be considered to its end. We note that the change in angular momentum of the fluid is due to the difference between the torque exerted on the fluid by the outer boundary condition and the torque that the fluid exerts on the steady inner sphere. Therefore, a stationary state is reached when

(3.1)\begin{equation} \Delta \varGamma \equiv \varGamma(1) - \varGamma(\eta) = 0, \end{equation}

where $\varGamma (1)$ and $\varGamma (\eta )$ are the torques about the $z$-axis exerted on the outer and inner boundary surfaces, respectively. The torque about the $z$-axis exerted on a layer located at some radius $r$ can be written

(3.2)\begin{equation} \varGamma(r)=\int_{\partial r} \boldsymbol{\sigma} \boldsymbol{\cdot} {\boldsymbol{e}}_z\,\textrm{d}S , \end{equation}

where $\boldsymbol {\sigma }$ is the local stress applied on the sphere of radius $r$. We expand $\boldsymbol {\sigma }$ on the spherical harmonics basis, hence the torque about the $z$-axis reads

(3.3)\begin{equation} \varGamma(r)= - \int_{\partial r} r \sin \theta \rho(r) \sum_{l} r^2 t^l \frac{\partial Y_l}{\partial\theta}\,\textrm{d}\varOmega, \end{equation}

where

(3.4)\begin{equation} t^l=r \frac{\partial}{\partial r} \left(\frac{w^l}{r} \right). \end{equation}

Using Legendre polynomials recurrence relations, (3.3) actually reads

(3.5)\begin{equation} \varGamma(r)= 4 \sqrt{\frac{\rm \pi}{3}} r^3 \rho(r) t^1(r) . \end{equation}

Applying boundary conditions (2.13), (3.1) can finally be written

(3.6)\begin{equation} \Delta \varGamma = -4 \sqrt{\frac{\rm \pi}{3}} \left( A \rho_s +\left. \eta^3 \frac{\partial w^1}{\partial r}\right|_{r=\eta}\right) , \end{equation}

where the non-dimensional surface density $\rho _s=\rho (1)=1$ for the considered incompressible model. We monitor the evolution of the relative difference between the torques $\Delta \varGamma /\varGamma (1)$ for various Ekman numbers and show the result in figure 2. It shows that on a time scale $O(E^{-1})$ the stress is completely communicated to the interior flow, and a steady state is reached.

Figure 2.

Relative difference between the torque applied on the outer sphere and the torque exerted by the fluid on the stationary inner sphere, as a function of the reduced time $E t$. We find the steady state to be reached on a $O(E^{-1})$ time scale, that is on a viscous time scale, in the asymptotic regime of small Ekman numbers.

Of course in the astrophysical or geophysical context, it may become necessary to relax the constant rotation of the inner sphere. The torque exerted by the fluid on this boundary can make the inner sphere spin-down. According to the angular momentum theorem

(3.7)\begin{equation} \varGamma_*(\eta R)=\frac{\textrm{d} L^{core}_z}{\textrm{d}t_*} = I_c \frac{\textrm{d} \varOmega_c}{\textrm{d}t_*}, \end{equation}

where $\varGamma _*$ is the dimensional torque, $I_c$ is the moment of inertia of the core assumed to be in solid-body rotation and $L^{core}_z$ is its total angular momentum. Starred quantities are dimensional.

The core spin-down then engenders an additional Euler force $\dot {\varOmega _c}{\boldsymbol {e}}_z \times {\boldsymbol {r}}$ since our frame is attached to the core. This force can, however, be neglected if the time scale associated with the core spin-down is larger than the (viscous) time scale during which the angular momentum is redistributed. In that case, the flow can reach a quasi-steady state. In the present work, for simplicity, we shall use this assumption ($\dot {\varOmega _c}=0$). Astrophysically, we justify its adoption by the Roche approximation, often used for massive stars, where the whole mass of the star is assumed to be in the core.

Another important time scale is that of the rise of the Stewartson layer. This layer appears as the result of the equatorial singularity of the inner Ekman boundary layer. It is well known that within an equatorial band of latitude $O(E^{1/5})$, the thickness of the Ekman boundary layer is $\delta _E = O(E^{2/5})$ (Roberts & Stewartson Reference Roberts and Stewartson1963). Hence, this singular equatorial viscous boundary layer is expected to be fully developed, and to initiate the development of the Stewartson layer, on the $O(\delta _E^2 /E)=O(E^{-1/5})$ time scale, that is much shorter than the $O(E^{-1})$ viscous time scale. Likewise, the central Stewartson layer of thickness $O(E^{1/3})$ is fully developed on the $O(E^{-1/3})$ time scale, again, much shorter than the viscous time scale. We therefore expect the Stewartson layer to start developing within a few revolutions of the inner shell and to be fully developed on a time scale that is much shorter than the viscous one on which we expect the shear flow inside it to evolve towards the steady state. We verify this in figure 3 showing the amplitude of the streamfunction $\psi _h/\psi _{h,st}$, taken at cylindrical radius $s=\eta$ and $z=1/2$, and normalised by its value at the steady state, as a function of the reduced time $Et$ for various Ekman numbers.

Figure 3.

Amplitude of the streamfunction taken at $s=\eta$ and $z=1/2$ normalised by its value at the steady state $\psi _h/\psi _{h,st}$, as a function of the reduced time $Et$ for various Ekman numbers. We find the Stewartson to appear within a few inner sphere revolutions and to evolve in a viscous time scale.

The analysis of the transient phase preceding the settling of a steady state thus underlines that the entire meridional flow evolves on a viscous time scale of order $E^{-1}$. In the astrophysical context characterised by extremely small Ekman numbers (typically $E<10^{-10}$), such a scaling implies that a steady state of the radiative envelope of massive stars is not reached during their lifetime.

3.2. The steady flow

Even if it is not reached during the lifetime of the star, the steady state is worth studying since it owns many simple features and its structure is very similar to that of the transient flow.

It is well known that the linear and stationary vorticity equation of an inviscid incompressible flow verifies Taylor–Proudman theorem (Proudman Reference Proudman1916; Taylor Reference Taylor1921), namely

(3.8)\begin{equation} ({\boldsymbol{e}}_z \boldsymbol{\cdot} \boldsymbol{\nabla} ) {\boldsymbol{u}} = \boldsymbol{0}. \end{equation}

For quasi-inviscid interior flows, that is for sufficiently small Ekman numbers, we thus expect a quasi-geostrophic solution for the velocity field. Figure 4 shows the normalised angular velocity as a function of the cylindrical radial coordinate $s$ for $E=10^{-7}$, $\eta =0.35$ and for various radial distances $r$. Figure 5 shows a 2-D view of the meridional circulation as well as the differential rotation in the rotating frame for the same model. The resulting flow is characterised by a nested Stewartson shear layer along the tangent cylinder $\mathcal {C}$ where the meridional circulation is essentially concentrated. This narrow shearing region separates two regions of quasi-geostrophic flow: the volume inside and outside $\mathcal {C}$. Inside $\mathcal {C}$, the no-slip conditions on the inner core impose an almost rigid rotation, while outside $\mathcal {C}$ a columnar differential rotation appears as a consequence of the surface stress.

Figure 4.

Normalised angular velocity as a function of the cylindrical radial coordinate $s$ for $E=10^{-7}$, $\eta =0.35$, $A=0.01$ and for various $r$.

Figure 5.

Meridional view of the streamfunction $\psi$ $(a)$ and of the differential rotation in the rotating frame of reference $\delta \varOmega = (\varOmega -\varOmega _c)/2\varOmega _c$ $(b)$ for $E=10^{-7}$, $\eta =0.35$ and $A=0.01$. Three regions can be distinguished: a Stewartson layer on the outer boundary of a tangent cylinder $\mathcal {C}$ of radius $\eta$ separating two regions of weak amplitude circulation. The region inside the cylinder $\mathcal {C}$ corotates with the inner sphere and the region outside $\mathcal {C}$ rotates differentially with a columnar profile.

3.2.1. Flow outside the core tangent cylinder $\mathcal {C}$

We now seek for an analytical solution for both the primary and secondary flows outside $\mathcal {C}$ ($s > \eta$). In the limit of small Ekman numbers, the divergenceless velocity satisfies the Taylor–Proudman theorem, we thus seek a geostrophic solution for both primary and secondary flows outside the tangent cylinder $\mathcal {C}$. Such a solution usually does not satisfy the viscous boundary conditions, however. We thus decompose the dynamical variables as

(3.9a,b)\begin{equation} {\boldsymbol{u}} = \bar{{\boldsymbol{u}}} + \tilde{{\boldsymbol{u}}}, \quad \textrm{and} \quad p = \bar{p}+ \tilde{p} , \end{equation}

where overlined variable correspond to the interior geostrophic fields, and tilded variables are their boundary layer corrections. Introducing the $O(1)$ stretched boundary layer coordinate $\zeta = (1-r)/\sqrt {E}$, and using the geostrophic equilibrium relation ${\boldsymbol {e}}_z \times \bar {{\boldsymbol {u}}} = - \boldsymbol {\nabla } \bar {p}$, the equation of motions, or Ekman boundary layer equation, reads

(3.10)\begin{equation} {\boldsymbol{e}}_z \times \tilde{{\boldsymbol{u}}} = - \frac{\partial \tilde{p}}{\partial \zeta} \boldsymbol{n} + \frac{\partial^2 \tilde{{\boldsymbol{u}}}}{\partial \zeta ^2}, \end{equation}

where $\boldsymbol {n} \equiv {\boldsymbol {e}}_r$ is the outwardly directed normal vector at the outer shell. Integrating the azimuthal component of (3.10) from $0$ to $\infty$ yields, at the lowest order,

(3.11)\begin{equation} \cos \theta \tilde{Q} = \left. - \frac{\partial \tilde{u}_{\phi}}{\partial \zeta}\right|_{\zeta=0} \quad \textrm{with}\ \tilde{Q}=\int_0^{\infty} \tilde{u}_{\theta}\,\textrm{d} \zeta. \end{equation}

Here, $\tilde {Q}$ is the outer Ekman layer $\theta$-directed volume flux.

Using the boundary condition for the prescribed tangential stress on the outer shell (2.5), (3.11) can be rewritten, at the lowest order

(3.12)\begin{equation} \cos \theta \tilde{Q} = \sqrt{E} \left[\tau(\theta) - \left. \frac{\partial}{\partial r} \left( \frac{F(s)}{r} \right)\right|_{r=1} \right], \end{equation}

where $F(s)= \bar {u}_{\phi }$ is the geostrophic solution for the azimuthal velocity in the interior. To ensure mass conservation, the divergence of $\tilde {{\boldsymbol {u}}}$ in the boundary layer generates further interior motion by establishing a weak amplitude secondary normal flow $\tilde {u}_r$ known as Ekman pumping. This radial velocity can be determined by integrating the mass conservation equation in the outer Ekman boundary layer

(3.13)\begin{equation} \left. \begin{aligned} \tilde{u}_r & = -\frac{\sqrt{E}}{\sin \theta} \frac{\partial \sin \theta \tilde{Q}}{\partial \theta}\\ & = -\frac{E}{\sin \theta} \frac{\partial}{\partial \theta} \frac{\sin \theta}{\cos \theta} \left[ \tau(\theta) - \frac{\partial}{\partial r}\left. \left(\frac{F(s)}{r}\right)\right|_{r=1} \right]. \end{aligned} \right\} \end{equation}

This equation shows that, near the outer shell, and in the Ekman layer, for the boundary condition on the tangential stress to be ensured, there must exist a $O(\sqrt {E})$ tangential flow that, in turn, generates a $O(E)$ circulation in the interior. We now express the geostrophic azimuthal velocity $F(s)$ as a function of the prescribed tangential surface stress $\tau (\theta )$ driving differential rotation. To do so, we ensure that the no-penetration boundary condition $u_r=0$ at $r=1$ is enforced, namely

(3.14)\begin{equation} s u_s + z u_z + \tilde{u}_r = 0 \quad \textrm{at}\ r=1, \end{equation}

where $u_s$ and $u_z$ are obtained from the azimuthal component of the momentum equation and from the mass conservation equation, respectively. Equation (3.14) finally yields a third-order ordinary differential equation for $F(s)$

(3.15)\begin{equation} -\frac{1-s^2}{s^2} \frac{\partial }{\partial s} s\left( \nabla^2 - \frac{1}{s^2}\right)F \!+\! s\left( \nabla^2 - \frac{1}{s^2}\right)F + \frac{\partial }{\partial s}\left(s^2 \frac{\partial }{\partial s}\frac{F}{s}\right) + sq\frac{\partial }{\partial s} \frac{F}{s} = -K(1+q(s)), \end{equation}

where $K=A \sqrt {3/4{\rm \pi} }$. The solutions of (3.15) can be expressed as an integral functional of the prescribed stress (Friedlander Reference Friedlander1976), that is

(3.16)\begin{equation} F(s)=-s \int \frac{\textrm{d}s}{s^3 \sqrt{1-s^2}} \int s\,\textrm{d}s \int \frac{\partial}{\partial s} \left( \frac{-K s^2 }{\sqrt{1-s^2}}\right)\,\textrm{d}s. \end{equation}

The solution for the geostrophic azimuthal velocity outside the tangent cylinder $\mathcal {C}$, avoiding singularities of the vorticity at $s = 1$, and accounting for the no-slip boundary condition at the inner shell finally reads

(3.17)\begin{equation} F(s)=-\frac{K}{3} \left( s \ln s - \frac{1}{s} + \alpha(\eta) s \right) , \end{equation}

with

(3.18)\begin{equation} \alpha(\eta)= \frac{1}{\eta^2} - \ln \eta, \end{equation}

and for which

(3.19a,b)\begin{equation} u_s = E\left(\varDelta - \frac{1}{s^2} \right)F(s)=- \frac{2KE}{3s} , \quad \textrm{and} \quad u_z=- \int_0^{z} \frac{1}{s}\frac{\partial s u_s}{\partial s}\,\textrm{d}z = 0. \end{equation}

We compare the angular velocity profiles from full numerical solutions to the analytical expression $\delta \varOmega = F(s)/s$, for various Ekman numbers and two inner core radii ($\eta =0.1$ and $0.35$) in figure 6. The analytical asymptotic solution reproduces rather well our numerical results when the Ekman number is below $10^{-9}$. Our boundary layer analysis also predicts the secondary meridional flow amplitude to be $O(E)$ outside the tangent cylinder $\mathcal {C}$. In figure 7 we show the streamfunction, defined as

(3.20a,b)\begin{equation} \frac{\partial \psi}{\partial r} = r \sin \theta u_{\theta}, \quad \frac{\partial \psi}{\partial \theta} = -r^2 \sin \theta u_r , \end{equation}

as a function of the cylindrical radial coordinate $s$ for various Ekman numbers. We note that $\psi$ indeed scales as $E$ outside the Stewartson layer, which appears as oscillations of $\psi$ near $s=\eta$. Furthermore, using (3.19a,b), the streamfunction associated with the geostrophic solution for the velocity field reads

(3.21)\begin{equation} \bar{\psi}(s,z)=-\frac{2KE}{3}z, \end{equation}

implying that streamlines are $z=\textrm {Cst}$ lines in a meridional plane. Hence, the outer Ekman boundary layer expels fluid in the $s$-direction, towards the Stewartson shear layer. The streamfunction $\psi$ as well as the corresponding streamlines, from (3.21) and from full numerical solutions are represented in figure 8, where we have masked the Stewartson layer. We see that the analytical and numerical solutions nicely match. Interestingly, the flow (3.19a,b) reconnects with the Stewartson layer and partly sources the upwelling flow in this region (see below § 3.2.3).

Figure 6.

Angular velocity in the rotating frame of reference as a function of the cylindrical radial coordinate $s$ for various Ekman numbers, $r=0.7$, $\eta =0.1$ $(a)$ and $\eta =0.35$ $(b)$, and $A=0.01$. The black dashed line corresponds to the analytical solution $\delta \varOmega = F(s)/s$ (see (3.17)).

Figure 7.

Value of $\psi E^{-1}$ along a meridian at $r=0.7$ as a function of the cylindrical radial coordinate $s$ for various Ekman numbers, $\eta =0.35$ and $A=0.01$. The predicted $O(E)$ scaling of the secondary flow amplitude is verified outside the Stewartson layer, which is located at $s=0.35$.

Figure 8.

Meridional view of the streamfunction and streamlines from full numerical simulations $(a)$, and from the boundary layer analysis $(b)$, for $E=10^{-8}$, $\eta =0.1$ and $A=0.01$. The tangent cylinder $\mathcal {C}$ as well as the adjacent reconnecting layer are masked.

The foregoing solution shows that the stress-driven spin-down flow is rather different from the spherical Taylor–Couette flow outside the tangent cylinder. For such a flow, the amplitude of the meridional circulation in this region is vanishing exponentially away from the Stewartson layer (Proudman Reference Proudman1956; Dormy, Cardin & Jault Reference Dormy, Cardin and Jault1998). In our case there is a residual flow, exactly perpendicular to the rotation axis, directed to it, and which scales as the Ekman number.

3.2.2. Flow inside the tangent cylinder $\mathcal {C}$

Let us now investigate the geostrophic flow inside the tangent cylinder $\mathcal {C}$. Writing $u_{\phi } = \ell /(r \sin \theta )$ where $\ell$ is the specific angular momentum, the azimuthal component of the momentum equation and vorticity equation respectively read (e.g. Goldstein Reference Goldstein1938, Reference Goldstein1965; Proudman Reference Proudman1956; Stewartson Reference Stewartson1966)

(3.22)\begin{equation} \left. \begin{gathered} \frac{\partial \psi}{\partial r}\cos \theta - \frac{1}{r}\frac{\partial \psi}{\partial \theta}\sin \theta = E D^2 \ell,\\ -\frac{\partial \ell}{\partial r}\cos \theta + \frac{1}{r}\frac{\partial \ell}{\partial \theta}\sin \theta = E D^4 \psi, \end{gathered} \right\} \end{equation}

where

(3.23)\begin{equation} D^2 = \frac{\partial^2}{\partial r^2}+ \frac{\sin \theta}{r^2} \frac{\partial}{\partial \theta} \left(\frac{1}{\sin \theta}\frac{\partial}{\partial \theta} \right). \end{equation}

The boundary conditions (2.5) in terms of streamfunction and specific angular momentum in turn read

(3.24a)\begin{gather} \psi = \frac{\partial \psi}{\partial r} = \ell = 0 \quad \textrm{at} \ r=\eta, \end{gather}

(3.24b)\begin{gather} \left. \begin{gathered} \psi=\frac{\partial^2 \psi}{\partial r^2} - 2 \frac{\partial \psi}{\partial r}=0,\\ \frac{\partial \ell}{\partial r}-2\ell=-K \sin^2 \theta, \end{gathered} \right\} \quad \textrm{at}\ r=1. \end{gather}

In the Ekman boundary layers, (3.22) may be rewritten

(3.25)\begin{equation} \left. \begin{gathered} \frac{\partial \psi}{\partial r}\cos \theta = E \frac{\partial^2 \ell}{\partial r^2},\\ -\frac{\partial \ell}{\partial r}\cos \theta = E \frac{\partial^4 \psi}{\partial r^4}. \end{gathered} \right\} \end{equation}

We integrate the first equation over $r$ to obtain the fourth-order differential equation on $\psi$

(3.26)\begin{equation} \frac{\partial^4 \psi}{\partial r^4}= - \frac{\cos^2 \theta}{E^2} (\psi-\bar{\psi}), \end{equation}

where $\bar {\psi }$ is a function of $\theta$ only that can be determined by the boundary conditions (Proudman Reference Proudman1956). Let us first focus on the outer Ekman layer localised near $r=1$. The solution of (3.26) satisfying the boundary conditions (3.24) as well as the condition $\psi (\xi \to \infty ) = \bar {\psi }$, where $\xi =(1-r)\sqrt {\cos \theta / 2E}$, reads

(3.27)\begin{equation} \psi=\bar{\psi} \left( 1- \textrm{e}^{-\xi} \cos \xi \right). \end{equation}

From the first equation of (3.25), we may further write

(3.28)\begin{equation} \ell - \bar{\ell}= \bar{\psi}\alpha \textrm{e}^{-\xi}\sin(\xi-{\rm \pi}/4), \end{equation}

where $\alpha =\sqrt {\cos \theta /(2E)} \gg 1$ in the asymptotic limit of small Ekman numbers. The boundary condition (3.24b) then implies

(3.29)\begin{equation} \left. \frac{\partial \ell}{\partial r}\right|_{r=1}- 2\ell(1)= -2\alpha^2\bar{\psi}- 2\bar{\ell} +2\alpha\bar{\psi} = -K\sin^2\theta, \end{equation}

where $\bar {\ell }$ is also a function of $\theta$ only that can be determined by the boundary conditions. Hence, the azimuthal and meridional components of the quasi-geostrophic velocity outside the Ekman and Stewartson shear layers, are not independent and satisfy the outer Ekman jump condition at the lowest order

(3.30)\begin{equation} \bar{\psi}= \frac{E}{\cos\theta} (K\sin^2\theta-2\bar{\ell}) . \end{equation}

In a similar way, the study of the inner Ekman layer with boundary conditions (3.24a) yields

(3.31)\begin{equation} \psi= \bar{\psi}(1-\textrm{e}^{-\xi'}(\cos \xi' + \sin \xi') ), \end{equation}

where $\xi '=(r-\eta )\alpha$, and the Ekman jump condition at the inner shell reads

(3.32)\begin{equation} \bar{\ell}=-2 \bar{\psi} \alpha . \end{equation}

Hence (3.32) implies $\bar {\psi }= O(\bar {\ell }\sqrt {E})$ and thus (3.30) implies $\bar {\psi }= O(E)$ and $\bar {\ell }= O(\sqrt {E})$. We may therefore rewrite (3.30) as

(3.33)\begin{equation} \bar{\psi}= \frac{E}{\cos\theta}K\sin^2\theta + O(E^{3/2}). \end{equation}

Finally, the geostrophic solution for the velocity inside the tangent cylinder $\mathcal {C}$ ($s < \eta$) reads, at the lowest order,

(3.34a,b)\begin{equation} \bar{\psi}(s) = \frac{EK s^2}{\sqrt{1-s^2}} \quad \textrm{and} \quad \delta\varOmega = \frac{\bar{\ell}(s)}{s^2} = - \frac{ \sqrt{2E}K}{\sqrt{1-s^2}} \left(1-\frac{s^2}{\eta^2}\right)^{1/4}, \end{equation}

and the meridional components of the geostrophic velocity at the lowest order

(3.35a,b)\begin{equation} u_s = \frac{1}{s} \frac{\partial \bar{\psi}}{\partial z} = 0 \quad \textrm{and} \quad u_z = -\frac{1}{s} \frac{\partial \bar{\psi}}{\partial s} = -\frac{E K}{\sqrt{1-s^2}} \left( 1 + \frac{1}{1-s^2}\right). \end{equation}

We compare these analytical $z$-directed velocity and angular velocity profiles with full numerical solutions in figure 9, for various $E$. We find our analytical expression to be in good agreement with the numerical solutions. Remarkably, the meridional flow inside the tangent cylinder is parallel to the rotation axis and directed towards the inner core, which is quite different from the outer-$\mathcal {C}$ meridional flow. Once inside the inner Ekman layer, the flow heads towards the equator where the Ekman layer thickens and eventually changes scale at the equatorial singularity (e.g. Roberts & Stewartson Reference Roberts and Stewartson1963; Stewartson Reference Stewartson1966; Hollerbach Reference Hollerbach1994; Dormy et al. Reference Dormy, Cardin and Jault1998; Marcotte, Dormy & Soward Reference Marcotte, Dormy and Soward2016). The fluid then returns to the outer Ekman layer following the Stewartson layer. Figure 10 shows the meridional view of two arbitrarily selected streamlines, one inside and one outside the tangent cylinder $\mathcal {C}$. It illustrates the different shape of the secondary flow in the two regions, as well as the reconnecting shear layer redirecting the $s$-direction flow outside $\mathcal {C}$ to the Stewartson layer where it flows parallel to the rotation axis towards the outer boundary.

Figure 9.

The $z$-directed velocity $u_z E^{-1}$ $(a)$ and angular velocity in the rotating frame of reference $(b)$ inside the tangent cylinder $\mathcal {C}$ as a function of the cylindrical radial coordinate $s$, for $E=10^{-7}$ and $E=10^{-8}$, $r=0.7$ and $\eta =0.35$. The black dashed line corresponds to the analytical geostrophic solutions (3.34a,b) and (3.35a,b).

Figure 10.

Meridional view of two arbitrarily selected streamlines, inside and outside the tangent cylinder $\mathcal {C}$. The arrows indicate the direction of propagation of the fluid along the streamlines.

3.2.3. Flow in the nested Stewartson shear layer

In figures 7 and 9 we see that the meridional flow strongly deviates from the analytical solution as one nears the Stewartson layer. From the numerical solution, the amplitude of the meridional circulation turns out to be much stronger in this layer than outside it. This is obvious if we compute the total meridional kinetic energy of the secondary flow

(3.36)\begin{equation} E_{k, tot}= \int V_m^2\,\textrm{d}V , \end{equation}

where $V_m^2= u_r^2 + u_{\theta }^2$. In figure 11(a) we show $E_{k, tot}$ as a function of the Ekman number. We find $E_{k, tot}$ to scale as $E$ in the asymptotic regime $E\to 0$, whereas the contribution of the flows outside the Stewartson layer remains $O(E^2)$. Hence, the Stewartson layer deserves some investigation.

Figure 11.

$(a)$ Total meridional kinetic energy of the stationary flow as a function of the Ekman number, for $A=0.01$ and $\eta =0.1$. The black dashed line corresponds to the fit in the asymptotic regime ($E \lesssim 10^{-9}$) and yields $E_{k, tot} \propto E^{1.074}$. $(b)$ $E_{k, tot} E^{-1.074}$ as a function of the reduced time $Et$ for various Ekman numbers, $A=0.01$ and $\eta =0.35$.

We first recall that in the Taylor–Couette flow the Stewartson layer is a nested shear layer that can be split into three layers of width $E^{2/7}$, $E^{1/3}$ and $E^{1/4}$ (Stewartson Reference Stewartson1966). The $E^{1/3}$-layer is the central layer, while the $E^{2/7}$-layer is on the inner side and the $E^{1/4}$ on the outer side. We now wish to find out if these three nested layers are present in our models. Let us start by noticing that, according to (3.13), the Ekman pumping may be written

(3.37)\begin{equation} \tilde{u}_r = - \sqrt{E}[\boldsymbol{\nabla} \times (\tilde{Q} {\boldsymbol{e}}_{\phi})]\boldsymbol{\cdot}{\boldsymbol{e}}_r. \end{equation}

Hence, introducing the adimensional Cartesian stretched coordinate of order $O(1)$ in a shear layer parallel to the rotation axis, outside the tangent cylinder, and of thickness $O(E^{\gamma })$, $\xi =E^{-\gamma }(s-\eta )$, we get

(3.38)\begin{equation} \tilde{u}_r =O\left[ \left(E \boldsymbol{\nabla} \times \left(E^{-\gamma}\frac{\partial F}{\partial \xi} {\boldsymbol{e}}_{\phi}\right) \right)\right] \boldsymbol{\cdot} {\boldsymbol{e}}_r = O(E^{1-2\gamma}) , \end{equation}

where $F=\bar {u}_{\phi }$. We conclude that, in such layers, the radial velocity resulting from the Ekman pumping/suction is of order $O(E^{1-2\gamma })$. The mass conservation equation, in turn, reads

(3.39)\begin{equation} E^{-\gamma} \frac{\partial u_s}{\partial \xi} + \frac{\partial u_z}{\partial z}= 0, \end{equation}

implying that in such a shear layer $u_z = O(E^{-\gamma }u_s)$, that is $u_z \gg u_s$ in the asymptotic regime of small Ekman numbers.

Simple manipulations of the momentum and vorticity equations allow us to eliminate the velocity $\boldsymbol {u}$ and lead to the following sixth-order partial differential equation for the pressure (Greenspan Reference Greenspan1968)

(3.40)\begin{equation} \frac{\partial^2 p}{\partial z^2} + E^{2-6\gamma} \frac{\partial^6 p}{\partial \xi^6} = 0, \end{equation}

representing a balance between Coriolis and viscous forces. The equality of these two terms thus requires $\gamma =1/3$, and indicates a shear layer of thickness $O(E^{1/3})$ as expected. Hence, in this central ageostrophic layer, we expect

(3.41a,b)\begin{equation} u_s = O(E^{2/3}), \quad \textrm{and} \quad u_z = O(E^{1/3}). \end{equation}

Let us now recall that the quasi-geostrophic $E^{1/4}$ layer is associated with the equilibrium between the Ekman pumping flow, and the internal friction of the $O(1)$ azimuthal flow (e.g. Stewartson Reference Stewartson1966; Barcilon Reference Barcilon1968). This equilibrium is, in fact, exactly that which is in place in the entire region outside the tangent cylinder $\mathcal {C}$, and the $E ^{1/4}$-layer is therefore not needed. Indeed, the azimuthal component of the momentum equation in the boundary layer of thickness $O(E^{\gamma })$ reads

(3.42)\begin{equation} u_s = E^{1-2\gamma} \frac{\partial^2 u_{\phi}}{\partial \xi^2}, \end{equation}

that is $u_s=O(E^{1-2\gamma })$ if we assume $u_{\phi }=O(1)$. In addition, (3.39) further implies $u_z=O(E^{1-3\gamma })$. However, the Ekman pumping still demands $u_z=O(E^{1-2\gamma })$ (see (3.37)). Hence, $\gamma =0$ and the equilibrium between the Ekman pumping flow and the internal friction of the $O(1)$ azimuthal flow does indeed occur in a $O(1)$ ‘layer’ that is the entire $s>\eta$ region. Note that if we had considered no-slip boundary conditions at the outer shell, that is the case of a spherical Taylor–Couette flow, the $z$-directed velocity resulting from the Ekman pumping would be of order $O(\sqrt {E} \boldsymbol {\nabla } \times \bar {u}_{\phi } {\boldsymbol {e}}_{\phi } )= O(E^{1/2 - \gamma })$. The matching of the $u_z$-flux and the Ekman pumping implies that $1/2 - \gamma =1-3\gamma$, and thus $\gamma =1/4$ as expected.

Let us finally consider the case of a quasi-geostrophic $z$-directed free shear layer of thickness $O(E^{\gamma '})$ located inside the tangent cylinder $\mathcal {C}$ and near $s= \eta$. It follows that in such a region, $\ell$ is a function of $s$ only. Note that in the case of a spherical Taylor–Couette flow, $\gamma '=2/7$ and $u_{\phi }(s=\eta ) = O(E^{1/28})$ (e.g. Stewartson Reference Stewartson1966; Marcotte et al. Reference Marcotte, Dormy and Soward2016). In this region, the azimuthal component of the momentum equation (3.22) integrated over $z$ reads

(3.43)\begin{equation} \psi = E z \frac{\partial^2 \ell }{\partial s^2} + f(s) , \end{equation}

where $f(s)$ is determined by applying the inner Ekman layer jump condition (3.32) at $z=0$ and $s - \eta \ll 1$ (Stewartson Reference Stewartson1966), yielding

(3.44)\begin{equation} f(s) \simeq -\sqrt{\frac{E}{2}} \left[2 \left(1-\frac{s}{\eta}\right) \right]^{-1/4} \ell. \end{equation}

In addition, at $z=\sqrt {1-\eta ^2}$, that is at $r=1$ and $s=\eta$, (3.30) reads

(3.45)\begin{equation} \psi = \frac{E K \eta^2}{\sqrt{1-\eta^2}}, \end{equation}

and (3.43) finally becomes

(3.46)\begin{equation} E^{1-2\gamma'} \frac{\sqrt{1-\eta^2}}{\eta^4} \frac{\textrm{d}^2 \ell}{\textrm{d}x^2}- E^{1/2 - \gamma'/4} (2x)^{-1/4} \ell =\frac{EK\eta^2}{\sqrt{1-\eta}}, \end{equation}

where we have introduced the stretched $O(1)$ shear layer coordinate $x=(1-s/\eta )/E^{\gamma '}$. At this point, it is interesting to determine the correct balancing in (3.46). We note that the only balance implying the existence of an asymptotically narrow shear layer is that of the first two terms, the third one being of higher order. This balance enforces $\gamma ' = 2/7$, that is a standard inner Stewartson shear layer of thickness $O(E^{2/7})$. Hence, the second-order ordinary differential equation for the specific angular momentum at the lowest order and in the $E^{2/7}$ shear layer reads

(3.47)\begin{equation} \frac{\textrm{d}^2 \ell}{\textrm{d}x^2}- \beta x^{-1/4}\ell=0, \end{equation}

where

(3.48)\begin{equation} \beta = \frac{\eta^2}{2^{3/4}\sqrt{1-\eta^2}}. \end{equation}

This equation can be transformed into a Bessel equation, and the solution asymptotically matching the solution (3.34a,b), that is decaying to zero as $x \to + \infty$, reads (see also Stewartson Reference Stewartson1966; Moore & Saffman Reference Moore and Saffman1968; Marcotte et al. Reference Marcotte, Dormy and Soward2016)

(3.49)\begin{equation} \ell(x) = \ell(0)C \sqrt{x}K_{4/7}(\tfrac{8}{7}\sqrt{\beta}x^{7/8}), \end{equation}

where $K_{\nu }(z)$ is the modified Bessel function of second kind and $C$ is a constant to be determined. We note that

(3.50)\begin{equation} \lim_{x\to 0} \sqrt{x}K_{4/7}\left(\frac{8}{7}\sqrt{\beta}x^{7/8}\right)= \frac{ \varGamma\left( 4/7\right)}{2\beta^{2/7}} \left(\frac{7}{4}\right)^{4/7} + O(x), \end{equation}

and thus, at the lowest order

(3.51)\begin{equation} C= \frac{2\beta^{2/7} }{\varGamma\left( 4/7\right)} \left(\frac{4}{7}\right)^{4/7} . \end{equation}

Finally, $\textrm {d} \ell /\textrm {d}s$ is rendered continuous across $s= \eta$ in the two quasi-geostrophic regions, that is

(3.52)\begin{equation} \left. \frac{\textrm{d} \ell}{\textrm{d}s}\right|_{s=\eta^-} = \left. \frac{\textrm{d}\ell}{\textrm{d}s}\right|_{s=\eta^+}, \end{equation}

which yields the specific angular momentum at $s=\eta$

(3.53)\begin{equation} \ell(0) = -\frac{K\eta^2 E^{2/7}}{3\beta^{4/7}} \left( \frac{2}{\eta^2} +1 \right) \left(\frac{7}{4}\right)^{1/7} \frac{\varGamma(4/7)}{\varGamma(3/7)} = O(E^{2/7}). \end{equation}

We verify the shear layer thickness as well as the analytical azimuthal velocity profile (3.49) comparing it with full numerical solutions in figure 12(a). The lowest-order solution for the azimuthal velocity in the entire $s \leq \eta$ domain is shown in figure 12(b). Finally, using (3.43), the corresponding streamfunction reads

(3.54)\begin{equation} \psi(x)= -\frac{2K\eta^2 E^{5/7} }{3\beta^{2/7}\varGamma(3/7)} \left( \frac{2}{\eta^2} +1 \right) \left(\frac{4}{7}\right)^{3/7} x^{1/4} K_{4/7} \left( \frac{8}{7}\sqrt{\beta}x^{7/8}\right) \left[ \frac{\beta z}{\eta^2} - \frac{1}{2^{3/4}} \right] , \end{equation}

and thus

(3.55)\begin{equation} \left. \begin{gathered} u_s =-\frac{B_s E^{5/7}}{1-xE^{2/7}} x^{1/4} K_{4/7} \left( \frac{8}{7}\sqrt{\beta}x^{7/8}\right),\\ u_z = \frac{B_z E^{3/7}}{x^{3/4}(1-xE^{2/7})}\left(K_{4/7}\left( \frac{8}{7}\sqrt{\beta}x^{7/8}\right) + 4 \sqrt{\beta}x^{7/8}K_{3/7}\left( \frac{8}{7}\sqrt{\beta}x^{7/8}\right)\right), \end{gathered} \right\} \end{equation}

in the $E^{2/7}$-layer, where

(3.56)\begin{equation} B_s = \frac{2K \beta^{5/7} }{3 \eta \varGamma(3/7)} \left( \frac{2}{\eta^2} +1 \right) \left(\frac{4}{7}\right)^{3/7} \end{equation}

and

(3.57)\begin{equation} B_z = \frac{K\left( 2 + \eta^2 \right) }{6 \eta^2 \beta^{2/7}\varGamma(3/7)} \left(\frac{4}{7}\right)^{3/7} \left( \frac{\beta z}{\eta^2} - \frac{1}{2^{3/4}} \right). \end{equation}

Figure 12.

$(a)$ Value of $u_{\phi } E^{-2/7}$ as a function of the stretched cylindrical radial coordinate $(s-\eta )E^{-2/7}$ for various Ekman numbers, $\eta =0.35$, $z=0.7$ and $A=0.01$. The black dashed line corresponds to the analytical geostrophic solution (3.49). $(b)$ Angular velocity in the rotating frame of reference as a function of the cylindrical radial coordinate $s$. The dashed lines correspond to the analytical solutions (3.34a,b) and (3.49).

We verify the validity of the solution (3.54) in the $E^{2/7}$-layer away from the $E^{1/3}$-layer in figure 13. We note that $\psi$ and its derivatives remain discontinuous across $s=\eta$, hence, contrary to the differential rotation (3.49), the solution (3.54) is not valid in the close vicinity of $s=\eta$, and it is the ageostrophic $E^{1/3}$-layer that smooths them out. Furthermore, since $3/7>1/3$, the $z$-component of the velocity is expected to be maximum in the $E^{1/3}$ layer.

Figure 13.

Value of $\psi E^{-5/7}$ as a function of the cylindrical radial coordinate $s$ for various Ekman numbers, $\eta =0.35$, $z=0.7$ and $A=0.01$. The dashed lines correspond to the analytical geostrophic solutions (3.54).

We further verify the scalings for the thickness of the two Stewartson shear layers, as well as that of the amplitude of their associated $s$- and $z$-directed velocity numerically, in figures 14 and 15. Although the amplitude of the $z$-directed velocity is, as expected, maximum in the central Stewartson layer, we find a slight discrepancy with the expected power laws for the velocity near $s=\eta$. This discrepancy is, in fact, not surprising as these two Stewartson shear layers are nested. The velocity in the region where the two coexist, that is the $E^{1/3}$-layer, is thus expected to follow both (3.41a,b) and (3.55). In practice, we measure

(3.58a,b)\begin{equation} u_s \simeq O(E^{0.7}), \quad \textrm{and} \quad u_z \simeq O(E^{0.37}), \end{equation}

which lies in between the velocity scalings in the two layers. For simplicity's sake, we may consider the amplitude scaling of meridional velocity in the $E^{1/3}$-layer to be that of (3.41a,b), that is $u_s=O(E^{2/3})$ and $u_z=O(E^{1/3})$.

Figure 14.

Value of $u_s E^{-2/3}$ $(a)$ and $u_z E^{-1/3}$ $(b)$ as a function of the stretched cylindrical radial coordinate $(s-\eta )E^{-1/3}$ for various Ekman numbers, $\eta =0.35$, $z=0.7$ and $A=0.01$.

Figure 15.

Value of $u_s E^{-5/7}$ $(a)$ and $u_z E^{-3/7}$ $(b)$ as a function of the stretched cylindrical radial coordinate $(s-\eta )E^{-2/7}$ for various Ekman numbers, $\eta =0.35$, $z=0.7$ and $A=0.01$. The black dashed lines correspond to the solution (3.55).

Ultimately, we find our results to be consistent with the scaling of the kinetic energy $E_{k, tot}\propto E$, which suggests that the main contribution comes from a $E^{1/3}$-velocity field spread over a free shear layer of width $E^{1/3}$. This $E^{1/3}$ amplitude scaling of the velocity field in the Stewartson layer is actually verified throughout the transient, as can be seen in figure 11(b). This implies that despite a steady state seemingly not being reached in the lifetime of massive stars, the time relevant for the advective transport of angular momentum and of chemicals within their radiative envelope scales as $E^{-2/3}$, since the whole meridional velocity grows on a viscous time scale (e.g. figure 3).

4.1. Description

Stable stratification introduces a restoring buoyancy force that acts to inhibit vertical motions and may partially or entirely eliminate the control of the flow exercised by the Ekman layers. Indeed, the tendency for viscous layers to generate secondary circulation is counteracted by the stable stratification. Thus, we may wonder to what extent the flow inside the stably stratified radiative envelope of a massive star losing angular momentum can remain properly modelled with a geostrophic solution. To investigate this question we introduce a radially stable stratification and use the Boussinesq approximation. We still neglect the centrifugal acceleration. The linearised and dimensionless vorticity, mass conservation and heat equations read

(4.1)\begin{equation} \left. \begin{gathered} \frac{\partial}{\partial t} \left(\boldsymbol{\nabla} \times {\boldsymbol{u}} \right)+\boldsymbol{\nabla} \times ({\boldsymbol{e}}_z \times {\boldsymbol{u}} - \delta T {\boldsymbol{r}}) = E \Delta \boldsymbol{\nabla} \times {\boldsymbol{u}},\\ \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}} = 0,\\ \frac{\partial \delta T}{\partial t} + Pr S r u_r = E \Delta \delta T, \end{gathered} \right\} \end{equation}

where we have used $\Delta T (2\varOmega _c)^2/ \mathcal {N}^2$ as the temperature scale; $\Delta T$ is the temperature difference between outer and inner shells, and $\mathcal {N}$ is the Brunt–Väisälä frequency defined as

(4.2)\begin{equation} \mathcal{N}^2= \frac{\alpha g \Delta T}{R}, \end{equation}

where $\alpha$ is the coefficient of thermal expansion and $\boldsymbol {g}= -g {\boldsymbol {r}}$ is the gravitational field. We have assumed the equilibrium temperature gradient to be produced by a uniform distribution of heat source (Chandrasekhar Reference Chandrasekhar1961), which implies that $\boldsymbol {\nabla } T_{eq, *} = (\Delta T / R^2) r {\boldsymbol {e}}_r$; $\delta T$ is the temperature perturbation from such thermal equilibrium, and $Pr= \nu /\kappa _T$ is the Prandtl number, which only enters the equations of motion combined with $S= \mathcal {N}^2/(2 \varOmega _c)^2$, as actually noted by Barcilon & Pedlosky (Reference Barcilon and Pedlosky1967). Here, $\kappa _T$ is the heat diffusivity that we have assumed to be constant. In what follows, we use the dimensionless parameter introduced by Garaud (Reference Garaud2002), namely

(4.3)\begin{equation} \lambda= Pr\frac{\mathcal{N}^2}{(2\varOmega_c)^2}, \end{equation}

as the parameter characterising the effective strength of thermal stratification.

We complete these equations with the boundary conditions (2.5) for the velocity field, and we impose $\delta T=0$ at the boundaries. The initial conditions are ${\boldsymbol {u}}=\boldsymbol {0}$, and $\delta T=0$. Projecting (4.1) on spherical harmonics yields the following system of equations for radial parts

(4.4)\begin{equation} \left. \begin{gathered} \frac{\partial w^l}{\partial t} - \left[ A_{l-1}^{l} r^{l-1}\frac{\partial}{\partial r}\left(\frac{r u^{l-1}}{r^{l-1}}\right) + A_{l+1}^{l} r^{-l-2}\frac{\partial}{\partial r}\left(r^{l+3} u^{l+1}\right) \right]= E \Delta_l w^l,\\ \frac{\partial}{\partial t}\left( \Delta_l r u^l \right) +B_{l-1}^{l} r^{l-1}\frac{\partial}{\partial r}\left(\frac{w^{l-1}}{r^{l-1}}\right) + \frac{B_{l+1}^{l}}{r^{l+2}} \frac{\partial}{\partial r}\left(r^{l+2} w^{l+1}\right) + \varLambda T^l = E\Delta_l \Delta_l r u^l,\\ \frac{\partial T^l}{\partial t}+ \lambda r u^l = E \Delta_l T^l, \end{gathered} \right\} \end{equation}

where

(4.5)\begin{equation} \delta T(r,\theta)= \sum_{l=0}^{+\infty} T^l(r) Y_l (\theta). \end{equation}

Figure 16 shows the differential rotation in the rotating frame of reference and the meridional circulation for $E=10^{-7}$ and various $\lambda$ once the steady state is settled. We see that the flow departs from the quasi-geostrophic equilibrium as $\lambda$ increases. At high $\lambda$ the angular velocity profile thus becomes shellular, namely it only depends on the radial distance.

Figure 16.

Meridional view of the streamfunction $\psi$ (a,c,e,g) and of the differential rotation in the rotating frame of reference $\delta \varOmega = (\varOmega -\varOmega _c)/ (2\varOmega _c$) (b,d,f,h) for $E=10^{-7}$, $\eta =0.35$, $A=0.01$ and various $\lambda$.

4.2. Horizontal boundary layers

In the case of a thermally homogeneous fluid, we have seen that horizontal boundary layers, and specifically Ekman layers play a crucial role on the interior dynamics of the flow. Hence, one may wonder how such layers impact the interior flow in a thermally stratified fluid. We first consider the condition of existence of such horizontal boundary layers and determine their thickness scaling with the relevant adimensional parameters.

Inside boundary layers, horizontal length scales are much larger than the transverse one. Radial derivatives are therefore expected to prevail. We thus write the simplified equations, assuming $\partial _r \gg \partial _{\theta }$ in horizontal boundary layers

(4.6)\begin{equation} \left. \begin{gathered} - u_{\phi} \sin \theta \simeq - \frac{\partial p}{\partial r} + \delta T r + E \frac{\partial^2 u_r}{\partial r^2},\\ - u_{\phi} \cos \theta \simeq - \frac{1}{r}\frac{\partial p}{\partial \theta} + E \frac{\partial^2 u_{\theta}}{\partial r^2},\\ u_{\theta} \cos \theta + u_r \sin \theta \simeq E \frac{\partial^2 u_{\phi}}{\partial r^2},\\ \frac{\partial u_r}{\partial r} + \frac{1}{r}\frac{\partial u_{\theta}}{\partial \theta} \simeq 0 ,\\ r u_r \simeq \frac{E}{\lambda} \frac{\partial^2 \delta T}{\partial r^2}. \end{gathered} \right\} \end{equation}

Combining these equations, only keeping the highest-order derivatives for each coefficient ($E^2$, $\lambda$ and 1) yields the general horizontal boundary layer equation (see also Raze, Lignières & Mimoun Reference Raze, Lignières and Mimoun2017)

(4.7)\begin{equation} E^2 r\frac{\partial^6 u_{\theta}}{\partial r^6} + \lambda \frac{\partial^2 u_{\theta}}{\partial \theta^2} - r \cos^2 \theta \frac{\partial^2 u_{\theta}}{\partial r^2} \simeq 0. \end{equation}

Introducing the $O(1)$ stretched radial coordinate $\zeta = (1-r)/\sqrt {E}$, (4.7) can be re-written, in the outer boundary layer and in the asymptotic regime of small Ekman numbers,

(4.8)\begin{equation} \frac{1}{E} \frac{\partial^6 u_{\theta}}{\partial \zeta^6} + \lambda \frac{\partial^2 u_{\theta}}{\partial \theta^2} - \frac{1}{E}\cos^2 \theta \frac{\partial^2 u_{\theta}}{\partial \zeta^2} \simeq 0. \end{equation}

We see that if $\lambda \ll 1/E$, we get the actual Ekman layer where the Coriolis force is balanced by the viscous shear. As remarked by Barcilon & Pedlosky (Reference Barcilon and Pedlosky1967), we note that, whenever Ekman layers are present, their structure is independent of the stratification strength because the forces in balance in such layers essentially involve horizontal motions. Similarly, if $\lambda \gg 1$, introducing the stretched radial coordinate $\gamma = (1-r)\sqrt {\lambda }$, (4.7) can be re-written, in the outer boundary layer,

(4.9)\begin{equation} E^2 \lambda^2 \frac{\partial^6 u_{\theta}}{\partial \gamma^6} + \frac{\partial^2 u_{\theta}}{\partial \theta^2} - \cos^2\theta \frac{\partial^2 u_{\theta}}{\partial \gamma^2} \simeq 0. \end{equation}

Hence, if $1 \ll \lambda ^2 \ll 1/E^2$, the viscous term drops out and buoyancy balances the Coriolis force. This is typical of a thermal boundary layer of width $\delta _{\lambda }= O(1/\sqrt {\lambda })$. We note that in the parameter regime $\lambda \ll 1$, only Ekman layers are present at the boundaries, while in the regime $1 \ll \lambda \ll 1/E$ both thermal and Ekman boundary layers coexist. In this latter case thermal boundary layers are always much thicker than Ekman layers.

4.3. Massive stars interior flows

From the foregoing discussion, it turns out that the dynamics of the flows may be quite different whether $\lambda \gg 1$ or $\lambda \ll 1$. Let us now place the case of rapidly rotating massive stars.

We note that, in the envelope of massive stars, the radiative viscosity and radiative heat diffusion largely dominate diffusion of collisional origin (Espinosa Lara & Rieutord Reference Espinosa Lara and Rieutord2013). These two quantities read

(4.10a,b)\begin{equation} \nu_{rad} = \frac{4aT^4}{15 c \kappa_R}, \quad \kappa_{rad} = \frac{4 a c T^3 }{3 \kappa_R c_p}, \end{equation}

where $T$ is the temperature, $\kappa _R$ is the Rosseland mean opacity, $a = 4 \sigma /c$ is the radiation density constant with $\sigma$ the Stefan–Boltzmann constant, $c$ the speed of light and $c_p$ is the specific heat capacity at constant pressure. In such a radiation-dominated system, the radiative Prandtl number reads

(4.11)\begin{equation} Pr_{rad} = \frac{3}{10}\frac{c_s^2}{c^2}, \end{equation}

where $c_s$ is the adiabatic sound speed, showing that, naturally, $Pr_{rad}\ll 1$.

However, the differential rotation of the radiative envelope is expected to drive some small-scale turbulence. Zahn (Reference Zahn1992) proposes that the turbulent vertical kinematic viscosity associated with marginal shear stability reads as

(4.12)\begin{equation} \nu_t = \frac{Ri_c \kappa_T}{3} \left( \frac{s}{\mathcal{N}} \frac{\textrm{d} \varOmega}{\textrm{d}s}\right)^2, \end{equation}

where $Ri_c \simeq 1/4$ is the critical Richardson number. In figure 17, we plot the radiative and turbulent Prandtl numbers radial profiles at the equator of a 15 solar mass star as given by a 2-D ESTER model (Rieutord et al. Reference Rieutord, Espinosa Lara and Putigny2016) for various rotation rates defined as the ratio $\omega$ between the equatorial angular velocity and the equatorial Keplerian angular velocity. The associated radial profiles of $\lambda$ are also shown. Although the turbulent Prandtl number is a few orders of magnitude larger than the radiative one in fast rotating stars, turbulent $\lambda$ is just one order of magnitude larger.

Figure 17.

The $Pr$ and $\lambda$ radial profiles at the equator measured from ESTER 2-D stellar models of a $15\ M_{\odot }$ star, for various angular velocity ratios $\omega$. The radiative Prandtl number $Pr_{rad}$ and the corresponding $\lambda$-parameter are represented in full lines, and their turbulent counterparts are represented in dashed lines.

We find that the turbulent $\lambda$-parameter is roughly independent of $\omega$, and never exceeds $10^{-3}$ in the considered models. The radiative $\lambda$-parameter, on the other hand, increases for decreasing $\omega$ to the point it may locally exceed $10^{-4}$ when $\omega \lesssim 0.1$. As massive stars are often considered to be fast rotators, we conclude that the thermal stratification regime of their radiative envelope corresponds to $\lambda \ll 1$.

4.4. Asymptotically weak thermal stratification regime

We now study the case of weak temperature stratification, relevant to the radiative envelope of rotating massive stars. In other contexts the $\lambda \ll 1$ regime may also be reached if the Brunt–Väisälä frequency is just small. We have seen that in this regime, only Ekman layers characterised by an $O(E)$ radial velocity at their edge, are present at the boundaries. We thus assume $\lambda \ll 1$ for the weak stratification, but also suppose that $E\ll \lambda$ as expected in stars.

4.4.1. The steady flow

We first focus on the steady flow. The radial component of the steady momentum equation reads

(4.13)\begin{equation} - u_{\phi} \sin \theta = -\frac{\partial p}{\partial r} + \delta T r + E \left[ \Delta u_r -\frac{2}{r^2}\left( u_r + \frac{1}{\sin \theta} \frac{\partial \sin \theta u_{\theta}}{\partial \theta}\right) \right], \end{equation}

where differential rotation is driven by the $O(1)$ surface stress, that is $u_{\phi }=O(1)$. In the asymptotic regime of small Ekman number, the temperature deviation from equilibrium is therefore, at most, $O(1)$. If that is the case, namely if $\delta T = O(1)$, the heat equation (4.1) implies that the interior radial velocity $u_r$ is $O(E/\lambda )$, namely larger than the $O(E)$ Ekman pumping (3.13), in the considered stratification regime. Because of the boundary conditions imposed by Ekman layers, an $O(E/\lambda )$ radial velocity must vanish at $r=\eta$ and $r=1$. Furthermore, the radial component of the vorticity equation reads

(4.14)\begin{equation} \frac{\partial u_r}{\partial z} = O(E), \end{equation}

implying that $u_r$ is $z$-independent up to $O(E)$ corrections (see Barcilon & Pedlosky Reference Barcilon and Pedlosky1967). Hence, a consistent solution for the interior flow is $u_r = O(E)$ and consequently $\delta T = O(\lambda )$, which is completed by the classical Ekman boundary layer with $\tilde {u}_r=O(E)$ and $\tilde {u}_{\theta }=O(\sqrt {E})$.

Let us now write the variables of the interior flow, as a power expansion of the small parameter $\lambda$, in the asymptotic regime of small Ekman numbers. They read

(4.15)\begin{equation} \left. \begin{gathered} u_r = E(u_{r,0} + \lambda u_{r,1} + \ldots),\\ u_{\phi} = u_{\phi,0} + \lambda u_{\phi,1} + \ldots,\\ \delta T = \lambda \delta T_1 + \lambda^2 \delta T_2 + \ldots,\\ p = p_0 + \lambda p_1 + \ldots. \end{gathered} \right\} \end{equation}

Injecting (4.15) into the momentum equation then yields, at zeroth and first orders

(4.16a,b)\begin{equation} u_{\phi,0} = \frac{\partial p_0}{\partial s}, \quad \textrm{and} \quad 0=\frac{\partial p_0}{\partial z}, \end{equation}

and

(4.17a,b)\begin{equation} u_{\phi,1} = \frac{\partial p_1}{\partial s} - s \delta T_1, \quad \textrm{and} \quad z \delta T_1=\frac{\partial p_1}{\partial z}. \end{equation}

Hence, at the lowest order, and identically to the homogeneous case, the Taylor–Proudman theorem is satisfied and the $O(1)$ geostrophic flow dynamics is consequently entirely controlled by the Ekman pumping. This is indeed verified in figure 16. On the other hand, the $O(\lambda )$ flow is affected by temperature deviation from equilibrium produced by the $O(E)$ radial velocity. We conclude that the thermal stratification regime in the radiative envelope of rotating massive stars has a negligible influence on the primary and secondary quasi-geostrophic and stationary flows.

4.4.2. The transient flow

As in the unstratified case, we now wish to determine the scaling of the transient time with the relevant non-dimensional parameters, that is, the Ekman number and the $\lambda$-parameter, in the $E\ll \lambda \ll 1$ limit. We write the radial component of the vorticity equation

(4.18)\begin{equation} \frac{\partial}{\partial t}\left( \frac{1}{r \sin \theta}\frac{\partial}{\partial \theta} \sin \theta u_{\phi}\right) -\frac{\partial u_r}{\partial z} - \frac{\sin \theta}{r} u_{\theta} = \frac{E}{r \sin \theta}\frac{\partial}{\partial \theta}\sin \theta \left( \nabla^2 -\frac{1}{r^2 \sin^2 \theta}\right) u_{\phi}, \end{equation}

as well as the $z$-derivative of the heat equation

(4.19)\begin{equation} \frac{\partial}{\partial t} \left( \frac{\cos \theta}{\lambda} \frac{\partial }{\partial z} \frac{\delta T}{ z}\right)+ \frac{\partial u_r}{\partial z} = \cos \theta \frac{E}{\lambda} \frac{\partial}{\partial z}\frac{1}{z}\Delta \delta T + \frac{\sin^2 \theta}{r \cos \theta} u_r. \end{equation}

Adding (4.18) and (4.19), and using the $\phi$-component of the momentum equation, we get

(4.20)\begin{align} &\frac{\partial}{\partial t} \left( \frac{z}{\lambda} \frac{\partial }{\partial z} \frac{ \delta T}{z} + \textrm{cotan}\theta \frac{\partial}{\partial \theta}\tan \theta u_{\phi} \right) \nonumber\\ &\quad = E \textrm{cotan} \theta \frac{\partial }{\partial \theta} \tan \theta \left( \nabla^2 -\frac{1}{r^2 \sin^2 \theta}\right) u_{\phi} + \frac{E}{\lambda}z \frac{\partial}{\partial z}\frac{1}{z}\Delta \delta T. \end{align}

Hence, if $\delta T$ remains $O(\lambda )$ during the transient, we expect the quasi-geostrophic steady state to be reached on the $O(E^{-1})$ viscous time scale in the limit $\lambda \ll 1$.

We verify this conclusion with our numerical solution and measure the transient time scale $\tau _t$ as the time required for the relative difference between the torques at the boundaries to be less than $0.01\,\%$. We show the scaled values $E\tau _t$ as a function of $\lambda$ for $E=10^{-6}$ in figure 18. We find that indeed, as for the thermally homogeneous case, the transient time scale is $O(E^{-1})$ in the weakly stratified limit.

Figure 18.

The scaled transient time scale $E\tau _t$ given by the numerical solution as a function of $\lambda$ for $E=10^{-6}$. The steady state is reached on a $O(E^{-1})$ time scale, both in the strongly and weakly stratified limits.

The thermal stratification regime in the interior of slowly rotating stars may, on the other hand, lie in the $\lambda \gg 1$ limit. In this regime, it is well known that the Eddington–Sweet time scale associated with the angular momentum redistribution by meridional circulation is so long that it is unlikely that the system would ever relax to a steady state. Hence, the study of the time-dependent interior dynamics of slowly rotating stars, which may largely depend on initial conditions and on the damping rate of the baroclinic modes, is outside the scope of this work. However, for the sake of completeness and to appreciate some effects of a strong thermal stratification in the stress-driven barotropic spin-down flow, we give a short account of the $\lambda \gg 1$ regime in appendix B.

5.1. The background

The next step towards a realistic model of the radiative envelope of a massive star is to include the strong density variations of the fluid between the convective core and the stellar surface. Typically density varies by a factor of $10^9$ in the envelope of a star of 15 solar masses. To take this density distribution into account in a simple way we assume that the gas can be described by a polytropic equation of state as usually done in stellar physics (e.g. Maeder Reference Maeder2009). Hence the hydrostatic background state verifies

(5.1a,b)\begin{equation} -\frac{\textrm{d}p_{*}}{\textrm{d}r_{*}}-\frac{GM}{r_*^2}\rho_{*}=0 \quad \textrm{and} \quad p_{*}=\kappa\rho_{*}^{\varGamma}, \end{equation}

where $p_{*}$, $\rho _*$ and $r_*$ are the dimensional pressure, density and radial coordinate, respectively. As previously, we also use the Roche approximation and assume that the core gathers all the mass $M$ of the star. In (5.1a,b) $G$ is the gravitational constant, $\kappa$ is a constant related to the thermal conditions at the core boundary and $\varGamma =1+1/n$, where $n$ is the so-called polytropic index. In the following we shall set $n=3$, which is typical for radiative envelopes. We note that for such an index the fluid is stably stratified (e.g. Dintrans & Rieutord Reference Dintrans and Rieutord2001; Rieutord & Dintrans Reference Rieutord and Dintrans2002, for instance), but as concluded from the previous section the small value of the Prandtl number allows us to neglect buoyancy effects. The following results therefore assume no buoyancy force.

Equation (5.1a,b) can be easily solved and gives the density profile

(5.2)\begin{equation} \rho_{*}(r_{*})=\rho_c(1+B j(r_{*}))^n, \end{equation}

where

(5.3a,b)\begin{equation} B= \frac{\rho_s^{1/n}-1}{1/\eta - 1}, \quad j(r_{*})=\frac{1}{\eta} - \frac{R}{r_{*}}. \end{equation}

Here $\rho _c = \rho _{*}(\eta R)$ is the density at the core boundary and $\rho _s=\rho _{*}(R)/\rho _c$ is the adimensional surface density. Note that such a background is also used in numerical simulations of convection in stellar envelope (Raynaud et al. Reference Raynaud, Rieutord, Petitdemange, Gastine and Putigny2018) or planetary atmospheres (Gastine & Wicht Reference Gastine and Wicht2012). Hence, in this section, we solve (2.1) with the background polytropic density profile (5.2), the boundary conditions (2.5) and with the initial condition ${\boldsymbol {u}}= \boldsymbol {0}$.

5.2. The transient phase

Again, we first focus on the transient phase preceding the settling of a steady state. In particular, we aim at determining whether the scaling of the governing time scales are modified by density variations of the background. We measure the transient time as the time for which relative difference between the torque exerted by the fluid on the stationary inner sphere and the torque applied on the outer sphere, is less than $0.01 \%$. We monitor the evolution of the relative difference between the torques $\Delta \varGamma / \varGamma (1)$ given by (3.6), for various surface density ratios. We show the resulting steady-state times in figure 19 for the polytropic indexes $n=3$ and $n=3/2$. The latter index corresponds to an isentropic monatomic gas, hence to a neutral thermal stratification. We find the steady-state time to depend on the density ratio between the core and the surface. For the polytropic index $n=3$, relevant to the envelope of massive stars, we find $E \tau _t \propto \rho _s^{0.3}$ for $10^{-3} \leq \rho _s$, and $E \tau _t \propto \rho _s^{0.09}$ for $\rho _s < 10^{-4}$. Of course, these power laws quantitatively depend on the chosen density profile, as can be seen with the $n=3/2$ case. Unlike thermal stratification, density stratification thus (mildly) influences the time scale required to reach a steady state in a stellar envelope.

Figure 19.

Scaled steady-state time scale $E\tau _t$ as a function of the surface density $\rho _s$, for $E=10^{-6}$ models with various values of $\rho _s$, $n=3/2$ and $n=3$. The dashed and dotted lines correspond to the steady-state time scalings with $\rho _s$.

5.3. The steady flow

We now focus on the steady flow. We may observe that if viscosity, nonlinearity and buoyancy are neglected then the steady flow obeys a Taylor–Proudman theorem applied to the momentum $\rho {\boldsymbol {u}}$ instead of the velocity, namely

(5.4)\begin{equation} ({\boldsymbol{e}}_z \boldsymbol{\cdot} \boldsymbol{\nabla} ) \rho {\boldsymbol{u}} = \boldsymbol{0} . \end{equation}

Unfortunately, we cannot reiterate the boundary layer analysis of the incompressible case since the viscous force (2.3) now includes terms that depend on $r=\sqrt {s^2+z^2}$ and make the partial differential equation not separable. We may, however, observe that since the density variations do not add any new length scale, the viscous balance in the shear layers remains similar and we can still expect the presence of a Stewartson layer along the tangent cylinder.

We now revert to numerical solutions to make progress. We thus solve (2.1) with boundary condition (2.5). As already mentioned, we neglect the effects of buoyancy. We first focus on the differential rotation and the meridional circulation. It is convenient to introduce the streamfunction $\chi$ associated with the meridional momentum, namely

(5.5a,b)\begin{equation} \frac{\partial \chi}{\partial r} = r \sin \theta \rho u_{\theta}, \quad \frac{\partial \chi}{\partial \theta} = -r^2 \sin \theta \rho u_r. \end{equation}

Figure 20 shows the differential rotation and this new streamfunction in the rotating frame of reference for $E=10^{-7}$.

Figure 20.

Meridional view of the streamfunction $\chi =\rho \psi$ (a) and of the differential rotation in the rotating frame of reference $\delta \varOmega = (\varOmega -\varOmega _c)/2\varOmega _c$ (b) for $E=10^{-7}$, $\eta =0.35$, $\rho _s=10^{-4}$ and $A=0.01$.

We note that as for the incompressible flow, the amplitude of the secondary flow is dominated by the Stewartson layer located at the edge of the tangent cylinder $\mathcal {C}$. The differential rotation, on the other hand, and as expected, is no longer cylindrical and very $z$-dependent, with a maximum value that is close to the outer shell and the tangent cylinder. Figure 21(a) shows the differential rotation as a function of the cylindrical radial coordinate for $E=10^{-7}$ at different radii $r$, illustrating the $z$-dependence of $u_{\varphi }$, while figure 21(b) shows that $q_{\varphi }=\rho (r)u_{\varphi }$ verifies the Taylor–Proudman constraint.

Figure 21.

Normalised angular velocity (a) and normalised angular velocity multiplied by the density $\rho$ (b), and as a function of the cylindrical radial coordinate $s$ for $E=10^{-7}$, $\eta =0.35$, $\rho _s=10^{-4}$, $A=0.01$ and various radii.

As far as meridional circulation is concerned, this is still an ${O}(E)$ flow outside the Stewartson layer as shown by figures 22(a) and 23(a). We note that outside the tangent cylinder, streamlines are no longer straight lines, even for the ${\boldsymbol {q}}$ momentum field. In the inner cylinder and near the rotation axis we still get streamlines parallel to the rotation axis for the momentum ${\boldsymbol {q}}$. This latter feature comes from the fact that the Ekman layer has the same structure as in the constant density case and drives an Ekman circulation which has a unique component along the $z$-axis. From, (3.35a,b) we deduce that

(5.6)\begin{equation} q_z = -\frac{E K\rho_s}{\sqrt{1-s^2}} \left( 1 + \frac{1}{1-s^2}\right) + O(E^2) \end{equation}

in this region.

Figure 22.

$(a)$ Value of $\chi E^{-1}$ as a function of the cylindrical radial coordinate $s$ for various Ekman numbers, $\eta =0.35$, $r=0.7$, $\rho _s=10^{-4}$ and $A=0.01$. The predicted $O(E)$ scaling of the secondary flow is verified away from the Stewartson nested layers. $(b)$ Meridional view of the streamfunction $\chi E^{-1}$ outside $\mathcal {C}$, for the $E=10^{-9}$ model.

Figure 23.

$(a)$ Value of $u_z E^{-1}$ as a function of the cylindrical radial coordinate $s$ for various Ekman numbers, $\eta =0.35$, $r=0.7$, $\rho _s=10^{-4}$ and $A=0.01$. $(b)$ Meridional view of the streamfunction $\chi$ inside $\mathcal {C}$ and away from the $O(E^{2/7})$ Stewartson layer, for the $E=10^{-9}$ model.

Finally, the Stewartson layer seems to conserve its structure if we focus on the momentum, as shown by figures 24 and 25. Hence, the $E^{1/3}$ scale is still the dominating scale.

Figure 24.

Value of $q_s E^{-2/3}$ (a) and $q_z E^{-1/3}$ (b) as a function of the stretched cylindrical radial coordinate $(s-\eta )E^{-1/3}$ for various Ekman numbers, $\eta =0.35$, $z=0.7$, $\rho _s=10^{-4}$ and $A=0.01$.

Figure 25.

Value of $q_s E^{-5/7}$ (a) and $q_z E^{-3/7}$ (b) as a function of the stretched cylindrical radial coordinate $(s-\eta )E^{-2/7}$ for various Ekman numbers, $\eta =0.35$, $z=0.7$, $\rho _s=10^{-4}$ and $A=0.01$.