2.1. Half-plane: inviscid

A half-plane, however, can permit a starting vorticity distribution which spatially grows towards the boundary. Consider a boundary at $y=0$ and ${\it\eta}>0$ so that this is an inflow boundary when the fluid domain is $y>0$ . Taking the Fourier transform in  $x$ of (2.3) forces

(2.5) $$\begin{eqnarray}{\rm\Delta}{\it\psi}(x,y,t)=f(y-{\it\eta}t)\text{e}^{\text{i}k(x-yt+{\it\eta}t^{2}/2)},\end{eqnarray}$$

where $f$ is an arbitrary function. If ${\it\psi}$ is to be a modal disturbance which depends exponentially on $t$ (and $x$ ), i.e. ${\it\psi}(x,y,t)=\hat{{\it\psi}}_{k}(y)\exp (\text{i}kx+{\it\sigma}t)$ , then

(2.6) $$\begin{eqnarray}f({\it\xi})=A\text{e}^{-({\it\sigma}{\it\xi}+\text{i}k{\it\xi}^{2}/2)/{\it\eta}}\end{eqnarray}$$

is the only possibility ( $A$ is an arbitrary normalisation) and therefore

(2.7) $$\begin{eqnarray}\left(\frac{\text{d}^{2}}{\text{d}y^{2}}-k^{2}\right)\hat{{\it\psi}}_{k}(y)=A\text{e}^{-({\it\sigma}y+\text{i}ky^{2}/2)/{\it\eta}}\end{eqnarray}$$

(no modal solution exists for ${\it\eta}=0$ ). We will show that this equation, where the initial vortical distribution decays exponentially as $y\rightarrow \infty$ , has growing modal disturbances (i.e. $\text{Re}({\it\sigma})>0$ ). No instability is possible if the boundary is an outflow boundary – i.e. ${\it\eta}<0$ – as $A$ is forced to be zero by boundedness. This suggests that the instability exists providing the boundary is not a zero-vorticity boundary which would force $A=0$ . In what follows, we adopt a non-slip boundary condition, as Ilin & Morgulis (Reference Ilin and Morgulis2013) originally did, but will revisit this issue in § 2.1.4 ((2.7) is a third-order differential equation for ${\it\psi}$ , integrated once to include an arbitrary constant – hence three boundary conditions are needed to specify a unique solution).

To confirm there is instability, (2.7) must be solved to derive the dispersion relation. Using a Green’s function approach, setting $A$ to 1 (w.l.o.g.) and imposing the two (usual inviscid) boundary conditions that $\hat{{\it\psi}}_{k}(0)=0$ and $\lim _{y\rightarrow \infty }\hat{{\it\psi}}_{k}(y)=0$ , (2.7) has the solution

(2.8) $$\begin{eqnarray}\hat{{\it\psi}}_{k}(y)=-\int _{0}^{y}\frac{\text{e}^{-ky}}{k}\sinh (k\bar{y})\text{e}^{-({\it\sigma}\bar{y}+\text{i}k\bar{y}^{2}/2)/{\it\eta}}\,\text{d}\bar{y}-\int _{y}^{\infty }\frac{\text{e}^{-k\bar{y}}}{k}\sinh (ky)\text{e}^{-({\it\sigma}\bar{y}+\text{i}k\bar{y}^{2}/2)/{\it\eta}}\,\text{d}\bar{y}.\end{eqnarray}$$

A third boundary condition needs to be applied to give the dispersion relation and, with no-slip at the influx boundary ( $\text{d}\hat{{\it\psi}}_{k}/\text{d}y=0$ at $y=0$ ), this is the relation

(2.9) $$\begin{eqnarray}\int _{0}^{\infty }\text{e}^{-k\bar{y}-({\it\sigma}\bar{y}+\text{i}k\bar{y}^{2}/2)/{\it\eta}}\,\text{d}\bar{y}=0.\end{eqnarray}$$

There is no solution to this for $k=0$ , indicating the absence of a 1D instability, but there is instability in 2D ( $k\neq 0$ ). The expression (2.9) is valid for finite ${\it\eta}$ , but it is now useful to consider small ${\it\eta}$ , where this instability can be understood as a critical layer instability. The (special) case where this critical layer is at the (influx) boundary (so it is in fact a boundary layer) was discussed by Ilin & Morgulis (Reference Ilin and Morgulis2013) in their inviscid Taylor–Couette set-up. The growth rates for such modes are $\text{Re}({\it\sigma})=O(\sqrt{{\it\eta}})$ . The other (generic) case when the critical layer is distinct from the inflow boundary has not been discussed before. The growth rate here is smaller – $\text{Re}({\it\sigma})=O({\it\eta}\log (1/{\it\eta}))$ (see below) – but is still larger than the default $O({\it\eta})$ which would be expected.

The asymptotic form of the dispersion relation (2.9) can be derived as ${\it\eta}\rightarrow 0^{+}$ using standard steepest descent/saddle point ideas. Here we need two parts of the integrand to contribute at the same leading order and precisely cancel. This can happen in two ways, since there is just one saddle point at $\bar{y}_{s}=\text{i}{\it\sigma}/k$ : the contribution from an end point (clearly $\bar{y}=0$ ) cancels the contribution from the saddle point (the interior critical layer case) or the saddle point and the end point are effectively one and the same asymptotically (the boundary layer case). In the former case, the leading contributions from the end point at $\bar{y}=0$ (first term on the left-hand side in (2.10)) and from the saddle point (second term on the left-hand side) must satisfy

(2.10) $$\begin{eqnarray}\frac{{\it\eta}}{{\it\sigma}}+\sqrt{\frac{2{\rm\pi}{\it\eta}}{\text{i}k}}\text{e}^{-\text{i}{\it\sigma}-\text{i}{\it\sigma}^{2}/(2k{\it\eta})}=0.\end{eqnarray}$$

Now, for the saddle point to be asymptotically separated from the endpoint $\bar{y}=0$ , $|{\it\sigma}/k|=O(1)$ . This, together with the fact that the magnitude of the saddle point contribution $O(\sqrt{{\it\eta}}\exp ({\it\sigma}_{r}{\it\sigma}_{i}/(k{\it\eta})))$ (where ${\it\sigma}={\it\sigma}_{r}+\text{i}{\it\sigma}_{i}$ ) has to be $O({\it\eta})$ requires either ${\it\sigma}_{r}=O(1)$ and then ${\it\sigma}_{i}\ll 1$ , or ${\it\sigma}_{r}\ll 1$ and ${\it\sigma}_{i}=O(1)$ . The former case is inconsistent because the saddle point is in the wrong part of the $\bar{y}$ -complex plane, leaving the contribution from the end point $\bar{y}=0$ unbalanced. The latter situation, however, does yield solutions. Defining

(2.11) $$\begin{eqnarray}{\it\sigma}=\text{i}{\it\sigma}_{i}+{\it\delta}({\it\eta})\tilde{{\it\sigma}}_{r},\end{eqnarray}$$

with ${\it\sigma}_{i}$ and $\tilde{{\it\sigma}}_{r}$ both $O(1)$ quantities and ${\it\delta}({\it\eta})\rightarrow 0$ as ${\it\eta}\rightarrow 0$ , then (2.10) requires

(2.12) $$\begin{eqnarray}{\it\delta}({\it\eta})\tilde{{\it\sigma}}_{r}=-\frac{k{\it\eta}}{2{\it\sigma}_{i}}\left(\log \frac{1}{{\it\eta}}+\log \frac{2{\rm\pi}{\it\sigma}_{i}^{2}}{k}\right)-k{\it\eta}+O({\it\eta}^{2}\log 1/{\it\eta}),\end{eqnarray}$$

and

(2.13) $$\begin{eqnarray}\text{Arg}[\text{e}^{\text{i}({\it\sigma}_{i}^{2}/(2k{\it\eta})+{\rm\pi}/4)}]=O({\it\eta}\log 1/{\it\eta}),\end{eqnarray}$$

so there is instability with asymptotic growth rate $O({\it\eta}\log 1/{\it\eta})$ for a discrete set of frequencies

(2.14) $$\begin{eqnarray}{\it\sigma}_{i}=-\sqrt{{\textstyle \frac{1}{2}}{\rm\pi}k{\it\eta}(8n-1)},\end{eqnarray}$$

where $n$ is an integer of $O(1/{\it\eta})$ . The form of the corresponding eigenfunction is discussed in § 2.1.2.

The other situation – the boundary layer case – arises when the saddle point is within $O(\sqrt{{\it\eta}})$ of the endpoint at $\bar{y}=0$ , i.e. $|{\it\sigma}|=O(\sqrt{{\it\eta}})$ . At this point, the endpoint and the saddle point contributions cannot be considered separately. Instead $\bar{y}$ and ${\it\sigma}$ must be rescaled so that (2.9) becomes

(2.15) $$\begin{eqnarray}\int _{0}^{\infty }\text{e}^{-\hat{{\it\sigma}}z-\text{i}z^{2}}\,\text{d}z=0,\end{eqnarray}$$

where $z\,:=\,\,\bar{y}/\sqrt{2{\it\eta}/k}$ and

(2.16) $$\begin{eqnarray}{\it\sigma}=\sqrt{k{\it\eta}/2}\,(\hat{{\it\sigma}}_{r}+\text{i}\hat{{\it\sigma}}_{i})+O({\it\eta}).\end{eqnarray}$$

This must be solved numerically, but a good estimate for the eigenvalues can be found by treating the contributions from the end point and saddle point as if they are separated. This means taking the integer $n$ to be $1\ll n\ll 1/{\it\eta}$ in the frequency expression (2.14) for the internal critical layer mode and calculating the corresponding growth using (2.12). This leads to the asymptotic form

(2.17a,b ) $$\begin{eqnarray}\hat{{\it\sigma}}_{i}\sim -\sqrt{{\rm\pi}(8n-1)}\quad \hat{{\it\sigma}}_{r}\sim \frac{1}{|\hat{{\it\sigma}}_{i}|}\log ({\rm\pi}\hat{{\it\sigma}}_{i}^{2})\quad \text{as integer}~n\rightarrow \infty ,\end{eqnarray}$$

which performs very well even for the first eigenvalue since $\hat{{\it\sigma}}_{i}$ is already ${<}\!-4.7$ then: see table 1. Figure 1 shows a typical critical layer eigenfunction and three boundary layer eigenfunctions for $k=1$ and ${\it\eta}=10^{-3}$ .

Figure 1. Unstable boundary layer eigenfunctions ((a): blue solid line is the most unstable, the red dashed line is the third most unstable and the black dash-dot line is the fifth most unstable mode) and an unstable critical layer eigenfunction ((b) with the critical layer shown as a red dashed line) for ${\it\eta}=10^{-3}$ and $k=1$ . In all, the real part of $u={\it\psi}_{y}$ is shown.

Table 1. The 10 most unstable eigenvalues from the dispersion relation (2.15) and asymptotic estimates from (2.17).

2.1.1. Inviscid asymptotics $0<{\it\eta}\ll 1$ for the boundary layer instability

Here a modal stream function solution of the form ${\it\Psi}(x,y,t)={\it\psi}(y)\text{e}^{\text{i}kx+{\it\sigma}t}$ is sought with a boundary layer of thickness $\sqrt{{\it\eta}}$ , as first uncovered by Ilin & Morgulis (Reference Ilin and Morgulis2013) (see their § 3.2). Defining the new boundary layer variable ${\it\xi}\,:=\,\,y\sqrt{k/(2{\it\eta})}$ and splitting the stream function into an expansion of large-scale parts and a boundary layer corrections (hatted variables)

(2.18) $$\begin{eqnarray}{\it\psi}\,:=\,\,({\it\psi}_{0}(y)+\hat{{\it\psi}}_{0}({\it\xi}))+\sqrt{{\it\eta}}({\it\psi}_{1}(y)+\hat{{\it\psi}}_{1}({\it\xi}))+\cdots \,,\end{eqnarray}$$

we look for an instability with vanishing frequency ( $=k\,\times$ speed of the inflow boundary) to leading order,

(2.19) $$\begin{eqnarray}{\it\sigma}\,:=\,\,\sqrt{{\it\eta}k/2}\,\hat{{\it\sigma}}+O({\it\eta}).\end{eqnarray}$$

The governing equation

(2.20) $$\begin{eqnarray}({\it\sigma}+\text{i}ky+{\it\eta}\partial _{y})({\it\psi}_{yy}-k^{2}{\it\psi})=0\end{eqnarray}$$

is then simplified to

(2.21) $$\begin{eqnarray}{\it\psi}_{0\,yy}-k^{2}{\it\psi}_{0}=0\end{eqnarray}$$

for the leading flow and

(2.22) $$\begin{eqnarray}(\hat{{\it\sigma}}+2\text{i}{\it\xi}+\partial _{{\it\xi}})\hat{{\it\psi}}_{0{\it\xi}{\it\xi}}=0\end{eqnarray}$$

for its boundary layer correction. The interesting observation here is that there is no large-scale solution ${\it\psi}_{0}$ which can handle both the ${\it\eta}=0$ boundary conditions that $v(x,0,t)=0$ and $\lim _{y\rightarrow \infty }v(x,y,t)=0$ . Crucially, this means that the boundary layer correction $\hat{{\it\psi}}_{0}$ must be $O(1)$ so as to contribute at leading order to fix up the $v(x,0,t)=0$ boundary condition, rather than the usual $O(\sqrt{{\it\eta}})$ to satisfy the extra ${\it\eta}\neq 0$ tangential boundary condition $u(x,0,t)=0$ . As a consequence, there is an $O(1/\sqrt{{\it\eta}})$ tangential velocity $u$ in the boundary layer which must vanish at ${\it\xi}=0$ . Since the large-scale flow cannot contribute at this order, this non-slip condition is solely on the boundary layer flow and is sufficient to determine the growth rate. Integrating (2.22) twice and incorporating the fact that $\lim _{{\it\xi}\rightarrow \infty }\hat{{\it\psi}}_{0\,{\it\xi}}=0$ , leads to

(2.23) $$\begin{eqnarray}\hat{{\it\psi}}_{0\,{\it\xi}}=-A\int _{{\it\xi}}^{\infty }\text{e}^{-\hat{{\it\sigma}}z-\text{i}z^{2}}\,\text{d}z,\end{eqnarray}$$

where $A$ is an arbitrary constant. Imposing $u=0$ at $y=0$ then forces $\hat{{\it\psi}}_{0{\it\xi}}(0)=0$ , which is precisely condition (2.15).

2.1.2. Inviscid asymptotics $0<{\it\eta}\ll 1$ for the critical layer instability

An interior critical layer instability is constructed by looking for a mode of $O(1)$ frequency. Adopting the expansion (2.11), where again both ${\it\sigma}_{i}$ and $\tilde{{\it\sigma}}_{r}$ are $O(1)$ , the critical layer is centred on $y_{c}\,:=\,\,-\!{\it\sigma}_{i}/k$ (so ${\it\sigma}_{i}<0$ ) and, as in the boundary layer case, has thickness $O(\sqrt{{\it\eta}})$ . Defining the critical layer variable

(2.24) $$\begin{eqnarray}{\it\xi}\,:=\,\,\frac{y-y_{c}}{\sqrt{{\it\eta}}}\end{eqnarray}$$

the (leading order) stream function in the critical layer, $\hat{{\it\psi}}$ , satisfies

(2.25) $$\begin{eqnarray}({\it\delta}\tilde{{\it\sigma}}_{r}/\sqrt{{\it\eta}}+\text{i}k{\it\xi}+\partial _{{\it\xi}})\hat{{\it\psi}}_{{\it\xi}{\it\xi}}=0,\end{eqnarray}$$

which can be integrated once to give

(2.26) $$\begin{eqnarray}\hat{{\it\psi}}_{{\it\xi}{\it\xi}}=\sqrt{{\it\eta}}\text{e}^{-{\it\delta}\tilde{{\it\sigma}}_{r}{\it\xi}/\sqrt{{\it\eta}}-\text{i}k{\it\xi}^{2}/2}\end{eqnarray}$$

(choosing the normalisation of the mode here for convenience later) and then two further times to give

(2.27) $$\begin{eqnarray}\displaystyle \hat{{\it\psi}}_{{\it\xi}}\! & = & \displaystyle \!\sqrt{{\it\eta}}\left(\int _{0}^{{\it\xi}}\text{e}^{-{\it\delta}\tilde{{\it\sigma}}_{r}x/\sqrt{{\it\eta}}-\text{i}kx^{2}/2}\,\text{d}x+A\right)\sim \sqrt{{\it\eta}}\left(\pm \sqrt{\frac{{\rm\pi}}{2\text{i}k}}+A\right)+\text{h.o.t.}\quad \text{as }{\it\xi}\rightarrow \pm \infty \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(2.28) $$\begin{eqnarray}\displaystyle \hat{{\it\psi}}\! & = & \displaystyle \!\sqrt{{\it\eta}}\left(\int _{0}^{{\it\xi}}({\it\xi}-x)\text{e}^{-{\it\delta}\tilde{{\it\sigma}}_{r}x/\sqrt{{\it\eta}}-\text{i}kx^{2}/2}\,\text{d}x+A{\it\xi}\right)+B\nonumber\\ \displaystyle \! & {\sim} & \displaystyle \!\sqrt{{\it\eta}}{\it\xi}\left[\pm \sqrt{\frac{{\rm\pi}}{2\text{i}k}}+A\right]+B+\text{h.o.t.}\quad \text{as }{\it\xi}\rightarrow \pm \infty ,\end{eqnarray}$$

where $A$ and $B$ are $O(1)$ constants.

Outside the critical layer, the stream function consists of a WKB-type solution and simple exponentials which satisfy Laplace’s equation,

(2.29) $$\begin{eqnarray}\displaystyle & \displaystyle y_{c}<y:\quad {\it\psi}_{+}=C_{+}\text{e}^{-ky}+\frac{E_{+}}{k^{2}(y-y_{c})^{2}}\text{e}^{-[{\it\delta}\tilde{{\it\sigma}}_{r}(y-y_{c})+\text{i}k(y-y_{c})^{2}/2]/{\it\eta}} & \displaystyle\end{eqnarray}$$

(2.30) $$\begin{eqnarray}\displaystyle & \displaystyle 0<y<y_{c}:\quad {\it\psi}_{-}=C_{-}\text{e}^{-ky}+D_{-}\text{e}^{ky}+\frac{E_{-}}{k^{2}(y-y_{c})^{2}}\text{e}^{-[{\it\delta}\tilde{{\it\sigma}}_{r}(y-y_{c})+\text{i}k(y-y_{c})^{2}/2]/{\it\eta}}, & \displaystyle\end{eqnarray}$$

where $C_{\pm }$ , $D_{-}$ and $E_{\pm }$ are constants whose order of magnitude will be set by matching to the critical layer stream function. The fact that ${\it\psi}_{+}$ and $\text{d}{\it\psi}_{+}/\text{d}y$ must vanish as $y\rightarrow \infty$ is imposed by including only the decaying exponential for $y>y_{c}$ : the WKB mode vanishes as $y\rightarrow \infty$ when ${\it\eta}>0$ provided $\tilde{{\it\sigma}}_{r}>0$ .

The governing equation (2.20) is third order and so ${\it\psi},{\it\psi}_{y}$ and ${\it\psi}_{yy}$ must be continuous everywhere. The oscillatory form of $\hat{{\it\psi}}_{{\it\xi}{\it\xi}}$ in the critical layer means it must match entirely to the outer WKB-type solution

(2.31) $$\begin{eqnarray}{\it\psi}^{WKB}\,:=\,\,\frac{-{\it\eta}^{3/2}}{k^{2}(y-y_{c})^{2}}\text{e}^{-[{\it\delta}\tilde{{\it\sigma}}_{r}(y-y_{c})+\text{i}k(y-y_{c})^{2}/2]/{\it\eta}}\end{eqnarray}$$

either side of the critical layer so $E_{\pm }=-{\it\eta}^{3/2}$ . This means that the outer WKB solution only contributes at $O(\sqrt{{\it\eta}})$ to the tangential velocity $u$ near the critical layer, whereas the critical layer stream function $\hat{{\it\psi}}$ forces an $O(1)$ tangential flow (see (2.27)). As a result there must be a large-scale flow of $O(1)$ in both $u$ and $v$ outside the critical layer: this explains the scaling of the integration constants in (2.28) and means that $C_{\pm }$ and $D_{-}$ are all $O(1)$ . There are then six (complex) conditions to be satisfied at leading order by the five remaining constants and complex eigenvalue ${\it\sigma}$ . The first four are matching conditions on $u$ ( ${\it\psi}_{y}$ ) and $v$ ( $-\text{i}k{\it\psi}$ ) as $y\rightarrow y_{c}^{+}$ ,

(2.32a,b ) $$\begin{eqnarray}-kC_{+}\text{e}^{-ky_{c}}=\sqrt{\frac{{\rm\pi}}{2\text{i}k}}+A\quad \text{and}\quad C_{+}\text{e}^{-ky_{c}}=B,\end{eqnarray}$$

and $y\rightarrow y_{c}^{-}$ ,

(2.33a,b ) $$\begin{eqnarray}-kC_{-}\text{e}^{-ky_{c}}+kD_{-}\text{e}^{ky_{c}}=-\sqrt{\frac{{\rm\pi}}{2\text{i}k}}+A\quad \text{and}\quad C_{-}\text{e}^{-ky_{c}}+D_{-}\text{e}^{ky_{c}}=B,\end{eqnarray}$$

or, eliminating $A$ and $B$ , simply two jump conditions across the critical layer

(2.34a,b ) $$\begin{eqnarray}[u]_{-}^{+}=\sqrt{\frac{2{\rm\pi}}{\text{i}k}}\quad \text{and}\quad [v]_{-}^{+}=0.\end{eqnarray}$$

The two remaining (boundary) conditions, $v(x,0,t)=0$ and $u(x,0,t)=0$ , require

(2.35a,b ) $$\begin{eqnarray}C_{-}+D_{-}=0\quad \text{and}\quad -kC_{-}+kD_{-}-\frac{\text{i}\sqrt{{\it\eta}}}{ky_{c}}\text{e}^{[{\it\delta}\tilde{{\it\sigma}}_{r}y_{c}-\text{i}ky_{c}^{2}/2]/{\it\eta}}=0,\end{eqnarray}$$

where, while the outer WKB-type solution contributes at $O(\sqrt{{\it\eta}})$ to $u$ near the critical layer, it must contribute at $O(1)$ to $u$ at the inflow boundary. The resulting dispersion relation is

(2.36) $$\begin{eqnarray}\frac{\text{i}\sqrt{{\it\eta}}}{ky_{c}}\text{e}^{[{\it\delta}\tilde{{\it\sigma}}y_{c}-\text{i}ky_{c}^{2}/2]/{\it\eta}}+\sqrt{\frac{2{\rm\pi}}{\text{i}k}}\text{e}^{-ky_{c}}=0,\end{eqnarray}$$

which leads to the same leading expressions for $\tilde{{\it\sigma}}_{r}$ and ${\it\sigma}_{i}$ given in (2.12) and (2.14). The necessary change in the scaling of the WKB solution contribution gives the dominant $O({\it\eta}\log 1/{\it\eta})$ contribution to $\tilde{{\it\sigma}}_{r}$ and the exact numerical counterbalancing of the large-scale tangential flow gives the subdominant $O({\it\eta})$ contribution.

2.1.3. The need for shear and inflow

To emphasise that shear is a key ingredient of the instability, the above problem can be solved for the shearless flow $\boldsymbol{U}=\hat{\boldsymbol{x}}+{\it\eta}\hat{\boldsymbol{y}}$ , leading to the requirement that there must exist $\text{Re}({\it\sigma})>0$ with

(2.37) $$\begin{eqnarray}\int _{0}^{\infty }\text{e}^{-k\bar{y}-{\it\sigma}\bar{y}/{\it\eta}}\,\text{d}\bar{y}=0,\end{eqnarray}$$

which is never satisfied. The presence of an inflow boundary is also crucial: converting the above inflow boundary to an outflow boundary ( ${\it\eta}\rightarrow -{\it\eta}$ ) removes the instability. This is why the instability discussed here is not relevant to the considerable literature on suction boundary layers, where suction is always a stabilising effect (e.g. Joslin Reference Joslin1998).

2.1.4. The third boundary condition

Imposing the parametrised boundary condition $(1-{\it\beta})u+{\it\beta}\,\text{d}u/\text{d}y=0$ as the second boundary condition at $y=0$ gives the modified dispersion relation

(2.38) $$\begin{eqnarray}\int _{0}^{\infty }\text{e}^{-k\bar{y}-({\it\sigma}\bar{y}+\text{i}k\bar{y}^{2}/2)/{\it\eta}}\,\text{d}\bar{y}=\frac{{\it\beta}}{1-{\it\beta}}\end{eqnarray}$$

( ${\it\beta}=0$ recovers non-slip and ${\it\beta}=1$ a no-vorticity or equivalently stress-free condition as $v=0$ along $y=0$ ). For the boundary layer instabilities, the left-hand side is $O(\sqrt{{\it\eta}})$ , so this dispersion relation can only be assumed similar to the non-slip relation (2.15) if ${\it\beta}\ll \sqrt{{\it\eta}}$ . In fact, numerical computations indicate that the boundary layer instability is suppressed by ${\it\beta}\approx 2.6\sqrt{{\it\eta}/k}\ll 1$ . For the critical layer instabilities, the condition (2.10) becomes a possible balance between three different terms

(2.39) $$\begin{eqnarray}\frac{{\it\eta}}{{\it\sigma}}+\sqrt{\frac{2{\rm\pi}{\it\eta}}{\text{i}k}}\text{e}^{-\text{i}{\it\sigma}-\text{i}{\it\sigma}^{2}/(2k{\it\eta})}=\frac{{\it\beta}}{1-{\it\beta}},\end{eqnarray}$$

where ${\it\sigma}$ now has an $O(1)$ frequency as in (2.11). If ${\it\beta}\gg {\it\eta}$ , the dominant balance is between the second and third terms (as opposed to the first and second for non-slip) and now leads to damped eigenvalues. It is then clear that, for instability to occur, there should be no restriction on the normal derivative of the tangential velocity at the inflow boundary. In practice this means that the instability only really occurs for non-slip boundary conditions when $0\leqslant {\it\eta}\ll 1$ , which incidentally is the one choice which, in concert with the no-normal flow condition, means no disturbance kinetic energy is being advected into the domain through the inflow boundary.

2.2. Half-plane: viscous

The base state (2.1) is unchanged (by design) when viscosity is introduced, but the linearised disturbance equation (2.2) now includes diffusion of vorticity:

(2.40) $$\begin{eqnarray}\left(\frac{\partial }{\partial t}+y\frac{\partial }{\partial x}+{\it\eta}\frac{\partial }{\partial y}\right){\it\omega}=\frac{1}{Re}{\rm\nabla}^{2}{\it\omega}.\end{eqnarray}$$

Looking for a modal solution ${\it\omega}(x,y,t)=\hat{{\it\omega}}(y)\exp (\text{i}kx+{\it\sigma}t)$ leads to the equation

(2.41) $$\begin{eqnarray}\frac{\text{d}^{2}\hat{{\it\omega}}}{\text{d}y^{2}}-Re\,{\it\eta}\frac{\text{d}\hat{{\it\omega}}}{\text{d}y}-[Re({\it\sigma}+\text{i}ky)+k^{2}]\hat{{\it\omega}}=0,\end{eqnarray}$$

which has the bounded solution (as $y\rightarrow \infty$ )

(2.42) $$\begin{eqnarray}\hat{{\it\omega}}=\text{e}^{Re\,{\it\eta}\,y/2}\text{Ai}\left[\frac{Re{\it\sigma}+k^{2}+{\textstyle \frac{1}{4}}Re^{2}{\it\eta}^{2}}{(\text{i}kRe)^{2/3}}+(\text{i}kRe)^{1/3}y\right],\end{eqnarray}$$

where $\text{Ai}(z)$ is the Airy function bounded as $|z|\rightarrow \infty$ with $|\!\arg (z)|<{\rm\pi}$ . The dispersion relation is then

(2.43) $$\begin{eqnarray}\int _{0}^{\infty }\text{e}^{(Re\,{\it\eta}/2-k)\bar{y}}\text{Ai}\left[\frac{Re{\it\sigma}+k^{2}+{\textstyle \frac{1}{4}}Re^{2}{\it\eta}^{2}}{(\text{i}kRe)^{2/3}}+(\text{i}kRe)^{1/3}\bar{y}\right]\text{d}\bar{y}=0\end{eqnarray}$$

and instability ( $\text{Re}({\it\sigma})$ crosses through 0 to become positive) occurs at a critical inflow ${\it\eta}_{crit}(k,Re)$ . This integral is actually the same as that treated by Gallet et al. (Reference Gallet, Doering and Spiegel2010) (see their expression (39)) after the transformations

(2.44a,b ) $$\begin{eqnarray}N_{g}^{2/3}\,:=\,\,2^{1/3}\frac{Re\,{\it\eta}-2k}{(kRe)^{1/3}},\quad a_{g}=-\!\left(\frac{4Re}{N_{g}k^{2}}\right)^{1/3}({\it\sigma}+k{\it\eta})\end{eqnarray}$$

( $N_{g}$ and $a_{g}$ from Gallet et al. (Reference Gallet, Doering and Spiegel2010)). For the $k$ of interest ( ${\leqslant}O(1)$ ), $k{\it\eta}\ll |{\it\sigma}|$ and we can reuse their critical value of $N_{g}$ , which has $a_{g}$ passing through the imaginary axis to give

(2.45a,b ) $$\begin{eqnarray}{\it\eta}_{crit}=\left(\frac{1}{2}N_{g}^{2}k\right)^{1/3}Re^{-2/3},\quad {\it\sigma}_{crit}=-\!\left(\frac{N_{g}k^{2}}{4Re}\right)^{1/3}a_{g},\end{eqnarray}$$

where $N_{g}=4.57557$ and $a_{g}=5.63551\text{i}$ (consistent with Gallet et al. (Reference Gallet, Doering and Spiegel2010), who quote ‘4.58’ and ‘5.62i’ for $N_{g}$ and $a_{g}$ , respectively). This threshold ${\it\eta}_{crit}$ tends to zero as $k\rightarrow 0$ , albeit with the unstable eigenfunction extending a distance $O((kRe)^{-1/3})$ in the $y$ direction. In terms of connecting this analysis to other problems, there are two notable cases: $k=O(1)$ , which is the interesting case in rotating flow where the wavenumber is forced to be an integer by periodicity, and $k=O(Re^{-1})$ , which is gives the most unstable disturbance in a domain bounded in $y$ (i.e. a channel, see Nicoud & Angilella Reference Nicoud and Angilella1997). In the former case, ${\it\eta}_{crit}=O(Re^{-2/3})$ and the implication from the scaling of the critical frequency is the growth rate away from criticality (i.e. $|{\it\eta}-{\it\eta}_{crit}|=O({\it\eta}_{crit})$ ) will be $O(Re^{-1/3})$ or $O(\sqrt{{\it\eta}})$ , as before. For $k=O(Re^{-1})$ , ${\it\eta}_{crit}=O(Re^{-1})$ , which is consistent with the numerical findings in Nicoud & Angilella (Reference Nicoud and Angilella1997) that the threshold ‘cross-flow’ Reynolds number for inflow instability in their plane Couette flow is independent of the shear Reynolds number.

For any given $k$ , further boundary layer-type instabilities exist as ${\it\eta}$ increases with the first six thresholds listed in table 2 for $k=1$ . Within this ‘boundary layer’ scaling of $\tilde{{\it\lambda}}\,:=\,\,-Re^{1/3}{\it\sigma}_{i}$ and $\tilde{{\it\eta}}_{crit}\,:=\,\,Re^{2/3}{\it\eta}_{crit}$ for $k=O(1)$ , asymptotic predictions for $\tilde{{\it\lambda}}\rightarrow \infty$ can be derived following the same route as in the inviscid case. This proves a little more involved, leading to two coupled relations,

(2.46) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{{\it\eta}}_{crit}=\tilde{{\it\lambda}}\left(\frac{3k}{2}\log (2{\rm\pi}\tilde{{\it\lambda}}/k)-\frac{k}{2}\log \left[\frac{2}{3k\log (2{\rm\pi}\tilde{{\it\lambda}}/k)}\right]\right)^{-1/3}, & \displaystyle\end{eqnarray}$$

(2.47) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\tilde{{\it\lambda}}^{2}}{\tilde{{\it\eta}}_{crit}}=\frac{1}{2}{\rm\pi}(8n-1)\quad n=1,2,3,\ldots , & \displaystyle\end{eqnarray}$$

$\tilde{{\it\lambda}}=O(10)$

$\log (\tilde{{\it\lambda}})$

Table 2. The six lowest inflow thresholds for instability in the dispersion relation (2.43) with $k=1$ as $Re\rightarrow \infty$ and asymptotic estimates from (2.46) and (2.47).

As in the inviscid case, there are also interior critical layer modes excited for even higher inflows: if the critical layer is at $y=-{\it\sigma}_{i}/k=O(1)$ , then

(2.48) $$\begin{eqnarray}{\it\eta}_{crit}=|{\it\sigma}_{i}|\left(\frac{2}{kRe\log Re}\right)^{1/3}\gg O(Re^{-2/3}).\end{eqnarray}$$

2.3. Inflow and suction together: inviscid plane Couette flow with suction

The half-plane system exhibits boundary inflow instability for all ${\it\eta}>0$ . This can be seen by a simple rescaling of space

(2.49a−c ) $$\begin{eqnarray}{\it\kappa}={\it\eta}^{1+{\it\beta}}k,\quad Y=y/{\it\eta}^{1+{\it\beta}}\quad \text{and}\quad {\it\epsilon}=1/{\it\eta}^{{\it\beta}},\end{eqnarray}$$

where ${\it\beta}>0$ , so that the (2.20) becomes

(2.50) $$\begin{eqnarray}({\it\sigma}+\text{i}{\it\kappa}Y+{\it\epsilon}\partial _{Y})({\it\psi}_{YY}-{\it\kappa}^{2}{\it\psi})=0,\end{eqnarray}$$

which is just the original equation with a new small number ${\it\epsilon}$ as ${\it\eta}\rightarrow \infty$ . However, this is rather artificial, since the applied pressure gradient also has to be increased with ${\it\eta}$ to maintain the constant shear in $y$ . Resorting to a constant pressure gradient instead now means the shear field decreases as ${\it\eta}$ increases, again making it difficult to draw general conclusions for large ${\it\eta}$ . Introducing another boundary is then the only alternative; and this must be an outflow boundary if the resulting base flow is to be steady and spatially non-developing. The simplest modification is to add an outflow boundary at $y=1$ which is moving at $\hat{\boldsymbol{x}}$ so that the constant-vorticity basic flow (2.1) is still a solution (Nicoud & Angilella Reference Nicoud and Angilella1997). The equivalent expression to (2.8) is then

(2.51) $$\begin{eqnarray}\displaystyle \hat{{\it\psi}}_{k}(y)\! & = & \displaystyle \!-\int _{0}^{y}\frac{\sinh (k\bar{y})}{k\sinh k}\sinh k(1-y)\text{e}^{-({\it\sigma}\bar{y}+\text{i}k\bar{y}^{2}/2)/{\it\eta}}\,\text{d}\bar{y}\nonumber\\ \displaystyle & & \displaystyle -\,\int _{y}^{1}\frac{\sinh k(1-\bar{y})}{k\sinh k}\sinh (ky)\text{e}^{-({\it\sigma}\bar{y}+\text{i}k\bar{y}^{2}/2)/{\it\eta}}\,\text{d}\bar{y},\end{eqnarray}$$

with the dispersion relation

(2.52) $$\begin{eqnarray}\int _{0}^{1}\sinh k(1-\bar{y})\text{e}^{-({\it\sigma}\bar{y}+\text{i}k\bar{y}^{2}/2)/{\it\eta}}\,\text{d}\bar{y}=0,\end{eqnarray}$$

when a non-slip condition is applied at the influx boundary $y=0$ for ${\it\eta}>0$ . This is essentially the same as the half-plane dispersion relation and will have unstable eigenvalues as there is an inflow boundary. The simple rescaling (2.49a−c ), however, is disallowed and instability is lost if ${\it\eta}$ becomes too large (see figure 12 of Nicoud & Angilella Reference Nicoud and Angilella1997). This could be interpreted as the stabilising influence of the newly introduced suction boundary ultimately overpowering the destabilising inflow boundary. Further evidence for this comes from the pressure-gradient-free version of this flow, which is also linearly unstable in the presence of viscosity (Doering et al. Reference Doering, Spiegel and Worthing2000). In this case, the base state varies exponentially in the cross-stream direction (see equation (2.13) of Doering et al. Reference Doering, Spiegel and Worthing2000) and possesses an area of linear instability in the $(\tan {\it\theta},Re)$ plane (Doering et al. (Reference Doering, Spiegel and Worthing2000, figure 3), where $\tan {\it\theta}$ is their proxy for ${\it\eta}$ ). The lower boundary of this instability region, ${\it\theta}_{lower}(Re)$ , plausibly scales like $Re^{-1}$ , which is the viscous threshold for the instability as described in § 2.2, whereas the upper boundary has ${\it\theta}_{upper}(Re)=\cot ^{-1}54,370$ , which appears to be suction ultimately stabilising the flow (Hocking Reference Hocking1974).

3.1. Basic state

There are two usual scenarios: (1) the swirl field ${\it\Omega}(r)$ is set by the boundary conditions as in Taylor–Couette flow or (2) the swirl field is determined by an imposed body (gravitational) force as in the astrophysical context. In the former, the swirl is determined by the azimuthal momentum equation given a radial flow with the radial momentum equation determining the pressure field. For example Ilin & Morgulis (Reference Ilin and Morgulis2013) discuss the inviscid, irrotational, axisymmetric Taylor–Couette-like flow

(3.11) $$\begin{eqnarray}{\it\Omega}_{inviscid}\,:=\,\,\frac{1}{r^{2}},\end{eqnarray}$$

which is the only possibility with axisymmetric radial flow (recall ${\it\Omega}(1)=1$ by non-dimensionalisation). Gallet et al. (Reference Gallet, Doering and Spiegel2010) consider the viscous Taylor–Couette situation, where possible base flows form a one-parameter family

(3.12) $$\begin{eqnarray}{\it\Omega}_{TC}=\frac{1-A}{r^{2}}+Ar^{-{\it\eta}Re}\end{eqnarray}$$

with $A$ a constant set by the motion of the outer cylinder, which is generally rotational ( $Z=A(2-{\it\eta}Re)r^{-{\it\eta}Re}$ ). In the (latter) astrophysical context, the acting gravitational (body) force sets the swirl field (via the radial momentum equation), which then sets the radial flow by the need to balance ensuing azimuthal viscous stresses. The focus here is on the latter situation, and we consider the general set of profiles

(3.13) $$\begin{eqnarray}{\it\Omega}=r^{{\it\alpha}}\end{eqnarray}$$

in order to understand how robust the boundary inflow instability is. Profiles with $-2<{\it\alpha}<0$ have angular momentum increasing with radius and are thus Rayleigh-stable (Rayleigh Reference Lord Rayleigh1917). The gradient of the vorticity $\text{d}Z/\text{d}r={\it\alpha}({\it\alpha}+2){\it\Omega}/r$ also does not change sign across the domain $r\in [1,a]$ , so that the flow is inviscidly stable (for disturbances which vanish at $r=1$ and $a$ ) by a rotating flow analogue of Rayleigh’s inflexion-point theorem (Rayleigh Reference Lord Rayleigh1880). Particularly interesting choices for ${\it\alpha}$ are ${\it\alpha}=-3/2$ , which is the Keplerian profile for a thin accretion disk due to the radial force balance

(3.14) $$\begin{eqnarray}-r{\it\Omega}^{2}=-\frac{GM}{r^{2}},\end{eqnarray}$$

(where $G$ is the gravitational constant and $M$ the mass of the central object), and ${\it\alpha}=-1$ , which is used to model spiral galaxies. Since we work with deviations away from the basic state, the exact body force required to maintain the underlying rotation profile is not explicitly needed in what follows. In contrast to the Taylor–Couette situation, the azimuthal component of the Navier–Stokes equations forces the existence of a small radial flow

(3.15) $$\begin{eqnarray}U=\frac{1}{Re}\left(\frac{\partial }{\partial r}\log \left|\frac{\partial }{\partial r}(r^{2}{\it\Omega})\right|-\frac{1}{r}\right)=\frac{{\it\alpha}}{Re}\frac{1}{r},\end{eqnarray}$$

which is an inflow if $\text{d}{\it\Omega}/\text{d}r<0$ ( ${\it\alpha}<0$ ). Studying the consequences of this small accretion flow is the motivation for this work.

3.2. 2D inviscid swirling flow with $\text{d}Z/\text{d}r=0$

We study the simplest case of vanishing vorticity gradient in the base flow first ( $\text{d}Z/\text{d}r=0$ ), so ${\it\Omega}=1/r^{2}$ or $1$ , to show how the analysis mirrors that in the half-plane case. The case of uniform rotation ${\it\Omega}=1$ initially looks uninteresting because the base flow needs an azimuthal as well as a radial body force to maintain it (since ${\it\alpha}=0$ in (3.15)). But it is worth considering, as then only the enforced radial inflow creates shear in the base flow and the question is whether this is enough to generate instability. The inviscid, linearised governing equation (3.8) is just

(3.16) $$\begin{eqnarray}\left(\frac{\partial }{\partial t}+{\it\Omega}(r)\frac{\partial }{\partial {\it\theta}}-\frac{{\it\eta}}{r}\frac{\partial }{\partial r}\right){\it\omega}=0,\end{eqnarray}$$

which is the ‘curved’ analogue of (2.2) and again third order rather than the usual second order. As before, the solution can be written down using characteristics as

(3.17) $$\begin{eqnarray}{\it\omega}(r,{\it\theta},t)={\it\omega}_{0}\left(r^{2}+2{\it\eta}t,{\it\theta}+\frac{1}{{\it\eta}}\int _{\sqrt{r^{2}+2{\it\eta}t}}^{r}{\it\Omega}(r^{\ast })r^{\ast }\,\text{d}r^{\ast }\right),\end{eqnarray}$$

where ${\it\omega}_{0}(r,{\it\theta})$ is the initial vorticity. After a discrete Fourier transform, the $\text{e}^{\text{i}m{\it\theta}}$ component is

(3.18) $$\begin{eqnarray}\left(\frac{1}{r}\frac{\partial }{\partial r}r\frac{\partial }{\partial r}-\frac{m^{2}}{r^{2}}\right){\it\psi}_{m}(r,t)=g(r^{2}+2{\it\eta}t)\exp \left(\frac{\text{i}m}{{\it\eta}}\int _{\sqrt{r^{2}+2{\it\eta}t}}^{r}{\it\Omega}(r^{\ast })r^{\ast }\,\text{d}r^{\ast }\right),\end{eqnarray}$$

where $g$ is an arbitrary function. For modal growth, ${\it\psi}_{m}(r,t)$ should depend on $t$ only through an $\exp ({\it\sigma}t)$ factor for some complex constant ${\it\sigma}$ , which forces

(3.19) $$\begin{eqnarray}g({\it\xi})=B\exp \left(\frac{{\it\sigma}{\it\xi}}{2{\it\eta}}-\frac{\text{i}m}{{\it\eta}}\int _{\sqrt{{\it\xi}}}^{a}{\it\Omega}(r^{\ast })r^{\ast }\,\text{d}r^{\ast }\right)\end{eqnarray}$$

( $B$ is some constant which can be set to 1 providing the influx boundary conditions do not force $B=0$ ) and, letting ${\it\psi}_{m}(r,t)=\hat{{\it\psi}}_{m}(r)\exp ({\it\sigma}t)$ , then

(3.20) $$\begin{eqnarray}\left(\frac{1}{r}\frac{\partial }{\partial r}r\frac{\partial }{\partial r}-\frac{m^{2}}{r^{2}}\right)\hat{{\it\psi}}_{m}(r)=\exp \left(\frac{{\it\sigma}r^{2}}{2{\it\eta}}-\frac{\text{i}m}{{\it\eta}}\int _{r}^{a}{\it\Omega}(r^{\ast })r^{\ast }\,\text{d}r^{\ast }\right).\end{eqnarray}$$

This has the solution

(3.21) $$\begin{eqnarray}\displaystyle \hat{{\it\psi}}_{m}(r)\! & = & \displaystyle \!\int _{1}^{r}G_{+}(r,\bar{r};a,m)\exp \left(\frac{{\it\sigma}\bar{r}^{2}}{2{\it\eta}}-\frac{\text{i}m}{{\it\eta}}\int _{\bar{r}}^{a}{\it\Omega}(r^{\ast })r^{\ast }\,\text{d}r^{\ast }\right)\,\text{d}\bar{r}\nonumber\\ \displaystyle & & \displaystyle +\,\int _{r}^{a}G_{-}(r,\bar{r};a,m)\exp \left(\frac{{\it\sigma}\bar{r}^{2}}{2{\it\eta}}-\frac{\text{i}m}{{\it\eta}}\int _{\bar{r}}^{a}{\it\Omega}(r^{\ast })r^{\ast }\,\text{d}r^{\ast }\right)\,\text{d}\bar{r},\end{eqnarray}$$

where

(3.22) $$\begin{eqnarray}\displaystyle & \displaystyle G_{-}(r,\bar{r};a,m)\,:=\,\,\frac{\bar{r}}{2m(a^{2m}-1)}\left(\bar{r}^{m}-\frac{a^{2m}}{\bar{r}^{m}}\right)\left(r^{m}-\frac{1}{r^{m}}\right) & \displaystyle\end{eqnarray}$$

(3.23) $$\begin{eqnarray}\displaystyle & \displaystyle G_{+}(r,\bar{r};a,m)\,:=\,\,\frac{\bar{r}}{2m(a^{2m}-1)}\left(\bar{r}^{m}-\frac{1}{\bar{r}^{m}}\right)\left(r^{m}-\frac{a^{2m}}{r^{m}}\right) & \displaystyle\end{eqnarray}$$

$\hat{{\it\psi}}_{m}(1)=0$

$\hat{{\it\psi}}_{m}(a)=0$

$\text{d}\hat{{\it\psi}}_{m}(a)/\text{d}r=0$

(3.24) $$\begin{eqnarray}\int _{1}^{a}(\bar{r}^{m+1}-\bar{r}^{1-m})\exp \left(\frac{{\it\sigma}\bar{r}^{2}}{2{\it\eta}}-\frac{\text{i}m}{{\it\eta}}\int _{\bar{r}}^{a}{\it\Omega}(r^{\ast })r^{\ast }\,\text{d}r^{\ast }\right)\,\text{d}\bar{r}=0.\end{eqnarray}$$

For ${\it\Omega}=1/r^{2}$ this is exactly expression (5.3) from Ilin & Morgulis (Reference Ilin and Morgulis2013), which they establish has unstable eigenvalues. The energy budget for an infinitesimal 2D disturbance,

(3.25) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\left\langle \frac{1}{2}\boldsymbol{u}^{2}\right\rangle ={\it\eta}\left\langle \frac{v^{2}-u^{2}}{r^{2}}\right\rangle -\left.{\it\eta}\int _{0}^{2{\rm\pi}}\frac{1}{2}v^{2}\right|_{r=1}\text{d}{\it\theta}-\left\langle ruv\frac{\text{d}{\it\Omega}}{\text{d}r}\right\rangle ,\end{eqnarray}$$

where $\langle ~\rangle \,:=\,\,\int _{0}^{2{\rm\pi}}\int _{1}^{a}r\,\text{d}r\,\text{d}{\it\theta}$ , $u=0$ at $r=1$ and $u=v=0$ at $r=a$ , clearly shows that enhanced growth rates of $O(\sqrt{{\it\eta}})$ are possible only if $\text{d}{\it\Omega}/\text{d}r$ is non-zero somewhere in the domain, emphasising the importance of azimuthal shear. In particular, for ${\it\Omega}=1$ , the last term drops, so that any instability can have a growth rate of only $O({\it\eta})$ at best. In fact, the dispersion relation appears only to have stable solutions, indicating that boundary inflow and just the shear due to the radial inflow are insufficient to drive any instability.

The ${\it\eta}\rightarrow 0$ asymptotic analysis for ${\it\Omega}=1/r^{2}$ mirrors that in the half-plane case with the saddle point at ${\it\sigma}+\text{i}m{\it\Omega}(r_{s})=0$ . For the boundary layer scenario, adopting the scalings

(3.26) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\sigma}=-\text{i}m{\it\Omega}(a)+\sqrt{{\it\eta}{\it\chi}/a}\,(\hat{{\it\sigma}}_{r}+\text{i}\hat{{\it\sigma}}_{i})+O({\it\eta}), & \displaystyle\end{eqnarray}$$

(3.27) $$\begin{eqnarray}\displaystyle & \displaystyle z\,:=\,\,\sqrt{a{\it\chi}}\,\frac{(a-r)}{\sqrt{{\it\eta}}}, & \displaystyle\end{eqnarray}$$

${\it\chi}\,:=\,\,-\!(m{\it\Omega}^{\prime }(a))/2$

${\it\Omega}^{\prime }(r)<0$

$z$

$\hat{{\it\sigma}}$

$O\left(\sqrt{{\it\eta}(-m{\it\Omega}^{\prime }(a)/(2a))}\right)$

For the critical layer scenario, the appropriate eigenvalue scaling is again (2.11), which, given ${\it\sigma}+\text{i}m{\it\Omega}(r_{s})=0$ , leads to the expansion

(3.28) $$\begin{eqnarray}r_{s}\,:=\,\,R_{s}+\frac{\text{i}{\it\delta}({\it\eta})\tilde{{\it\sigma}}_{r}}{m{\it\Omega}^{\prime }(R_{s})}+\cdots\end{eqnarray}$$

for the saddle point position, where ${\it\sigma}_{i}=-m{\it\Omega}(R_{s})$ defines $R_{s}$ . The balance between the leading contribution from the saddle point at $r=r_{s}$ (first term) and the endpoint $r=a$ (second term) is then

(3.29) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\sqrt{-\text{i}{\it\eta}\,{\rm\pi}}(r_{s}^{m+1}-r_{s}^{1-m})}{\sqrt{-{\textstyle \frac{1}{2}}mr_{s}{\it\Omega}^{\prime }(r_{s})}}\exp \left(\displaystyle \frac{{\it\sigma}r_{s}^{2}}{2{\it\eta}}-\frac{\text{i}m}{{\it\eta}}\int _{r_{s}}^{a}{\it\Omega}(r^{\ast })r^{\ast }\,\text{d}r^{\ast }\right) & \displaystyle \nonumber\\ \displaystyle & \quad -~\displaystyle \frac{\text{i}{\it\eta}\,(a^{m}-a^{-m})}{({\it\sigma}_{i}+m{\it\Omega}(a))}\exp \left(\displaystyle \frac{{\it\sigma}a^{2}}{2{\it\eta}}\right)=0, & \displaystyle\end{eqnarray}$$

which leads to the frequency condition

(3.30) $$\begin{eqnarray}\text{Arg}\left\{\exp \left[\text{i}\left(\frac{{\it\sigma}_{i}(R_{s}^{2}-a^{2})}{2{\it\eta}}-\frac{m}{{\it\eta}}\int _{R_{s}}^{a}{\it\Omega}(r^{\ast })r^{\ast }\,\text{d}r^{\ast }+\frac{{\rm\pi}}{4}\right)\right]\right\}=0\end{eqnarray}$$

and growth rate

(3.31) $$\begin{eqnarray}\displaystyle {\it\delta}\tilde{{\it\sigma}}_{r}=\frac{{\it\eta}}{(a^{2}-R_{s}^{2})}\left[\log \left(\frac{1}{{\it\eta}}\right)+2\log \left(\frac{(R_{s}^{m+1}-R_{s}^{1-m})(-{\it\sigma}_{i}-m{\it\Omega}(a))\sqrt{2{\rm\pi}}}{(a^{m}-a^{-m})\sqrt{-mR_{s}{\it\Omega}^{\prime }(R_{s})}}\right)+o(1)\right]. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Figure 2 shows that this performs well for numerical calculations with ${\it\eta}=10^{-3}$ and $10^{-4}$ . The most unstable boundary layer eigenfunction for ${\it\eta}=10^{-3}$ is shown in figure 3 along with a critical layer mode with ${\it\sigma}_{i}\approx -0.5$ which is typical.

Figure 2. Inviscid 2D instabilities for ${\it\alpha}=-2$ , $m=1$ and $a=2$ with scaled growth rate $\hat{{\it\sigma}}_{r}={\it\sigma}_{r}/\sqrt{-m{\it\eta}{\it\Omega}^{\prime }(a)/(2a)}$ plotted against frequency ${\it\sigma}_{i}$ . The right (left) vertical dashed line is ${\it\sigma}_{i}=-m{\it\Omega}(a)$ ( ${\it\sigma}_{i}=-m{\it\Omega}(1)$ ). The eigenvalues from a full 2D eigenvalue calculation are shown as (the upper) black dots for ${\it\eta}=10^{-3}$ ( $N=1000$ ) and (the lower) blue dots for ${\it\eta}=10^{-4}$ (since $N=4000$ here, the 15 most unstable eigenvalues are marked with a normal-sized dot and the rest by smaller blue dots). The red squares are the 20 most unstable modes from the boundary layer asymptotics (§ 3.2) with the dash-dot red line tracing the path of the rest (eigenvalues calculated from the boundary layer equation (2.22) are indistinquishable from the asymptotics at this scale). The green dashed lines are the critical layer asymptotic expression (3.31) with ${\it\eta}=10^{-3}$ and $10^{-4}$ plugged in and appropriately rescaled: $\hat{{\it\sigma}}_{r}={\it\delta}_{1}\tilde{{\it\sigma}}_{r}\sqrt{(2a)/(-m{\it\eta}{\it\Omega}^{\prime }(a))}$ and ${\it\sigma}_{i}=-m{\it\Omega}(R_{s})$ with $R_{s}$ taken over $[1,a]$ . Notice that, at ${\it\eta}=10^{-4}$ , even $N=4000$ struggles to fully resolve the critical layer eigenvalues near the inner radius – see the hump in the numerical data which breaks the otherwise excellent correspondence with the asymptotic prediction. This hump is delayed to lower ${\it\sigma}_{i}$ if $N$ is increased until it eventually disappears – e.g. the ${\it\eta}=10^{-3}$ curve.

Figure 3. The most unstable boundary layer eigenfunction (a) and an unstable critical layer eigenfunction with ${\it\sigma}_{i}\approx -0.5$ (b) for ${\it\eta}=10^{-3}$ , $m=1$ , $a=2$ and ${\it\alpha}=-2$ . In both, the green dashed line is $\text{Re}(u)$ and the blue solid line is $\text{Re}(v)$ . The critical layer position is shown as a red dashed line for the critical layer eigenfunction.

3.3. 2D inviscid swirling flow with $\text{d}Z/\text{d}r\neq 0$

The asymptotic analysis can easily be extended to treat more general rotation profiles ${\it\Omega}(r)$ where the gradient of vorticity is non-zero. The idea here is to assume some small viscosity to induce a radial inflow (via (3.15)) and then to show that the inviscid instability mechanism is robust to the exact form of the azimuthal flow (e.g. whether it has a vorticity gradient or not). So, taking the linear, inviscid version of perturbation equation (3.8) and inserting the usual ansatz ${\it\psi}(r,{\it\theta},t)={\it\psi}(r)\exp (\text{i}m{\it\theta}+{\it\sigma}t)$ gives

(3.32) $$\begin{eqnarray}\left({\it\sigma}+\text{i}m{\it\Omega}(r)-\frac{{\it\eta}}{r}\frac{\text{d}}{\text{d}r}\right)\mathscr{L}{\it\psi}=\frac{\text{i}m{\it\psi}}{r}\frac{\text{d}Z}{\text{d}r},\end{eqnarray}$$

where

(3.33) $$\begin{eqnarray}\mathscr{L}\,:=\,\,\frac{\text{d}^{2}}{\text{d}r^{2}}+\frac{1}{r}\frac{\text{d}}{\text{d}r}-\frac{m^{2}}{r^{2}}.\end{eqnarray}$$

The ${\it\eta}=0$ problem

(3.34) $$\begin{eqnarray}({\it\sigma}+\text{i}m{\it\Omega}(r))\mathscr{L}{\it\psi}=\frac{\text{i}m{\it\psi}}{r}\frac{\text{d}Z}{\text{d}r}\end{eqnarray}$$

subject to the two boundary conditions ${\it\psi}(1)={\it\psi}(a)=0$ seems to have no discrete modal solutions for the profiles ${\it\Omega}=r^{{\it\alpha}}$ studied here. This is easily proved for the special choices ${\it\alpha}=0$ and $-2$ (so $\text{d}Z/\text{d}r=0$ ) (Ilin & Morgulis Reference Ilin and Morgulis2013), but is only suggested by numerical evidence generally. This observation is significant because it forces a non-standard singular perturbation analysis in which the additional flow component for ${\it\eta}\neq 0$ has to contribute at leading order to help satisfy one of these two boundary conditions rather than just ‘fixing up’ the third boundary condition.

The analysis of the boundary layer instability is exactly the same as the ${\it\Omega}=1/r^{2}$ case because the presence of a vorticity gradient does not affect the boundary layer problem at leading order. So the growth rate $\text{Re}({\it\sigma})$ scales like $\sqrt{((-m{\it\Omega}^{\prime }(a){\it\eta})/2)/a}$ with azimuthal wavenumber $m$ , local shear at the boundary $-{\it\Omega}^{\prime }(a)>0$ , and radial flow ${\it\eta}$ . The key point is that the boundary layer instability only depends on the shear at the inflow boundary.

The analysis for the critical layer instability also proceeds in a similar fashion, with one key difference: the problem for the large-scale flow is now complicated by a non-vanishing vorticity gradient, but this only has a secondary effect: $\text{d}Z/\text{d}r$ affects the growth rate at $O({\it\eta})$ rather than at the leading $O({\it\eta}\log 1/{\it\eta})$ level. The starting point is the ansatz

(3.35) $$\begin{eqnarray}{\it\sigma}=\text{i}{\it\sigma}_{i}+{\it\delta}\tilde{{\it\sigma}}_{r}\sqrt{{\it\chi}/R_{s}},\end{eqnarray}$$

where ${\it\sigma}_{i}=-m{\it\Omega}(R_{s})$ and $R_{s}$ is the position of the critical layer, ${\it\chi}\,:=\,\,-m{\it\Omega}^{\prime }(R_{s})/2$ and ${\it\delta}\ll \sqrt{{\it\eta}}$ . In the critical layer

(3.36) $$\begin{eqnarray}(-{\it\delta}\tilde{{\it\sigma}}_{r}/\sqrt{{\it\eta}}+2\text{i}{\it\xi}+\partial _{{\it\xi}})\hat{{\it\psi}}_{{\it\xi}{\it\xi}}=0\end{eqnarray}$$

to leading order, where ${\it\xi}:=\sqrt{R_{s}{\it\chi}/{\it\eta}}(r-R_{s})$ . This has the arbitrarily normalised solution

(3.37) $$\begin{eqnarray}\hat{{\it\psi}}_{{\it\xi}{\it\xi}}=\text{e}^{{\it\delta}\tilde{{\it\sigma}}{\it\xi}/\sqrt{{\it\eta}}-\text{i}{\it\xi}^{2}}.\end{eqnarray}$$

Matching this to ${\it\psi}_{rr}$ outside the critical layer requires the WKB solution

(3.38) $$\begin{eqnarray}\displaystyle {\it\psi}^{WKB}=\frac{-R_{s}{\it\chi}\,{\it\eta}}{m^{2}r^{2}[{\it\Omega}(r)-{\it\Omega}(R_{s})]^{2}}\exp \left[\frac{\text{i}m}{{\it\eta}}\int _{R_{s}}^{r}\bar{r}({\it\Omega}(\bar{r})-{\it\Omega}(R_{s}))\,\text{d}\bar{r}+\sqrt{\frac{{\it\chi}}{R_{s}}}\frac{{\it\delta}\tilde{{\it\sigma}}(r^{2}-R_{s}^{2})}{2{\it\eta}}\right] & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

to exist there. There are two other ‘outer’ large-scale solutions either side of the critical layer which solve

(3.39) $$\begin{eqnarray}[{\it\Omega}(r)-{\it\Omega}(R_{s})]\mathscr{L}{\it\psi}=\frac{{\it\psi}_{0}}{r}\frac{\text{d}Z}{\text{d}r}\end{eqnarray}$$

and together accommodate the jump conditions

(3.40a,b ) $$\begin{eqnarray}[u]_{-}^{+}=0\quad \text{and}\quad [v]_{-}^{+}=\sqrt{\frac{{\rm\pi}R_{s}{\it\chi}}{\text{i}{\it\eta}}}\end{eqnarray}$$

across the critical layer (found by integrating (3.37) twice) and the no-normal velocity boundary conditions at $r=1$ and $r=a$ . Then, as in § 2.1.2, balancing the large-scale component of the azimuthal velocity with that from the WKB solution at $r=a$ furnishes the growth rate and dispersion relation for the frequency. Only the higher-order $O({\it\eta})$ part of the growth rate depends on $\text{d}Z/\text{d}r$ , or indeed ${\it\Omega}(r)$ , whereas the leading $O({\it\eta}\log 1/{\it\eta})$ part does not (note ${\it\chi}$ , which contains ${\it\Omega}^{\prime }(R_{s})$ , needs to be non-zero, but otherwise scales out).

3.4. Viscous asymptotics for $0<{\it\eta}\ll 1$ for the boundary layer instability

The inviscid asymptotics can be generalised to include viscosity, which we do here just for the more dangerous boundary layer instability, as this defines the viscous threshold for instability. In the absence of any other physics, the radial inflow is set by the viscosity present via (3.15). However, to extract the scaling law for the viscous instability threshold for more generally maintained radial inflows, we ignore this connection and assume ${\it\eta}$ can vary independently of $Re$ . It is a simple matter to reinstate this connection later to establish stability or instability for a purely hydrodynamic situation. On this basis, the linearised disturbance equation for ${\it\psi}={\it\psi}(r)\text{e}^{\text{i}m{\it\theta}+{\it\sigma}t}$ is

(3.41) $$\begin{eqnarray}\left({\it\sigma}+\text{i}m{\it\Omega}(r)-\frac{{\it\eta}}{r}\frac{\text{d}}{\text{d}r}\right)\mathscr{L}{\it\psi}=\frac{\text{i}m}{r}\frac{\text{d}Z}{\text{d}r}{\it\psi}+\frac{1}{Re}\mathscr{L}^{2}{\it\psi}.\end{eqnarray}$$

Setting ${\it\sigma}+\text{i}m{\it\Omega}(a)={\it\epsilon}\hat{{\it\sigma}}$ , with $\hat{{\it\sigma}}=O(1)$ , and where ${\it\epsilon}$ is the width of the boundary layer at $r=a$ , and balancing the frequency, inflow and viscous terms (Gallet et al. Reference Gallet, Doering and Spiegel2010)

(3.42) $$\begin{eqnarray}\hat{{\it\sigma}}/{\it\epsilon}\sim {\it\eta}/{\it\epsilon}^{3}\sim \frac{1}{{\it\epsilon}^{4}Re}\end{eqnarray}$$

leads to ${\it\epsilon}=Re^{-1/3}$ and ${\it\eta}=\hat{{\it\eta}}{\it\epsilon}^{2}$ , where $\hat{{\it\eta}}=O(1)$ . Defining the boundary layer variable

(3.43) $$\begin{eqnarray}y\,:=\,\,\frac{\sqrt{{\it\chi}}}{\sqrt{\hat{{\it\eta}}/a}}\left(\frac{a-r}{{\it\epsilon}}\right)\end{eqnarray}$$

with, recall, ${\it\chi}\!\,:=\,\,\!-m{\it\Omega}^{\prime }(a)/2$ and adopting the standard boundary layer decomposition ${\it\psi}={\it\psi}_{0}+\hat{{\it\psi}}_{0}+{\it\epsilon}({\it\psi}_{1}+\hat{{\it\psi}}_{1})+\cdots$ , leads to the boundary layer equation at $r=a$ (the boundary layer at $r=1$ is very weak and can be ignored at leading order)

(3.44) $$\begin{eqnarray}\left[\frac{\text{d}^{2}}{\text{d}y^{2}}-N\frac{\text{d}}{\text{d}y}-N(\hat{{\it\sigma}}_{0}+2\text{i}y)\right]\frac{\text{d}^{2}\hat{{\it\psi}}_{0}}{\text{d}y^{2}}=0\end{eqnarray}$$

with

(3.45a,b ) $$\begin{eqnarray}N:=\frac{(\hat{{\it\eta}}/a)^{3/2}}{{\it\chi}^{1/2}}\quad \text{and}\quad \hat{{\it\sigma}}_{0}:=\hat{{\it\sigma}}/(N{\it\chi}^{2})^{1/3}.\end{eqnarray}$$

This is exactly equation (36) in Gallet et al. (Reference Gallet, Doering and Spiegel2010) if their $\text{`}a\text{'}:=-\hat{{\it\sigma}}_{0}$ . The solution which vanishes for $y\rightarrow \infty$ is

(3.46) $$\begin{eqnarray}\frac{\text{d}^{2}\hat{{\it\psi}}_{0}}{\text{d}y^{2}}=\text{e}^{Ny/2}\text{Ai}\left[\frac{N^{1/3}\hat{{\it\sigma}}_{0}}{(2\text{i})^{2/3}}+\frac{N^{4/3}}{4(2\text{i})^{2/3}}+(2\text{i})^{1/3}N^{1/3}y\right],\end{eqnarray}$$

where Ai is the Airy function. Imposing the further non-slip condition $\text{d}\hat{{\it\psi}}_{0}/\text{d}y=0$ at $y=0$ defines the dispersion relation

(3.47) $$\begin{eqnarray}\int _{0}^{\infty }\text{e}^{Ny/2}\text{Ai}\left[\frac{N^{1/3}\hat{{\it\sigma}}_{0}}{(2\text{i})^{2/3}}+\frac{N^{4/3}}{4(2\text{i})^{2/3}}+(2\text{i})^{1/3}N^{1/3}y\right]\text{d}y=0.\end{eqnarray}$$

The onset of instability is found for $N=N_{c}$ where $\text{Re}(\hat{{\it\sigma}}_{0})$ passes through zero (in the positive direction). Numerically, it is better to solve (3.44) directly as an eigenvalue problem rather than treat this integral equation (see appendix A). We find $N_{c}=4.57557$ and $\hat{{\it\sigma}}_{0}=-5.63551\text{i}$ , which is consistent with Gallet et al. (Reference Gallet, Doering and Spiegel2010): the threshold radial flow rate for instability is

(3.48) $$\begin{eqnarray}{\it\eta}_{crit}=aN_{c}^{2/3}{\it\chi}^{1/3}Re^{-2/3}.\end{eqnarray}$$

Since the viscously induced radial flow in an accretion disk is $O(Re^{-1})\ll {\it\eta}_{crit}$ , the instability is not expected to be triggered unless it is substantially subcritical. Whether this is the case will be considered below in § 6 after we examine the effect of introducing some extra physics to the instability.

5.1. 3D extension of the boundary layer instability

Working with the primitive variables, the appropriate scalings to capture the boundary layer instability are

(5.7a,b ) $$\begin{eqnarray}{\it\sigma}=-\text{i}m{\it\Omega}(a)+\hat{{\it\sigma}}\sqrt{\frac{-m{\it\eta}{\it\Omega}^{\prime }(a)}{2a}}+\cdots \,,\quad {\it\xi}:=(a-r)\sqrt{\frac{-ma{\it\Omega}^{\prime }(a)}{2{\it\eta}}},\end{eqnarray}$$

(5.8) $$\begin{eqnarray}(u,v,w,p)=\left(\sqrt{m{\it\eta}}\,\hat{u} ,v,{\hat{w}}\sqrt{\frac{2{\it\Omega}(a)}{Z(a)}},\sqrt{\frac{{\it\eta}}{m}}\,\hat{p}\right)\end{eqnarray}$$

in the boundary layer ${\it\xi}=O(1)$ , where

(5.9) $$\begin{eqnarray}{\it\mu}:=\frac{k}{m}\end{eqnarray}$$

measures the degree of three-dimensionality. With these, the equation set (5.1)–(5.4) reduces to

(5.10) $$\begin{eqnarray}\displaystyle & \displaystyle 0=2{\it\Omega}v+\sqrt{-\frac{1}{2}a{\it\Omega}^{\prime }(a)}\frac{\text{d}\hat{p}}{\text{d}{\it\xi}}, & \displaystyle\end{eqnarray}$$

(5.11) $$\begin{eqnarray}\displaystyle & \displaystyle \sqrt{\frac{-{\it\Omega}^{\prime }(a)}{2a}}\left(\hat{{\it\sigma}}+2\text{i}{\it\xi}+\frac{\text{d}}{\text{d}{\it\xi}}\right)v=-Z(a)\hat{u} -\frac{\text{i}}{a}\hat{p}, & \displaystyle\end{eqnarray}$$

(5.12) $$\begin{eqnarray}\displaystyle & \displaystyle \sqrt{\frac{2{\it\Omega}(a)}{Z(a)}}\sqrt{\frac{-{\it\Omega}^{\prime }(a)}{2a}}\left(\hat{{\it\sigma}}+2\text{i}{\it\xi}+\frac{\text{d}}{\text{d}{\it\xi}}\right){\hat{w}}=-\text{i}{\it\mu}\hat{p}, & \displaystyle\end{eqnarray}$$

(5.13) $$\begin{eqnarray}\displaystyle & \displaystyle -\sqrt{-\frac{1}{2}a{\it\Omega}^{\prime }(a)}\frac{\text{d}\hat{u} }{\text{d}{\it\xi}}+\frac{\text{i}v}{a}+\text{i}{\it\mu}\sqrt{\frac{2{\it\Omega}(a)}{Z(a)}}{\hat{w}}=0 & \displaystyle\end{eqnarray}$$

$\sqrt{m{\it\eta}}$

$\hat{u}$

$\hat{p}$

(5.14) $$\begin{eqnarray}\displaystyle & \displaystyle \left(\hat{{\it\sigma}}+2\text{i}{\it\xi}+\frac{\text{d}}{\text{d}{\it\xi}}\right)\frac{\text{d}v}{\text{d}{\it\xi}}=-\text{i}{\it\gamma}{\hat{w}}, & \displaystyle\end{eqnarray}$$

(5.15) $$\begin{eqnarray}\displaystyle & \displaystyle \left(\hat{{\it\sigma}}+2\text{i}{\it\xi}+\frac{\text{d}}{\text{d}{\it\xi}}\right)\frac{\text{d}{\hat{w}}}{\text{d}{\it\xi}}=\text{i}{\it\gamma}v-2\text{i}{\hat{w}}, & \displaystyle\end{eqnarray}$$

(5.16) $$\begin{eqnarray}{\it\gamma}\,:=\,\,2a{\it\mu}\frac{\sqrt{2({\it\alpha}+2)}}{(-{\it\alpha})}.\end{eqnarray}$$

This reduced system only involves $m$ through ${\it\mu}$ and is valid for $m\ll 1/{\it\eta}$ . So, for ${\it\mu}$ held fixed, the growth rate scales with $\sqrt{m}$ (see the rescaling in (5.7)) across the astrophysical regime of $1\ll m\ll 1/{\it\eta}$ .

In the particular case treated by Ilin & Morgulis (Reference Ilin and Morgulis2013) ( ${\it\alpha}=-2$ so ${\it\gamma}=0$ ) adding an axial wavenumber does nothing to the 2D instability provided ${\it\eta}k$ stays small. In the general case of interest, $-2\leqslant {\it\alpha}<0$ and ${\it\mu}>0$ (clearly the stability problem is symmetric under the transformation $({\it\mu},v,{\hat{w}})\rightarrow (-{\it\mu},v,-{\hat{w}})$ , so only ${\it\mu}>0$ needs be considered), a solution is available in terms of Kummer functions. Equations (5.14) and (5.15) can be combined to a single equation for $v$

(5.17) $$\begin{eqnarray}(\hat{\mathscr{G}}^{2}+2\text{i}\hat{\mathscr{G}}-{\it\gamma}^{2})v=0,\end{eqnarray}$$

where

(5.18) $$\begin{eqnarray}\hat{\mathscr{G}}:=\frac{\text{d}^{2}}{\text{d}{\it\xi}^{2}}+(\hat{{\it\sigma}}+2\text{i}{\it\xi})\frac{\text{d}}{\text{d}{\it\xi}}.\end{eqnarray}$$

This can be factored as $(\hat{\mathscr{G}}+4\text{i}{\it\Gamma}_{1})(\hat{\mathscr{G}}+4\text{i}{\it\Gamma}_{2})v=0$ , where ${\it\Gamma}_{1}:=(1-\sqrt{1-{\it\gamma}^{2}})/4$ and ${\it\Gamma}_{2}:=(1+\sqrt{1-{\it\gamma}^{2}})/4$ . Under the change of variable $z:=\text{i}(\hat{{\it\sigma}}+2\text{i}{\it\xi})^{2}/4$ , each factor is just Kummer’s differential equation,

(5.19) $$\begin{eqnarray}z\frac{\text{d}^{2}v}{\text{d}{\it\xi}^{2}}+\left(\frac{1}{2}-z\right)\frac{\text{d}v}{\text{d}{\it\xi}}-{\it\Gamma}_{j}v=0\quad j=1,2.\end{eqnarray}$$

The solution for $v$ (and therefore indirectly ${\hat{w}}$ ) which decays exponentially for ${\it\xi}\rightarrow \infty$ and $\hat{{\it\sigma}}_{r}>0$ is

(5.20) $$\begin{eqnarray}v=\text{e}^{z}(\mathscr{A}U({\textstyle \frac{1}{2}}-{\it\Gamma}_{1},{\textstyle \frac{1}{2}},-z)+\mathscr{B}U({\textstyle \frac{1}{2}}-{\it\Gamma}_{2},{\textstyle \frac{1}{2}},-z)),\end{eqnarray}$$

where $U$ is the multivalued Kummer’s function (see NIST Reference Olver, Lozier, Boisvert and Clark2010, § 13.2.25) and $\mathscr{A}$ and $\mathscr{B}$ are constants. If these are not both to vanish when the further boundary conditions $v={\hat{w}}=0$ are imposed at ${\it\xi}=0$ , then the product

(5.21) $$\begin{eqnarray}U({\textstyle \frac{1}{2}}-{\it\Gamma}_{1},{\textstyle \frac{1}{2}},-{\textstyle \frac{1}{4}}\text{i}\hat{{\it\sigma}}^{2})U({\textstyle \frac{1}{2}}-{\it\Gamma}_{2},{\textstyle \frac{1}{2}},-{\textstyle \frac{1}{4}}\text{i}\hat{{\it\sigma}}^{2})=0,\end{eqnarray}$$

which defines unstable ( $\hat{{\it\sigma}}_{r}>0$ ) eigenfrequencies $\hat{{\it\sigma}}$ . The eigenvalues for ${\it\mu}>0$ which smoothly connect with those at ${\it\mu}=0$ are given by $U(1/2-{\it\Gamma}_{1},1/2,-\text{i}\hat{{\it\sigma}}^{2}/4)=0$ , since ${\it\mu}\rightarrow 0$ implies ${\it\Gamma}_{1}\rightarrow 0$ and convergence to the 2D equation $\hat{\mathscr{G}}v=0$ . These eigenfunctions can be plotted using

(5.22) $$\begin{eqnarray}v=\left\{\begin{array}{@{}ll@{}}\displaystyle 2\frac{{\it\Gamma}({\textstyle \frac{1}{2}})}{{\it\Gamma}(1-{\it\Gamma}_{1})}\text{e}^{z}M\!\left(\frac{1}{2}-{\it\Gamma}_{1},\frac{1}{2},-z\right)\!-\text{e}^{z}U\!\left(\frac{1}{2}-{\it\Gamma}_{1},\frac{1}{2},-z\right)\quad & 0\leqslant {\it\xi}<-\displaystyle \frac{1}{2}(\hat{{\it\sigma}}_{r}+\hat{{\it\sigma}}_{i})\\ \text{e}^{z}U({\textstyle \frac{1}{2}}-{\it\Gamma}_{1},{\textstyle \frac{1}{2}},-z)\quad & -{\textstyle \frac{1}{2}}(\hat{{\it\sigma}}_{r}+\hat{{\it\sigma}}_{i})\leqslant {\it\xi}\end{array}\right.\end{eqnarray}$$

and ${\hat{w}}=2{\it\Gamma}_{1}v/{\it\gamma}$ (with $z:=\text{i}(\hat{{\it\sigma}}+2\text{i}{\it\xi})^{2}/4$ and ‘ ${\it\Gamma}(z)$ ’ indicating the Gamma function), which compensates for the branch cut along the negative real axis in $U$ routinely imposed by packages such as Maple ( $M$ is the single-valued Kummer function: § 13.1.2, Abramowitz & Stegun Reference Abramowitz and Stegun1964): see figure 4.

Figure 4. The leading growing $v$ eigenfunction calculated using (5.22) in Maple for $m=1$ , $k=0.3$ , ${\it\alpha}=-3/2$ , corresponding to $\hat{{\it\sigma}}=0.64021-4.59526\text{i}$ (real/imaginary part is red solid/blue dashed line).

Equations (5.14) and (5.15) can also be directly treated numerically (see appendix A). This boundary layer approach works well for the dominant instabilities as ${\it\mu}$ increases from zero, showing how they are systematically suppressed until their growth rates are comparable to the interior critical layer modes: see figure 5, which shows this for the leading instability. At this point (which is ${\it\mu}\approx 0.4$ for $a=2,{\it\eta}=10^{-4},{\it\alpha}=-1.5$ ) the full 3D eigenvalue (5.1)–(5.4) must be solved (see appendix A), which shows the complete suppression of all the instabilities by ${\it\mu}=1$ : see figure 6.

Figure 5. Leading inviscid 3D instability for ${\it\alpha}=-1.5$ , $m=1$ , $a=2$ with scaled growth rate $\hat{{\it\sigma}}_{r}={\it\sigma}_{r}/\sqrt{((-m{\it\eta}{\it\Omega}^{\prime }(a))/2)/a}$ plotted against ${\it\mu}$ , the degree of three-dimensionality, as computed using the full eigenvalue code – ${\it\eta}=10^{-4}$ uppermost (blue) line with triangles, ${\it\eta}=10^{-5}$ second uppermost (red) line with squares, ${\it\eta}=10^{-6}$ third uppermost (green) line with diamonds and ${\it\eta}=10^{-7}$ lowest (cyan) line with filled circles – where $N=1000$ proves sufficient (the markers on each line actually show the $N=500$ stability results for ${\it\mu}=0.4,0.45$ and 0.5 to demonstrate convergence). The leading instability is also shown for the boundary layer problem (5.14) and (5.15) as a dashed black line ( $N=4000$ ). The boundary layer scaling works well, with the full eigenvalue prediction smoothly converging to the boundary layer prediction as ${\it\eta}\rightarrow 0$ for ${\it\mu}\lesssim 0.3$ , but beyond this the $\sqrt{{\it\eta}}$ scaling is no longer correct. This plot demonstrates that increasing three-dimensionality acts to suppress the 2D instability.

Figure 6. Inviscid 3D instabilities for ${\it\alpha}=-1.5$ , $m=1$ , $a=2$ and ${\it\eta}=3\times 10^{-4}$ with scaled growth rate $\hat{{\it\sigma}}_{r}={\it\sigma}_{r}/\sqrt{-m{\it\eta}{\it\Omega}^{\prime }(a)/(2a)}$ plotted against frequency ${\it\sigma}_{i}$ for various values of ${\it\mu}$ . The right (left) vertical dashed line is ${\it\sigma}_{i}=-m{\it\Omega}(a)$ ( ${\it\sigma}_{i}=-m{\it\Omega}(1)$ ). The eigenvalues are calculated from a full 3D eigenvalue calculation with $N=2400$ for (in order downwards from the uppermost 2D eigenvalues): ${\it\mu}=0$ (blue dots); ${\it\mu}=0.35$ (red dots); ${\it\mu}=0.4$ (magenta dots); ${\it\mu}=0.43$ (green dots); ${\it\mu}=0.45$ (cyan dots); ${\it\mu}=0.5$ (black dots); ${\it\mu}=0.6$ (blue dots with squares); ${\it\mu}=0.75$ (red dots with diamonds) and ${\it\mu}=1$ (magenta dots with circles). This plot demonstrates that increasing three-dimensionality suppresses the 2D instability (in fact completely by ${\it\mu}=0.75$ here). Notice that numerical errors start to creep into the eigenvalues whose eigenfunctions have critical layers far from the inflow boundary (the developing corrugations in the curves) as ${\it\mu}$ increases (not unexpected as the numerical truncation $N$ is only $O(1/{\it\eta})$ ).

The general conclusion is that 3D disturbances are less unstable than 2D disturbances under boundary inflow. In fact, since the boundary layer instability depends only on the shear at the boundary, this could have been anticipated by using a Squires transformation to map a 3D disturbance to a more unstable 2D disturbance (Squires Reference Squires1933). Normally, such a transformation fails in a rotating system, due to curvature terms, but here these are marginalised by the dispersion relation being only sensitive to the shear at the boundary.

6.1. Weakly nonlinear analysis of 2D viscous boundary layer instability

The analysis revolves almost completely around the boundary layer and can be handled relatively straightforwardly: see appendix B. The growing instability is found to be supercritical, so that the branch of 2D solutions exists for ${\it\eta}\geqslant {\it\eta}_{crit}$ . Along this branch of solutions, the azimuthally averaged radial pressure drop, ${\rm\Delta}p$ , across the domain, which is a more natural control parameter for a laboratory experiment, changes. Furthermore, in an accretion disk, maintaining constant angular momentum of the flow is a more natural constraint under which to study flow bifurcations. These, however, just represent different perspectives of viewing the same bifurcation and the ensuing 2D branch of solutions. For example, solutions with constant ${\rm\Delta}p$ and varying ${\it\eta}$ are the same as those with constant ${\it\eta}$ and varying ${\rm\Delta}p$ , provided the rest of the boundary conditions are consistent (e.g. azimuthally asymmetric radial flow components vanish at the boundaries in both, and the azimuthally asymmetric pressure distribution on each boundary is unconstrained in both). To make this clear and to be able to make statements away from the vicinity of the bifurcation point, we now compute the fully nonlinear 2D solution branch.

6.2. Fully nonlinear 2D solutions

To study (non-helical) 2D bifurcations off the 1D steady state (3.1), a Reynolds number of $Re=10^{4}$ is chosen as a compromise, being hopefully large enough to be in the asymptotic regime but not too large that flows become too arduous to follow numerically. According to the asymptotics, ${\it\eta}_{crit}=aN_{c}^{2/3}{\it\chi}^{1/3}Re^{-2/3}\approx 0.0061\,(0.0087)$ for $m_{0}=1$ ( $m_{0}=3$ ), whereas actually the thresholds ${\it\eta}_{crit}=0.009635\,(0.0136)$ are found numerically at $Re=10^{4}$ ( ${\it\alpha}=-3/2$ and $a=2$ ). Figure 7 shows how the spectrum of the linear operator depends on ${\it\eta}$ for $Re=10^{5}$ . Instability is first possible at $Re\approx 5\times 10^{3}$ (not shown), with the second mode of instability appearing for $Re>10^{4}$ , and by $Re=10^{5}$ there are three unstable modes. Further computations for $Re=10^{6}$ (not shown) show that further modes of instability emerge (now eight) and the persistent feature that each unstable mode is only unstable for a finite range of ${\it\eta}$ . The lower threshold ${\it\eta}_{l}$ (which is the focus in this work and Gallet et al. (Reference Gallet, Doering and Spiegel2010)) must tend to zero as $Re\rightarrow \infty$ to be consistent with the inviscid analysis, whereas the upper threshold ${\it\eta}_{u}$ must tend to a finite limit to be consistent with the analysis of Ilin & Morgulis (Reference Ilin and Morgulis2015a ).

The bifurcation is a Hopf bifurcation and the oscillation can be made to look steady by going into a frame rotating at the phase speed $c_{{\it\theta}}\,:=\,\,{\it\sigma}_{i}/m_{0}$ at the bifurcation point. Subsequently, the 2D solution branch is traced out by using the representation

(6.1) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}u\\ v\\ p\end{array}\right]:=\mathop{\sum }_{n=0}^{N-1}\mathop{\sum }_{m=0}^{M}\left[\begin{array}{@{}c@{}}u_{mn}(t)(T_{n+2}({\it\xi})-T_{n}({\it\xi}))\\ v_{mn}(t)(T_{n+2}({\it\xi})-T_{n}({\it\xi}))\\ p_{mn}(t)T_{n}({\it\xi})\end{array}\right]\text{e}^{\text{i}mm_{0}({\it\theta}-c_{{\it\theta}}t)}+\text{c.c.},\end{eqnarray}$$

where $T_{n}(x)\,:=\,\,\cos (n\cos ^{-1}x)$ are Chebyshev polynomials, ${\it\xi}\,:=\,\,(2r-a-1)/(a-1)$ , $m_{0}$ indicates the rotational symmetry of the bifurcating eigenfunction and $c_{{\it\theta}}$ is the azimuthal phase speed, which is found as part of solution process. Figure 8 shows the new branch of 2D solutions arising from the one unstable mode present at $Re=10^{4}$ for $m_{0}=1$ and $m_{0}=3$ found by a Newton–Raphson root-finding algorithm with truncations varying from $(M,N)=(20,60)$ to $(10,150)$ (so $O(10^{4})$ degrees of freedom). Both bifurcations at the lower threshold ${\it\eta}_{l}$ are supercritical, consistent with the weakly nonlinear analysis of § 6.1, and the solution branches reconnect with the 1D solution (3.1) at the upper threshold ${\it\eta}_{u}$ : see figure 8.

Along the 2D solution branch, the ‘surplus’ azimuthally averaged radial pressure drop,

(6.2) $$\begin{eqnarray}{\it\delta}p:=({\rm\Delta}p)_{2D}-({\rm\Delta}p)_{1D}=\int _{1}^{a}2{\it\Omega}(r)\overline{v}+\frac{\overline{v^{2}-u^{2}}}{r}\,\text{d}r,\end{eqnarray}$$

(where $\overline{(\cdot )}:=1/2{\rm\pi}\int _{0}^{2{\rm\pi}}(\cdot )\,\text{d}{\it\theta}$ is just an azimuthal average), the surplus of 2D rotational energy compared to the 1D solution,

(6.3) $$\begin{eqnarray}{\it\delta}E:=2{\rm\pi}\int _{1}^{a}(r{\it\Omega}\overline{v}+1/2\overline{v}^{2})r\,\text{d}r\end{eqnarray}$$

and the disturbance angular momentum ${\it\delta}I$ are all positive quantities: see figure 9. The initial increase in the pressure drop and the total angular momentum indicates that if either were used as a control parameter instead of the azimuthally averaged radial velocity, the bifurcation would remain supercritical. For example, figure 11 shows a plot of ${\rm\Delta}p$ against the azimuthally averaged radial flow parameter ${\it\eta}$ for the 2D solution branch already plotted in figure 8. The weakly nonlinear analysis in appendix B takes ${\it\eta}$ as the control parameter and finds that the bifurcated 2D solution branch exists for ${\it\eta}>{\it\eta}_{crit}$ (marked as the point A in figure 11) with ${\rm\Delta}p$ larger for the 2D solutions than the equivalent (same ${\it\eta}$ ) 1D solution (see points B and C in figure 11). If ${\rm\Delta}p$ is the control parameter, the bifurcation is still supercritical, as 2D solutions can only exist for ${\rm\Delta}p\geqslant {\rm\Delta}p$ at A, but now the bifurcated 2D solutions give rise to two smaller radial inflow solutions (points C or D in figure 11) compared to the 1D solution E.

Figure 10. (a)  $\int _{0}^{2{\rm\pi}}r^{2}uv\,\text{d}{\it\theta}$ against $r$ for the bifurcating eigenfunction at ${\it\eta}=9.645\times 10^{-3}$ (left dot on figure 8(a)). (b) The mean angular velocity profile for the $m_{0}=1$ 2D solution at ${\it\eta}=3.67\times 10^{-2}$ $Re=10^{4}$ , $a=2$ and ${\it\alpha}=-3/2$ (solid blue line) and the profile ${\it\Omega}=r^{{\it\alpha}}$ corresponding to the 1D basic state (red dashed line). This plot clearly indicates that angular momentum has been transported outwards by the saturated instability.

Another interesting issue for accretion disks is whether this instability acts to transfer angular momentum $I$ outwards. Forming the integral $\int _{0}^{2{\rm\pi}}r^{2}$  (3.4)  $\text{d}{\it\theta}$ gives the conservation equation

(6.4) $$\begin{eqnarray}\left(\mathscr{I}(r):=\int _{0}^{2{\rm\pi}}r[r{\it\Omega}(r)+v]r\,\text{d}{\it\theta}\right)_{t}+J_{r}=0,\end{eqnarray}$$

where $I=\int _{1}^{a}\mathscr{I}(r)\,\text{d}r$ and

(6.5) $$\begin{eqnarray}J:=\int _{0}^{2{\rm\pi}}\left\{ruv-{\it\eta}v-\frac{1}{Re}r^{2}\frac{\text{d}}{\text{d}r}\left(\frac{v}{r}\right)\right\}r\,\text{d}{\it\theta}\end{eqnarray}$$

is the associated radial flux of angular momentum (the flux vanishes for the 1D basic state (3.1)). The first (Reynolds-stress) term in $J$ is computable using the bifurcating eigenfunction alone and can be used to determine the effect of the instability on the angular momentum distribution (the other two terms, which involve nonlinear aspects of the instability, subsequently balance the first term to give a finite amplitude 2D steady state: figure 9(b) actually shows that $r^{2}\overline{uv}\approx {\it\eta}r\overline{v}$ ). Figure 10 plots this term over the radius for the $m_{0}=1$ state just after the bifurcation (shown in figure 8(a) as the leftmost point at ${\it\eta}=9.645\times 10^{-3}$ ) and shows that it alternates in sign but is mostly positive, indicating outwards angular momentum transport. The angular velocity of the 2D state at its maximum amplitude (the middle $m_{0}=1$ point in figure 8(a) and the middle cross-section in figure 8(c)) provides more definitive evidence of this outward transport: see figure 10.

Figure 11. The azimuthally averaged radial pressure drop ${\rm\Delta}p$ is plotted again the azimuthally averaged radial flow ${\it\eta}$ for the undisturbed 1D basic state (blue solid line connecting points A, B and E; ${\rm\Delta}p_{1D}:=(1-1/a)+({\it\eta}^{2}(1-1/a^{2}))/2$ ) and the 2D solutions (red solid line connecting points A, C and D) which emanate out of the bifurcation for $m_{0}=1$ , $a=2$ , $Re=10^{4}$ and ${\it\alpha}=-3/2$ (as shown in figure 8). The bifurcation point is point A, and the plot shows that the bifurcation is supercritical for either ${\it\eta}$ or ${\rm\Delta}p$ being the control parameter. Note that, if ${\rm\Delta}p$ is fixed, the bifurcation leads to states with lower radial flow (compare states C and D with E).

Figure 12. The 3D solution branches continued out of the six Hopf bifurcations (see the inset, which magnifies the dashed box for detail) from the 2D branch (thick blue loop as in figure 8) found for $m_{0}=1$ at $Re=10^{4}$ , ${\it\alpha}=-3/2$ , $k=1$ and $a=2$ . All the bifurcations are supercritical. The ordinate is the disturbance energy (2D and 3D) and the abscissa the radial flow ${\it\eta}$ . The black dot indicates the particular 3D flow state shown in detail in figure 13. Typically a resolution of $(N,M,L)=(40,20,4)$ was used to follow these solutions, with curve segments only shown if they are robust under truncation changes.

6.3. Fully nonlinear 3D solutions

To briefly explore the possibility of 3D states, 3D bifurcations were sought off the 2D branch of $m_{0}=1$ solutions. Concentrating on 3D states with an axial wavenumber $k=1$ , six Hopf bifurcations were found between the initial 2D bifurcation point ${\it\eta}={\it\eta}_{2D}=0.0096$ and ${\it\eta}=0.0152$ . These were then traced using the fully nonlinear steady representation

(6.6) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}u\\ v\\ w\\ p\end{array}\right]:=\mathop{\sum }_{n=0}^{N-1}\mathop{\sum }_{m=-M}^{M}\mathop{\sum }_{l=0}^{L}\left[\begin{array}{@{}c@{}}u_{lmn}(t)(T_{n+2}({\it\xi})-T_{n}({\it\xi}))\\ v_{lmn}(t)(T_{n+2}({\it\xi})-T_{n}({\it\xi}))\\ w_{lmn}(t)(T_{n+2}({\it\xi})-T_{n}({\it\xi}))\\ p_{lmn}(t)T_{n}({\it\xi})\end{array}\right]\text{e}^{\text{i}(mm_{0}{\it\theta}^{\ast }+lkz^{\ast })}+\text{c.c.}\end{eqnarray}$$

by moving into a frame moving with the phase speed $c_{z}$ in $z$ (initially $c_{z}:={\it\sigma}_{i}/k$ at the bifurcation point, where ${\it\sigma}_{i}$ is the instability frequency). Since ${\it\theta}^{\ast }:={\it\theta}-c_{{\it\theta}}t$ already incorporates the phase speed of the initial 2D instability, these secondary 3D states are steady in a rotating and translating frame. Figure 12 shows that all the traced bifurcations are supercritical, i.e. there is no evidence of solution branches reaching ${\it\eta}<{\it\eta}_{2D}$ . Solution branches are shown as far as they are reproducible using different truncation levels. The maximum realistic resolution was $(N,M,L)=(40,20,4)$ , giving 58 963 degrees of freedom, since the branch continuation approach was a direct Newton–Raphson solver, albeit with multithreaded linear algebra software. A typical 3D flow is shown in figure 13, with the presence of vertical jets at the outflow boundary and the lack of any large-scale structure being particularly noteworthy.

Figure 13. A typical 3D state (shown as a dot on figure 12) calculated using $(N,M,L)=(40,28,4)$ (82 003 degrees of freedom) plotted at $z=0$ (a),  $z={\rm\pi}/2$  (b),  $z={\rm\pi}$ (c) and $z=3{\rm\pi}/2$ (d) (axial wavelength $=2{\rm\pi}$ ). The axial speed is shown by 10 contours going from $-0.025$ to 0.035 (the colour/shading of the zero contour is shown on the outside of the flow domain). The arrows indicate the speed and direction of the flow in the $(r,{\it\theta})$ plane, with the longest arrow representing a speed of 0.272.