1. Introduction
A predictive model of the transition from laminar flow to turbulence in viscous parallel shear flows has remained elusive since seminal empirical observations of turbulence transition in the late 19th century (Reynolds Reference Reynolds1883). One of the earliest and most thoroughly explored approaches to this problem examines the linear stability of the Navier–Stokes equations when subjected to small velocity perturbations (Drazin & Reid Reference Drazin and Reid2004). Evaluation of the temporal growth of the perturbations typically proceeds by solution of the Orr–Sommerfeld eigenvalue equation (see Appendix A), which depends on the Reynolds number of the flow, 
$Re = U_0L_0/\nu _0$, where 
$U_0$ and 
$L_0$ are characteristic flow speed and length scales, respectively, and 
$\nu _0$ is the kinematic viscosity of the fluid. Velocity perturbations are constrained to satisfy a no-slip condition at the fluid–solid interfaces. The eigenvalue problem conventionally incorporates this requirement into the homogeneous boundary condition:
\begin{equation} \tilde{\boldsymbol{u}}(\boldsymbol{x}_{wall},t) = 0, \end{equation}
where 
$\tilde {\boldsymbol {u}}$ is the perturbation velocity vector amplitude, and 
$\boldsymbol {x}_{wall}$ is the location of the fluid–solid interfaces.
The aforementioned analysis predicts that plane Couette flow is stable to all two-dimensional linear perturbations at all Reynolds numbers (Davey Reference Davey1973). This prediction is contradicted by empirical observations of transition to turbulence at channel Reynolds numbers as low as 
$Re \approx 360$ (Drazin & Reid Reference Drazin and Reid2004). The analysis method similarly finds no linear instability to account for the observed transition of Hagen–Poiseuille pipe flow to turbulence at a Reynolds number 
$Re \approx 2000$. For plane Poiseuille flow, linear stability analysis does identify a single unstable eigenmode, which first appears at a Reynolds number 
$Re \approx 5772$ (Orszag Reference Orszag1971). In practice, however, transition to turbulence is observed to occur at significantly lower Reynolds numbers 
$Re \approx 1000$. Moreover, the growth rate associated with the unstable eigenvalue 
$\hat {\omega }$ at 
$Re \approx 5772$ is relatively weak, e.g. its imaginary part is 
$\mathbb {I}[\hat {\omega }] \approx 0.0037$ for wavenumber 
$\alpha \approx 1.02$. The associated unstable eigenmode may therefore be insufficient to trigger turbulence.
This apparent inability of linear stability theory to accurately predict turbulence transition in a variety of viscous parallel shear flows has motivated the exploration of alternative frameworks to explain and possibly predict turbulence transition. These include consideration of finite-sized velocity perturbations (Orszag & Kells Reference Orszag and Kells1980), nonlinear transition processes (Schmid & Henningson Reference Schmid and Henningson2001), and transient growth mechanisms such as those associated with non-normality of the Orr–Sommerfeld eigenvalue equation (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993). While it is beyond the scope of this paper to comprehensively review the rich literature on the theory of turbulence transition, the reader is referred to an excellent recent review article (Avila, Barkley & Hof Reference Avila, Barkley and Hof2023).
The aforementioned body of work suggests that linear instability may not be a necessary condition for turbulence transition. Nonetheless, in this paper we show that linear stability analysis is sufficient to quantitatively predict the occurrence of flow instabilities in each of three canonical viscous parallel shear flows: plane Couette flow, plane Poiseuille flow and Hagen–Poiseuille pipe flow. We proceed by considering two alternative models for incipient flow perturbations in these canonical flows. In the first case, we use a more general ansatz for the perturbation boundary condition in the Orr–Sommerfeld equation, of which the no-slip condition is a limiting case. In the second case, we model the perturbation behaviour at the wall using a variation of Stokes’ second problem (Landau & Lifschitz Reference Landau and Lifschitz1987). While the latter model represents a more significant departure from the conventional perturbation no-slip condition, it leads to physically realistic, quantitative predictions of a Reynolds-number-dependent transition to linear instability that is consistent with empirical observations, i.e. 
$Re_{crit} \sim O(10^2\unicode{x2013}10^4)$.
Both of the alternative models for incipient perturbations exhibit near-wall behaviour that is, in many cases, empirically indistinguishable from the classical no-slip condition. The superior predictive capability of these new models – without the need to appeal to nonlinear or transient growth processes – suggests a re-examination of the no-slip condition as a sufficiently accurate model of the fluid–solid interface physics. More generally, these results hint at the possibility that more detailed study of flow physics at fluid–solid interfaces could lead to a better understanding of turbulence transition.
Section 2 explores the first alternative model for incipient perturbations. This model includes the conventional no-slip condition as a limiting case, enabling more direct comparison with prior results in the literature. While this feature is pedagogically useful, this perturbation model exhibits an associated trade-off, namely, that the predicted linear instability does not exhibit an explicit Reynolds number dependence. Section 3 examines the second alternative model, which represents a more significant departure from the form of no-slip condition conventionally applied to velocity perturbations in linear stability analysis. This model does successfully predict a Reynolds-number-dependent linear instability that is consistent with empirical observations of turbulence. Section 4 concludes the paper by discussing the implications of these findings for our understanding of the role of fluid–solid interface physics in prediction of turbulence transition.
2. Perturbation model I
2.1. Generalized boundary condition
Let us replace the streamwise component of the no-slip condition in (1.1) with a more general, homogeneous boundary condition that depends on the shear rate of the velocity perturbation at the wall:
\begin{equation} \tilde{u}(\boldsymbol{x}_{wall},t) \mp S\,\tilde{u}'(\boldsymbol{x}_{wall},t) = 0, \end{equation}
where 
$\tilde {u}$ is the streamwise component of the velocity perturbation amplitude, with 
$\tilde {u}' = \partial \tilde {u}/\partial y$ for the plane Couette and Poiseuille flows, and 
$\tilde {u}' = \partial \tilde {u}/\partial r$ for the Hagen–Poiseuille pipe flow. Here, 
$S$ is a characteristic slip length constant (normalized by the channel half-width or pipe radius), and the sign 
$\mp$ applies to the wall at 
$y$ (or 
$r$) 
$= 1$ and 
$y = -1$, respectively.
Physically, the boundary condition in (2.1) requires that any non-zero perturbation velocity at the wall is in the same direction as the corresponding shear exerted on the fluid by the wall. While it is similar in form to the Navier slip boundary condition (Lauga, Brenner & Stone Reference Lauga, Brenner and Stone2007), there are two key differences between this model and prior studies of the effect of wall slip on linear stability (e.g. Lauga & Cossu Reference Lauga and Cossu2005; Chai & Song Reference Chai and Song2019; Ceccacci et al. Reference Ceccacci, Calabretto, Thomas and Denier2022). First, in the present analysis, the boundary condition (2.1) applies only to velocity perturbations; the base flow is assumed to exhibit the conventional no-slip condition. Physically, this assumption is consistent with the ansatz that the distribution of velocity perturbations exhibits a zero mean, and therefore the velocity perturbations do not change the base flow. By contrast, in the aforementioned studies, the base flow is assumed to exhibit wall slip even in the absence of velocity perturbations. While the assumption of non-zero base flow slip is appropriate in the context of those prior studies (e.g. where the walls were assumed to exhibit hydrophobicity or other surface treatment), the goal of the present study is to examine wall-bounded flows more generally.
The second key distinction between the present perturbation model and prior studies is that the perturbation velocity at the wall is applied here in the same direction as the corresponding shear exerted on the fluid by the wall. In contrast, the previous studies apply perturbation wall slip in the opposite direction of the wall shear exerted on the fluid. Physically, the present model of wall slip assumes that the velocity perturbations are caused by wall shear perturbations on the adjacent fluid. This distinction is significant, and we confirmed that the stability predictions that follow do indeed depend on relative directions of the velocity perturbation and the wall shear on the fluid (see Appendix E).
This generalized boundary condition (2.1) reduces to the no-slip condition (1.1) when the slip parameter is 
$S = 0$. The no-slip condition is also satisfied by velocity perturbations with non-zero values of the slip parameter 
$S$ if the corresponding wall-normal gradient of streamwise flow speed is zero at the wall, i.e. 
$\tilde {u}'(\boldsymbol {x}_{wall},t) = 0$. While we presently focus on the case of real-valued slip parameter 
$S$, we have observed similar results for complex 
$S$. In the latter case, a phase difference exists between the perturbation wall slip and the perturbation wall shear.
2.2. Linear stability maps
The Orr–Sommerfeld eigenvalue equation was solved for plane Couette flow, plane Poiseuille flow and Hagen–Poiseuille pipe flow using Chebyshev collocation (Schmid & Henningson Reference Schmid and Henningson2001; Malik & Skote Reference Malik and Skote2019); see Appendix A for governing equations and data repository for Matlab implementations. The boundary condition in (2.1) was evaluated for values of the slip parameter 
$S = 0$ and ranging from 
$S = 10^{-7}$ to 
$S = 1$ in ten equally-spaced increments per decade. The Reynolds number was varied from 
$Re = 50$ to 
$Re = 10\,000$ in increments of 10.
For combinations of slip parameter 
$S$ and Reynolds number 
$Re$ spanning this range, figure 1 plots contours of the maximum eigenvalue imaginary part, i.e. 
$\mathbb {I}[\hat {\omega }]$, for Orr–Sommerfeld solutions corresponding to plane Couette flow, plane Poiseuille flow and Hagen–Poiseuille pipe flow, respectively. (Similar contour maps of the second-largest eigenvalue are provided in Appendix C for reference.)
Figure 1. Contour maps of the maximum Orr–Sommerfeld eigenvalue imaginary part 
$\mathbb {I}[\hat {\omega }]$ versus slip parameter 
$S$ and Reynolds number 
$Re$ for (a) plane Couette flow, (b) plane Poiseuille flow, and (c) Hagen–Poiseuille pipe flow. Blue contours indicate regions of linear stability, and red contours indicate regions of linear instability. Black contours indicate neutral stability boundaries. The Reynolds number 
$Re \approx 5772$ of the single unstable eigenmode in plane Poiseuille flow for 
$S = 0$ is indicated on the ordinate axis in (b). Vertical grey lines correspond to a discontinuous change in predicted hydrodynamic stability, reflecting the inability of the Chebyshev expansion (
$N = 400$) to resolve unstable eigenmodes for values of slip parameter 
$S$ below the grey line. Wavenumbers are 
$(\alpha, \beta ) = (1, 0)$ for the planar flows, and 
$(\alpha, n) = (1, 1)$ for Hagen–Poiseuille pipe flow. See Appendices A and B for details of each calculation.
For values of slip parameter 
$S$ approaching zero, the linear stability maps are each consistent with the results of conventional linear stability analyses using (1.1) as the boundary condition for velocity perturbations. Specifically, both the plane Couette flow and the Hagen–Poiseuille pipe flow indicate linear stability of the flow (i.e. 
$\mathbb {I}[\hat {\omega }] < 0$) for all Reynolds numbers investigated. The plane Poiseuille flow is also in agreement with previous studies for 
$S$ approaching zero, with a single unstable eigenmode appearing between 
$Re = 5770$ and 
$Re = 5780$ (Orszag Reference Orszag1971).
However, a striking discontinuity is observed in each of the linear stability maps as the slip parameter 
$S$ is increased above a critical threshold (i.e. vertical grey lines). For values of the slip parameter greater than this threshold 
$S_{crit}$, we observe a region of linear instability spanning up to four decades in 
$S$. These strongly unstable eigenmodes – with growth rates 
$\mathbb {I}[\hat {\omega }]$ exceeding 
$O(10^7)$ in some cases – appear at all Reynolds numbers investigated, making them potentially relevant to the process of turbulence transition.
The apparent lack of unstable eigenmodes for values of slip parameter 
$S < S_{crit}$ is an artefact of the limited spatial resolution of the Chebyshev expansion using a finite number of modes. To illustrate this, figure 2(a) plots contours of the maximum Orr–Sommerfeld eigenvalue imaginary part 
$\mathbb {I}[\hat {\omega }]$ versus slip parameter 
$S$ and number of Chebyshev modes 
$N$, for the case of plane Couette flow at 
$Re = 360$ and wavenumbers 
$(\alpha, \beta ) = (1, 0)$. Eigenvalue contours to the right of the grey curve are oriented vertically, indicating that the corresponding eigenvalues are insensitive to the number of Chebyshev modes used to solve the Orr–Sommerfeld equation. This suggests that the computed eigenmodes are not spurious artefacts of the numerical method (cf. Dawkins, Dunbar & Douglass Reference Dawkins, Dunbar and Douglass1998). By contrast, the location of the grey stability discontinuity shifts to lower values of the slip parameter 
$S$ as the number of Chebyshev modes used to compute the eigenmodes increases. This boundary reflects the finite spatial resolution of the Chebyshev representation of unstable eigenmodes for finite mode numbers 
$N$. As the number of Chebyshev modes is increased, unstable eigenmodes with smaller slip length scale 
$S$ are successfully resolved. The minimum number of Chebyshev modes, 
$N_{min}$, required to resolve unstable eigenmodes for a given slip parameter 
$S$ is well approximated by 
$N_{min} \approx 2.3 S^{-1/2}$ (i.e. the green dotted line).
Figure 2. Contour maps of the maximum Orr–Sommerfeld eigenvalue imaginary part 
$\mathbb {I}[\hat {\omega }]$ versus slip parameter 
$S$ and (a) number of Chebyshev modes 
$N$ for plane Couette flow at 
$Re = 360$ and wavenumbers 
$(\alpha, \beta ) = (1, 0)$, (b) Reynolds number 
$Re$ computed using 
$N=50$ Chebyshev modes. Blue contours indicate regions of linear stability, and red contours indicate regions of linear instability. The neutral stability boundary is indicated by a black line. Discontinuity in stability is indicated by a grey curve. The location of the stability discontinuity shifts to lower values of the slip parameter 
$S$ as 
$N$ increases, reflecting the ability to resolve unstable eigenmodes with decreasing slip length scale 
$S$ for increasing numbers of modes 
$N$ used in the Chebyshev expansion. The minimum number of Chebyshev modes 
$N_{min}$ required to resolve unstable eigenmodes for a given slip parameter 
$S$ is well approximated by 
$N_{min} \approx 2.3 S^{-1/2}$ (green dotted line). The stability discontinuity (grey dotted line) from figure 1(a), computed using 
$N = 400$ Chebyshev modes, is reproduced to illustrate the rightward shift in the threshold value of the slip parameter 
$S_{crit}$ corresponding to the stability discontinuity when fewer Chebyshev modes are used.
Figure 2(b) provides a complementary illustration of this trend, showing the linear stability map computed using 
$N = 50$ Chebyshev modes. For comparison with the stability map computed using 
$N = 400$ Chebyshev modes in figure 1(a), the stability discontinuity (i.e. the grey dotted line) from that plot is reproduced here. The results highlight the need to implement Chebyshev collocation with a sufficient number of collocation points to resolve the distinct shape of unstable eigenmode profiles in close proximity to the wall. Calculations using an insufficient number of Chebyshev modes erroneously predict the absence of any unstable modes. In light of the trend 
$N_{min} \sim S^{-1/2}$, the number of Chebyshev modes required to resolve unstable eigenmodes becomes prohibitive as the conventional no-slip condition is approached (i.e. 
$S \rightarrow 0$). To be sure, a similar resolution challenge will be faced by other numerical methods.
2.3. Asymptotic analysis of Couette flow
To further analyse the discovered unstable eigenmodes while circumventing the aforementioned numerical resolution challenges, we conducted an asymptotic analysis of plane Couette flow in the limit of strongly unstable eigenmodes, i.e. 
$\mathbb {I}[\hat {\omega }] \gg 1$ (see Appendix D). This analysis predicts that the unstable eigenmode growth rate scales with slip parameter 
$S$ and Reynolds number 
$Re$ as 
$\mathbb {I}[\hat {\omega }] \sim S^{-2}\,Re^{-1}$. Moreover, the unstable plane Couette eigenmode amplitude components have the predicted asymptotic form
\begin{align} \tilde{u}(y) &\approx
({\rm i}/\alpha)\left[(1/S)\, {\rm e}^{(y-1)/S} + (1/S)\,
{\rm e}^{-(y+1)/S} \vphantom{\frac{\alpha\cosh
\alpha(y+1)}{\sinh 2\alpha}}\right.\nonumber\\ &\quad -
\frac{\alpha\cosh \alpha(y+1)}{\sinh 2\alpha} -
\left.\frac{\alpha\cosh \alpha(y-1)}{\sinh 2\alpha}\right]
+ O(S),
\end{align}
\begin{align}
\tilde{v}(y) &\approx {\rm
e}^{(y-1)/S} - {\rm e}^{-(y+1)/S} - \frac{\sinh
\alpha(y+1)}{\sinh 2\alpha}\nonumber\\
&\quad - \frac{\sinh
\alpha(y-1)}{\sinh 2\alpha} + O(S) .
\end{align}
Figure 3 compares the predictions of the asymptotic analysis with the calculations of plane Couette flow using Chebyshev collocation with 
$N = 400$. The agreement is excellent for both the maximum unstable eigenvalues (figure 3a) and the shape of each component of the unstable eigenmode (figure 3b). These results further support the conclusion that the discovered unstable eigenmodes are not spurious numerical artefacts.
Figure 3. (a) Maximum unstable eigenvalues of plane Couette flow plotted versus slip parameter 
$S$. Thick coloured curves are computed using Chebyshev collocation (
$N = 400$) for 
$Re = 50$ (green), 
$Re = 360$ (blue) and 
$Re = 3000$ (magenta). Thin grey lines are corresponding analytical model prediction 
$\mathbb {I}[\hat {\omega }] \sim S^{-2}\,Re^{-1}$. (b) Streamwise (
$\tilde {u}$, cyan) and wall-normal (
$\tilde {v}$, magenta) components of unstable eigenmode amplitude profiles computed using Chebyshev collocation (
$N = 400$, 
$Re = 360$, 
$S = 1 \times 10^{-2}$ and 
$(\alpha, \beta ) = (1, 0)$). Corresponding eigenmode profiles predicted by asymptotic analyses in (2.2) and (2.3) are shown in superimposed thin grey curves.
Nonetheless, the predicted scaling of the eigenmode growth rate 
$\mathbb {I}[\hat {\omega }] \sim S^{-2}\,Re^{-1}$ is counterintuitive, as it suggests stronger instability at lower Reynolds numbers. Moreover, the analysis does not exhibit a critical Reynolds number below which instability is not predicted to occur. This prediction is contradicted by the known stability of viscous parallel shear flows at sufficiently low Reynolds number. We speculate that the apparent inconsistency can be resolved by consideration of the predicted shape of the unstable eigenmodes. For example, from (2.2), we see that the magnitude of the streamwise component of the unstable eigenmode velocity gradient at the wall goes as 
$\|\partial \tilde {u}/\partial y\| \sim \mathbb {I}[\hat {\omega }]\,Re \gg 1$. In a fluid with finite dynamic viscosity 
$\mu$, the corresponding wall shear required to create this unstable perturbation is therefore 
$\tilde {\tau }_{wall} = \mu (\partial \tilde {u}/\partial y) \gg 1$, which may not be physically realizable. Hence although the analysis predicts that the flow is linearly unstable to eigenmodes of the shape given in (2.2) and (2.3), that prediction of eigenmode growth is relevant only if those eigenmodes are actually present in the flow. In practice, the viscosity of real fluids may prevent the appearance of these unstable eigenmodes at low Reynolds numbers due to the large wall shear required to produce them. To be sure, a nonlinear analysis would be necessary to resolve this inconsistency conclusively.
In addition, not all values of slip parameter 
$S$ will be physically realizable at all Reynolds numbers. Reynolds-number-dependent slip is commonly observed in flow over complex interfaces, e.g. porous media (Beavers & Joseph Reference Beavers and Joseph1967; Wu & Mirbod Reference Wu and Mirbod2018; Guo et al. Reference Guo, Li, Nian, Xu, Liu and Wang2020). If a similar phenomenon is associated with flow perturbations over solid walls (i.e. as opposed to the base flow), then the Reynolds number dependence of the present stability predictions may be implicit in the slip parameter 
$S$.
3. Perturbation model II
3.1. Further departure from the no-slip condition
Perturbation model I was formulated with the intent of parametrizing modest departures from the no-slip condition, such that the conventional no-slip condition (1.1) is a limiting case. Here, we consider the implications of more significant departures from the classical no-slip condition as applied to velocity perturbations in linear stability analyses of viscous parallel shear flows. In the analysis that follows, we maintain the assumption that the base flow satisfies the no-slip condition.
Following Schmid & Henningson (Reference Schmid and Henningson2001), we consider velocity perturbations of the general Cartesian form
\begin{equation} \boldsymbol{u}(x,y,z,t) = \mathbb{R}[\tilde{\boldsymbol{u}}(y)\, {\rm e}^{{\rm i}(\alpha x + \beta z - \omega t)}], \end{equation}
where the Cartesian components of the perturbation velocity vector are 
$\boldsymbol {u} = u \hat {\boldsymbol {i}} + v \hat {\boldsymbol {j}} + w \hat {\boldsymbol {k}}$ in the streamwise (i.e. 
$\hat {\boldsymbol {i}}$), wall-normal (i.e. 
$\hat {\boldsymbol {j}}$) and transverse (i.e. 
$\hat {\boldsymbol {k}}$) directions. The temporal evolution of the perturbations is examined by treating the spatial wavenumbers as 
$\alpha,\beta \in \mathbb {R}$ and the frequency as 
$\omega = \alpha c$, where 
$c \in \mathbb {C}$ is the complex phase speed of the perturbation. We will consider the case 
$\beta = 0$ below without loss of generality.
If we permit non-zero values of the streamwise perturbation velocity magnitude at the wall, then the corresponding fluid particle trajectories 
$X(t)$ along the wall are given by
\begin{equation} \frac{\partial X}{\partial t} = \|\tilde{\boldsymbol{u}}_{wall}\|\cos (\alpha x - \omega_R t)\, {\rm e}^{\omega_I t}, \end{equation}
where 
$\omega _R = \mathbb {R}[\omega ]$, 
$\omega _I = \mathbb {I}[\omega ]$, and the phase of 
$\tilde {\boldsymbol {u}}_{wall}$ is omitted without loss of generality. Inspection of (3.2) shows that fluid particles at the wall are stationary relative to the wall – and therefore satisfy the no-slip condition – where
\begin{equation} \alpha x - \omega_R t = \frac{\rm \pi}{2} + {\rm \pi}j, \end{equation}or equivalently,
\begin{equation} x = \mathbb{R}[c]\,t + \frac{\dfrac{{\rm \pi} }{2} + {\rm \pi}j}{\alpha}, \end{equation}
where 
$j$ is an integer. For single eigenmodes, the discrete locations where the no-slip condition is satisfied are not fixed, but instead propagate with the phase speed 
$\mathbb {R}[c]$ of the associated eigenmode. Consequently, all fluid particles on the wall will eventually experience non-zero slip. By contrast, for paired eigenmodes with real part 
$\pm \omega _R$, the velocity of fluid particles at the wall is given by
\begin{equation} \frac{\partial X}{\partial t} = \|\tilde{\boldsymbol{u}}_{wall}\|\,[\cos(\alpha x - \omega_R t) + \cos(\alpha x + \omega_R t)]\, {\rm e}^{\omega_I t}. \end{equation}The term in brackets in (3.5) can be rewritten as
\begin{equation} \cos(\alpha x - \omega_R t) + \cos(\alpha x + \omega_R t) = 2\cos(\alpha x)\cos(\omega_R t). \end{equation} Hence the no-slip condition is satisfied at all times at the discrete, fixed locations 
$x = ({{\rm \pi} }/{2} + {\rm \pi}j)/\alpha$. Once fluid particles are advected to these locations along the wall, they asymptotically satisfy the no-slip condition thereafter. During any initial transient wall slip to reach these locations, the maximum fluid particle displacement along the wall is 
$O({\rm \pi} /\alpha )$, i.e. the distance between adjacent locations where 
${\partial X}/{\partial t} = 0$. The duration of the initial wall slip transient is 
$O({\rm \pi} /\|\tilde {\boldsymbol {u}}_{wall}\|\,\alpha )$.
For stable eigenmodes, i.e. 
$\omega _I < 0$, the wall slip decays exponentially per (3.5), rendering any perturbation-induced motion of fluid particles along the wall likely imperceptible in practice. For unstable eigenmodes, the exponentially increasing wall slip speed associated with 
$\omega _I > 0$ will advect particles to the stationary nodes on the wall more rapidly, thereby hastening their asymptotic approach to the no-slip condition in those cases. To be sure, there will always remain fluid particles in motion along some portion of the wall in order to satisfy continuity. However, the linearity of the eigensolutions implies that this motion can theoretically be arbitrarily small at the incipient stage of an instability. More specifically, if 
$\tilde {\boldsymbol {u}}$ is an eigenmode of the Orr–Sommerfeld equations, then any constant multiple 
$\gamma$ of the eigenmode (i.e. 
$\gamma \tilde {\boldsymbol {u}}$) is also an eigenmode, including those for which 
$\gamma \ll 1$. Together, these considerations lead to the interesting possibility that the wall slip associated with both stable and unstable velocity perturbations will be empirically indistinguishable from the conventional no-slip condition, even if the perturbations do not formally satisfy the no-slip boundary condition.
If the no-slip condition is not the only plausible model for the physics of velocity perturbations in viscous parallel shear flows, then exploration of other physically realistic boundary conditions for linear stability analysis is warranted. Perturbation model I was one such alternative. In the following subsections, we explore another boundary condition that not only produces linearly unstable eigenmodes, but also leads to quantitative predictions of an explicit critical Reynolds number range associated with linear instability.
3.2. Womersley model of velocity perturbations
The Orr–Sommerfeld equation is a fourth-order system, therefore we must specify four boundary conditions on the two walls of the plane Couette or Poiseuille flow. The first two boundary conditions prohibit flow through the walls:
$$\begin{gather} \tilde{v}(y = 1,t) =0, \end{gather}$$
$$\begin{gather}\tilde{v}(y ={-}1,t) = 0. \end{gather}$$In lieu of the conventional no-slip condition, we instead hypothesize that the fundamentally oscillatory nature of the streamwise velocity perturbations (cf. (3.1)) can be modelled with a boundary condition based on Stokes’ second problem (Landau & Lifschitz Reference Landau and Lifschitz1987) as
\begin{equation} \frac{\partial \tilde{u}}{\partial t} - \frac{1}{Re}\,\frac{\partial^2 \tilde{u}}{\partial y^2} = 0. \end{equation} Equation (3.9) provides the remaining two required boundary conditions at 
$y = \pm 1$. Note that (3.9) and (2.1) are only two of many physically motivated ansatzes that could be explored in general. Indeed, a goal of this paper is to demonstrate that there exist boundary conditions in addition to the no-slip condition that may also accurately represent incipient perturbations, and that these alternative boundary conditions may be associated with linear stability that is more consistent with empirical observations of turbulence transition.
Substituting the form of velocity perturbation in (3.1) into (3.9), and using 
$\tilde {u} = ({\rm i}/\alpha )\,\mathcal {D}\tilde {v}$ from the continuity equation (see Appendix A), the boundary condition can be expressed as
\begin{equation} \omega\,\frac{\partial \tilde{v}}{\partial y} - \frac{{\rm i}}{Re}\,\frac{\partial^3 \tilde{v}}{\partial y^3} = 0. \end{equation} We introduce the perturbation Womersley number 
$\widetilde {Wo}$ as a measure of the transient inertial force associated with the velocity perturbation 
$\boldsymbol {u}$ relative to its associated viscous force (cf. Womersley Reference Womersley1955):
\begin{equation} \widetilde{Wo} = \frac{1}{\alpha}\,\sqrt{\frac{\omega}{\nu_0}}, \end{equation}
where 
$\nu _0$ is the kinematic viscosity of the fluid. The frequency 
$\omega$ is complex in general; however, we make the simplifying assumption 
$\|\mathbb {I}[\omega ]\| \ll 1$, from which it follows that 
$\widetilde {Wo} \in \mathbb {R}$. This approximation holds identically for neutral stability, i.e. 
$\mathbb {I}[\hat {\omega }] = 0$. We will see in the next subsection that this approximation also captures large regions of parameter space away from the neutral stability boundaries.
Since the planar flows under present consideration are defined with 
$U_0 = 1$ and 
$L_0 = 1$, the boundary condition (3.10) can be written as
\begin{equation} \frac{\widetilde{Wo}^2\,\alpha^2}{Re}\,\frac{\partial \tilde{v}}{\partial y} - \frac{{\rm i}}{Re}\,\frac{\partial^3 \tilde{v}}{\partial y^3} = 0. \end{equation} Linear stability analysis proceeds by solving the Orr–Sommerfeld equation using the boundary conditions (3.7), (3.8) and (3.12) at 
$y = \pm 1$. Here, we use Chebyshev collocation to solve the Orr–Sommerfeld eigenvalue equation following Schmid & Henningson (Reference Schmid and Henningson2001) for plane Couette and Poiseuille flow, and following Malik & Skote (Reference Malik and Skote2019) for Hagen–Poiseuille pipe flow (see Appendices A and B).
3.3. Linear stability maps
Figure 4 plots contours of the maximum eigenvalue imaginary part, i.e. 
$\mathbb {I}[\hat {\omega }]$, for Orr–Sommerfeld solutions corresponding to plane Couette flow with 
$\widetilde {Wo} = 10$ and 
$(\alpha, \beta ) = (1, 0)$. For relatively low numbers of Chebyshev modes 
$N$ used to solve the Orr–Sommerfeld equation, we observe spurious unstable eigenmodes at increasing Reynolds numbers (i.e. the upper left corner of the contour plot). The spurious nature of these eigenvalues is evidenced by the dependence of their value on the number of Chebyshev modes used to solve the Orr–Sommerfeld equation (Dawkins et al. Reference Dawkins, Dunbar and Douglass1998). However, as the number of Chebyshev modes is increased, spurious eigenmodes disappear and the eigenvalues converge. This convergence is indicated by the horizontal alignment of the eigenvalue contours for large 
$N$. Note that the present achievement of a resolution-independent stability map is distinct from the behaviour observed in the case of perturbation model I, wherein the required minimum number of Chebyshev modes increased monotonically as 
$N_{min} \approx 2.3 S^{-1/2}$ (e.g. figure 2). For subsequent analyses, we use 
$N = 160$ to ensure that the eigenvalues are converged and not numerically spurious.
Figure 4. Contour map of the maximum Orr–Sommerfeld eigenvalue imaginary part 
$\mathbb {I}[\hat {\omega }]$ versus the number of Chebyshev modes 
$N$ and and Reynolds number 
$Re$ for plane Couette flow with 
$\widetilde {Wo} = 10$ and 
$(\alpha, \beta ) = (1, 0)$.
Intriguingly, a robust band of linear instability is observed in the range 
$100 < Re < 1000$. To explore this in greater detail, figure 5 plots a stability map of plane Couette flow versus the perturbation Womersley number 
$\widetilde {Wo}$ and Reynolds number 
$Re$ for 
$(\alpha, \beta ) = (1, 0)$. Several interesting features emerge. At Reynolds numbers below 
$Re \approx 30$, the flow is stable for all values of perturbation Womersley number 
$\widetilde {Wo}$. Similarly, if the perturbation Womersley number is below 
$\widetilde {Wo} \approx 5$, then the flow is stable for all values of Reynolds number. Both of these effects may be attributable to the dampening effect of the fluid viscosity on the base flow and on the incipient perturbations, respectively. The latter dependence of the flow stability on the nature of the incipient perturbations, as measured by perturbation Womersley number 
$\widetilde {Wo}$, may explain empirical observations of a range of different Reynolds numbers corresponding to turbulence transition in plane Couette flow experiments (Drazin & Reid Reference Drazin and Reid2004). Specifically, the present analysis predicts that plane Couette flow becomes linearly unstable for 
$Re > \widetilde {Wo}^2$, as indicated by the green dashed boundary in figure 5.
Figure 5. Contour map of the maximum Orr–Sommerfeld eigenvalue imaginary part 
$\mathbb {I}[\hat {\omega }]$ versus perturbation Womersley number 
$\widetilde {Wo}$ and Reynolds number 
$Re$ for plane Couette flow with 
$(\alpha, \beta ) = (1, 0)$. The region of linear instability (red contours) is well described by 
$\widetilde {Wo}^2 < Re < \widetilde {Wo}^3$, as indicated by the green and yellow dashed lines, respectively.
Interestingly, this analysis also predicts that linear instability can be avoided at sufficiently high Reynolds numbers if the flow satisfies 
$Re > \widetilde {Wo}^3$, i.e. the region above the yellow dashed boundary in figure 5. This condition may be difficult to achieve in practice for three reasons. First, an arbitrary incipient perturbation may comprise a spectrum of Womersley numbers 
$\widetilde {Wo}$ that spans a broad range of the abscissa in figure 5. Satisfying 
$Re > \widetilde {Wo}^3$ across that full spectrum may be difficult in the absence of a means to control the frequency spectrum of incipient perturbations a priori. Second, a flow generated from rest increases in Reynolds number from zero, therefore is likely to cross the region of linear instability from below. Third, the condition 
$Re > \widetilde {Wo}^3$ can be recast in terms of the governing physical parameters as
\begin{equation} \frac{U_0L_0}{\nu_0} > \frac{\omega^{3/2}}{\nu_0^{3/2}\alpha^3}. \end{equation}Hence any attempt to increase the Reynolds number by reducing the viscosity will lead to a concomitant and larger increase in the perturbation Womersley number, making it more difficult to satisfy the stability criterion in (3.13).
These challenges notwithstanding, the stability map in figure 5 is a testable prediction of this model that can be examined in experiments. Confirmation of the model predictions would suggest that this alternative to the no-slip condition may improve the fidelity of linear stability analysis as a predictive tool. Furthermore, if the regime 
$Re > \widetilde {Wo}^3$ is empirically accessible, then it may be possible to achieve stable, laminar flow and associated fluid dynamic performance gains at arbitrarily high Reynolds numbers.
Similar results were observed in analysis of plane Poiseuille flow as well as for Hagen–Poiseuille pipe flow (see Appendix A for pipe flow governing equations). As shown in figure 6, the region of linear instability has a similar shape for both Poiseuille flows as in the plane Couette flow. The neutral stability boundaries have a similar scaling but distinct prefactors, i.e. 
$0.5\,\widetilde {Wo}^2 < Re < 0.5\,\widetilde {Wo}^3$ and 
$4\,\widetilde {Wo}^2 < Re < 2.5\,\widetilde {Wo}^3$, respectively. The numerical discrepancy in prefactors of the two planar flows can be resolved by defining the Couette flow Reynolds number based on the maximum flow speed contrast in the domain (i.e. 
$U_0 = 2$) instead of the maximum flow speed, or by using the full-width distance between the walls as the characteristic length instead of the half-width.
Figure 6. Contour map of the maximum Orr–Sommerfeld eigenvalue imaginary part 
$\mathbb {I}[\hat {\omega }]$ versus perturbation Womersley number 
$\widetilde {Wo}$ and Reynolds number 
$Re$ for (a) plane Poiseuille flow with 
$(\alpha, \beta ) = (1, 0)$, and (b) Hagen–Poiseuille pipe flow with 
$(\alpha, n) = (1, 1)$. The region of linear instability (red contours) is well described by 
$\widetilde {Wo}^2 < Re < \widetilde {Wo}^3$, as indicated by the green and yellow dashed lines, respectively. See the main text for discussion of scaling prefactors.
A potential conceptual inelegance of the model described herein is its characterization of incipient perturbations using an independent parameter 
$\widetilde {Wo}$. The frequency 
$\omega$ used to define the perturbation Womersley number is not necessarily identical to the eigenvalue 
$\hat {\omega }$ of the Orr–Sommerfeld system. Indeed, the frequency 
$\omega$ is an input to the boundary condition (3.12), whereas the eigenvalue 
$\hat {\omega }$ is a solution output. It is unclear whether this distinction is significant, but it is noted here for the sake of completeness.
4. Conclusions
The present analyses are motivated by the premise that incipient velocity perturbations in viscous parallel shear flows may not satisfy the no-slip condition exactly. Two alternative models of the fluid–solid interface physics demonstrated that even small departures from the exact no-slip condition – at levels that may be empirically indistinguishable from the the classical no-slip condition – can significantly affect the predictions of linear stability analysis. Intriguingly, perturbation model II predicts a minimum critical Reynolds number for linear instability that is quantitatively consistent with empirical observations of turbulence transition, e.g. transition Reynolds numbers of 
$O(10^3\unicode{x2013}10^4)$ in pipe flow (Avila et al. Reference Avila, Barkley and Hof2023). Moreover, the model predicts an unexpected, second regime of high-Reynolds-number flow stability. Experimental validation of that prediction would not only provide compelling support for perturbation model II as a realistic model of the perturbation flow physics, it would also potentially have important implications in engineering applications.
More generally, the discovered linearly unstable eigenmodes in both perturbation models imply the presence of distinct, near-wall flow kinematics that may also be observable in experiments that confirm their existence. The transient nature of these near-wall fluid particle motions and their appearance as perturbations to the base flow may require development of more sensitive measurement techniques than are available currently. Nonetheless, the present theoretical framework can guide experimental efforts to test these concepts further.
Acknowledgements
The authors thank M. Malik for discussion of the pipe flow calculations and M.K. Fu for manuscript feedback.
Funding
J.O.D. was supported by the US National Science Foundation.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
Matlab implementations of the linear stability codes are openly available in the Caltech Data Repository at https://doi.org/10.22002/4sm2w-br988.
Appendix A. Orr–Sommerfeld equations
A.1. Plane Couette and plane Poiseuille flow
Following Schmid & Henningson (Reference Schmid and Henningson2001), we consider velocity perturbations of the general Cartesian form
\begin{equation} \boldsymbol{u}(x,y,z,t) = \mathbb{R}[\tilde{\boldsymbol{u}}(y)\, {\rm e}^{{\rm i}(\alpha x + \beta z - \omega t)}], \end{equation}
where the Cartesian components of the velocity vector are 
$\boldsymbol {u} = u \hat {\boldsymbol {i}} + v \hat {\boldsymbol {j}} + w \hat {\boldsymbol {k}}$ in the streamwise (i.e. 
$\hat {\boldsymbol {i}}$), wall-normal (i.e. 
$\hat {\boldsymbol {j}}$) and transverse (i.e. 
$\hat {\boldsymbol {k}}$) directions. The temporal evolution of the perturbations is examined by treating the spatial wavenumbers as 
$\alpha,\beta \in \mathbb {R}$ and the frequency as 
$\omega = \alpha c$, where 
$c \in \mathbb {C}$ is the complex phase speed of the perturbation. The stability of perturbations of a given spatial wavenumber is evaluated by solving the Orr–Sommerfeld equation
\begin{equation} \left[(-{\rm i}\omega + {\rm i} \alpha U)(\mathcal{D}^2 - k^2)-{\rm i}\alpha U''-\frac{1}{Re}\,(\mathcal{D}^2 - k^2)^2\right]\tilde{v} = 0, \end{equation}
where 
$k^2 = \alpha ^2 + \beta ^2$, and 
$Re = U_0L_0/\nu _0$ is the Reynolds number. Here, 
$U_0$ is taken as the maximum velocity of the base flow in the domain, and 
$L_0$ is the half-width of the distance between the solid walls of the Couette or Poiseuille flow. The operator 
$\mathcal {D}$ and the prime both denote a derivative with respect to the wall-normal 
$y$-coordinate direction.
For plane Couette flow, the base flow is given by 
$U(y) = y$ on the domain 
$-1 \leq y \leq 1$. For plane Poiseuille flow, the base flow is given by 
$U(y) = 1 - y^2$ on the domain 
$-1 \leq y \leq 1$. In perturbation model I, (A2) is solved using the boundary condition 
$\tilde {v} = 0$ at 
$y = \pm 1$ (i.e. no flow penetration at the solid walls), along with the generalized boundary condition for the streamwise perturbation component 
$\tilde {u}$ in (2.1) to replace the conventional no-slip condition. Given the form of velocity perturbation in (A1), and using 
$\tilde {u} = ({\rm i}/\alpha )\,\mathcal {D}\tilde {v}$ from the continuity equation, the generalized boundary condition for 
$\tilde {u}$ can be written in terms of the normal velocity 
$\tilde {v}$ as
\begin{equation} \frac{{\rm i}}{\alpha}\,\mathcal{D}\tilde{v} \mp \frac{{\rm i}}{\alpha}\,S\,\mathcal{D}^2\tilde{v} = 0 \quad \textrm{for}\ y = \pm1. \end{equation} In perturbation model II, (A2) is solved using boundary conditions (3.7), (3.8), and (3.12) at 
$y = \pm 1$.
A.2. Hagen–Poiseuille pipe flow
Following Malik & Skote (Reference Malik and Skote2019), we consider velocity perturbations in cylindrical coordinates of the form
\begin{equation} \boldsymbol{u}(x,r,\theta,t) = \mathbb{R}[\tilde{\boldsymbol{u}}(r)\, {\rm e}^{{\rm i}(\alpha x + n \theta - \omega t)}], \end{equation}
where 
$(x,r,\theta )$ are the axial, radial and azimuthal directions, respectively, and 
$n$ is the azimuthal wavenumber. The Orr–Sommerfeld equations can be expressed as
\begin{align} \left.\begin{array}{c@{}} {\rm i}(\alpha U-\omega)(\mathcal{D}\tilde{u}-{\rm i} \alpha\tilde{v}) ={-}{\rm i}\alpha U'\tilde{u}-(U''+U'\,\mathcal{D})\tilde{v}+\dfrac{1}{Re}\, (\varDelta(\mathcal{D}\tilde{u}-{\rm i} \alpha\tilde{v})-2 {\rm i} n r^{{-}2}\tilde{\eta}), \\ {\rm i}(\alpha U-\omega)[\mathcal{D}(r\tilde{w})-{\rm i}n\tilde{v}] ={-}{\rm i} \alpha r U'\tilde{w}+\dfrac{1}{Re}\, (\bar{\varDelta}[\mathcal{D}(r\tilde{w})-{\rm i}n\tilde{v}]),\\ {\rm i}(\alpha U-\omega)\tilde{\eta} ={-}{\rm i}nU'r^{{-}1}\tilde{v}+\dfrac{1}{Re}\,(\varDelta\eta + 2nr^{{-}2}(\alpha\tilde{v}+{\rm i}\,\mathcal{D}\tilde{u})), \end{array}\right\} \end{align}
where 
$\varDelta = [\mathcal {D}^2+r^{-1}\,\mathcal {D}-r^{-2}(d+1)]$, 
$\bar {\varDelta } = [\mathcal {D}^2-r^{-1}\,\mathcal {D}-r^{-2}(d-1)]$, 
$\tilde {\eta } = \textrm {i}(nr^{-1}\tilde {u} -\alpha \tilde {w})$, 
$d = n^2 + \alpha ^2 r^2$, and the operator 
$\mathcal {D}$ and the prime both denote a derivative with respect to the radial direction. The base flow for Hagen–Poiseuille pipe flow is given by 
$U(r) = 1-r^2$ on the domain 
$0 \leq r \leq 1$.
The Orr–Sommerfeld equations are solved following Malik & Skote (Reference Malik and Skote2019) by using surrogate analytic functions 
$(\phi, \varOmega )$ that facilitate more straightforward implementation of the perturbation models (see below) and regularity conditions on the axis of symmetry 
$r = 0$. The wall-normal velocity and vorticity are expressed in terms of these functions as
\begin{equation} (\tilde{v}, \tilde{\eta}) = \begin{cases} (r^l\phi, r^l\varOmega), & n \neq 0,\\ (r\phi, r\varOmega), & n = 0, \end{cases} \end{equation}
where 
$l = |n| - 1$. The no-penetration boundary condition is therefore given by 
$\phi = 0$ at 
$r = 1$.
Using continuity and the normal vorticity definition, the streamwise (i.e. 
$\tilde {u}$) velocity perturbation component is given by
\begin{equation} \tilde{u} = \begin{cases} \dfrac{{\rm i}r^{l+1}}{d}\,[\alpha(l+1)\phi + \alpha r\,\mathcal{D}\phi - n \varOmega], & n \neq 0,\\ \dfrac{{\rm i}}{\alpha}\,(2\phi + r\mathcal{D}\phi), & n = 0. \end{cases} \end{equation}A.2.1. Perturbation model I
From (2.1) and (A7), the generalized boundary condition for the streamwise velocity component at the wall (
$r = 1$) is
\begin{equation} \left. \begin{array}{c@{}} [(a_1 + a_3 + a_6)\mathcal{D} + a_4\,\mathcal{D}^2]\phi + [(a_2+a_7) + a_5\,\mathcal{D}]\varOmega = 0, \quad n \neq 0,\\[3pt] {[}(a_1^0 + a_2^0)\,\mathcal{D} + a_3^0\,\mathcal{D}^2 ] \phi = 0, \quad n = 0, \end{array}\right\}, \end{equation}where
\begin{equation} \left.\begin{array}{c@{}}
a_1 = \dfrac{{\rm i} \alpha}{d}, \\ a_2
={-}\dfrac{{\rm i}n}{d}, \\ a_3 ={-}\dfrac{S\,{\rm
i} \alpha(l+2)}{d}, \\ a_4 ={-}\dfrac{S\,{\rm i}
\alpha}{d}, \\ a_5 = \dfrac{S\,{\rm
i}n}{d}, \\ a_6 ={-}S \alpha\left[\dfrac{{\rm
i}(l+1)d-2\,{\rm i} \alpha^2}{d^2}\right], \\ a_7
= S n\left[\dfrac{{\rm i}(l+1)d-2\,{\rm i}
\alpha^2}{d^2}\right],\\ a_1^0 = \dfrac{{\rm
i}}{\alpha}, \\ a_2^0 ={-}3S\,\dfrac{{\rm
i}}{\alpha},\\ a_3^0 ={-}S\,\dfrac{{\rm
i}}{\alpha}. \end{array}\right\} \end{equation}
A.2.2. Perturbation model II
From (3.9) and (A7), the boundary condition for the streamwise velocity component at the wall (
$r = 1$) is
\begin{align} \left.\begin{array}{c@{}}
\begin{aligned} &
\left[a_1\,\mathcal{D}-\dfrac{1}{Re}\,(b_1\alpha+b_2\alpha(l+2)
+b_3\alpha(l+2))\,\mathcal{D}-\dfrac{1}{Re}\,
(b_2\alpha+b_3\alpha+b_4\alpha(l+3))\,\mathcal{D}^2
\right.\\ & \quad
\left.{}-\dfrac{1}{Re}\,(b_4\alpha)\,\mathcal{D}^3\right]\phi+
\left[a_2+\dfrac{1}{Re}\,(b_1n)+\dfrac{1}{Re}(b_2+b_3)
n\,\mathcal{D}+\dfrac{1}{Re}\,(b_4n)\,\mathcal{D}^2\right]\varOmega
= 0, \quad n \neq 0, \end{aligned}\\ {[}(a_1^0 +
a_2^0)\,\mathcal{D} + a_3^0\,\mathcal{D}^2 +
a_4^0\,\mathcal{D}^3]\phi = 0, \quad n = 0,
\end{array}\right\} \end{align}
where
\begin{align} \left.\begin{array}{c@{}} a_1
= \dfrac{(\alpha\,
\widetilde{Wo})^2}{Re}\left(\dfrac{\alpha}{n^2+\alpha^2}\right),\\
a_2 = \dfrac{-(\alpha\,
\widetilde{Wo})^2}{Re}\left(\dfrac{n}{n^2+\alpha^2}\right),\\
b_1 =\dfrac{{\rm i}(l+1)(2\alpha^2)+{\rm
i}(l+1)^2(n^2+\alpha^2)- 2{\rm
i}(l+3)\alpha^2}{(n^2+\alpha^2)^2} - \dfrac{4\alpha^2({\rm
i}(l+1)(n^2+\alpha^2)-2\,{\rm i}
\alpha^2)}{(n^2+\alpha^2)^3},\\ b_2 =\dfrac{{\rm
i}(l+1)(n^2+\alpha^2)-2{\rm i}
\alpha^2}{(n^2+\alpha^2)^2},\\ b_3 =\dfrac{{\rm
i}(l+2)(n^2+\alpha^2)-2{\rm i}
\alpha^2}{(n^2+\alpha^2)^2},\\ b_4 =\dfrac{{\rm
i}}{n^2+\alpha^2},\\ a_1^0
=\dfrac{(\alpha\,\widetilde{Wo})^2}{Re},\\ a_2^0 =
\dfrac{-3{\rm i}}{Re},\\ a_3^0 = \dfrac{-5{\rm
i}}{Re},\\ a_4^0= \dfrac{-{\rm i}}{Re}.
\end{array} \right\} \end{align}
Appendix B. Numerical solution
The Orr–Sommerfeld eigenvalue equations were solved using Chebyshev collocation as formulated in plane Cartesian coordinates (Schmid & Henningson Reference Schmid and Henningson2001) and cylindrical coordinates (Malik & Skote Reference Malik and Skote2019). Matlab implementations of the algorithm for each flow are provided in the Caltech Data Repository.
The Matlab codes were verified by comparing the computed eigenvalues for the case 
$S = 0$ (i.e. no-slip condition) to the tabulated eigenvalues for plane Couette and Poiseuille flows in Schmid & Henningson (Reference Schmid and Henningson2001), as well as the tabulated eigenvalues for Hagen–Poiseuille pipe flow in Malik & Skote (Reference Malik and Skote2019).
B.1. Perturbation model I
For each of the three flows (i.e. plane Couette flow, plane Poiseuille flow and Hagen–Poiseuille pipe flow), 400 collocation points were used to compute each of the 
$996\,Re \times 66\,S = 65\,736$ total cases in the 
$S$–
$Re$ parameter space of figure 1. Wavenumbers were fixed at 
$\alpha = 1$ and 
$\beta = 0$ for plane Couette flow and plane Poiseuille flow, and 
$\alpha = 1$ and 
$n = 1$ for Hagen–Poiseuille pipe flow. These wavenumbers typically corresponded to the most unstable eigenmodes, but qualitatively similar results were observed for other wavenumbers. The eigenvalue with maximum imaginary part was recorded for each parameter set 
$(S, Re)$ studied. If the eigenvalue with second-largest imaginary part was within 
$10^{-3}$ of the eigenvalue with largest imaginary part, then the eigenvalues were treated as a combined pair for the purposes of subsequent analysis.
B.2. Perturbation model II
For each of the three flows (i.e. plane Couette flow, plane Poiseuille flow and Hagen–Poiseuille pipe flow), 160 collocation points were used to compute each of the 
$64\,Re \times 100\,\widetilde {Wo} = 6400$ total cases in the 
$\widetilde {Wo}$–
$Re$ parameter space of figures 5 and 6. Wavenumbers were fixed at 
$\alpha = 1$ and 
$\beta = 0$ for plane Couette flow and plane Poiseuille flow, and 
$\alpha = 1$ and 
$n = 1$ for Hagen–Poiseuille pipe flow. These wavenumbers typically corresponded to the most unstable eigenmodes, but qualitatively similar results were observed for other wavenumbers. The eigenvalue with maximum imaginary part was recorded for each parameter set 
$(\widetilde {Wo}, Re)$ studied. If the eigenvalue with second-largest imaginary part was within 
$10^{-3}$ of the eigenvalue with largest imaginary part, then the eigenvalues were treated as a combined pair for the purposes of subsequent analysis.
Appendix C. Second-largest Orr–Sommerfeld eigenvalue imaginary part for perturbation model I
Figure 7 plots contours of the second-largest eigenvalue imaginary part for perturbation model I. These results demonstrate that the predicted linear instability is not necessarily limited to a single eigenvalue with positive imaginary part.
Figure 7. Contour maps of the second-largest Orr–Sommerfeld eigenvalue imaginary part 
$\mathbb {I}[\hat {\omega }]$ versus slip parameter 
$S$ and Reynolds number 
$Re$ for (a) plane Couette flow, (b) plane Poiseuille flow, and (c) Hagen–Poiseuille pipe flow. Blue contours indicate regions of linear stability, and red contours indicate regions of linear instability. Black contours indicate neutral stability boundaries. Vertical grey lines correspond to a discontinuous change in predicted hydrodynamic stability, reflecting the inability of the Chebyshev expansion (
$N = 400$) to resolve unstable eigenmodes for values of slip parameter 
$S$ below the grey line.
Appendix D. Asymptotic analysis of unstable Couette eigenmodes for perturbation model I
The Orr–Sommerfeld equation for plane Couette flow can be written as
\begin{equation} \frac{\partial}{\partial t}\nabla^{2}\tilde{v} - \frac{1}{Re}\,\nabla^{2}\nabla^{2}\tilde{v}+{\rm i} \alpha\,\nabla^{2}\tilde{v} = 0, \end{equation}
where 
$\alpha$ is the axial wavenumber. The temporal behaviour of linear velocity perturbations is given by 
$\tilde {v} \sim \textrm {e}^{-\textrm {i}\omega t}$, where the complex frequency is 
$\omega = \mathbb {R}[\omega ] + \textrm {i}\,\mathbb {I}[\omega ]$.
We consider the limiting case of large, unstable eigenvalues, i.e. 
$\mathbb {I}[\hat {\omega }] \gg 1$. To leading order, (D1) can therefore be written as
\begin{equation} \left(\mathbb{I}[\hat{\omega}]\,Re + {\alpha}^2 - \frac{\partial^2}{\partial y^2}\right)\nabla^{2}\tilde{v} = 0. \end{equation}
Equation (D2) has leading-order solutions of the form 
$\tilde {v} \sim \textrm {e}^{\pm \kappa _1 y},\, \textrm {e}^{\pm \alpha y}$, where 
$\kappa _1 = \sqrt {\mathbb {I}[\hat {\omega }]\,Re + {\alpha }^2}$ 
$\gg 1$. The solutions must satisfy a no-penetration condition at the walls, i.e. 
$\tilde {v}(y = \pm 1) = 0$. This boundary condition is satisfied for
\begin{equation} \tilde{v}(y) = {\rm e}^{\kappa_1 (y-1)} - \frac{\sinh \alpha(y+1)}{\sinh 2\alpha} + \kappa_2 \left({\rm e}^{-\kappa_1(y+1)} + \frac{\sinh \alpha(y-1)}{\sinh 2\alpha}\right), \end{equation}
since 
$\tilde {v}(y = \pm 1) \sim O(\textrm {e}^{-2\kappa _1}) \approx 0$. The free parameter 
$\kappa _2$ can be constrained by applying the boundary condition in (2.1) as formulated for plane Couette flow in (A3). Specifically, at 
$y = 1$,
\begin{equation} \frac{\kappa_1 - \alpha\,\dfrac{\cosh 2\alpha}{\sinh 2\alpha} + \dfrac{\alpha \kappa_2}{\sinh 2\alpha}}{\kappa_1^2-\alpha^2 } =S, \end{equation}
where terms of 
$O(\textrm {e}^{-2\kappa _1})$ or smaller have been neglected. Similarly, at 
$y = -1$,
\begin{equation} \frac{-\dfrac{\alpha}{\sinh \alpha} - \kappa_1 \kappa_2 + \alpha \kappa_2\,\dfrac{\cosh 2\alpha}{\sinh 2\alpha}}{\kappa_2 (\kappa_1^2-\alpha^2)} ={-}S. \end{equation} Combining (D4) and (D5) gives 
$\kappa _2 = \pm 1$. The negative solution 
$\kappa _2 = -1$ gives antisymmetric eigenmodes consistent with the opposite direction of travel of the walls at 
$y = \pm 1$. Substituting this solution for 
$\kappa _2$ into the system of equations gives the approximation 
$\kappa _1 \approx S^{-1}$, or equivalently,
\begin{equation} \mathbb{I}[\hat{\omega}] \sim S^{{-}2}\,Re^{{-}1}. \end{equation}This predicted scaling is consistent with the unstable eigenvalues computed using Chebyshev collocation with 
$N = 400$ in figure 3(a) for 
$Re = 50$, 360 and 3000. The preceding asymptotic analysis also predicts that eigenmode profiles corresponding to large, unstable eigenvalues have the approximate shape given in (2.2) and (2.3). This prediction is in excellent agreement with the profiles computed using Chebyshev collocation, as shown in figure 3(b).
Appendix E. Alternative boundary conditions for perturbation model I
To explore the dependence of the discovered unstable eigenmodes in perturbation model I on the form of the boundary condition ansatz in (2.1), figure 8 plots the streamwise component of the eigenmode with maximum imaginary part for 
$Re = 360$ and 
$S = 1 \times 10^{-1}$, using the following alternative boundary conditions:
\begin{equation}
\left.\begin{array}{cc@{}} {(\textit{a})} \quad
\tilde{u}(\boldsymbol{x}_{wall},t) - S
\tilde{u}'(\boldsymbol{x}_{wall},t) = 0, & y
={\pm} 1,\\[3pt] {(\textit{b})} \quad
\tilde{u}(\boldsymbol{x}_{wall},t) + S
\tilde{u}'(\boldsymbol{x}_{wall},t) = 0, & y
={\pm} 1,\\[3pt] {(\textit{c})} \quad
\tilde{u}(\boldsymbol{x}_{wall},t) \mp S
\tilde{u}'(\boldsymbol{x}_{wall},t) = 0, & y
={\pm} 1,\\[3pt] {(\textit{d})} \quad
\tilde{u}(\boldsymbol{x}_{wall},t) \pm S
\tilde{u}'(\boldsymbol{x}_{wall},t) = 0, & y
={\pm} 1.\\ \end{array}\right\}
\end{equation}
Figure 8. Streamwise (i.e. 
$\mathbb {R}[\tilde {u}]$) eigenmode components for eigenvalues with maximum imaginary part in plane Couette flow (
$Re = 360$, 
$S = 1 \times 10^{-1}$) with generalized boundary condition in (2.1) (i.e. (c)) and alternative boundary conditions in (E1) (i.e. (a,b,d)). Directional alignment of perturbation wall slip and associated wall shear is enforced on (a) top wall, (b) bottom wall, (c) both walls, (d) neither wall. Eigenmodes in (a–c) are unstable; eigenmode in (d) is stable.
The boundary condition in figure 8(a) constrains the streamwise component 
$\tilde {u}$ of any velocity perturbation at the top wall (
$y = 1$) to be oriented in the same direction as the associated shear exerted by the wall on the fluid due to the perturbation. Conversely, on the bottom wall (
$y = -1$), the boundary condition sets the direction of any velocity perturbation to be opposite to the direction of the associated shear exerted by the wall on the fluid. The boundary condition in figure 8(b) is the inverse of that in figure 8(a), constraining the streamwise component 
$\tilde {u}$ of any velocity perturbation at the bottom wall (
$y = -1$) to be oriented in the same direction as the associated shear exerted by the wall on the fluid due to the perturbation. The boundary condition on the top wall in figure 8(b) is identical to the boundary condition on the bottom wall in figure 8(a). Figure 8(c) illustrates a boundary condition that is symmetric with respect to its treatment of the top and bottom walls, as in (2.1). In this case, the streamwise component of velocity perturbations at both walls is constrained to be oriented in the same direction as the associated shear exerted by the wall on the fluid due to the perturbation. In figure 8(d), the boundary condition is also symmetric with respect to its treatment of the top and bottom walls. However, the streamwise component of velocity perturbations at both walls is constrained to be oriented in the direction opposite to the associated wall shear. This boundary condition is similar to that used in prior studies of the effect of wall slip on linear instability (e.g. Lauga & Cossu Reference Lauga and Cossu2005; Chai & Song Reference Chai and Song2019; Ceccacci et al. Reference Ceccacci, Calabretto, Thomas and Denier2022).
Figure 9 shows that the corresponding wall-normal component of velocity perturbations, while satisfying the no-penetration condition at the walls (i.e. 
$\tilde {v} = 0$), is biased towards walls for which the wall slip of the perturbation is aligned with the wall shear exerted on the fluid due to the velocity perturbations.
Figure 9. Wall-normal (i.e. 
$\mathbb {R}[\tilde {v}]$) eigenmode components for eigenvalues with maximum imaginary part in plane Couette flow (
$Re = 360$, 
$S = 1 \times 10^{-1}$) with generalized boundary condition in (2.1) (i.e. (c)) and alternative boundary conditions in (E1) (i.e. (a,b,d)). Large near-wall 
$\tilde {v}$-component perturbations arise near walls with boundary condition enforcing directional alignment of perturbation wall slip and associated wall shear. Specifically, directional alignment of perturbation wall slip and associated wall shear is enforced on (a) top wall, (b) bottom wall, (c) both walls, (d) neither wall. Eigenmodes in (a–c) are unstable; eigenmode in (d) is stable.
Notably, each boundary condition for which the velocity perturbation and wall shear are aligned on at least one wall (i.e. figures 9a–c) corresponds to an unstable eigenmode. In each of the present cases, the growth rate of the perturbation is 
$\mathbb {I}[\hat {\omega }] \approx 0.199$. By contrast, for the boundary condition in figure 9(d), all eigenmodes are stable, in agreement with previous studies (e.g. Lauga & Cossu Reference Lauga and Cossu2005; Chai & Song Reference Chai and Song2019; Ceccacci et al. Reference Ceccacci, Calabretto, Thomas and Denier2022). The least stable eigenmode has growth rate 
$\mathbb {I}[\hat {\omega }] \approx -0.175$.
Open access
