2.1. Formulation of the problem

Consider incompressible fluid flow with kinematic viscosity $\nu$ between two co-axial cylinders of dimensional inner and outer radii $r_i^*$ and $r_o^*$ . We choose the cylinder gap, $d^*=r_o^*-r_i^*$ , as the length scale. The arrangement of the cylinders is hence uniquely determined by the radius ratio $\eta =r_i^*/r_o^*\leqslant 1$ . The diffusive time scale, $(d^*)^2/\nu$ , is the conventional choice for describing Taylor–Couette flow. Using the associated velocity scale, $\nu /d^*$ , and the pressure scale, $\rho ^* \nu ^2/(d^*)^2$ , where $\rho ^*$ is the fluid density, the dimensionless velocity $\mathbf {v}$ and pressure $p$ fields are governed by the Navier–Stokes equations

(2.1) \begin{eqnarray} \partial _{t}\mathbf {v} + (\mathbf {v}\cdot \nabla )\mathbf {v} = -\nabla p + \nabla ^{2}\mathbf {v} \end{eqnarray}

and the divergence-free condition $ \nabla \cdot \mathbf {v} = 0$ . Due to the use of the viscous scale, the Reynolds number does not appear in the diffusion terms. Instead, the normalised angular velocities of the cylinders define the two conventional Reynolds numbers. In non-dimensional cylindrical coordinates $(r,\theta ,z)$ , the no-slip boundary conditions at the cylinder walls are prescribed as

(2.2) \begin{equation} \mathbf {v}(r_i,\theta ,z)=(0,R_i,0),\quad \mathbf {v}(r_o,\theta ,z)=(0,R_o,0), \end{equation}

where $r_i= {\eta }/({1-\eta })$ and $r_o={1}/({1-\eta })$ . The Reynolds numbers and the dimensional angular velocities of the inner and outer cylinders, $\Omega _i^*$ and $\Omega _o^*$ , are related through

(2.3) \begin{equation} R_i=\frac {\Omega _i^* r_i^* d^*}{\nu },\qquad R_o=\frac {\Omega _o^* r_o^* d^*}{\nu }. \end{equation}

As is well known, the laminar circular Couette solution can be written as $(u,v,w,p)=(0,R_o\,r\Omega (r),0,R_o^2\int ^r r_*\Omega ^2(r_*) \textrm{d}r_*)$ , where the scaled angular velocity has the analytic form

(2.4) \begin{eqnarray} \Omega (r)=A+B/r^2 \end{eqnarray}

with the constants

(2.5) \begin{eqnarray} A\equiv \frac {1-\eta ^2 a}{1+\eta },\qquad B\equiv \frac {a-1}{1+\eta }r_i^2, \qquad a\equiv \frac {R_i}{\eta R_o}. \end{eqnarray}

The Rayleigh line corresponds to $a=\eta ^{-2}$ . The Rayleigh stable region is thus written as $a\in [0,\eta ^{-2}]$ , and it includes the solid-body rotation $a=1$ .

2.2. Narrow-gap limit

Now we derive the reduced problem valid at the narrow-gap limit $\eta \rightarrow 1$ . We introduce a small perturbation parameter $\delta =\sqrt {1-\eta }$ and expand

(2.6) \begin{eqnarray} (u,v,w,p)=(u^+, \delta ^{-1}v^+,w^+,p^+)+\cdots , \,\, (R_o,R_i)=\delta ^{-1}(Re,aRe)+\cdots . \end{eqnarray}

Changing the coordinates as $y=r-r_i$ , $x=\delta ^{-1}\theta$ and only retaining the leading-order terms, we get

(2.7a) \begin{eqnarray} (\partial _t+v^+\partial _x+u^+\partial _y+w^+\partial _z)u^+-v^{+2}=-p^+_y+(\partial _y^2+\partial _z^2) u^+, \end{eqnarray}

(2.7b) \begin{eqnarray} (\partial _t+v^+\partial _x+u^+\partial _y+w^+\partial _z)v^+=(\partial _y^2+\partial _z^2) v^+, \end{eqnarray}

(2.7c) \begin{eqnarray} (\partial _t+v^+\partial _x+u^+\partial _y+w^+\partial _z)w^+=-p^+_z+(\partial _y^2+\partial _z^2) w^+, \end{eqnarray}

(2.7d) \begin{eqnarray} \partial _x v^+ +\partial _y u^+ +\partial _z w^+=0, \end{eqnarray}

and the boundary conditions

(2.8a) \begin{eqnarray} (u^+,v^+,w^+)=(0,Re,0)\qquad \text {at}\qquad y=1, \end{eqnarray}

(2.8b) \begin{eqnarray} (u^+,v^+,w^+)=(0,aRe,0)\qquad \text {at}\qquad y=0. \end{eqnarray}

The limiting process here is essentially that used by Hall (Reference Hall1982, Reference Hall1983, Reference Hall1988) for the Görtler vortex problem, and the final set of equations is often referred to as the boundary region equations. The laminar circular Couette flow now becomes a linear profile

(2.9) \begin{eqnarray} {(u^+,v^+,w^+,p^+)=(0,Re \, v_b(y),0,Re^2p_b(y))},\qquad v_b(y)=(1 - a)y+a. \end{eqnarray}

The last term on the left-hand side of (2.7a ) is the Görtler term, representing the deviation from the plane Couette problem.

Stability of the basic flow to infinitesimal perturbations can be analysed by writing

(2.10) \begin{eqnarray} (u^+,v^+,w^+,p^+)&=&(0,Re \, v_b(y),0,Re^2p_b(y))\nonumber \\ &+&(\hat {u}(y),\hat {v}(y),\hat {w}(y),Re\, \hat {p}(y))\exp (i\alpha (x-Re \, ct)+ikz)+\text {c.c.}\qquad\; \end{eqnarray}

and linearising (2.7). Here, c.c. denotes the complex conjugate. To ease notation, we shall drop the hats after substitution. Eliminating the pressure and the axial velocity component, we arrive at the coupled equations

(2.11) \begin{eqnarray} \mathcal {L}v+{{v_b'}} u=0,\qquad \mathcal {L}( u''-k^2u)+2v_b k^2v=0, \end{eqnarray}

involving the linear operator $\mathcal {L}=i\alpha (v_b-c)-Re^{-1}(\partial _y^2-k^2)$ . The prime denotes differentiation with respect to $y$ . The above equations and the boundary conditions $u=u'=v=0$ imposed at $y=0$ and 1 constitute an eigenvalue problem for the complex growth rate $c=c_r+ic_i$ . From the symmetry of the system, it is sufficient to consider positive scaled azimuthal wavenumber $\alpha$ and axial wavenumber $k$ . The flow is unstable if $c_i\gt 0$ .

The following mathematical facts are immediately apparent.

1. (i) Squire’s theorem does not hold due to the presence of the Görtler term.

2. (ii) If $k=0$ , the stability problem is equivalent to that for plane Couette flow. Therefore, the well-known proof by Romanov (Reference Romanov1973) can be applied, showing that no instability exists.

3. (iii) The parameter range $a\in [0,1]$ corresponds to the Rayleigh stable case. In this regime, the flow is stable if $\alpha =0$ , i.e. the flow is axisymmetric. As remarked in § 1, the Rayleigh criterion assumes the flow to be inviscid. However, even if we include the viscous terms, we can still show the absence of instability, as the proof by Synge (Reference Synge1938) for Taylor–Couette flow is still valid for the narrow-gap limit equations.

It follows that for the Rayleigh stable case $a\in [0,1]$ , instability is impossible if either $\alpha$ or $k$ is zero.

In the rest of § 2, we focus on the stationary inner cylinder case $a=0$ for simplicity. The eigenvalue problem (2.11) can be numerically solved by the Chebyshev collocation method. With increasing $Re$ from zero, an instability emerges at $(Re,\alpha ,k)=(2.831\times 10^7,0.978,1.350)$ , as reported by Deguchi (Reference Deguchi2017); see figure 2. Like many other shear flows, the size of the instability region in the wavenumber plane expands as the Reynolds number increases further. However, in the current eigenvalue problem, the flow becomes unstable only within a very narrow region, making detection of this instability through numerical computations extremely difficult.

Figure 2. Linear instability found in the narrow-gap Taylor–Couette flow problem (2.11) for the stationary inner cylinder case $a=0$ .The magenta bullet indicate the wavenumber pair $(\alpha ,k)=(0.978,1.350)$ at which the instability first emerges at $Re=2.831\times 10^7$ . The cyan dashed and black solid curves are the neutral curves for $Re=2.9 \times 10^7$ and $10^8$ , respectively. The crosses are found by the large- $Re$ asymptotic result (2.16). The dash-dotted line is $k = (5\pi/8\sqrt{2})\alpha$ found by the analysis of the Kummer’s function. In the inset, $K = k-(5\pi/8\sqrt{2})\alpha$ is used to illustrate the structure of the unstable regions clearly.

Figure 3. (a) Neutral curve of the narrow-gap Taylor–Couette flow problem (2.11) for $\alpha =0.978$ , $a=0$ and (b) the corresponding phase speed $c$ . The full numerical stability results are shown by the black solid curves. The dashed curves are the large- $Re$ asymptotic results (2.28), (2.30) and (2.45). The blue dotted line in panel (a) is the leading-order approximation $k=1.442$ obtained by (2.16). The numerical resolution is verified using up to 3000 Chebyshev modes.

Figure 4. Comparison of the eigenfunctions between the linearised Navier–Stokes computations (lines) and asymptotic predictions (symbols). (a) Comparison of results without scaling the wall-normal coordinate. The solid and dashed curves are the real and imaginary parts of the lower branch neutral eigenfunction at $(Re,\alpha ,k,a)=(10^{10},0.978,1.434,0)$ . The points are the leading-order core solution $u_0$ with $(\alpha ,k_0)=(0.978,1.442)$ ; see (2.15). Note that $u_0$ attains its maximum value of $0.1523$ at $y=0.3698$ , and $u$ is normalised to have the same property. (b) Comparison in terms of the viscous wall layer scaling. The points are $U_0-Y$ , where $U_0$ is the leading-order wall-layer solution (2.22) with $(\alpha ,k_0,k_1,c_1)=(0.978,1.442,-14.9,2.31)$ and $Y$ is the stretched coordinate $Y=Re^{1/3}y$ . Note that the eigenfunction is neutral, and thus the real part of the asymptotic solution approaches $U_d$ as $Y\rightarrow \infty$ .

Figure 3(a) shows the neutral curve for $\alpha = 0.978$ , with the corresponding phase speed depicted in figure 3(b). The two branches of the neutral curve emanating from the aforementioned critical point exhibit large-Reynolds-number asymptotic characteristics somewhat similar to TS waves. Indeed, we shall see in §§ 2.3 and 2.4 that the asymptotic structure of the perturbation near the wall bears similarities to that studied by Tollmien (Reference Tollmien1929); Schlichting (Reference Schlichting1933) and Lin (Reference Lin1945, Reference Lin1955). For this reason, we adopt the terms upper and lower branches in figure 2 so that they preserve this analogy. The red dashed curves in the figures are the asymptotic results to be derived in §§ 2.3 and 2.4. For the lower branch, a reasonably accurate asymptotic approximation is obtained for a fairly wide range of $Re$ . In contrast, the upper branch exhibits a wiggle, and only beyond this point does the asymptotic result provide a good approximation. These features were also observed in neutral curves of TS waves in channel flows and boundary layer flows (see Drazin & Reid (Reference Drazin and Reid1981), for example).

The phase velocity $c$ provides important insight into the structure of the eigenfunctions because a thin vortex layer, known as the critical layer, typically forms around the region where the base flow speed matches $c$ . The magnitude of $c$ drops with increasing $Re$ , as seen in figure 3(b), suggesting that if there is a critical layer, it should be sitting around $y=0$ , where the base-flow velocity is small. The radial eigenfunction $u$ for the lower branch, shown in figure 4, indicates that there is indeed some structure near the inner cylinder, while the perturbation does not diminish in the core region away from the wall surface. The nature of the interaction between the core flow and the near-wall layer will naturally be revealed by the method of matched asymptotic expansion. Our theory will also demonstrate that constant wavelengths are necessary to maintain the core structure. Thus, as seen in figure 3(a), the wavenumber does not decrease with increasing $Re$ , unlike TS waves.

Figure 5. Sketch of large- $Re$ asymptotic structure of (a) the lower branch type and (b) the upper branch type. A Stokes layer of thickness $O(Re^{-1/2})$ exists near $y=1$ , but it is omitted.

2.3. Lower branch asymptotic analysis

Let us start by analysing the lower branch, which has a relatively simple asymptotic structure. We fix $\alpha$ at a positive constant throughout the analysis. When $Re$ is asymptotically large, the majority of the fluid layer is occupied by an inviscid core region, while near the inner cylinder ( $y=0$ ), a thin viscous wall layer of thickness $O(\epsilon )$ emerges (see figure 5 a). Since we have set $a=0$ , within the boundary layer, the size of $v_b$ is $O(\epsilon )$ . Therefore, the advective effect of $O(\epsilon )$ and the viscous effect of $O(Re^{-1}\epsilon ^{-2})$ in the operator $\mathcal {L}$ in (2.11) balance when $\epsilon =Re^{-1/3}$ . This wall layer thickness is typical of critical layers, where the phase speed of the wave matches the base flow speed.

Using the small parameter $\epsilon$ , we expand $k=k_0+\epsilon k_1+\cdots$ and $c=\epsilon c_1 +\cdots$ . The overall structure of the asymptotic problem is as follows. First, the leading-order inviscid core problem determines $k_0$ . Then, at the next order, the core flow is coupled with the wall layer problem to yield the dispersion relation linking $c_1$ and $k_1$ . In the neutral case, $c_1$ and $k_1$ can be solved explicitly.

In the inviscid core region, we expand the flow field as

(2.12) \begin{eqnarray} u=u_0(y)+\epsilon u_1(y)+\cdots ,\qquad v=v_0(y)+\epsilon v_1(y)+\cdots . \end{eqnarray}

Substituting (2.12) into (2.11), retaining only the leading-order terms and eliminating $v_0$ , we have the single second-order ordinary differential equation

(2.13) \begin{eqnarray} \mathscr {L} u_0 =0, \qquad \mathscr {L}\equiv \partial _y^2+k_0^2 \left (\frac {2}{\alpha ^2 y}-1 \right ). \end{eqnarray}

A similar manipulation at the next order yields the inhomogeneous version of the equation

(2.14) \begin{eqnarray} \mathscr {L} u_1 = - \left \{\frac {4 k_0^2c_1}{\alpha ^2 y^2} +2k_1 k_0 \left (\frac {2}{\alpha ^2y}-1 \right )\right \} u_0. \end{eqnarray}

From the no-penetration condition at the wall, we impose $u_0(1)=u_1(1)=u_0(0)=0$ . However, it will turn out from the boundary layer analysis that $u_1(0)$ cannot be zero, so we denote this displacement velocity as $U_d=u_1(0)$ .

The no-slip conditions at the outer cylinder ( $y=1$ ) are satisfied through a Stokes layer of thickness $O(Re^{-1/2})$ . However, this layer is passive and does not need to be analysed.

The leading-order core system is an eigenvalue problem for the eigenvalue $k_0^2$ . The solution can be written down using Kummer’s function of the first kind $M$ :

(2.15) \begin{eqnarray} u_0(y)=ye^{-k_0 y}\frac {M(\gamma , 2, 2k_0 y)}{M(\gamma , 2, 0)},\qquad \gamma \equiv 1-\frac {k_0}{\alpha ^2}. \end{eqnarray}

The boundary condition at $y=1$ is fulfilled only when

(2.16) \begin{eqnarray} M(\gamma , 2, 2k_0)=0. \end{eqnarray}

The crosses in figure 2 are obtained by looking for the smallest positive $k_0$ that satisfies (2.16) for each $\gamma$ . For a small $\alpha$ , the approximation

(2.17) \begin{eqnarray} k_0\approx \frac {5\pi }{8\sqrt {2}}\alpha \end{eqnarray}

is available, thanks to the property of the Kummer’s function (Abramowitz & Stegun Reference Abramowitz and Stegun1965). Both results are useful for estimating the wavenumber region where instability might occur. For $\alpha =0.978$ used in figure 3(a), the prediction by (2.16) gives $k_0\approx 1.442,$ which is shown by the blue dotted line in the figure. Also, the leading-order core solution (2.15) agrees very well with the numerically obtained neutral eigenfunction, as depicted in figure 4(a). Note that to simplify the subsequent analysis, we have normalised the eigenfunction imposing $u_0'(0)=1$ .

Within the viscous wall layer, we use the stretched variable $Y=y/\epsilon$ , and expand the wall-normal and streamwise velocity components as

(2.18) \begin{eqnarray} u=\epsilon U_0(Y)+\cdots , \qquad v=V_0(Y)+\cdots , \end{eqnarray}

respectively. Here, the normalisation $u_0'(0)=1$ and the existence of the displacement velocity $u_1(0)=U_d$ imply that $U_0 \sim Y+U_d$ for large $Y$ to match the boundary layer solution to the core solutions. While at $Y=0$ , of course, $V_0,U_0$ and $U_0'$ must vanish to satisfy the no-slip conditions. Substituting (2.18) and the Taylor expansion $v_b=\lambda y+\cdots$ into (2.11), the leading-order equations are obtained as

(2.19) \begin{eqnarray} (i\alpha ({\lambda }Y-c_1)-\partial _Y^2)V_0+{\lambda }U_0=0, \end{eqnarray}

(2.20) \begin{eqnarray} (i\alpha ({\lambda }Y-c_1)-\partial _Y^2)U_0''=0. \end{eqnarray}

From (2.9), $\lambda$ is simply unity when $a=0$ ; however, for our later purposes, we will carry out the calculations while keeping the general $\lambda$ . The change of the variable $\xi =\sigma Y+\xi _0$ with the constants $\xi _0=-\sigma c_1/{\lambda }$ and $\sigma =(i\lambda \alpha )^{1/3}$ transforms (2.20) to Airy’s equation for $\partial _{\xi }^2U_0$ :

(2.21) \begin{eqnarray} (\xi -\partial _{\xi }^2)\partial _{\xi }^2U_0=0. \end{eqnarray}

After some algebra, the solution satisfying the boundary conditions can be found as

(2.22) \begin{eqnarray} \hspace {-3mm} U_0=\sigma ^{-1}\left [\xi +\frac {1}{\kappa (\xi _0)}\left (-\xi \kappa (\xi )-\text {Ai}'(\xi )+\text {Ai}'(\xi _0) \right ) \right ]. \end{eqnarray}

For large $Y$ , the magnitude of the functions $\kappa (\xi )=\int ^{\infty }_{\xi }\text {Ai}(\check {\xi })\textrm{d}\check {\xi }$ and $\text {Ai}'(\xi )$ rapidly diminish. Therefore, referring to the aforementioned matching condition, the displacement velocity is found as

(2.23) \begin{eqnarray} U_{d}= \frac {\text {Ai}'(\xi _0)}{\sigma \kappa (\xi _0)}-\frac {c_1}{{\lambda }}. \end{eqnarray}

This is the only information we need for the dispersion relation in the wall layer analysis. Hence, we do not solve (2.19) to find $V_0$ .

The dispersion relation readily follows from the solvability condition for the inhomogeneous equation (2.14). The adjoint solution is merely $u_0$ because the homogeneous equation (2.13) is self-adjoint. Multiplying (2.14) by $u_0$ and integrating from 0 to 1,

(2.24) \begin{eqnarray} \hspace {-3mm} \int ^1_0 u_0 \mathscr {L} u_1 dy= - \int ^1_0 \left \{\frac {4 k_0^2c_1}{\alpha ^2 y^2} +2k_1 k_0 \left (\frac {2}{\alpha ^2y}-1 \right )\right \} u_0^2 \textrm{d}y. \end{eqnarray}

The left-hand side becomes $U_d$ after integration by parts and using (2.13). Therefore, $k_1$ and $c_1$ are now linked by

(2.25) \begin{eqnarray} -\frac {1}{(\lambda \alpha )^{1/3}} \frac {i^{5/3}\text {Ai}'(\xi _0)}{ \kappa (\xi _0)}+{c_1} \left(I_1-\frac 1{{\lambda }}\right)+k_1I_2 =0, \end{eqnarray}

where

(2.26) \begin{eqnarray} I_1(\alpha )=\frac {4 k_0^2}{\alpha ^2} \int ^1_0 \frac {u_0^2}{y^2} \textrm{d}y,\qquad I_2(\alpha )=2k_0 \int ^1_0 \left (\frac {2}{\alpha ^2y}-1 \right )u_0^2 \textrm{d}y, \end{eqnarray}

and $\xi _0=i^{7/3}(\lambda \alpha )^{1/3} c_1/{{\lambda }}$ .

If the mode is neutral (i.e. $c_1$ is purely real), the imaginary part of (2.25) yields

(2.27) \begin{eqnarray} \Im \left \{ \frac {i^{5/3}\text {Ai}'(i^{7/3}s)}{ \kappa (i^{7/3}s)} \right \}=0,\qquad s=(\lambda \alpha )^{1/3} c_1/{{\lambda }}. \end{eqnarray}

A real root of this function can be found at $s=s_0\approx 2.2972$ (Miles Reference Miles1960). Setting $\lambda =1$ , the following approximation is obtained for the lower branch shown in figure 3(b):

(2.28) \begin{eqnarray} c=2.2972(0.978Re)^{-1/3}. \end{eqnarray}

The real part of (2.25) can now be solved for $k_1$ :

(2.29) \begin{eqnarray} k_1 =\frac {q_0 +s_0 (1-I_1)}{\alpha ^{1/3}I_2 },\qquad q_0=\frac {i^{5/3}\text {Ai}'(i^{7/3}s_0)}{ \kappa (i^{7/3}s_0)}\approx {-1.00.} \end{eqnarray}

The integrals (2.26) can be worked out numerically after substituting $u_0$ found in (2.15). For $\alpha =0.978$ , using $k_0=1.442$ estimated before, the right-hand side of (2.29) can be computed as approximately 14.9. Therefore, we have the asymptotic result

(2.30) \begin{eqnarray} k=1.442-14.9Re^{-1/3} \end{eqnarray}

used in figure 3(a). The wall layer solution (2.22) is compared with the numerical neutral solution in figure 4(b).

2.4. Upper branch asymptotic analysis

The analysis for the upper branch is more complicated, and thus the neutrality of the perturbation is assumed from the beginning to facilitate discussion. The critical difference to the lower branch is that the majority of the wall layer becomes inviscid. Viscosity acts within two thinner sublayers: the critical layer sitting in the middle of the inviscid zone and the Stokes layer adjacent to the wall (see figure 5 b).

The core structure is identical to the lower branch case, but the wall layer thickness, again denoted as $\epsilon$ , is different. Within the wall layer, the radial perturbation velocity is $O(\epsilon )$ to leading order, and it is merely a linear profile. Since the critical layer appears within the wall layer, the phase speed $c$ must be $O(\epsilon )$ . A displacement velocity induced across the critical layer affects the next-order component of the wall layer radial velocity, which is $O(\epsilon ^2)$ . A displacement velocity is also induced by the Stokes layer, and it must match to the contribution from the critical layer. The former velocity is proportional to the Stokes layer thickness $\Delta$ , so we require $\triangle =\epsilon ^2$ . Finally, by balancing the $O(c)$ advection term in (2.11) with the $O(Re^{-1}\triangle ^{-2})$ diffusion term within the Stokes layer, we arrive at the scaling $\epsilon =Re^{-1/5}$ and $\triangle =Re^{-2/5}$ . The critical layer thickness $\tilde {\epsilon }=Re^{-1/3}$ is indeed thinner than the wall layer thickness $\epsilon$ .

We begin the formal asymptotic analysis by writing $k=k_0+\epsilon k_1+\cdots$ and $c=\epsilon c_1 +\cdots$ . Within the wall layer, the perturbation velocity components are expanded as

(2.31) \begin{eqnarray} u=\epsilon U_0(Y)+\epsilon ^2 U_1(Y)+\cdots ,\qquad v=V_0(Y)+\epsilon V_1(Y)+\cdots , \end{eqnarray}

where $Y=y/\epsilon$ . Substituting them into (2.11), the leading-order equations are obtained as

(2.32a) \begin{eqnarray} i\alpha ( Y- c_1)V_0+U_0=0, \end{eqnarray}

(2.32b) \begin{eqnarray} i\alpha ( Y- c_1) U''_{0}=0, \end{eqnarray}

and the next-order equations are

(2.33a) \begin{eqnarray} i\alpha ( Y- c_1)V_1+U_1=0, \end{eqnarray}

(2.33b) \begin{eqnarray} i\alpha ( Y- c_1) U_1''+2Yk_0^2V_0=0. \end{eqnarray}

For $U_0$ , we impose $U_0(0)=0$ to satisfy the non-penetrating condition on the wall. Also, $U_0\sim Y+U_{d}$ as $Y\rightarrow \infty$ to match the normalised core solution. Then, the solution to the leading-order problem can be found as

(2.34) \begin{eqnarray} U_0=Y, \qquad V_0=-\frac {Y}{i\alpha ( Y- c_1)}, \end{eqnarray}

which yields $U_{d}=0$ . The singularity that appears in $V_0$ should be resolved within the critical layer around $Y=c_1$ . Given $V_0$ , we can integrate (2.33b ) twice to find

(2.35) \begin{eqnarray} U_1=-\frac {k_0^2}{\alpha ^2}\left \{(Y-c_1)^2 + 4c_1( Y- c_1)[\ln |Y- c_1|+\phi ] \right .\nonumber \\ \left . -2c_1^2[\ln |Y- c_1|+\psi ] \right \}. \end{eqnarray}

The constants $\phi$ and $\psi$ are unknown and may differ for $Y\gt c_1$ and $Y\lt c_1$ . To distinguish the constants in the former region from those in the latter region, we attach a subscript $+$ to the former and $-$ to the latter. The constants $\phi _+$ and $\psi _+$ must be purely real to match the core solution which is still inviscid at the corresponding order. The critical layer analysis reveals that

(2.36) \begin{eqnarray} \phi _- =\phi _+ -i\pi ,\qquad \psi _-=\psi _+ -i\pi . \end{eqnarray}

This ‘logarithmic phase-shift’ may not be surprising to readers familiar with the asymptotic analysis of shear flows. To find this shift, one needs to use the expansions

(2.37a) \begin{eqnarray} u&=&\epsilon c_1+\tilde {\epsilon } \zeta +O(\epsilon ^2\ln \epsilon )+\epsilon ^2 \tilde {U}_1(\zeta )+O(\epsilon \tilde {\epsilon } \ln \epsilon )+\epsilon \tilde {\epsilon } \tilde {U}_2(\zeta )+\cdots , \end{eqnarray}

(2.37b) \begin{eqnarray} v&=&\epsilon \tilde {\epsilon }^{-1} \tilde {V}_0(\zeta )+\tilde {V}_1(\zeta )+\cdots , \end{eqnarray}

valid when $\zeta =\tilde {\epsilon }^{-1}(y- \epsilon c_1)$ is $O(1)$ . Equations (2.35) and (2.36) imply that

(2.38) \begin{eqnarray} \Im (u)\rightarrow -\epsilon ^2 \frac {6\pi c_1^2 k_0^2}{\alpha ^2} \qquad \text {as}\qquad Y\rightarrow 0. \end{eqnarray}

Within the Stokes layer, we write $Z=y/\epsilon ^2$ and expand

(2.39) \begin{eqnarray} u=\epsilon ^2\check {U}_0(Z)+\cdots ,\qquad v=\check {V}_0(Z)+\cdots . \end{eqnarray}

The no-slip boundary conditions are written as $\check {V}_0(0)=\check {U}_0(0)=\check {U}'_0(0)=0$ . Substitution of (2.39) to (2.11) yields the leading-order equations

(2.40) \begin{eqnarray} (i\alpha c_1+\partial _{Z}^2)\check {V}_0=0,\qquad (i\alpha c_1+\partial _{Z}^2)\check {U}_0''=0. \end{eqnarray}

The solution of the second equation, which does not grow exponentially for large $Z$ , can be found as follows:

(2.41) \begin{eqnarray} \check {U}_0=Z+\frac {e^{-qZ}-1}{q}, \qquad q=(1-i)\sqrt {\frac {\alpha c_1}{2}}. \end{eqnarray}

The exponential part is negligible when $Z$ is large. Thus, at the outer edge of the Stokes layer,

(2.42) \begin{eqnarray} \Im (u)\rightarrow -\epsilon ^2 \frac {1}{(2\alpha c_1)^{1/2}} \qquad \text {as}\qquad Z\rightarrow \infty . \end{eqnarray}

Matching (2.38) and (2.42), the scaled phase speed is found as

(2.43) \begin{eqnarray} c_1= \frac {\alpha ^{3/5}}{2^{1/5} (6\pi k_0^2)^{2/5}}. \end{eqnarray}

Here, we can use (2.16) to find $k_0$ , as the core analysis is identical to that for the lower branch case, while the dispersion relation (2.29) simplifies to

(2.44) \begin{eqnarray} k_1 =-\frac {c_1I_1}{I_2} \end{eqnarray}

because $u_1(0)=U_{d}=0$ , as remarked just below (2.34).

Substitution of $\alpha =0.978$ and $k_0=1.442$ into (2.43) and (2.44) gives $c_1\approx 0.198$ and $k_1\approx 1.89$ . The asymptotic approximations shown in figure 3 are thus obtained as

(2.45) \begin{eqnarray} c=0.198Re^{-1/5},\qquad k=1.442-1.89Re^{-1/5}. \end{eqnarray}

3.1. Lower branch asymptotic analysis

We introduce a small parameter $\epsilon =R_o^{-1/3}$ that describes the thickness of the wall layer near the inner cylinder. Note that this $\epsilon$ is not the same as that used in § 2.3; we will clarify their relationship later. The rotation rate $a$ is taken to be $O(\epsilon )$ so that $a_1=a/\epsilon$ is of order unity. The azimuthal wavenumber $ m$ is held fixed, while the axial wavenumber and the complex phase speed are written in the forms $k=k_0+\epsilon k_1+\cdots$ and $c=\epsilon c_1+\cdots$ , respectively.

In the inviscid core region, we expand

(3.4) \begin{eqnarray} (u,v,w,p)=(u_0,v_0,w_0,p_0)+\epsilon (u_1,v_1,w_1,p_1)+\cdots . \end{eqnarray}

The terms appearing on the right-hand side are all functions of $r$ . Substitution of those expansions to (3.2) yields the leading-order equation for $u_0(r)$ ,

(3.5) \begin{equation} \mathscr {L}u_0=0, \end{equation}

and the next-order equation for $u_1(r)$ ,

(3.6) \begin{equation} \mathscr {L}u_1=k_1F_2 u_0'-\Big (c_1F_0+ k_1F_1+a_1F_3 \Big )u_0. \end{equation}

The differential operator $\mathscr {L}$ and the functions $F_0, F_1, F_2$ and $F_3$ depend on $m,k_0,\eta$ , as can be found in the Appendix. The boundary conditions are $u_0(r_i)=u_0(r_o)=u_1(r_o)=0$ and $u_1(r_i)=U_d$ , where $U_d$ is the displacement velocity to be found by the wall layer analysis.

Within the wall layer near the inner cylinder, the suitable expansions are

(3.7) \begin{eqnarray} (u,v,w,p)=(\epsilon U_0,V_0,W_0,\epsilon P_0)+\cdots , \end{eqnarray}

where the terms on the right-hand side are functions of $Y=\epsilon ^{-1}(r-r_i)$ . Within the layer where $Y$ is of order unity, the base-flow angular velocity (2.4) behaves like $\Omega =\epsilon (\lambda Y+b a_1 )+\cdots$ with the constants

(3.8) \begin{eqnarray} \lambda \equiv \Omega _0'(0)=\frac {2}{r_i(1+\eta )},\qquad b\equiv \Omega _1(0)={1-\eta }. \end{eqnarray}

After some manipulations, from the leading-order problem, we can find that $U_0(Y)$ satisfies

(3.9) \begin{eqnarray} (im (\lambda Y+ba_1-c_1)-\partial _Y^2)U_0''=0 \end{eqnarray}

which, after re-reading $m, (c_1-ba_1)$ as $\alpha$ , $c_1$ , coincides with (2.20). Therefore, we can reuse the narrow-gap result to obtain the displacement velocity

(3.10) \begin{eqnarray} U_d= -\frac {i^{5/3}\text {Ai}'(\xi _0)}{(\lambda m)^{1/3} \kappa (\xi _0)}-\frac {c_1-b a_1}{\lambda },\qquad \xi _0=\frac {i^{7/3}(\lambda m)^{1/3}(c_1-b a_1)}{\lambda }. \end{eqnarray}

To solve the inhomogeneous problem (3.6), we use the inner product $\langle f,g \rangle \equiv \int _{r_i}^{r_o} f g \textrm{d}r$ . The homogeneous problem (3.5) is not self-adjoint, so we introduce the adjoint problem

(3.11) \begin{equation} \mathscr {L}^\dagger u^\dagger =0, \qquad u^\dagger (r_i)=u^\dagger (r_o)=0, \end{equation}

where the operator $\mathscr {L}^\dagger$ is given in the Appendix. The solution $u^\dagger (r)$ of this problem with the normalisation $(\hat u^\dagger )'(r_i)=1$ can be uniquely determined. Taking the inner product of $u^\dagger$ and (3.6), we obtain the dispersion relation

(3.12) \begin{equation} U_d+{\langle u^\dagger ,F_0 u_0 \rangle }c_1-k_1\langle u^\dagger ,F_2 u_0'-F_1 u_0\rangle +a_1\langle u^\dagger ,F_3 u_0\rangle =0, \end{equation}

where $U_d(m,\eta ,c_1,a_1)$ is given in (3.10).

Figure 6. Asymptotic results obtained by numerically solving (3.5). (a) Dependence of the leading-order axial wavenumber $k_0$ on the azimuthal wavenumber $ m$ for $\eta =0.99$ (red), 0.8 (green), 5/7 (blue), 0.5 (pink), 0.1 (black), 0.01 (orange). (b) The red and green curves are the results for $\eta =0.99$ and 0.8, respectively, shown in the $\delta \, m$ – $k_0$ plane. The circles are the narrow-gap limit results assuming $(k_0,\delta \, m)=(\overline {k}_0,\alpha )$ .

Figure 7. Dependence of $k_1$ on $a_1$ for the neutral modes computed by (3.12). (a) Results for wide-gap configurations. The azimuthal wavenumber is fixed at $ m=1$ . (b) Results for narrow-gap configurations. The azimuthal wavenumber is varied as $ m=0.978/\delta$ to observe the convergence to the narrow-gap results with $\alpha =0.978$ (circles). In the narrow-gap results, we set $(\delta ^{1/3}k_1,\delta ^{1/3}a_1)=(\overline {k}_1,\overline {a}_1)$ .

3.2. Numerical analysis

The wide-gap version of the asymptotic results requires numerical work to extract useful information. First, the leading-order problem, (3.5) and the associated boundary conditions, needs to be solved. The two boundary conditions can be satisfied only when specific values of $k_0(m,\eta )$ are taken. The admissible values of $k_0$ for various $\eta$ are summarised in figure 6(a). Only integer values of $ m$ are physically relevant, though real values of $ m$ can be used in numerical computation without issue. If the perturbations are assumed to be neutral, the next step is to find $c_1(m,\eta ,a_1)$ from the imaginary part of (3.12). Then the real part of (3.12) can be used to compute $k_1(m,\eta ,a_1)$ after some numerical integrations. The dependence of $k_1$ on $a_1$ is perfectly linear, as shown in figure 7(a) for the case $ m=1$ .

Figure 8. The scaled growth rate. In panels (a) and (b), the wide gap $\eta =5/7$ is used with $ m=2$ . The red, green and blue curves are the linearised Navier–Stokes results for $R_o= 10^9$ , $10^{10}$ and $10^{11}$ , respectively. The symbols are the asymptotic results. (a) $a=0$ ; (b) $k=k_0+R_o^{-1/3}k_1$ with $(k_0,k_1)=(1.651,-46.0)$ . The circles are computed by the wide-gap formula (3.12). Panels (c) and (d) are similar results but for the narrow gap $\eta =0.99$ with $ m=0.978/\delta$ . (c) $a=0$ with $(k_0,a_1)=(1.442,0)$ ; (d) $k=k_0+R_o^{-1/3}k_1$ with $(k_0,k_1)=(1.15,-10)$ . The asymptotic computation uses the narrow-gap limit formula (3.13) with $\alpha =0.978$ .

Now, let us explore the connection between the numerical analysis presented in this section and the narrow-gap limit analysis discussed in § 2.3. Recall that in the latter analysis, we required $Re = \delta \, R_o = \sqrt {1 - \eta } \, R_o$ to be large to define the small parameter $\overline {\epsilon }=Re^{-1/3}$ . Here, the overline is added to distinguish from $\epsilon =R_o^{-1/3}$ defined in § 3.1; likewise, we use $\overline {k}_0, \overline {k}_1,\overline {c}$ and $\overline {c}_1$ to denote the corresponding quantities that appeared in § 2.3. The two wavenumber expansions $k=\overline {k}_0+\overline {\epsilon } \overline {k}_1+\ldots$ and $k=k_0+\epsilon k_1+\ldots$ imply the equivalence $\overline {k}_0=k_0$ and $\overline {k}_1=\delta ^{1/3}k_1$ . Of course, to observe the convergence to the narrow-gap result, we must set $\alpha =\delta \, m$ because the wavenumber $\alpha$ is defined on the variable $x$ . Also, from the growth rate equality $i\alpha Re \overline {c}=imR_o c$ , we expect $\overline {c}_1=\delta ^{-5/3}c_1$ to hold. Moreover, we write $\overline {a}_1=\delta ^{1/3}a_1$ because, with this and the above rescaling, in the limit $\eta \rightarrow 1$ , the dispersion relation (3.12) reduces to

(3.13) \begin{eqnarray} -\frac {1}{ \alpha ^{1/3}} \frac {i^{5/3}\text {Ai}'(\overline {\xi }_0)}{ \kappa (\overline {\xi }_0)}+\overline {c}_1 (I_1-1)+\overline {k}_1I_2 +\overline {a}_1\left (1-\frac {I_1}{2}\right )=0, \end{eqnarray}

where $\overline {\xi }_0=i^{7/3}\alpha ^{1/3}(\overline {c}_1-\overline {a}_1)$ . This is precisely a generalisation of (2.25) for non-zero rotation rates. The convergence of $k_0$ to the narrow-gap result as $\eta \rightarrow 1$ is graphically evident in figure 6(b). Figure 7(b) depicts the similar excellent convergence to the narrow-gap result for the $a_1$ – $k_1$ relation. The narrow-gap result in this figure was produced by (3.13) for $ m=0.978/\delta$ , so the value of $\overline {k}_1$ at $a_1=0$ corresponds to the lower branch asymptotic curve seen in figure 3(a).

Figure 9. Comparison of the eigenfunctions between the leading-order asymptotic predictions (symbols) and the linearised Navier–Stokes results (lines) with $R_o=10^{10}, \ 10^{11},\ 10^{12}$ and $10^{13}$ . The parameters are $\eta =5/7$ , $a=0$ , $ m=2$ and $k= 1.65$ , where the mode is unstable. Red, real part; green, imaginary part. (a) Comparison in terms of the inviscid core scaling. (b) Comparison in terms of the viscous wall layer scaling.

The dispersion relation (3.12) can also be used for non-neutral modes. Figure 8 presents the scaled growth rate computed by the asymptotic analysis. Panels (a) and (b) are the results for $\eta =5/7$ , $ m=2$ , which are used by Deguchi (Reference Deguchi2017). The value of $k_0$ is found as 1.651 from (3.5), and the variation of $c_{1i}=\Im (c_1)$ on $k_1$ for $a_1=0$ can be found by (3.12) as the circles in panel (a). The solid curves represent the full Navier–Stokes computations based on (3.2), which, as expected, converge to the asymptotic results as $Re$ increases. To compare the results, we plotted $R_o^{1/3}c_i$ against $R_o^{1/3}(k-k_0)$ ; recall the scaling described just above (2.12). The flow becomes unstable when $k_1 \approx R_o^{1/3}(k-k_0)$ is smaller than $-20$ , with the growth rate reaching its maximum around $k_1=-46.0$ . Fixing $k_1$ at this value and varying $a_1$ results in panel (b). An increase in the rotation rate has a stabilising effect. To observe the convergence in the Navier–Stokes results, we set $k=k_0+R_o^{-1/3}k_1$ .

Figures 8(c) and 8(d) compare the linearised Navier–Stokes results for $\eta =0.99$ (curves) with the narrow-gap asymptotic results (circles). Again, we choose $\alpha =0.978$ used in § 2 for the asymptotic computation. Thus, $k_0$ is found as $1.442$ by (2.16), and in panel (c), where $c_{1i}$ is plotted against $k_1$ for $a_1=0$ , the neutral point $k_1\approx -14.9$ recovers the approximation (2.30). In the Navier–Stokes computations, the choice $\eta =0.99$ implies that the suitable azimuthal wavenumber is $m=0.978/\delta =9.78$ . The asymptotic convergence of the numerical data $(\delta ^{-5/3}R_o^{1/3}c_i,\delta ^{1/3}R_o^{1/3}(k-k_0))$ towards $(\overline {c}_{1i},\overline {k}_1)$ computed by (3.13) is reasonable, considering that there are now two factors causing errors: the finiteness of $\delta$ and $Re^{-1}$ . Even if we choose the closest integer value of $ m=10$ , the result does not differ significantly. Panel (d) shows the similar comparison, but now the rotation rate is varied for fixed $k_1=-10$ . Overall, we can observe the similar trend as the wide-gap case. Figure 9 illustrates a comparison between the eigenfunction of the linearised Navier–Stokes equations (3.2) and the solution of the asymptotic problem. All the above results unequivocally demonstrate the convergence of the Navier–Stokes computation towards the asymptotic results.

Instability occurs for all values of $\eta \in (0,1)$ , if $R_o$ is large enough. This is demonstrated in figure 6, where the case of a stationary inner cylinder case ( $a=0$ ) is shown. Of particular interest from an application point of view is the extent to which the D17 mode can penetrate into the Rayleigh stable region, i.e. how much we can increase the rotation rate $a$ before the flow completely stabilises. Figure 7 shows that the lower branch type instability always exists as long as $a$ is $O(R_o^{-1/3})$ . The upper branch type survives a bit larger $a$ of $O(R_o^{-1/5})$ . When $a$ reaches $O(1)$ , the stability is determined by the inviscid problem similar to (3.5), but with a critical layer located away from the walls (Drazin & Reid Reference Drazin and Reid1981; Billant & Gallaire Reference Billant and Gallaire2005; Caillol & Maslowe Reference Caillol and Maslowe2007). This type of instability may appear when the wavenumbers are fixed and the Reynolds number is increased. For the narrow-gap limit case, taking long wavelengths $k=O(Re^{-1})$ and $\alpha =O(Re^{-1})$ leads to another large- $Re$ asymptotic limit of instability for $a=O(1)$ , However, this limit, similar to that considered by Cowley & Smith (Reference Cowley and Smith1985), cannot be realised for wide gaps, since $m$ must be an asymptotically small non-zero number. Continuing the neutral curve for general $a$ using (3.2) is a computationally challenging task, so we do not pursue further investigation.

4.1. Asymptotic suction boundary layer

The generalised version of the instability theory for the asymptotic suction boundary layer, incorporating the effect of wall curvature, can be naturally derived from Taylor–Couette flow by writing the boundary conditions as follows:

(4.1) \begin{equation} \mathbf {v}(r_i,\theta ,z)=(-S,R_i,0),\quad \mathbf {v}(r_o,\theta ,z)=(-\eta S,R_o,0). \end{equation}

Here, $S$ represents the non-dimensional suction velocity on the inner cylinder, which introduces a uniform radial cross-flow; the same flow configuration was studied by Gallet et al. (Reference Gallet, Doering and Spiegel2010). The base flow solution can be found to be $(u,v,w)=(-(r_i/r)S,R_o r \Omega ,0)$ , with the function $\Omega (r)$ in the azimuthal component modified as

(4.2) \begin{eqnarray} \Omega =Ar^{-r_iS}+B/r^2,\qquad A=\frac {1}{r_o^{1-r_iS}}\frac {1-\eta ^2 a}{1-\eta ^{2-r_iS}},\qquad B=r_i\frac {\eta a-\eta ^{1-r_iS}}{1-\eta ^{2-r_iS}}. \end{eqnarray}

The associated stability problem is governed by (3.2), but the linear operator must include the cross-flow effect as $\mathcal {L}=im (\Omega -c)-R_o^{-1}(\partial _r^2+ (r_iS/r)\partial_r+r^{-1}\partial _r-r^{-2}m^2-k^2)$ .

If we take large values of $S$ , the base flow develops a boundary layer near the inner cylinder. For simplicity, we set $a=0$ . Then, it is easy to check that as $S\rightarrow \infty$ , the base flow behaves as

(4.3) \begin{eqnarray} v\sim \frac {R_or_o}{r}+O(S^{-1}) \qquad \text {if } r-r_i=O(1), \end{eqnarray}

(4.4) \begin{eqnarray} v\sim \frac {R_o}{\eta }(1-e^{-y})+O(S^{-1})\qquad \text {if} \ y=S(r-r_i)=O(1). \end{eqnarray}

This observation suggests that the boundary layer thickness is $O(S^{-1})$ .

To clarify the connection with the standard boundary layer analysis, we apply the following manipulations.

1. (i) Replace the length scale with one based on boundary layer thickness, using $y$ and the wavenumbers $\beta \equiv S^{-1}k$ and $\alpha \equiv m/Sr_i$ .

2. (ii) Denote $(v,u,w)$ as $(u_x,u_y,u_z)$ to express the equations in the familiar form used the boundary layer analysis. Also, rewrite $\eta p$ as $p$ .

3. (iii) Redefine the Reynolds number as $Re\equiv R_o/\eta S$ , and write $u_{xb}\equiv r_i\eta \Omega$ and $C=r_i\eta c$ . This Reynolds number is naturally defined by the free stream velocity at the edge of the boundary layer and the boundary layer thickness.

These transformations recast the stability problem into

(4.5a) \begin{align} \mathcal {L}u_x+\left(\!2\frac {u_{xb}}{Sr_i} +\frac {r}{r_i} \frac {du_{xb}}{dy}\! \right)u_y-\frac {r_i}{r^2 S Re}u_y+\frac {1}{ S Re}\left(\!\frac {u_x}{Sr^2}-\frac {2i r_i\alpha u_y}{r^2}\!\right)=-i\alpha \frac {r_i}{r} p, \end{align}

(4.5b) \begin{align} \mathcal {L}u_y-2\frac {u_{xb}}{Sr_i}u_x +\frac {r_i}{r^2 S Re}u_y+\frac {1}{ S Re}\left(\frac {u_y}{Sr^2} +\frac {2 i r_i\alpha u_x}{r^2}\right)=- \frac {dp}{dy}, \end{align}

(4.5c) \begin{align} \mathcal {L}u_z=-i\beta { p}, \end{align}

(4.5d) \begin{align} i\alpha \frac {r_i}{r}u_x+\frac {du_y}{dy}+\frac {u_y}{Sr}+i\beta u_z=0, \end{align}

where $\mathcal {L}=i\alpha (u_{xb}-C)-Re^{-1}(\partial _y^2+(r_i/r)\partial_y+S^{-1}r^{-1}\partial _y-(r_i/r)^2\alpha ^2-\beta ^2)$ . No asy- mptotic approximations have been made at this stage, provided that $\Omega$ defined in (4.2) is used. Therefore, together with the no-slip conditions on the cylinder walls, the linearised problem remains valid throughout the entire annular region.

If $S$ is large, the stability problem may be simplified to

(4.6a) \begin{eqnarray} \mathcal {L}{u_x}+u'_{xb}{u_y}=-i{\alpha } {p},\qquad \mathcal {L}{u_y}=-{p}'+g u_{xb}{u_x}, \end{eqnarray}

(4.6b) \begin{eqnarray} \mathcal {L}{u_z}=-i{\beta }{p},\qquad i{\alpha } {u_x}+{u_y}'+i{\beta }\,{u_z}=0, \end{eqnarray}

where $\mathcal {L}={\mathop {\textrm {i}}\nolimits } \alpha (u_{xb}-c)-Re^{-1}(\partial _{y}^2+\partial _{y}-\alpha ^2- \beta ^2)$ , $u_{xb}=1-e^{-y}$ and $g=2/Sr_i$ . The primes denote derivatives with respect to $y$ . At $y=0$ , the no-slip conditions $u_x=u_y=u_z=0$ are imposed, and we assume the perturbations vanish as $y\rightarrow \infty$ . The above boundary layer problem, formulated near the inner cylinder wall, constitutes an ‘intermediate reduced system’, as the magnitudes of the terms in (4.6) are not strictly balanced. To find the correct reduced system, a detailed asymptotic analysis in § 4.2 is required. Nonetheless, the same reduced problems are obtained whether the asymptotic analysis begins with (4.6) or (4.5), indicating that the former equations, despite their simplicity, include all the essential terms for the asymptotic reduction.

It should also be noted that the intermediate reduced system, excluding the Görtler term proportional to $g$ , has the similar structure as the local stability problem of general boundary layer flows. The constant $g$ is the parameter that represents the ratio of the boundary layer thickness to the radius of curvature of the wall. If the wavelengths are assumed to be comparable to the boundary layer thickness, the Görtler term is often omitted in the local stability studies. However, when the streamwise wavelength becomes sufficiently long, this term begins to significantly influence stability, as we shall see shortly.

In the context of the current parallel boundary layer problem, when the Görtler term is absent, (4.6) reduces to the Orr–Sommerfeld and Squire equations for the asymptotic suction boundary layer, where the onset of the instability is well known to occur at $(\alpha ,\beta ,Re)=(0.15547,0,27189)$ ; see Hocking (Reference Hocking1975); Drazin & Reid (Reference Drazin and Reid1981). The mechanism of this instability is essentially of the TS wave type. Note also that Görtler vortex type instabilities appear when $G=gRe^2$ is an $O(1)$ negative number, akin to Hall’s (Reference Hall1983) work on growing boundary layers. In this case, we need to recover several terms we omitted from (4.5). As noted earlier, those terms do not affect our analysis that follows.

Figure 10. Linear stability of the asymptotic suction boundary layer. The streamwise wavenumber is fixed at $\alpha =0.0036$ . (a) The red solid curves represent the neutral curves for the convex wall case $g=2\times 10^{-4}$ , with stable/unstable regions indicated by the labels. These results are computed using the intermediate reduced equations (4.6) within the boundary layer. The green circles show the same results, but they are computed using the full linearised Navier–Stokes equations (4.5d ) with $S=50.25$ . The black dashed curve corresponds to the flat wall case $g=0$ . The blue dash-dotted curve is the asymptotic result given by (4.19) for $g=2\times 10^{-4}$ . (b) Phase speed of the neutral modes for $g=2\times 10^{-4}$ .

Figure 10 shows the neutral spanwise wavenumber $\beta$ and phase speed $c$ as the Reynolds number is varied, with the streamwise wavenumber fixed at $\alpha =0.0036$ . The red solid curves are computed by the intermediate boundary layer approximation (4.6) with $g=2\times 10^{-4}$ ; the flow becomes unstable in the region between them. The neutral curves closely match the green circles, which are the neutral points computed using the full linearised Navier–Stokes equations (4.5d ) with $S=50.25$ . In the linearised Navier–Stokes computation, the radius ratio of the cylinders are chosen as $\eta =0.995$ , resulting in $g=2(1-\eta )/S\eta \approx 2\times 10^{-4}$ . The value of the azimuthal wavenumber is $ m=36$ , which corresponds to $\alpha \approx 0.0036$ . The Reynolds number, as mentioned earlier, has been rescaled as $Re=R_o/\eta S$ .

In the figure, the black dashed curve represents the result for the pure asymptotic suction boundary layer case ( $g=0$ ). Comparing the black and red curves reveals that even small variations in $g$ can cause significant changes in the shape of the neutral curve. At the large-Reynolds-number regime, the effect of wall curvature generates a new ‘D17’ instability in a narrow band of $\beta$ . The asymptotic behaviour of the D17 modes can be explained by the lower and upper branch theories, analogues to those discussed in § 2, with the lower branch elaborated in § 4.2.

For spanwise-independent perturbations, we can theoretically show that the stability remains unchanged from the $g=0$ case. Therefore, the unstable region at $g=2\times 10^{-4}$ represents a mixture of the TS mode and the D17 mode. Varying $\beta$ from zero along the red curve illustrates a smooth transition from TS waves to D17 modes, with the Reynolds number reaching the minimum value $Re=3.01\times 10^7$ at $\beta =0.0467$ . Note that this Reynolds number is not a critical value, as it is not optimised across all wavenumber pairs. Clearly, for any value of $g$ , the critical Reynolds number cannot fall below 27 819, which occurs at $\beta =0$ . We have confirmed that the critical Reynolds number remains at this value for $g=2\times 10^{-4}$ . In general, for relatively large $\alpha$ , small values of $g$ have no significant impact on stability. Numerical computations were necessary to verify this, as Squire’s theorem does not hold for $g\neq 0$ .

4.2. Lower-branch asymptotic analysis

Here, we conduct a lower-branch-type matched asymptotic expansion analysis for the stability problem (4.5). The key assumptions in the asymptotic analysis, in addition to the largeness of $S$ , are $Re\gg 1$ and $g\gg Re^{-2}$ . The second condition is satisfied at the high-Reynolds-number limit if $g\gt 0$ is held constant.

The analysis proceeds along the same lines as in § 2.3. An inviscid main layer emerges at which $y=O(1)$ , while the no-slip conditions are satisfied via the viscous wall layer of thickness $O(\epsilon )$ with $\epsilon =(gRe^2)^{-1/6}\ll 1$ . Those regions are somewhat similar to the main deck and lower deck of the triple-deck theory, respectively, but a structure like the upper deck does not emerge. Using the small parameter $\epsilon$ , we expand the spanwise wavenumber and the phase speed as $\beta =\beta _0+\epsilon \beta _1+\cdots$ and $C=\epsilon C_1+\cdots$ , respectively. Also, we rescale the streamwise wavenumber as $\alpha =g^{1/2}\alpha _0$ .

In the main layer, we write

(4.7) \begin{eqnarray} \left [ \begin{array}{c} u_x\\ u_y\\ u_z\\ p \end{array} \right ] = \left [ \begin{array}{c} u_{x0}(y)+\epsilon u_{x1}(y)+\cdots \\ g^{1/2}(u_{y0}(y)+\epsilon u_{y1}(y)+\cdots )\\ g^{1/2}(u_{z0}(y)+\epsilon u_{z1}(y)+\cdots )\\ g(p_0(y)+\epsilon p_1(y)+\cdots ) \end{array} \right ]. \end{eqnarray}

Substituting (4.7) into (4.5) (or (4.6)), the leading-order equations can be combined to yield

(4.8) \begin{eqnarray} \mathscr {L}u_{y0}=0,\qquad \mathscr {L}=\partial _y^2-\beta _0^2 -\frac {u''_{xb}}{u_{xb}} +\frac {\beta _0^2}{\alpha _0^2}\frac {u'_{xb}}{u_{xb}}. \end{eqnarray}

While, as expected, the next-order equations yield the inhomogeneous equation

(4.9) \begin{eqnarray} \mathscr {L}u_{y1}=\frac {C_1}{u_{xb}}\left (u''_{y0}-\beta _0^2u_{y0}-\frac {\beta _0^2}{\alpha _0^2}\frac {u'_{xb}}{u_{xb}}u_{y0}\right )+2\beta _0\beta _1\left (1-\frac {1}{\alpha _0^2}\frac {u'_{xb}}{u_{xb}}\right )u_{y0}. \end{eqnarray}

The non-trivial solution of (4.8) satisfying $u_{y0}=0$ at $y=0$ and $u_{y0}\rightarrow 0$ as $y\rightarrow \infty$ can be found as

(4.10) \begin{eqnarray} u_{y0}=\lambda ^{-1}u_{xb} e^{-\beta _0y} \end{eqnarray}

only when

(4.11) \begin{eqnarray} \beta _0=2\alpha _0^2. \end{eqnarray}

The constant $\lambda$ represents the value of $u'_{xb}$ at $y=0$ , and the factor $\lambda ^{-1}$ in (4.10) normalises $u'_{y0}=1$ at $y=0$ . Note that when the sign of $g$ is negative (i.e. in the case of a concave wall), constructing a non-trivial solution becomes impossible.

The solvability condition of the inhomogeneous equation (4.9) then determines the dispersion relation. The far-field condition $u_{y1}\rightarrow 0$ as $y\rightarrow \infty$ remains unchanged from the leading-order problem. Near the wall, however, the effect of the displacement velocity from the viscous wall layer must be accounted, such that $u_{y1}=U_{yd}$ at $y=0$ . The constant $U_{yd}$ can be found by the wall layer analysis.

The viscous wall layer expansions are

(4.12) \begin{eqnarray} \left [ \begin{array}{c} u_x\\ u_y\\ u_z\\ p \end{array} \right ] = \left [ \begin{array}{c} U_{x0}(Y)+\cdots \\ g^{1/2}(\epsilon U_{y0}(Y)+\cdots )\\ g^{1/2}(U_{z0}(Y)+\cdots )\\ g(\epsilon P_0(Y)+\cdots ) \end{array} \right ]. \end{eqnarray}

The no-slip conditions are $U_{x0}=U_{y0}=U'_{y0}=0$ at $Y=0$ . The leading-order wall-normal component $U_{y0}$ satisfies the equation

(4.13) \begin{eqnarray} (i\alpha (\lambda Y-C_1)-\partial _Y^2)U_{y0}''=0, \end{eqnarray}

which is the same equation as (2.20). Therefore, the displacement velocity $U_{yd}=\lim _{Y\rightarrow \infty } U_{y0}-Y$ is obtained as

(4.14) \begin{eqnarray} U_{yd}= -\frac {i^{5/3}\text {Ai}'(\xi _0)}{(\lambda \alpha _0)^{1/3} \kappa (\xi _0)}-\frac {C_1}{\lambda },\qquad \xi _0=\frac {i^{7/3}(\lambda \alpha _0)^{1/3}C_1}{\lambda }. \end{eqnarray}

Multiplying (4.10) by the inhomogeneous equation (4.9) and applying integration by parts, we obtain a dispersion relation similar to (2.25):

(4.15) \begin{eqnarray} -\frac {1}{(\lambda \alpha _0)^{1/3}} \frac {i^{5/3}\text {Ai}'(\xi _0)}{ \kappa (\xi _0)}+{C_1} \left (I_1-\frac 1{{\lambda }}\right )+\beta _1I_2 =0, \end{eqnarray}

where

(4.16) \begin{eqnarray} I_1=\int _0^{\infty }\frac {4\beta _0u'_{xb}-u''_{xb}}{\lambda ^2}e^{-2\beta _0y}\textrm{d}y,\qquad I_2=\int _0^{\infty }\frac {4u_{xb}u'_{xb}-2\beta _0u_{xb}^2}{\lambda ^2}e^{-2\beta _0y}\textrm{d}y,\,\,\,\, \end{eqnarray}

and $\xi _0=i^{7/3}(\lambda \alpha _0)^{1/3} C_1/{{\lambda }}$ . For the neutral disturbances, the result simplifies to

(4.17) \begin{eqnarray} C_1=\frac {\lambda s_0}{(\lambda \alpha _0)^{1/3}},\qquad \beta _1 =\frac {q_0+s_0 (1-\lambda I_1)}{(\lambda \alpha _0)^{1/3}I_2 } , \end{eqnarray}

where $q_0$ is the number defined in (2.29), and recall that $s_0\approx 2.2972$ .

So far, we considered a general $u_{xb}(y)$ , allowing the results to be applicable to other boundary layer flows. The specific outcomes for the asymptotic suction boundary layer are obtained by setting $u_{xb}=1-e^{-y}$ and $\lambda =1$ . For example, the integrals in (4.17) can be solved analytically, resulting in

(4.18) \begin{eqnarray} C_1=\frac {s_0}{\alpha _0^{1/3}},\qquad \beta _1 =\frac {8\alpha _0^4+6\alpha _0^2+1}{\alpha _0^{1/3} } \left ( q_0 - \frac {4s_0\alpha _0^2}{4\alpha _0^2+1} \right ). \end{eqnarray}

For the parameters used in figure 10, the leading-order wavenumbers are obtained as $(\alpha _0,\beta _0)=(0.255,0.130)$ , and we have the approximations

(4.19) \begin{eqnarray} \beta =0.130-3.31 g^{-1/6} Re^{-1/3},\qquad C=3.62 g^{-1/6}Re^{-1/3}, \end{eqnarray}

plotted by the blue dash-dotted curves. When the Reynolds number is large, these analytical results provide an excellent approximation for the ‘lower branch’, which appears as the upper branch in figure 10(a). Figure 11 examines the complex phase velocity $C=C_r+iC_i$ near this branch. As noted in § 3, predicting the asymptotic behaviour of the non-neutral modes requires using the full dispersion relation (4.15). As the Reynolds number increases, the numerical data $((\beta -\beta _0)/\epsilon ,C/\epsilon )$ converges to the asymptotic result $(\beta _1,C_1)$ . The points where $C_i$ becomes zero correspond to (4.19). Figure 12 compares the eigenfunctions at those neutral points. Once again, we observe clear convergence of the numerical results to the leading-order asymptotic solution.

Figure 11. Asymptotic convergence of the phase speed for $\alpha =0.0036$ , $g=2\times 10^{-4}$ . The curves are computed by (4.6). The points are the asymptotic result (4.15). (a) Real part; (b) imaginary part.

Figure 12. Comparison of the eigenfunctions at the lower-branch neutral point at $(Re,g,\alpha ,\beta )=(10^{12},$ $2\times 10^{-4},0.0036,0.130-3.30\epsilon )$ . The solid and dashed curves are the real and imaginary parts of $u_y$ computed by (4.6). The symbols are the corresponding asymptotic predictions. (a) Comparison in the core scaling. The symbols are the leading-order solution $u_{y0}$ found in (4.10). (b) Comparison in the wall-layer scaling. The symbols are $U_{y0} - Y$ , where $U_{y0}$ is obtained in the same manner as in the derivation of (2.22).