Linear stability of slip pipe flow

The classic pipe flow with no-slip boundary condition has been proved linearly stable to axisymmetric perturbations (Herron, 1991, 2017), and numerical studies suggest that the flow is linearly stable to any perturbations at arbitrary Reynolds numbers (Meseguer & Trefethen, 2003). The recent work of Chen et al. (2019) presented a rigorous proof of the linear stability of the flow to general perturbations at high Reynolds number regime. Therefore, transition to turbulence in pipe flow is subcritical via finite-amplitude perturbations (see e.g. Eckhardt et al. (2007); Avila et al. (2011)).

However, velocity slip of viscous fluid can occur on super-hydrophobic surfaces (Voronov et al., 2008; Rothstein, 2010), for which slip boundary condition instead of the classic no-slip condition should be adopted for the momentum equations, and the slip boundary condition can potentially influence the stability of the flow. A simplified and widely used slip boundary condition is the Navier slip boundary condition, which has been shown to apply to many flow problems and frequently adopted for linear stability studies (Vinogradova, 1999; Lauga & Cossu, 2005; Min & Kim, 2005; Gan & Wu, 2006; Ren et al., 2008; Ghosh et al., 2014; Seo & Mani, 2016; Chattopadhyay et al., 2017, to list a few). For pipe geometry, although many studies have investigated the linear stability of immiscible and miscible multi-fluid flows with either no-slip or Navier slip boundary condition (Hu & Joseph, 1989; Joseph, 1997; Li & Renardy, 1999; Selvam et al., 2007; Sahu, 2016; Chattopadhyay et al., 2017, etc.), much fewer studies were dedicated to the linear stability of single-phase pipe flow with slip boundary condition. Průša (2009) investigated this problem and showed that, subject to Navier slip boundary condition, pipe flow becomes less stable compared to the no-slip case, however, the destabilization effect is constrained to small Reynolds numbers and is not sufficient to render the flow linearly unstable. Their results indicated that the stability property of pipe flow is not qualitatively affected by the slip boundary condition, regardless of the slip length. For its counterpart in plane geometry, i.e. channel flow, on the contrary, Min & Kim (2005); Lauga & Cossu (2005) reported a stabilizing effect of velocity slip on the linear stability.

Usually, slip length is assumed homogeneous and isotropic, i.e. independent of position and direction at the wall in stability analysis. However, anisotropy in the effective slip length can be incurred by anisotropy in the texture pattern on superhydrophobic surfaces, such as parallel periodic slats, grooves and grates (Lecoq et al., 2004; Bazant & Vinogradova, 2008; Ng & Wang, 2009; Belyaev & Vinogradova, 2010; Asmolov & Vinogradova, 2012; Pralits et al., 2017). For example, Ng & Wang (2009) reported a ratio of down to about 0.25 between the transverse slip length (in the direction perpendicular to the slats) and longitudinal slip length (parallel to the slats). The linear stability of channel flow with anisotropic slip caused by parallel micro-graves was analyzed by Pralits et al. (2017) using the the tensorial formulation of slip boundary condition proposed by Bazant & Vinogradova (2008). Their results showed possibilities of linear instability using special alignment of the micro-graves. Recently, Chai & Song (2019) studied the linear stability of single-phase channel flow subject to anisotropy in slip length by considering streamwise and azimuthal slip separately as the limiting cases, which can potentially be realized or approximated by using specially designed surface texture, e.g. specially aligned micro-grates/graves, according to Bazant & Vinogradova (2008). Their results showed that streamwise slip mainly stabilizes the flow (with increased critical Reynolds number), although it surprisingly destabilizes the flow slightly in a small Reynolds number range, and that azimuthal slip can greatly destabilize the flow and reduce the critical Reynolds number given sufficiently large slip length. The critical Reynolds number can be reduced to a few hundred with a dimensionless azimuthal slip length of $\mathcal{O}(0.1)$, in contrast to $Re_{cr}=5772$ for the no-slip case. Their study also indicated that Squire’s theorem (Squire, 1933) ceases to apply when the wall normal velocity and vorticity are coupled via the slip boundary condition, such that the leading instability becomes three dimensional (3-D) rather than two dimensional (2-D) when slip length is sufficiently large, in agreement with Pralits et al. (2017). The stability of 3-D perturbations was not considered by Min & Kim (2005); Lauga & Cossu (2005) in which Squire’s theorem was seemingly assumed.

Differing from channel flow, linear instability is absent at arbitrary Reynolds numbers in classic pipe flow. This raises the question of whether the anisotropy in slip length can also cause linear instability in pipe flow. To our knowledge, this problem has not been studied in pipe geometry. The pseudospectrum analysis of classic pipe flow of Schmid & Henningson (1994); Meseguer & Trefethen (2003) suggests that, despite the linear stability, at sufficiently large Reynolds numbers, a small perturbation to the linear operator associated with the governing equation can possibly change the stability of the system. The slip boundary condition can be thought of as a perturbation to the linear operator with no-slip boundary condition. However, Průša (2009) showed that homogeneous and isotropic slip does not change the spectrum qualitatively no matter how large the slip length (i.e. operator perturbation) is. Following Chai & Song (2019), in this work, we still consider anisotropic slip length in the limiting cases and explore the possibility of linear instability for pipe flow. Aside from the critical Reynolds number as focused on by Chai & Song (2019), here we also investigate the effects of the slip on the spectrum and on the scaling of the leading eigenvalues with Reynolds number. Besides numerical calculations, we also perform analytical studies on the eigenvalues and eigenvectors of the 3-D yet streamwise-independent modes, and discuss about their structure as well as their dependence on the slip length on a theoretical basis, which to our knowledge have not been reported in the literature.

The nondimensional incompressible Navier-Stokes equations read

where $\bm{u}$ denotes velocity and $p$ denotes pressure. For pipe geometry, cylindrical coordinates $(r,\theta,x)$ are considered, where $r$, $\theta$ and $x$ denote the radial, azimuthal, and streamwise coordinates, respectively. Velocity components $u_{r}$, $u_{\theta}$ and $u_{x}$ are normalized by $2U_{b}$ where $U_{b}$ is the bulk speed (the average of the streamwise velocity on the pipe cross-section), length by pipe radius $R$ and time by $R/U_{b}$. The Reynolds number is defined as $Re=U_{b}R/\nu$ where $\nu$ is the kinematic viscosity of the fluid. In order to eliminate the pressure and impose the incompressibility condition, we adopt the velocity-vorticity formulation of Schmid & Henningson (1994), with which the governing equations of disturbances reduce to only two equations about the wall normal velocity $u_{r}$ and wall normal vorticity $\eta$. With a Fourier-Fourier-Chebyshev collocation discretization, considering perturbations of the form of $\{u_{r},\eta\}=\{\hat{u}_{r}(r),\hat{\eta}(r)\}\text{e}^{-i(\alpha x+n\theta)}$, the governing equations in the Fourier spectral space read

where

$\tau=\dfrac{t}{Re}$ is the scaled time, and unknowns are

The real number $\alpha$ is the axial wave number and $n$, which is an integer, is the azimuthal wavenumber. The base flow is denoted as $U$, $k^{2}=\alpha^{2}+\dfrac{n^{2}}{r^{2}}$, $i=\sqrt{-1}$ and the prime denotes the derivative with respect to $r$. The operators $\Gamma$ and $\phi$ are defined as $\Gamma=\dfrac{1}{r^{2}}-\dfrac{1}{r}\dfrac{\text{d}}{\text{d}r}\left(\dfrac{1}{k^{2}r}\dfrac{\text{d}}{\text{d}r}\right)$ and $\phi=k^{4}r^{2}-\dfrac{1}{r}\dfrac{\text{d}}{\text{d}r}\left(k^{2}r^{3}\dfrac{\text{d}}{\text{d}r}\right)$. The other two velocity components $\hat{u}_{x}$ and $\hat{u}_{\theta}$ can be calculated as

We use the Robin-type Navier slip boundary condition at the pipe wall for streamwise and azimuthal velocities separately, i.e.

where $l_{x}\geqslant 0$ and $l_{\theta}\geqslant 0$ are streamwise and azimuthal slip lengths, respectively, and are independent of each other. In spectral space, these boundary conditions apply identically to $\hat{u}_{x}$ and $\hat{u}_{\theta}$ given the homogeneity of the slip length. We use the no-penetration condition for the wall-normal velocity component at the wall, i.e. $u_{r}(1,\theta,x,t)=0$. Lauga & Cossu (2005); Chai & Song (2019) considered the same boundary conditions for slip channel flow. Note that in the isotropic slip case considered by Průša (2009), $l_{x}$ and $l_{\theta}$ are related as $l_{\theta}=\dfrac{l_{x}}{1+l_{x}}$, which gives $l_{\theta}\approx l_{x}$ for small slip lengths. With boundary condition (7), given that we impose the same volume flux as in the no-slip case, i.e.

the velocity profile of the constant-volume-flux base flow reads

where $\hat{\bm{x}}$ represents the unit vector in the streamwise direction. Note that the base flow is independent of $l_{\theta}$. Converting to the $(\hat{\Omega},\hat{\Phi})$ system, the boundary condition (7) reads

and

It can be seen that $\hat{\Omega}$ and $\hat{\Phi}$, i.e. $\hat{u}_{r}$ and $\hat{\eta}$, are coupled via the slip boundary condition.

In order to avoid the singularity at the pipe center, i.e. $r=0$, the domain [0, 1] is extended to [-1, 1] and an even number of Chebyshev grid points over [-1, 1] are used such that there is no grid point at $r=0$. This extension also allows us to use the Chebyshev collocation method for the discretization in the radial direction and the resulted redundancy is circumvented by setting proper parity conditions on $\hat{\Phi}$ and $\hat{\Omega}$ with respect to $r$ (Trefethen, 2000; Meseguer & Trefethen, 2003). In this way, no explicit boundary condition is imposed at the pipe center.

To determine whether a mode $(\alpha,n)$ is linearly stable or not, one only needs to calculate the eigenvalues of the operator $-M^{-1}L$ and check if any eigenvalue has a positive real part, $\lambda_{r}$, which determines the asymptotic growth/decay rate of the corresponding eigenvector as $t\to\infty$.

We consider the case of $l_{x}\neq 0$ and $l_{\theta}=0$ as the limiting case of streamwise slip being significant and azimuthal slip being negligible.

The effect of the slip on the spectrum is investigated for $Re=3000$ and is shown in Figure 1 for the modes $(\alpha,n)=(0,1)$ and $(0.5,1)$. Firstly, panel (a) shows that the eigenvalues of the $(\alpha,n)=(0,1)$ mode visually all fall on the $\lambda_{i}=0$ line ($\lambda_{i}$ denotes the imaginary part of the eigenvalue) and in the left half-plane, which suggests that the eigenvalues are all real and negative. Meseguer & Trefethen (2003) reported the same finding for the no-slip case in a large Reynolds number range up to $10^{7}$. In fact, the eigenvalues being real and negative can be rigorously proved, see our proof in Section 5.1. Secondly, as $l_{x}$ increases, it can be observed that there are two groups of eigenvalue, one of which stays constant and the other of which shifts to the right, see the two insets in panel (a). Specifically, as $l_{x}$ is increased to 0.5, the left eigenvalue in the left inset has moved from the circle to the triangle and finally to the square while the right eigenvalue stays constant. Nontheless, the rightmost eigenvalue increases as $l_{x}$ increases (see the right inset) which indicates that the flow becomes less stable. In Section 5.2, we will show that the former group corresponds to disturbances with $\Phi\not\equiv 0$, i.e. $u_{r}\not\equiv 0$ and the latter group, on the contrary, is associated with disturbances with $\Phi\equiv 0$, i.e. $u_{r}\equiv 0$ and, the rightmost eigenvalue belongs to the latter group (see Figure 13).

Panel (b) shows the case for the mode $(\alpha,n)=(0.5,1)$. The slip does not qualitatively change the shape of the spectrum. As $l_{x}$ increases, the eigenvalues overall move to the right. Besides a horizontal shift, there is a shift in the vertical direction either, and meanwhile the spectrum is compressed in the vertical direction, see the comparison between the $l_{x}=0.5$ and the other two cases. Using the term of Schmid & Henningson (1994); Meseguer & Trefethen (2003), the horizontal branch of the spectrum (the part with $\lambda_{r}\lesssim-600$) corresponds to mean modes, the upper branch corresponds to wall modes and the lower branch to center modes. Note that the speed of a wave is given by $\frac{-\lambda_{i}}{\alpha Re}$ in our formulation. It has been known that the wave speed of the mean modes follows the mean velocity of the ‘two-dimensional’ axial base flow, i.e. $\int_{0}^{1}U_{x}(r)\text{d}r$ in pipe flow (see e.g. Drazin & Reid (1981)), which gives $\frac{2}{3}$ in the no-slip case (Schmid & Henningson, 1994). In our case, the wave speed of the mean modes is decreased by the slip, reducing to 0.5559 for $l_{x}=0.5$ ($\frac{833.868}{0.5\times 3000}$, see the eigenvalue in Table 1) which is very close to $\frac{5}{9}$ given by $\int_{0}^{1}U_{x}(r)\text{d}r$ with the base flow shown in (9). The wall modes, which are located close to the wall, move at lower speed than the center modes, which are located close to the pipe center and move at speeds close to the centerline velocity. Since we fix the volume flux of the flow while the slip length is varied, the speed of the base flow close to the wall increases as $l_{x}$ increases, whereas the speed near the pipe center decreases, i.e. the velocity profile becomes flatter, see the base flow given by (9). Therefore, it can be expected that as $l_{x}$ increases, the speed of the wall-modes increases and that of the center modes decreases, and all three types of modes move at closer speeds. This is exactly what the compression in the vertical direction of the spectrum reveals. The other noticeable effect is that the slip brings the adjacent eigenvalues associated with the mean modes closer as the slip length increases, causing a seemingly degeneracy of the spectrum, see Figure 1(b).

Figure 2 shows the maximum of the real part of the eigenvalue, $\max{\lambda_{r}}$, as a function of the streamwise wavenumber, $\alpha$, for $Re=3000$ and $10^{4}$. For each $Re$, slip lengths $l_{x}=0.1$ and 1.0, and azimuthal wavenumbers $n=0$, 1, 2, 3 and 4 are considered. The trend shown in the figure suggests that, for both Reynolds numbers, $\alpha=0$ is nearly the least stable mode, i.e. the slowest decaying mode given that all $\max{\lambda_{r}}$’s are negative, regardless of the slip length. At small $\alpha$, where $\max{\lambda_{r}}$ is largest, the results suggest that $n=1$ is always the least stable one. At larger $\alpha$, however, $n=1$ is still the least stable when $l_{x}$ is small, see the case of $l_{x}=0.1$ in panel (a, c), but is not in a range of $\alpha$ around $\alpha=1$, see the case of $l_{x}=1.0$ in panel (b, d). Nevertheless, in this range, $\max{\lambda_{r}}$ is much smaller than that in the small $\alpha$ regime. Therefore, as we are most interested in the least stable mode, in the following, we will focus on the $n=1$ modes. In fact, for the $\alpha=0$ modes, we can rigorously prove that $n=1$ is the least stable azimuthal wavenumber, see Appendix B.

Figure 3 shows $\max{\lambda_{r}}$ as a function of $\alpha$ of the $n=1$ modes for $Re=3000$ (a) and $10^{4}$ (b). For each $Re$, overall $\max{\lambda_{r}}$ increases as $l_{x}$ increases, i.e. the $n=1$ modes decay more slowly as $l_{x}$ increases. The insets show the close-up of the small $\alpha$ region, in which the dependence of $\max{\lambda_{r}}$ on $\alpha$ is not monotonic, with the maximum appears at some small but finite $\alpha$ instead of $\alpha=0$. Nevertheless, the difference between the peak value and the value for $\alpha=0$ is very small, i.e. $\alpha=0$ is nearly the least stable mode, as aforementioned. In fact, the dependence on $l_{x}$ is not fully monotonic either, see the very small region around $\alpha=0.03$ for the $l_{x}=0.005$ (the thin black line) and $l_{x}=0.05$ (the blue line) cases as shown in the inset in (a) and around $\alpha=0.01$ in the inset in (b). However, for $\alpha=0$ and in most range of $\alpha$, our results show a monotonic increase of $\max{\lambda_{r}}$ as $l_{x}$ increases.

Figure 4 illustrates the dependence of $\max{\lambda_{r}}$ of the $n=1$ modes on $l_{x}$ in a broader range of $l_{x}$. For each $Re$, $\alpha=0$, 0.1, 0.5, 1 and 2 are shown. The trend shows that as $l_{x}$ keeps increasing, $\max{\lambda_{r}}$ seems to asymptotically approach a plateau with a negative value, i.e. all the modes shown in the figure appear to be linearly stable, for both Reynolds numbers.

The above results suggest that, with streamwise slip, the flow is linearly stable to any perturbations, regardless of the slip length. In order to show evidences in a broader parameter regime, we numerically searched for the global maximum of $\max{\lambda_{r}}$ over $\alpha$ and $n$ and explored a wider range of $l_{x}$ up to 10 and of $Re$ up to $10^{6}$. Practically, based on our analysis, we only need to search in a small range of $\alpha$ immediately above zero (see the insets in Figure 3) while setting $n=1$. Specifically, the range of [0, 1.2] is searched at $Re=100$, and the range is decreased as $Re^{-1}$ for higher Reynolds numbers. Then we plotted the global maximum of $\max{\lambda_{r}}$, still denoted as $\max{\lambda_{r}}$, as a function of $l_{x}$, for a few Reynolds numbers ranging from $10^{2}$ to $10^{6}$ in Figure 5(a).

It is interesting to note that our data for high Reynolds numbers all collapse over the whole $l_{x}$ range investigated, see the cases with $Re$ above $1\times 10^{4}$ in Figure 5, suggesting that the maximum eigenvalue of the system is independent of $Re$. At lower Reynolds numbers, e.g. $Re=100$ and $10^{3}$ in the figure, there is almost a collapse for small $l_{x}$ ($\lesssim 0.1$) but a small deviation from the high Reynolds number cases can be seen, see the inset that shows the zoom-in at $l_{x}=0.005$. As $l_{x}$ increases further, the maximum eigenvalue for $Re=100$ and $10^{3}$ approaches to that of $\alpha=0$ modes, which is strictly $Re$-independent (see the proof in Section 5.1). Besides, the figure also shows that the global maximum of $\max{\lambda_{r}}$ is slightly larger than the maximum of the $\alpha=0$ modes over the whole $l_{x}$ range and the difference is most significant at small $l_{x}$. We did not explore further larger $l_{x}$ considering that the range we investigated is already much larger than the slip length that can be encountered in applications ($\lesssim 0.1$ in set-ups with characteristic length of one millimeter or larger, because so far the maximum slip length achieved in experiments is $\mathcal{O}(100)$ micron, see Voronov et al. (2008); Lee et al. (2008); Lee & Jim (2009)). Nevertheless, the $S$-shaped trend as $l_{x}$ increases suggests that the flow stays stable no matter how large the slip length is. In fact, as $l_{x}\to\infty$, full slip boundary condition is recovered, and the velocity profile of the base flow will be completely flat and no mean shear exists, in which case linear stability can be expected for any perturbations. Panel (b) shows the product of $Re$ and the $\alpha$ at which $\max{\lambda_{r}}$ maximizes globally, denoted as $\alpha_{\max{\lambda_{r}}}$. Interestingly, it seems that this product is a constant when $l_{x}$ is small ($\lesssim 0.1$) for all the $Re$’s investigated and approaches a constant as $Re$ is sufficiently high ($\gtrsim 10^{4}$) if $l_{x}\gtrsim 0.1$. This indicates that $\alpha_{\max{\lambda_{r}}}$ scales as $Re^{-1}$ for either not very large $l_{x}$ or in high Reynolds number regime. It should be noted that we observed a non-monotonic dependence of $\alpha_{\max{\lambda_{r}}}\cdot Re$ on the slip length, which minimizes at around $l_{x}=0.1$.

That the global $\max{\lambda_{r}}$ is $Re$-independent, as our results suggest, indicates that the slowest exponential decay rate (referred to as decay rate for simplicity hereafter) of perturbations scales as $Re^{-1}$ given that the scaled time $\tau=\frac{t}{Re}$ is used in our formulation, see Eqs. (2). The same scaling was observed by the calculation of Meseguer & Trefethen (2003) for the $(\alpha,n)=(0,1)$ mode of the classic pipe flow. Therefore, our results suggest that, as $Re\to\infty$, the decay rate of perturbations asymptotically approaches zero and stays negative, i.e. the flow is linearly stable at arbitrary Reynolds number. The $Re^{-1}$-scaling of the slowest decay rate can be rigorously proved for the $\alpha=0$ modes, see Section 5.1.

In a word, in the pure streamwise slip case, we did not observe any linear instability in the large ranges of $l_{x}$ and $Re$ that we considered, and based on the data shown in Figure 5, we propose that streamwise slip destabilizes the flow but does not cause linear instability, regardless of the slip length and Reynolds number. A similar destabilizing effect was reported by Průša (2009) for the isotropic slip case.

We consider the case of $l_{\theta}\neq 0$ and $l_{x}=0$ as the limiting case of azimuthal slip being significant and streamwise slip being negligible.

The effect of azimuthal slip on the spectrum is investigated for $Re=3000$ and is shown in Figure 6 for the modes $(\alpha,n)=(0,1)$ and $(0.5,1)$. Similar to the streamwise slip case, the eigenvalues of the $\alpha=0$ mode also fall on the $\lambda_{i}=0$ line and in the left half-plane, see panel (a). This suggests that the eigenvalues of streamwise-independent modes are all real and negative. We will show a rigorous proof of this observation in Section 5.1. As $l_{\theta}$ increases, similar to the streamwise slip case, we also observed two groups of eigenvalues. One group stays constant as the azimuthal slip length changes and the other shifts to the right, see the inset in panel (a). As we will theoretically show in Section 5.2 and 5.3, the former group is associated with the disturbances with $\Phi\equiv 0$ and is independent of $l_{\theta}$, and the latter group is associated with the disturbances with $\Phi\not\equiv 0$. The rightmost eigenvalue belongs to the former group for $l_{\theta}<1$ and can only be overtaken by the latter group if $l_{\theta}>1$ (the two groups precisely overlap when $l_{\theta}=1$), i.e. the rightmost eigenvalue can only increase with $l_{\theta}$ if $l_{\theta}>1$. For the $\alpha=0.5$ and $n=1$ mode, the mean-mode branch overall does not show either a vertical or horizontal shift, but adjacent eigenvalues are brought closer by the increasing slip length, and for $l_{\theta}=0.5$ there is almost an eigenvalue degeneracy (see the inset in panel (b)). The center-mode branch is nearly unchanged as $l_{\theta}$ increases. However, the wall-mode branch is significantly affected. As $l_{\theta}$ increases, the wall mode overtakes the center mode and becomes the least stable perturbation, and as $l_{\theta}$ is sufficiently large, the rightmost eigenvalue appears in the right half-plane, indicating the onset of a linear instability. In contrast to the streamwise slip case, the wave speed of the mean modes does not change because the speed follows $\int_{0}^{1}U_{x}(r)\text{d}r$ as aforementioned and the base flow $\bm{U}(r)$ is not affected by the azimuthal slip. The speed of the center modes is also not affected, whereas the wave speed of the wall modes is considerably decreased by the slip. This is reasonable because the slip boundary condition should mostly affect the flow close to the wall and should not affect significantly the flow far from the wall.

Figure 7 shows $\max{\lambda_{r}}$ maximized over $\alpha$ (over [0, 2] in practice), still denoted as $\max{\lambda_{r}}$, as a function of $l_{\theta}$ for $n=0$, 1, 2, 3 and 4 at $Re=3000$. Overall, $\max{\lambda_{r}}$ increases monotonically as $l_{\theta}$ increases, while the $n=0$ case seems to stay constant until it starts to increase at around $l_{\theta}=0.4$. In the small $l_{\theta}$ regime, all modes are linearly stable. As $l_{\theta}$ is increased to around 0.1, $\max{\lambda_{r}}$ of the $n=1$ mode becomes positive, indicating a linear instability. As $l_{\theta}$ increases further, $n=2$ and 3 also become unstable. In the whole range of $l_{\theta}$ investigated, $n=1$ is the least stable/most unstable one, which is also the case for other Reynolds numbers we investigated. Therefore, in the following, we mainly discuss about $n=1$ modes.

Figure 8 shows $\max{\lambda_{r}}$ of modes $\alpha=0.1$, 0.5, 1.0 and 2.0 for $Re=3000$ and $n=1$ as a function of $l_{\theta}$. The results show that when $l_{\theta}$ is small, overall $\max{\lambda_{r}}$ decreases as $\alpha$ increases. As $l_{\theta}$ is increased, some moderate $\alpha$ turns to be the least stable/most unstable mode, see the crossover of $\alpha=0.1$ (cyan thin line) and 0.5 (red dashed line) cases in the figure. Panel (b) shows the small $l_{\theta}$ range, in which it appears that $\max{\lambda_{r}}$ first stays nearly unchanged and then starts to increase, and the trend shows that the larger $\alpha$, the later $\max{\lambda_{r}}$ starts to increase as $l_{\theta}$ is increased. The same behavior is also observed for $\alpha=0$ modes and we will show a rigorous proof of this behavior in Section 5.3. Interestingly, the case of $\alpha=2$ seems to stay unchanged up to $l_{\theta}=2.0$.

The dependence of $\max{\lambda_{r}}$ on $\alpha$ is more comprehensively shown in Figure 9. The smallest $l_{\theta}=0.005$ shows a monotonic decrease with increasing $\alpha$, which completely collapses onto the curve for $l_{\theta}=0$, i.e. the classic no-slip case. However, as $l_{\theta}$ increases, $\max{\lambda_{r}}$ significantly increases in the region of $0<\alpha\lesssim 1$ such that a bump appears in the curves, see those for $l_{\theta}=0.05$, 0.1 and 0.5. At certain point, the bump reaches $\max{\lambda_{r}=0}$ and the flow starts to become linearly unstable if $l_{\theta}$ increases further, see the cases of $l_{\theta}=0.1$ and 0.5. As observed in Figure 8 for the $\alpha=2$ case, the results suggest that $\max{\lambda_{r}}$ of sufficiently large $\alpha$ seems unaffected by azimuthal slip in the $l_{\theta}$ range investigated, see the collapse of all curves above $\alpha\simeq 1.2$ in Figure 9. It should be noted that as $l_{\theta}$ becomes larger, the bump widens up, i.e. $\max{\lambda_{r}}$ is affected by the slip in a wider range of $\alpha$.

Figure 10 shows the velocity field of the most unstable perturbation of mode $(\alpha,n)=(0.383,1)$ for $Re=3000$ with $l_{\theta}=0.1$. Panel (a) shows the in-plane velocity field in the $r-\theta$ pipe cross-section and (b) shows the pattern of the streamwise velocity in the $x-r$ cross-section. The patterns shown suggest that the flow manifests with a pair of helical waves. The flow structures are mostly located near the wall ($r\gtrsim 0.5$), indicating that the most unstable mode is a wall-mode, which can also been seen in Figure 6(b).

Obviously, azimuthal slip can cause linear instability given sufficiently large slip length. We can search in the $l_{\theta}$-$\alpha$ plane to obtain the neutral stability curve for given $Re$ and $n$. Figure 11(a) shows the neutral stability curves for $Re=3000$ and $n=1$, 2 and 3 ($n=4$ and higher are all stable, see Figure 7). It can be seen that, as $n$ increases, the unstable region shifts to the right and upward. However, as the results in Figure 7 show, $n=1$ is the most unstable based on the eigenvalue maximized over $\alpha$ and therefore, we only investigate the $n=1$ case in the following. Panel (b) shows the neutral stability curves for $n=1$ and $Re=3000$, 5000 and 7000. As $Re$ increases, the neutral stability curve moves towards the smaller $l_{\theta}$ region, indicating that, for a given $l_{\theta}$, the flow becomes more unstable as $Re$ increases, as expected. The data show that the wavelength of the unstable modes is comparable or significantly larger than the pipe diameter, whereas very long waves ($\alpha\to 0$) and short waves ($\alpha\gg 1.0$) are generally stable. That the flow is always stable to perturbations with $\alpha=0$, regardless of the value of $l_{\theta}$ and $Re$, can be rigorously proved (see Section 5.1).

Further, for each $l_{\theta}$, a critical Reynolds number $Re_{cr}$ can be determined by searching the first appearance of a positive $\max{\lambda_{r}}$ in the $l_{\theta}$-$\alpha$ plane by varying $Re$. Figure 11(c) shows $Re_{cr}$ as a function of $l_{\theta}$. As shown, $Re_{cr}$ is a few hundred if $l_{\theta}$ is large ($l_{\theta}\gtrsim 0.3$), but the trend suggests that it does not reduce to zero if $l_{\theta}\to\infty$. Since the classic pipe flow is linearly stable for arbitrary Reynolds number, there is an explosive increase in $Re_{cr}$ as $l_{\theta}$ decreases, which can be expected because the classic pipe flow will be recovered if $l_{\theta}\to 0$. We also explored the limit of $l_{\theta}\to\infty$, in which case the boundary condition for the azimuthal velocity becomes the full slip condition of

The neutral stability curve for $n=1$ in the $Re-\alpha$ plane is shown in Figure 12, which shows that the unstable modes are still long waves with $\alpha$ approximately between 0 and 0.8. The critical Reynolds number (the nose of the curve) appears approximately at $Re_{\text{cr}}=260$.

We can rigorously prove the linear stability of the base flow to perturbations with $\alpha=0$. In the following, we do not consider the $(\alpha,n)=(0,0)$ mode, which should be strictly stable as it is purely dissipative and there can be no energy production mechanism associated with it. In fact, the stability of the classic pipe flow to streamwise independent perturbations has already been proved by Joseph & Hung (1971) using an energy analysis. Nevertheless, here we also account for the effect of the velocity slip and perform analytical studies on the eigenvalues and eigenvectors of $\alpha=0$ modes.