Inertial enhancement of the polymer diffusive instability

Abstract Beneitez et al. (Phys. Rev. Fluids, vol. 8, 2023, L101901) have recently discovered a new linear ‘polymer diffusive instability’ (PDI) in inertialess rectilinear viscoelastic shear flow using the finitely extensible nonlinear elastic constitutive model of Peterlin (FENE-P) when polymer stress diffusion is present. Here, we examine the impact of inertia on the PDI for both plane Couette and plane Poiseuille flows under varying Weissenberg number ${W}$, polymer stress diffusivity $\varepsilon$, solvent-to-total viscosity ratio $\beta$ and Reynolds number ${Re}$, considering the FENE-P and simpler Oldroyd-B constitutive relations. Both the prevalence of the instability in parameter space and the associated growth rates are found to significantly increase with ${Re}$. For instance, as $Re$ increases with $\beta$ fixed, the instability emerges at progressively lower values of $W$ and $\varepsilon$ than in the inertialess limit, and the associated growth rates increase linearly with $Re$ when all other parameters are fixed. For finite $Re$, it is also demonstrated that the Schmidt number $Sc=1/(\varepsilon Re)$ collapses curves of neutral stability obtained across various $Re$ and $\varepsilon$. The observed strengthening of PDI with inertia and the fact that stress diffusion is always present in time-stepping algorithms, either implicitly as part of the scheme or explicitly as a stabilizer, implies that the instability is likely operative in computational work using the popular Oldroyd-B and FENE-P constitutive models. The fundamental question now is whether PDI is physical and observable in experiments, or is instead an artifact of the constitutive models that must be suppressed.


Introduction
The addition of polymers to a Newtonian solvent can induce dramatically different flow behaviours compared to those observed in the Newtonian fluid alone (Datta et al. 2022;Sánchez et al. 2022).In industrial processes, for instance, viscous polymer melts are susceptible to instabilities which constrain the maximum extrusion rate (Petrie & Denn 1976), while polymer additives are used in oil pipelines to reduce turbulent wall drag (Virk 1975).Two particularly important viscoelastic phenomena are the existence of 'elastic turbulence' (ET), a chaotic flow state sustained in the absence of inertia (Groisman & Steinberg 2000;Steinberg 2021), and 'elasto-inertial turbulence' (EIT), an inherently two-dimensional state arising when both inertia and elasticity are present (Samanta et al. 2013;Sid et al. 2018;2 Choueiri et al. 2021).While the initial pathway to ET in curvilinear geometries is understood (Larson et al. 1990;Pakdel & McKinley 1996;Shaqfeh 1996;Datta et al. 2022), relatively little is known about what happens in rectilinear viscoelastic flows.
Initial progress in characterizing ET in rectilinear situations arose through consideration of Kolmogorov flow over a 2-torus, where Boffetta et al. (2005) found a linear instability leading to ET (Berti & Boffetta 2010).Garg et al. (2018) subsequently discovered a centremode instability in viscoelastic pipe flow at finite Reynolds number , which was later also identified in plane Poiseuille flow (PPF, hereafter referred to as channel flow) (Khalid et al. 2021a) but notably not in plane Couette flow (PCF).Interestingly, this instability could only be traced down to  = 0 in channel flow (Khalid et al. 2021b;Buza et al. 2022b).The finite-amplitude state resulting from this instability is an 'arrowhead' travelling wave (Page et al. 2020;Buza et al. 2022a;Morozov 2022) which has been observed in channel flow EIT (Dubief et al. 2022) and, in retrospect, ET in 2D Kolmogorov flow (Berti & Boffetta 2010).In channel flow, efforts have began to establish a dynamical link between the arrowhead solution and both ET (Lellep et al. 2023) and EIT (Beneitez et al. 2023a).In the latter case in two-dimensions, there does not appear to be a simple dynamical pathway between these arrowhead solutions, where the dynamics is concentrated near the midplane, and EIT (Beneitez et al. 2023a), which seems more dependent on a near-wall mechanism (Shekar et al. 2019(Shekar et al. , 2021;;Dubief et al. 2022).
The very recent discovery of a new wall-mode "Polymer Diffusive Instability" (PDI) in plane Couette flow at  = 0 (Beneitez et al. 2023b), however, has added another intriguing possibility for the origin of ET.This instability is dependent on the existence of small but non-vanishing polymer stress diffusion which is invariably present in any timestepping algorithm, whether added explicitly to stabilise a numerical scheme like a spectral method (see e.g.Dubief et al. (2023)) or arising implicitly such as through a finite difference formulation (see e.g.Zhang et al. (2015); Pimenta & Alves (2017)).The PDI wall mode is primarily confined to a boundary layer of thickness

√
, where  is the (small) diffusion coefficient, traveling at the wall speed with a streamwise wavelength on the order of the boundary layer thickness.The instability is robust to the choice of boundary conditions applied to the polymer conformation equation, and has growth rates which remain  (1) as  → 0. Direct numerical simulations have demonstrated that PDI can lead to a sustained three-dimensional turbulent state, thus providing a potential mechanism for the origin of an ET-like state in the popular finitely-extensible nonlinear elastic constitutive model of Peterlin (FENE-P) (Beneitez et al. 2023b).
While PDI has the potential to be a viscoelastic instability of significant importance, there is an important caveat: the instability emerges at small length scales which can, depending on the parameters, approach the order of the polymer gyration radius (Beneitez et al. 2023b), violating the continuum approximation.There is thus a question of whether the instability is physical or actually an artifact of the Oldroyd-B and FENE-P models.Either possibility has important implications: if the instability is a physical phenomenon then it provides a pathway to ET and EIT, albeit one which will likely be challenging to establish experimentally due to the small length scales involved; or, it is an artificial feature of the popular FENE-P model, which can compromise its predictions.It thus appears important to now establish the prevalence of PDI across a much wider region of parameter space than was considered in the initial study of Beneitez et al. (2023b).Therefore, we here map out the regions where PDI is operative at finite , considering both plane Couette flow ('PCF') and the more experimentally-relevant channel flow scenarios.Both the prevalence of PDI and the associated growth rates are found to significantly increase at finite  and are relatively insensitive to the bulk flow geometry.PDI is therefore a candidate to trigger both ET and EIT in simulations using the FENE-P model.

Formulation
We consider the following dimensionless equations governing the flow of an incompressible viscoelastic fluid: where u = (, , ) and  denote the velocity and pressure fields, respectively, and  denotes the polymeric contribution to the stress tensor.Following Beneitez et al. (2023b), the equations have been nondimensionalized using the channel half-width  and a characteristic flow speed  0 , taken to be the wall speed in plane Couette flow or the centerline velocity of the base channel flow.In (2.1a), the Reynolds number  :=  0 /  describes the ratio of inertial to viscous forces (with   =   +   denoting the total kinematic viscosity comprised of solvent and polymer components), and  :=   /  denotes the ratio of solvent to total viscosity.The polymeric stress tensor  may be described in terms of the polymer orientation through the conformation tensor c as in (2.2).We emphasize that the inclusion of a polymer stress diffusion term ∇ 2 c (associated with diffusivity  := () −1 , where  denotes the Schmidt number) is the crucial ingredient for the polymer diffusive instability (PDI) identified by Beneitez et al. (2023b).In the inertialess limit,  = Sc −1 when the governing equations are non-dimensionalized using viscous scales.
To close equations (2.1-2.2), the FENE-P constitutive model is used to relate  and c: where I is the identity matrix,  denotes the maximum extensibility of the polymer chains and  :=  0 /, the Weissenberg number, describes the ratio of time-scales for polymer relaxation () to the flow.In the limit  → ∞, the simpler Oldroyd-B model is obtained.Inspection of (2.1-2.3)reveals five parameters of interest governing the flow dynamics: , , , , .
We analyze the linear stability of (2.1-2.3) by perturbing them about their base state: u = U + u * ,  =  +  * ,  = T +  * , c = C + c * .Our coordinate system (, , ) is aligned with the streamwise, wall-normal, and spanwise directions of the channel, respectively.The base state  (),    (),    (),   (),    () satisfies where primes indicates derivatives in the wall-normal () direction.We use pressure gradients    = 0 and    = −2 for plane Couette and channel flow, respectively.While we here compute the base flow accounting for the presence of a finite diffusivity  in (2.4b), it is worth noting that the presence of diffusion in the base flow is not required to induce PDI; the instability still arises if  = 0 in (2.4b) as the inclusion of  ≠ 0 does not significantly change the base state except in the limit of large  = O (1).
Normal mode solutions of the perturbed flow are sought using the ansatz  * (, , ) = φ ()   ( − ) , where real-valued  denotes the streamwise wavenumber and  =   +   is a complex wavespeed, with instability arising if   > 0. The perturbed state is governed by These parameters are chosen to roughly lie on the corresponding neutral curves plotted in Figure 2a.Results using bases of  = 300 and 400 Chebyshev polynomials are depicted by blue and red markers, respectively, to demonstrate convergence of the discrete unstable PDI modes highlighted in the zoomed inset above each panel.
the following system of seven equations for ( ũ, ṽ, p, c , c , c , c ): By expanding the variables in (2.5) in terms of a basis of Chebyshev polynomials, an eigenvalue problem is obtained which may be solved using standard Python libraries (Beneitez et al. 2023b).A basis of  = 300 polynomials was used here to target the eigenvalues associated with PDI via an inverse iterations algorithm, which yielded sufficient convergence.Representative eigenvalue spectra are shown in Figure 1 for plane Couette (PCF) and channel (PPF) geometries, illustrating the results using bases of both  = 300 and  = 400 Chebyshev polynomials to demonstrate convergence.The unstable PDI mode emerges with characteristic wavespeeds in the vicinity of |  | ≈ 1 and   ≈ 0 for PCF and PPF geometries, respectively.
Following Beneitez et al. (2023b), we choose boundary conditions that set  = 0 at the walls in the governing equations to match the typical configuration used in direct numerical simulations (DNS).While  ∼ O 10 6 is characteristic of polymer diffusion in a typical solvent like water (yielding  = 1/() ≈ 10 −9 − 10 −6 depending on ), in DNS much larger values of  ≈ O (10 −3 ) are often employed to achieve numerical stability (see e.g.§2.2 of Dubief et al. (2023)).Thus, it is commonplace to set  = 0 at the walls in order to recover the idealised limit  → 0 (Sureshkumar et al. 1997;Samanta et al. 2013;Dubief et al. 2022).

Results of linear stability analysis
The results of our linear stability analysis are presented for both the Oldroyd-B ( §3.1) and FENE-P ( §3.2) constitutive relations, delineating the trajectory of neutral curves (characterized by growth rate  :=   = 0) associated with PDI in the four-dimensional parameter space spanned by (, , , ).For a given point in parameter space, we perform a sweep over wavenumbers  to ensure that we are tracking the most unstable mode.The results for plane-Couette (PCF) and channel (PPF) geometries are overlaid on the same axes using blue and red colouring, respectively.To promote a collapse of the curves for both geometries, the Weissenberg () axis is scaled by the respective wall shear rates:  ′ wall = {1 (PCF), 2 (PPF)}.As described in §1, the PDI is a wall mode confined to a boundary layer of thickness O ( √ ) and so its behaviour will primarily be influenced by the wall shear.However, we will demonstrate that this collapse of geometries fails in certain scenarios (such as for large ), when the boundary layer grows sufficiently large to be influenced by the non-uniform shear profile away from the wall.In the case of finite , we also demonstrate that curves associated with variable  and  may be collapsed based on  = 1/().

Oldroyd-B
The trajectory of neutral curves in the -- volume are presented in Figure 2 for various  ∈ [0.7, 0.98].Figure 2a considers the - plane at a fixed  = 10 −5 , while Figure 2b considers the - plane at three fixed  ∈ {0, 1000, 5000}.Regions of stability and instability are found to the left and right of each curve, respectively.For a given , increasing  > 0 is found to promote instability at progressively lower values of both  (Figure 2a) and  (Figure 2b).Notably, in the inertialess limit ( = 0), Beneitez et al. (2023b) found that the neutral curves tracked a roughly fixed  ∼ /(1 − ) for  ≲ 10 −2 (see their Figure 2b) leading to the disappearance of PDI at finite  in the limit  → 1.Here, for  > 0, our Figure 2b demonstrates that while for small  = O (10 −7 ) the neutral curves do indeed still track a constant , they then begin to significantly deviate to lower  beyond  ≳ 10 −6 , with the  = 5000 case deviating at lower  than the  = 1000 case.Inertial effects thus play a significant role in promoting a greater prevalence of PDI in the parameter space.In Figure 2d, we demonstrate that the  = {1000, 5000} curves plotted in Figure 2b may be collapsed for both geometries based on the Schmidt number  = 1/().
The streamwise wavenumbers  associated with the neutral curves in Figure 2b are plotted in Figure 2c as a function of .At  = 0,  follows the 1/ √  scaling reported by Beneitez et al. (2023b) for all , whereas for  > 0,  deviates significantly from this scaling to plateau to a roughly constant value for  ≳ 10 −5 , with the unstable modes at higher  being more tightly confined to the wall (as indicated by a larger ).This deviation from the 1/ √  scaling for  > 0 corresponds to the previously noted deviation of the neutral curves away from a constant  in Figure 2b, where the curves begin turning to lower  in the vicinity of  ≈ 10 −6 − 10 −5 .The behaviour of  in Figure 2c also explains the mismatch in the collapse of the PCF and PPF curves in Figure 2b at higher  for  = 0 (see solid lines).Specifically, at higher ,  becomes sufficiently small at  = 0 such that the instability is no longer strongly confined to the wall (see square eigenfunction, Figure 2e) and so the local wall shear,  ′ wall , does not accurately describe the non-uniform shear profile influencing the instability.Conversely, the curves for  = {1000, 5000} do remain collapsed at high , as the wavenumbers  in Figure 2c plateau to sufficiently large values such that the instability remains confined to the wall (see triangle eigenfunction, Figure 2e).In Figure 2c, it is also worth noting that the prefactor of the PPF scaling (1/5) is roughly twice that of the PCF scaling (1/8), thus indicating that PDI is more tightly confined to the wall in the channel geometry.

FENE-P
We now consider how a finite polymer extensibility  modifies the behaviour of the instability as compared to the Oldroyd-B case ( → ∞) presented in §3.1.Focusing first on  = 200, Figure 3 illustrates the much richer behaviour of the neutral curves for the FENE-P case, presented in an analogous manner to Figures 2a-d.In the - plane (Figure 3a), a finite  introduces two notable differences compared to Oldroyd-B.First, the neutral curves for a given  now have a left and right branch and so the range of instability is bounded by an upper value of .Second, there is now a critical value of  (≈ 0.865 for both geometries) above which the neutral curves no longer intersect the  = 0 axis, at fixed  = 10 −5 .Therefore, inertial effects are once again demonstrated to promote PDI, generating instability at finite  for ultra-dilute polymer solutions ( → 1) that would otherwise remain stable at  = 0.
Neutral curves in the - plane for the inertialess case  = 0 are shown in Figure 3b, exhibiting notable differences to the Oldroyd-B  = 0 curves presented in Figure 2b.While for Oldroyd-B the plane Couette and channel curves adopt roughly the same shape, the FENE-P curves display dramatically different behaviours at large .Specifically, the plane Couette curves have an inverted 'U'-shape, highlighting that the instability ceases to exist at large  for certain  (e.g.see the  = 0.86 curve which reaches a maximum value at  ≈ 6 × 10 −3 ).Conversely, the channel curves are roughly 'U'-shaped, and the range of instability only increases with increasing .
At finite , as for the Oldroyd-B curves in Figure 2d, the FENE-P curves shown in Figure 3c at  = {1000, 5000} are again found to collapse for both geometries based on the Schmidt number  = 1/().Comparing Figures 2d and 3c reveals two main differences in the structure of the FENE-P and Oldroyd-B curves.First, for  ≈ 0.925, in both geometries a pinch-off phenomenon occurs for the FENE-P case where the neutral curves form a "bubble" of instability within the  −  plane in an otherwise stable region.Second, in the limit of  → 0 ( → ∞), the neutral curves behave differently as compared to the  = 0 case reported by Beneitez et al. (2023b).In particular, there now appears to be a critical value of  ≈ 0.865 (note the similarity to the critical  in Figure 3a corresponding to lift-off from the  = 0 axis) at which the two branches of the neutral curve asymptotically approach each other in the limit of small .Curves associated with  greater than this critical value are thus observed to turn back at higher  (see e.g. the  = 0.9 curves, Figures 3c), suggesting that, at finite , PDI will not exist in the limit  → 0 for all .We note another possibility, however, which is that an 'hourglass'-like pinch-off behaviour occurs, in which the critical curves touch at some finite  (appearing here to be around  ≈ 10 4 ), but then separate again for lower .The concave-up neutral curves (e.g. = 0.9) seen here, might then be reflected as concave-down branches at much lower  and instabilities for such  could then still exist in the  → 0 limit.As  → 0 isolates the instability to an increasingly thin boundary layer at the wall, significantly increased computational power is required to resolve the neutral curves for  > 10 4 and so we do not further consider this behaviour here.It is also worth noting in Figures 3c that while increasing  does not significantly influence the horizontal position of the neutral curves along the  axis, the increase in inertial effects does induce a significant downward shift of the curves in  (by a factor proportional to ), thus increasing the prevalence of PDI in parameter space.
The streamwise wavenumbers  associated with the FENE-P channel (PPF) neutral curves at  = 1000 in Figure 3c are presented in Figure 3d as a function of , to compare with the scaling of the Oldroyd-B curves in Figure 2c.As for Oldroyd-B, the left branches follow the 1/

√
scaling reported by Beneitez et al. (2023b) for  ≲ 10 −5 , when the neutral curves are roughly independent of  in Figure 3c.In contrast, the right branches follow this  scaling for the entire range of  considered here, thus explaining the slight mismatch in the collapse of the right branches of the two geometries at low  in Figure 3c (see e.g. the right branches of the  = 0.7 and 0.8 curves); in this regime,  has become sufficiently small such that the instability is no longer confined to the wall and thus the wall shear  ′ wall is not entirely suitable for scaling the  axis.In Figure 3d, we also note that for  ≳ 0.864, the pinchoff value seen in Figure 3c at  = O (10 4 ), the neutral curves turn back to higher  before the 1/

√
scaling regime is reached.
In Figure 4, we further consider the behaviour of the neutral curves presented in Figure 3 by varying parameters that were previously held fixed.While the - plane was considered with a fixed  = 10 −5 in Figure 3a, in Figure 4a we focus on the  = 0.87 neutral curve and demonstrate its dependence on variable , where it is found that decreasing  pushes the region of PDI to progressively higher .In Figure 4b, we collapse the curves from Figure 4a based on the Schmidt number  = 1/().A near perfect collapse is observed for the plane Couette (PCF) curves, while some channel (PPF) curves diverge to  → ∞ due to their intersection with the  = 0 axis in Figure 4a.Additionally, having only considered a fixed polymer extensibility  = 200 in Figure 3, in Figures 4c,d we consider the influence of variable  on the neutral curves in the - and - planes, respectively.In Figure 4c, decreasing  induces a turn-back of the neutral curves at progressively higher values of , decreasing the extent of PDI in the - plane.Figure 4d, at fixed  = 1000, displays the same qualitative features as those reported in the inertialess limit by Beneitez et al. (2023b) (see their Figure 2c), again demonstrating that decreasing  decreases the prevalence of PDI in the - plane.

Growth Rates
While we have thus far considered the behaviour of the neutral curves associated with PDI, it is also informative to consider how the growth rate of PDI evolves with  in regions of instability.Starting with each  curve in Figure 2a (Oldroyd-B), we first fix the value of  at which the neutral curve intersects the  = 0 axis, and where the growth rate is thus zero by definition.At this fixed , we then increase  incrementally and track the growth rate of the most unstable PDI mode, as shown in Figure 5a.In Figure 5b, we repeat this process for the FENE-P curves presented in Figure 3a.As some FENE-P curves do not intersect the  = 0 axis, in these cases we fix  to correspond to the minimum  of the neutral curve, and then increase  from there.In Figures 5c-d, we again repeat this process for the neutral curves plotted in Figures 4c-d, respectively.
In all cases, for neutral curves that intersect the  = 0 axis, the growth rate is observed to grow linearly with  as one moves away from the neutral curve, emphasizing the intensification of PDI due to the presence of inertia.Notably, the streamwise wavenumber  remains virtually constant during this scaling, until  ≳ 10 3 at which point the most unstable  beings to vary and the linear scaling breaks down.For the subset of FENE-P curves that do not intersect the  = 0 axis (e.g. = 0.87, 0.88, 0.90 in Figure 4b), it is also down when the streamwise wavenumber  approaches O (1), as the instability is no longer confined to the wall and thus feels a non-monotonic shear profile in the channel's interior, as occurs for large  at  = 0 (but notably not at higher , see Figure 2c), and small  at high  in FENE-P fluids.The finite extensibility of the polymer chains () is also found to have a significant impact on the prevalence of PDI, as compared to that predicted by the Oldroyd-B model.Considering  = 200, we found that for sufficiently high , the instability is suppressed at  = 0 and only appears at progressively larger  (Figure 3a).Similarly, beyond a critical value of , the instability may also be suppressed at small values of  (Figure 3c).Given that PDI emerges as a wall mode, we have also confirmed its presence in cylindrical pipe flow as well as Taylor-Couette flow.
The lengthscale associated with PDI raises some intriguing questions.As indicated by Beneitez et al. (2023b), PDI can emerge at a lengthscale roughly on the order of the polymer gyration radius, where the continuum assumptions of the FENE-P model may not hold.It is thus possible that PDI is an unphysical feature of the widely-used FENE-P model, which would have significant implications for computational studies of ET, EIT and polymer drag reduction using this model.Until now, such studies may have been unknowingly influenced by PDI, due to the ubiquity of stress diffusion in numerical schemes, either introduced explicitly as a regularisation term or arising implicitly through the discretization scheme.Assessing the relevance of PDI to real viscoelastic fluids is now a key challenge to be confronted.

Figure 1 :
Figure1: Eigenvalue spectra, plotted in terms of the real   and imaginary   components of the complex wavespeed, obtained by solving the system (2.5a-g) using the Oldroyd-B constitutive model at  = 1000,  = 10 −5 and  = 0.9 for a) plane Couette flow ( = 44.7, = 37) and b) channel flow ( = 22.1,  = 52).These parameters are chosen to roughly lie on the corresponding neutral curves plotted in Figure2a.Results using bases of  = 300 and 400 Chebyshev polynomials are depicted by blue and red markers, respectively, to demonstrate convergence of the discrete unstable PDI modes highlighted in the zoomed inset above each panel.

Figure 2 :
Figure 2: Curves of neutral stability for plane Couette ('PCF', blue) and channel ('PPF', red) geometries, using the Oldroyd-B constitutive relation for five values of  ∈ [0.7, 0.98].a) The - plane for fixed  = 10 −5 , noting that the PCF and PPF curves are virtually indistinguishable.b) The - plane at three fixed  = {0, 1000, 5000}.In panels a-b, the  axis is scaled by the wall shear rate:  ′ wall = {1 (PCF), 2 (PPF)}.c) The streamwise wavenumber  of the PDI, as a function of , along each of the neutral curves in panel b.d) A collapse of the  = {1000, 5000} curves for both geometries from panel b) based on the inverse Schmidt number 1/ = .By plotting 1/ rather than , small  occurs at the bottom of panel d) facilitating an easier comparison with panel b).e) Colormap of the trace of the polymer conformation tensor tr(c) (red and blue denote positive and negative values, respectively), with contours of the streamfunction superimposed, for PPF eigenfunctions in the upper half channel with  = 0.9, at locations indicated by the square and triangular markers in panels b-c.One wavelength  = 2/ of each eigenfunction is shown.

Figure 3 :
Figure 3: Curves of neutral stability using the FENE-P constitutive relation with a fixed extensibility  = 200, presented analogously to the Oldroyd-B curves in Figure 2 for various .Curves are shown in a) the - plane with a fixed  = 10 −5 , b) the - plane for fixed  = 0, and c) the 1/ =  versus  plane, which collapses curves obtained at  = 1000 and 5000 for both geometries.By plotting 1/ rather than , small  occurs at the bottom of panel c) facilitating an easier comparison with panel b).Panel d) illustrates the dependence of the optimal streamwise wavenumber  on  for the PPF curves in panel c) at  = 1000.Left and right branches of the curves in panel c) are distinguished here using solid and dashed lines, respectively.Comparison with the Oldroyd-B curves in Figure 2c reveals that the left branches behave similarly to the single Oldroyd-B branch, deviating from the 1/ √  scaling at large , while the right branches retain this scaling to the highest  considered.

Figure 4
Figure 4: a) Dependence of the  = 0.87 neutral curve presented in Figure 3a on variable , with fixed  = 200.b) A collapse of the curves in panel a) based on the inverse Schmidt number 1/ = , plotted analogously to Figures 2d and 3c such that small  appears at the bottom of the plot.The plane Couette curves (PCF) exhibit a perfect collapse and are virtually indistinguishable, whereas the  = 10 −3 and 10 −4 channel curves intersect the  = 0 axis in panel a) and thus diverge to infinite  in panel b).c) Dependence of the  = 0.87 neutral curve presented in Figure 3a on variable , with fixed  = 10 −5 .d) Neutral curves in the - plane for fixed  = 10 −5 and  = 1000.