1. Introduction
The addition of polymers to a Newtonian solvent can induce dramatically different flow behaviours compared with those observed in the Newtonian fluid alone (Datta et al. Reference Datta2022; Sánchez et al. Reference Sánchez, Jovanović, Kumar, Morozov, Shankar, Subramanian and Wilson2022). In industrial processes, for instance, viscous polymer melts are susceptible to instabilities which constrain the maximum extrusion rate (Petrie & Denn Reference Petrie and Denn1976), while polymer additives are used in oil pipelines to reduce turbulent wall drag (Virk Reference Virk1975). Two particularly important viscoelastic phenomena are the existence of ‘elastic turbulence’ (ET), a chaotic flow state sustained in the absence of inertia (Groisman & Steinberg Reference Groisman and Steinberg2000; Steinberg Reference Steinberg2021), and ‘elasto-inertial turbulence’ (EIT), an inherently two-dimensional state arising when both inertia and elasticity are present (Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013; Sid, Terrapon & Dubief Reference Sid, Terrapon and Dubief2018; Choueiri et al. Reference Choueiri, Lopez, Varshney, Sankar and Hof2021). While the initial pathway to ET in curvilinear geometries is understood (Larson, Shaqfeh & Muller Reference Larson, Shaqfeh and Muller1990; Pakdel & McKinley Reference Pakdel and McKinley1996; Shaqfeh Reference Shaqfeh1996; Datta et al. Reference Datta2022), 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 two-torus, where Boffetta et al. (Reference Boffetta, Celani, Mazzino, Puliafito and Vergassola2005) found a linear instability leading to ET (Berti & Boffetta Reference Berti and Boffetta2010). Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) subsequently discovered a centre-mode instability in viscoelastic pipe flow at finite Reynolds number 
$Re$, which was later also identified in plane Poiseuille flow (PPF, hereafter referred to as channel flow) (Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021a) but notably not in plane Couette flow (PCF). Interestingly, this instability could only be traced down to 
$Re=0$ in channel flow (Khalid, Shankar & Subramanian Reference Khalid, Shankar and Subramanian2021b; Buza, Page & Kerswell Reference Buza, Page and Kerswell2022b). The finite-amplitude state resulting from this instability is an ‘arrowhead’ travelling wave (Page, Dubief & Kerswell Reference Page, Dubief and Kerswell2020; Buza et al. Reference Buza, Beneitez, Page and Kerswell2022a; Morozov Reference Morozov2022) which has been observed in channel flow EIT (Dubief et al. Reference Dubief, Page, Kerswell, Terrapon and Steinberg2022) and, in retrospect, ET in two-dimensional Kolmogorov flow (Berti & Boffetta Reference Berti and Boffetta2010). In channel flow, efforts have begun to establish a dynamical link between the arrowhead solution and both ET (Lellep, Linkmann & Morozov Reference Lellep, Linkmann and Morozov2023) and EIT (Beneitez et al. Reference Beneitez, Page, Dubief and Kerswell2023a). 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. Reference Beneitez, Page, Dubief and Kerswell2023a), which seems more dependent on a near-wall mechanism (Shekar et al. Reference Shekar, McMullen, Wang, McKeon and Graham2019, Reference Shekar, McMullen, McKeon and Graham2021; Dubief et al. Reference Dubief, Page, Kerswell, Terrapon and Steinberg2022).
The very recent discovery of a new wall-mode ‘polymer diffusive instability’ (PDI) in plane Couette flow at 
$Re=0$ (Beneitez, Page & Kerswell Reference Beneitez, Page and Kerswell2023b), 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 time-stepping algorithm, whether added explicitly to stabilize a numerical scheme like a spectral method (see e.g. Dubief, Terrapon & Hof Reference Dubief, Terrapon and Hof2023) or arising implicitly such as through a finite difference formulation (see e.g. Zhang et al. Reference Zhang, Jiang, Liang and Cheng2015; Pimenta & Alves Reference Pimenta and Alves2017). The PDI wall mode is primarily confined to a boundary layer of thickness 
$\sqrt {\varepsilon }$, where 
$\varepsilon$ is the (small) diffusion coefficient, travelling at the wall speed with a streamwise wavelength of 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 
$O(1)$ as 
$\varepsilon \rightarrow 0$. Direct numerical simulations (DNS) 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. Reference Beneitez, Page and Kerswell2023b).
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. Reference Beneitez, Page and Kerswell2023b), 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. (Reference Beneitez, Page and Kerswell2023b). Therefore, we here map out the regions where PDI is operative at finite 
$Re$, considering both plane Couette flow 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 
$Re$ and are relatively insensitive to the bulk flow geometry. Therefore, PDI is a candidate to trigger both ET and EIT in simulations using the FENE-P model.
2. Formulation
We consider the following dimensionless equations governing the flow of an incompressible viscoelastic fluid:
\begin{gather} Re\left(\frac{\partial\boldsymbol{u}}{\partial t}+\left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\boldsymbol{u}\right)={-}\boldsymbol{\nabla} p+\beta\nabla^{2}\boldsymbol{u}+\left(1-\beta\right)\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau}, \quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}
\begin{gather} \frac{\partial\boldsymbol{\boldsymbol{c}}}{\partial t}+\left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right) \boldsymbol{c}+\boldsymbol{\tau}=\boldsymbol{c}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u} +\left(\boldsymbol{\nabla}\boldsymbol{u}\right)^{T} \boldsymbol{\cdot}\boldsymbol{\boldsymbol{c}} +\varepsilon\nabla^{2}\boldsymbol{\boldsymbol{c}}, \end{gather}
where 
$\boldsymbol {u}=(u,v,w)$ and 
$p$ denote the velocity and pressure fields, respectively, and 
$\boldsymbol {\tau }$ denotes the polymeric contribution to the stress tensor. Following Beneitez et al. (Reference Beneitez, Page and Kerswell2023b), the equations have been non-dimensionalized using the channel half-width 
$H$ and a characteristic flow speed 
$U_{0}$, taken to be the wall speed in plane Couette flow or the centreline velocity of the base channel flow. In (2.1a), the Reynolds number 
$Re:=U_{0} H/\nu _T$ describes the ratio of inertial to viscous forces (with 
$\nu _T = \nu _S + \nu _P$ denoting the total kinematic viscosity comprised of solvent and polymer components), and 
$\beta := \nu _S/\nu _T$ denotes the ratio of solvent to total viscosity. The polymeric stress tensor 
$\boldsymbol {\tau }$ may be described in terms of the polymer orientation through the conformation tensor 
$\boldsymbol {c}$ as in (2.2). We emphasize that the inclusion of a polymer stress diffusion term 
$\varepsilon \nabla ^{2}\boldsymbol {\boldsymbol {c}}$ (associated with diffusivity 
$\varepsilon := (Re\,Sc)^{-1}$, where 
$Sc$ denotes the Schmidt number) is the crucial ingredient for the PDI identified by Beneitez et al. (Reference Beneitez, Page and Kerswell2023b). In the inertialess limit, 
$\varepsilon =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 
$\boldsymbol {\tau }$ and 
$\boldsymbol {c}$:
\begin{equation} \boldsymbol{\tau}:=\frac{f\left(\mathrm{tr}\,\boldsymbol{c}\right)\boldsymbol{c}-\boldsymbol{\mathsf{I}}}{W},\quad f\left(s\right):=\left(1-\frac{s-3}{L^{2}}\right)^{{-}1}, \end{equation}
where 
$\boldsymbol{\mathsf{I}}$ is the identity matrix, 
$L$ denotes the maximum extensibility of the polymer chains and 
$W:=U_{0} \lambda /H$, the Weissenberg number, describes the ratio of time scales for polymer relaxation (
$\lambda$) to the flow. In the limit 
$L\rightarrow \infty$, the simpler Oldroyd-B model is obtained. Inspection of (2.1)–(2.3) reveals five parameters of interest governing the flow dynamics: 
$Re$, 
$W$, 
$\beta$, 
$\varepsilon$, 
$L$.
We analyse the linear stability of (2.1)–(2.3) by perturbing them about their base state: 
$\boldsymbol {u}=\boldsymbol {U}+\boldsymbol {u}^{*}$, 
$p=P+p^{*}$, 
$\boldsymbol {\tau }=\boldsymbol {T}+\boldsymbol {\tau }^{*}$, 
$\boldsymbol {c}=\boldsymbol {C}+\boldsymbol {c}^*$. Our coordinate system (
$x$, 
$y$, 
$z$) is aligned with the streamwise, wall-normal and spanwise directions of the channel, respectively. The base state 
$(U(\kern0.7pt y),C_{xx}(\kern0.7pt y),C_{yy}(\kern0.7pt y),C_{zz}(\kern0.7pt y),C_{xy}(\kern0.7pt y))$ satisfies
$$\begin{gather} -\partial_{x}P+\beta U''+\left(1-\beta\right)\partial_{y}T_{xy}=0, \end{gather}$$
$$\begin{gather}f\left(\mathrm{tr}\boldsymbol{C}\right)\boldsymbol{C}-\varepsilon W\boldsymbol{C}''-\begin{pmatrix}2WU'C_{xy}+1 & WU'C_{yy} & 0\\ WU'C_{yy} & 1 & 0\\ 0 & 0 & 1 \end{pmatrix}=0, \end{gather}$$
where primes indicate derivatives in the wall-normal (
$y$) direction. We use pressure gradients 
$\partial _{x} P = 0$ and 
$\partial _{x} P = -2$ for plane Couette and channel flow, respectively. While we here compute the base flow accounting for the presence of a finite diffusivity 
$\varepsilon$ 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 
$\varepsilon =0$ in (2.4b) as the inclusion of 
$\varepsilon \neq 0$ does not significantly change the base state except in the limit of large 
$\varepsilon = O(1)$.
Normal mode solutions of the perturbed flow are sought using the ansatz 
$\phi ^{*}(x,y,t)=\tilde {\phi }(\kern0.7pt y){\rm e}^{{\rm i}k(x-ct)}$, where real-valued 
$k$ denotes the streamwise wavenumber and 
$c=c_{r}+{\rm i}c_{i}$ is a complex wave speed, with instability arising if 
$c_i> 0$. The perturbed state is governed by the following system of seven equations for 
$(\tilde {u},\tilde {v},\tilde {p},\tilde {c}_{xx},\tilde {c}_{yy},\tilde {c}_{zz},\tilde {c}_{xy})$:
$$\begin{gather} {\rm i}k\tilde{u}+\tilde{v}' = 0, \end{gather}$$
$$\begin{gather}Re\left(-{\rm i}kc\tilde{u}+\tilde{v}U'+{\rm i}kU\tilde{u}\right)+{\rm i}k\tilde{p}-\beta\left({-}k^{2}\tilde{u}+\tilde{u}''\right)-\left(1-\beta\right)\left({\rm i}k\tilde{\tau}_{xx}+\tilde{\tau}'_{xy}\right)=0, \end{gather}$$
$$\begin{gather}Re\left(-{\rm i}kc\tilde{v}+{\rm i}kU\tilde{v}\right)+ \tilde{p}'-\beta\left({-}k^{2}\tilde{v}+\tilde{v}''\right)-\left(1-\beta\right)\left({\rm i}k\tilde{\tau}_{xy}+\tilde{\tau}'_{yy}\right)=0, \end{gather}$$
$$\begin{gather}\left[\varepsilon k^{2}+{\rm i}k\left(U-c\right)\right]\tilde{c}_{xx}+\tilde{v}C'_{xx}+\tilde{\tau}_{xx}-\varepsilon\tilde{c}''_{xx}-2\left({\rm i}kC_{xx}\tilde{u}+C_{xy}\tilde{u}'+\tilde{c}_{xy}U'\right)=0, \end{gather}$$
$$\begin{gather}\left[\varepsilon k^{2}+{\rm i}k\left(U-c\right)\right]\tilde{c}_{yy}+\tilde{v}C'_{yy}+\tilde{\tau}_{yy}-\varepsilon\tilde{c}''_{yy} -2\left({\rm i}kC_{xy}\tilde{v}+C_{yy}\tilde{v}'\right)=0, \end{gather}$$
$$\begin{gather}\left[\varepsilon k^{2}+{\rm i}k\left(U-c\right)\right]\tilde{c}_{zz}+\tilde{v}C'_{zz}+\tilde{\tau}_{zz}-\varepsilon\tilde{c}''_{zz} =0, \end{gather}$$
$$\begin{gather}\left[\varepsilon k^{2}+{\rm i}k\left(U-c\right)\right]\tilde{c}_{xy}+\tilde{v}C'_{xy}+\tilde{\tau}_{xy}-\varepsilon\tilde{c}''_{xy} -{\rm i}kC_{xx}\tilde{v}-U'\tilde{c}_{yy}-C_{yy}\tilde{u}'=0. \end{gather}$$ 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. Reference Beneitez, Page and Kerswell2023b). A basis of 
$N=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 
$N=300$ and 
$N=400$ Chebyshev polynomials to demonstrate convergence. The unstable PDI mode emerges with characteristic wave speeds in the vicinity of 
$|c_{r}|\approx 1$ and 
$c_r \approx 0$ for PCF and PPF geometries, respectively.
Figure 1. Eigenvalue spectra, plotted in terms of the real 
$c_r$ and imaginary 
$c_i$ components of the complex wave speed, obtained by solving the system (2.5a)–(2.5g) using the Oldroyd-B constitutive model at 
$Re=1000$, 
$\varepsilon =10^{-5}$ and 
$\beta =0.9$ for (a) plane Couette flow (
$W=44.7$, 
$k=37$) and (b) channel flow (
$W=22.1$, 
$k=52$). These parameters are chosen to roughly lie on the corresponding neutral curves plotted in figure 2(a). Results using bases of 
$N=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.
Following Beneitez et al. (Reference Beneitez, Page and Kerswell2023b), we choose boundary conditions that set 
$\varepsilon =0$ at the walls in the governing equations to match the typical configuration used in DNS. While 
$Sc\sim O(10^{6})$ is characteristic of polymer diffusion in a typical solvent like water (yielding 
$\varepsilon = 1/(Re\,Sc) \approx 10^{-9}-10^{-6}$ depending on 
$Re$), in DNS much larger values of 
$\varepsilon \approx O(10^{-3})$ are often employed to achieve numerical stability (see e.g. § 2.2 of Dubief et al. Reference Dubief, Terrapon and Hof2023). Thus, it is commonplace to set 
$\varepsilon =0$ at the walls in order to recover the idealized limit 
$\varepsilon \to 0$ (Sureshkumar, Beris & Handler Reference Sureshkumar, Beris and Handler1997; Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013; Dubief et al. Reference Dubief, Page, Kerswell, Terrapon and Steinberg2022).
3. 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 
$\sigma := kc_i=0$) associated with PDI in the four-dimensional parameter space spanned by 
$(Re, W, \beta,\varepsilon )$. For a given point in parameter space, we perform a sweep over wavenumbers 
$k$ 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 (
$W$) axis is scaled by the respective wall shear rates: 
$U_{{wall}}'=\{ 1\ (\textrm {PCF}),\ 2\ (\textrm {PPF})\}$. As described in § 1, the PDI is a wall mode confined to a boundary layer of thickness 
$O(\sqrt {\varepsilon })$ 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 
$\varepsilon$), 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 
$Re$, we also demonstrate that curves associated with variable 
$Re$ and 
$\varepsilon$ may be collapsed based on 
$Sc = 1/(\varepsilon Re)$.
3.1. The Oldroyd-B case
The trajectory of neutral curves in the 
$Re$–
$W$–
$\varepsilon$ volume are presented in figure 2 for various 
$\beta \in [0.7,0.98]$. Figure 2(a) considers the 
$Re$–
$W$ plane at a fixed 
$\varepsilon = 10^{-5}$, while figure 2(b) considers the 
$\varepsilon$–
$W$ plane at three fixed 
$Re\in \{ 0,1000,5000\}$. Regions of stability and instability are found to the left and right of each curve, respectively. For a given 
$\beta$, increasing 
$Re>0$ is found to promote instability at progressively lower values of both 
$W$ (figure 2a) and 
$\varepsilon$ (figure 2b). Notably, in the inertialess limit (
$Re=0$), Beneitez et al. (Reference Beneitez, Page and Kerswell2023b) found that the neutral curves tracked a roughly fixed 
$W \sim \beta / (1-\beta )$ for 
$\varepsilon \lesssim 10^{-2}$ (see their figure 2b) leading to the disappearance of PDI at finite 
$W$ in the limit 
$\beta \rightarrow 1$. Here, for 
$Re>0$, our figure 2(b) demonstrates that while for small 
$\varepsilon =O(10^{-7})$ the neutral curves do indeed still track a constant 
$W$, they then begin to significantly deviate to lower 
$W$ beyond 
$\varepsilon \gtrsim 10^{-6}$, with the 
$Re=5000$ case deviating at lower 
$\varepsilon$ than the 
$Re=1000$ case. Inertial effects thus play a significant role in promoting a greater prevalence of PDI in the parameter space. In figure 2(d), we demonstrate that the 
$Re=\{ 1000,5000\}$ curves plotted in figure 2(b) may be collapsed for both geometries based on the Schmidt number 
$Sc = 1/(\varepsilon Re)$.
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 
$\beta \in [0.7,0.98]$. (a) The 
$Re$–
$W$ plane for fixed 
$\varepsilon = 10^{-5}$, noting that the PCF and PPF curves are virtually indistinguishable. (b) The 
$\varepsilon$–
$W$ plane at three fixed 
$Re=\{ 0,1000,5000\}$. In panels (a,b), the 
$W$ axis is scaled by the wall shear rate: 
$U_{{wall}}'=\{ 1\ (\textrm {PCF}),\ 2\ (\textrm {PPF})\}$. (c) The streamwise wavenumber 
$k$ of the PDI, as a function of 
$\varepsilon$, along each of the neutral curves in panel (b). (d) A collapse of the 
$Re=\{ 1000,5000\}$ curves for both geometries from panel (b) based on the inverse Schmidt number 
$1/Sc = \varepsilon Re$. By plotting 
$1/Sc$ rather than 
$Sc$, small 
$\varepsilon$ occurs at the bottom of panel (d) facilitating an easier comparison with panel (b). (e,f) Colourmaps of the trace of the polymer conformation tensor 
$\mathrm {tr}(\boldsymbol {c})$ (red and blue denote positive and negative values, respectively), with contours of the stream function superimposed, for PPF eigenfunctions in the upper half-channel with 
$\beta =0.9$, at locations indicated by the square and triangular markers in panels (b,c). One wavelength 
$\lambda = 2{\rm \pi} / k$ of each eigenfunction is shown.
The streamwise wavenumbers 
$k$ associated with the neutral curves in figure 2(b) are plotted in figure 2(c) as a function of 
$\varepsilon$. At 
$Re=0$, 
$k$ follows the 
$1/\sqrt {\varepsilon }$ scaling reported by Beneitez et al. (Reference Beneitez, Page and Kerswell2023b) for all 
$\varepsilon$, whereas for 
$Re>0$, 
$k$ deviates significantly from this scaling to plateau to a roughly constant value for 
$\varepsilon \gtrsim 10^{-5}$, with the unstable modes at higher 
$Re$ being more tightly confined to the wall (as indicated by a larger 
$k$). This deviation from the 
$1/\sqrt {\varepsilon }$ scaling for 
$Re>0$ corresponds to the previously noted deviation of the neutral curves away from a constant 
$W$ in figure 2(b), where the curves begin turning to lower 
$W$ in the vicinity of 
$\varepsilon \approx 10^{-6}-10^{-5}$. The behaviour of 
$k$ in figure 2(c) also explains the mismatch in the collapse of the PCF and PPF curves in figure 2(b) at higher 
$\varepsilon$ for 
$Re=0$ (see solid lines). Specifically, at higher 
$\varepsilon$, 
$k$ becomes sufficiently small at 
$Re=0$ such that the instability is no longer strongly confined to the wall (see square eigenfunction, figure 2e) and so the local wall shear, 
$U'_{wall}$, does not accurately describe the non-uniform shear profile influencing the instability. Conversely, the curves for 
$Re=\{ 1000,5000\}$ do remain collapsed at high 
$\varepsilon$, as the wavenumbers 
$k$ in figure 2(c) plateau to sufficiently large values such that the instability remains confined to the wall (see triangle eigenfunction, figure 2f). In figure 2(c), 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.
3.2. The FENE-P case
We now consider how a finite polymer extensibility 
$L$ modifies the behaviour of the instability as compared with the Oldroyd-B case (
$L\rightarrow \infty$) presented in § 3.1. Focusing first on 
$L=200$, figure 3 illustrates the much richer behaviour of the neutral curves for the FENE-P case, presented in an analogous manner to figures 2(a–d). In the 
$Re$–
$W$ plane (figure 3a), a finite 
$L$ introduces two notable differences compared with Oldroyd-B. First, the neutral curves for a given 
$\beta$ now have a left-hand and right-hand branch and so the range of instability is bounded by an upper value of 
$W$. Second, there is now a critical value of 
$\beta$ (
$\approx 0.865$ for both geometries) above which the neutral curves no longer intersect the 
$Re=0$ axis, at fixed 
$\varepsilon =10^{-5}$. Therefore, inertial effects are once again demonstrated to promote PDI, generating instability at finite 
$Re$ for ultradilute polymer solutions (
$\beta \rightarrow 1$) that would otherwise remain stable at 
$Re=0$.
Figure 3. Curves of neutral stability using the FENE-P constitutive relation with a fixed extensibility 
$L=200$, presented analogously to the Oldroyd-B curves in figure 2 for various 
$\beta$. Curves are shown in (a) the 
$Re$–
$W$ plane with a fixed 
$\varepsilon =10^{-5}$, (b) the 
$\varepsilon$–
$W$ plane for fixed 
$Re=0$ and (c) the 
$1/Sc = \varepsilon Re$ vs 
$W$ plane, which collapses curves obtained at 
$Re=1000$ and 
$5000$ for both geometries. By plotting 
$1/Sc$ rather than 
$Sc$, small 
$\varepsilon$ occurs at the bottom of panel (c), facilitating an easier comparison with panel (b). Panel (d) illustrates the dependence of the optimal streamwise wavenumber 
$k$ on 
$\varepsilon$ for the PPF curves in panel (c) at 
$Re=1000$. Left-hand and right-hand branches of the curves in panel (c) are distinguished in panel (d) using solid and dashed lines, respectively. Comparison with the Oldroyd-B curves in figure 2(c) reveals that the left-hand branches behave similarly to the single Oldroyd-B branch, deviating from the 
$1/\sqrt {\varepsilon }$ scaling at large 
$\varepsilon$, while the right-hand branches retain this scaling to the highest 
$\varepsilon$ considered.
Neutral curves in the 
$\varepsilon$–
$W$ plane for the inertialess case 
$Re=0$ are shown in figure 3(b), exhibiting notable differences to the Oldroyd-B 
$Re=0$ curves presented in figure 2(b). 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 
$\varepsilon$. Specifically, the plane Couette curves have an inverted ‘U’-shape, highlighting that the instability ceases to exist at large 
$\varepsilon$ for certain 
$\beta$ (e.g. see the 
$\beta =0.86$ curve which reaches a maximum value at 
$\varepsilon \approx 6\times 10^{-3}$). Conversely, the channel curves are roughly ‘U’-shaped, and the range of instability only increases with increasing 
$\varepsilon$.
At finite 
$Re$, as for the Oldroyd-B curves in figure 2(d), the FENE-P curves shown in figure 3(c) at 
$Re=\{ 1000,5000\}$ are again found to collapse for both geometries based on the Schmidt number 
$Sc=1/(\varepsilon Re)$. Comparing figures 2(d) and 3(c) reveals two main differences in the structure of the FENE-P and Oldroyd-B curves. First, for 
$\beta \approx 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 
$\varepsilon -W$ plane in an otherwise stable region. Second, in the limit of 
$\varepsilon \rightarrow 0$ (
$Sc \rightarrow \infty$), the neutral curves behave differently as compared with the 
$Re=0$ case reported by Beneitez et al. (Reference Beneitez, Page and Kerswell2023b). In particular, there now appears to be a critical value of 
$\beta \approx 0.865$ (note the similarity to the critical 
$\beta$ in figure 3a corresponding to lift-off from the 
$Re=0$ axis) at which the two branches of the neutral curve asymptotically approach each other in the limit of small 
$\varepsilon$. Curves associated with 
$\beta$ greater than this critical value are thus observed to turn back at higher 
$\varepsilon$ (see e.g. the 
$\beta =0.9$ curves, figure 3c), suggesting that, at finite 
$Re$, PDI will not exist in the limit 
$\varepsilon \rightarrow 0$ for all 
$\beta$. We note another possibility, however, which is that an ‘hourglass’-like pinch-off behaviour occurs, in which the critical curves touch at some finite 
$\varepsilon$ (appearing here to be around 
$Sc \approx 10^{4}$), but then separate again for lower 
$\varepsilon$. The concave-up neutral curves (e.g. 
$\beta =0.9$) seen here, might then be reflected as concave-down branches at much lower 
$\varepsilon$ and instabilities for such 
$\beta$ could then still exist in the 
$\varepsilon \rightarrow 0$ limit. As 
$\varepsilon \rightarrow 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 
$Sc > 10^{4}$ and so we do not further consider this behaviour here. It is also worth noting in figure 3(c) that while increasing 
$Re$ does not significantly influence the horizontal position of the neutral curves along the 
$W$ axis, the increase in inertial effects does induce a significant downward shift of the curves in 
$\varepsilon$ (by a factor proportional to 
$Re$), thus increasing the prevalence of PDI in parameter space.
Figure 4. (a) Dependence of the 
$\beta =0.87$ neutral curve presented in figure 3(a) on variable 
$\varepsilon$, with fixed 
$L=200$. (b) A collapse of the curves in panel (a) based on the inverse Schmidt number 
$1/Sc = \varepsilon Re$, plotted analogously to figures 2(d) and 3(c) such that small 
$\varepsilon$ appears at the bottom of the plot. The PCF curves exhibit a perfect collapse and are virtually indistinguishable, whereas the 
$\varepsilon = 10^{-3}$ and 
$10^{-4}$ channel curves intersect the 
$Re=0$ axis in panel (a) and thus diverge to infinite 
$Sc$ in panel (b). (c) Dependence of the 
$\beta =0.87$ neutral curve presented in figure 3(a) on variable 
$L$, with fixed 
$\varepsilon = 10^{-5}$. (d) Neutral curves in the 
$\beta$–
$W$ plane for fixed 
$\varepsilon =10^{-5}$ and 
$Re=1000$.
The streamwise wavenumbers 
$k$ associated with the FENE-P channel (PPF) neutral curves at 
$Re=1000$ in figure 3(c) are presented in figure 3(d) as a function of 
$\varepsilon$, to compare with the scaling of the Oldroyd-B curves in figure 2(c). As for Oldroyd-B, the left-hand branches follow the 
$1/\sqrt {\varepsilon }$ scaling reported by Beneitez et al. (Reference Beneitez, Page and Kerswell2023b) for 
$\varepsilon \lesssim 10^{-5}$, when the neutral curves are roughly independent of 
$W$ in figure 3(c). In contrast, the right-hand branches follow this 
$k$ scaling for the entire range of 
$\varepsilon$ considered here, thus explaining the slight mismatch in the collapse of the right-hand branches of the two geometries at low 
$\varepsilon$ in figure 3(c) (see e.g. the right-hand branches of the 
$\beta =0.7$ and 0.8 curves); in this regime, 
$k$ has become sufficiently small such that the instability is no longer confined to the wall and thus the wall shear 
$U_{{wall}}'$ is not entirely suitable for scaling the 
$W$ axis. In figure 3(d), we also note that for 
$\beta \gtrsim 0.864$, the pinchoff value seen in figure 3(c) at 
$Sc = O(10^{4})$, the neutral curves turn back to higher 
$\varepsilon$ before the 
$1/\sqrt {\varepsilon }$ 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 
$Re$–
$W$ plane was considered with a fixed 
$\varepsilon = 10^{-5}$ in figure 3(a), in figure 4(a) we focus on the 
$\beta =0.87$ neutral curve and demonstrate its dependence on variable 
$\varepsilon$, where it is found that decreasing 
$\varepsilon$ pushes the region of PDI to progressively higher 
$Re$. In figure 4(b), we collapse the curves from figure 4(a) based on the Schmidt number 
$Sc = 1/(\varepsilon Re)$. A near perfect collapse is observed for the PCF curves, while some channel (PPF) curves diverge to 
$Sc \rightarrow \infty$ due to their intersection with the 
$Re=0$ axis in figure 4(a). Additionally, having only considered a fixed polymer extensibility 
$L=200$ in figure 3, in figures 4(c,d) we consider the influence of variable 
$L$ on the neutral curves in the 
$Re$–
$W$ and 
$\beta$–
$W$ planes, respectively. In figure 4(c), decreasing 
$L$ induces a turn-back of the neutral curves at progressively higher values of 
$Re$, decreasing the extent of PDI in the 
$Re$–
$W$ plane. Figure 4(d), at fixed 
$Re=1000$, displays the same qualitative features as those reported in the inertialess limit by Beneitez et al. (Reference Beneitez, Page and Kerswell2023b) (see their figure 2c), again demonstrating that decreasing 
$L$ decreases the prevalence of PDI in the 
$\beta$–
$W$ plane.
3.3. 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 
$Re$ in regions of instability. Starting with each 
$\beta$ curve in figure 2(a) (Oldroyd-B), we first fix the value of 
$W$ at which the neutral curve intersects the 
$Re=0$ axis, and where the growth rate is thus zero by definition. At this fixed 
$W$, we then increase 
$Re$ incrementally and track the growth rate of the most unstable PDI mode, as shown in figure 5(a). In figure 5(b), we repeat this process for the FENE-P curves presented in figure 3(a). As some FENE-P curves do not intersect the 
$Re=0$ axis, in these cases we fix 
$W$ to correspond to the minimum 
$Re$ of the neutral curve, and then increase 
$Re$ from there. In figures 5(c,d), we again repeat this process for the neutral curves plotted in figures 4(b,c), respectively.
Figure 5. Growth rates 
$\sigma := kc_{i}$ of the most unstable mode as a function of Reynolds number 
$Re$ obtained at a fixed 
$W$ corresponding to either the intersection of the neutral curve with the 
$Re=0$ axis or the minimum of the neutral curve in 
$Re$ for cases in which there is no intersection with 
$Re=0$. Growth rates are shown for neutral curves spanning various 
$\beta$ in (a) figure 2(a) (Oldroyd-B, fixed 
$\varepsilon =10^{-5}$) and (b) figure 3(a) (fixed 
$L=200$, 
$\varepsilon =10^{-5}$), and neutral curves at fixed 
$\beta =0.87$ spanning (c) various 
$\varepsilon$ in figure 4(b) (fixed 
$L=200$) and (d) various 
$L$ in figure 4(c) (fixed 
$\varepsilon =10^{-5}$).
In all cases, for neutral curves that intersect the 
$Re=0$ axis, the growth rate is observed to grow linearly with 
$Re$ as one moves away from the neutral curve, emphasizing the intensification of PDI due to the presence of inertia. Notably, the streamwise wavenumber 
$k$ remains virtually constant during this scaling, until 
$Re \gtrsim 10^3$ at which point the most unstable 
$k$ beings to vary and the linear scaling breaks down. For the subset of FENE-P curves that do not intersect the 
$Re=0$ axis (e.g. 
$\beta =0.87$, 0.88, 0.90 in figure 4b), it is also intriguing that the growth rates display a dramatic increase with 
$Re$ to quickly join the linear 
$Re$ scaling of the curves that do intersect the 
$Re=0$ axis. We note that the relative vertical translation of the various curves in figure 5 is due to differences in slope between the neutral curves at their intersection with the 
$Re=0$ axis; at a fixed 
$W$, these differences in slope will govern the rate at which one moves away from the neutral curve as 
$Re$ is increased, and hence the observed difference in the growth rate magnitudes.
4. Conclusions
In this study, we have demonstrated that the PDI is active in both plane Couette and channel flows with or without inertial effects, and that the instability intensifies with increasing Reynolds number 
$Re$. Through exploration of a variety of dimensionless parameters, we have found that PDI is operational across large regions of the parameter space including those relevant to many prior experiments (Choueiri, Lopez & Hof Reference Choueiri, Lopez and Hof2018; Qin et al. Reference Qin, Salipante, Hudson and Arratia2019; Choueiri et al. Reference Choueiri, Lopez, Varshney, Sankar and Hof2021; Jha & Steinberg Reference Jha and Steinberg2021). In particular, increasing 
$Re$ enhances the prevalence of the instability, promoting instability at progressively smaller values of both 
$W$ and 
$\varepsilon$ than in the inertialess limit. Our results therefore significantly extend the conclusion of Beneitez et al. (Reference Beneitez, Page and Kerswell2023b) that PDI could also present a possible transition mechanism to EIT as well as ET in FENE-P fluids.
The eigenfunction for PDI is a wall mode, confined to a boundary layer of thickness 
$\sqrt {\varepsilon }$. As a result, the neutral curves for plane Couette and channel flow are found to nearly overlap in most regions of parameter space when 
$W$ is scaled by the wall-shear rate. This collapse breaks down when the streamwise wavenumber 
$k$ 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 
$\varepsilon$ at 
$Re=0$ (but notably not at higher 
$Re$, see figure 2c), and small 
$\varepsilon$ at high 
$W$ in FENE-P fluids. The finite extensibility of the polymer chains (
$L$) is also found to have a significant impact on the prevalence of PDI, as compared with that predicted by the Oldroyd-B model. Considering 
$L=200$, we found that for sufficiently high 
$\beta$, the instability is suppressed at 
$Re=0$ and only appears at progressively larger 
$Re$ (figure 3a). Similarly, beyond a critical value of 
$\beta$, the instability may also be suppressed at small values of 
$\varepsilon$ (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 length scale associated with PDI raises some intriguing questions. As indicated by Beneitez et al. (Reference Beneitez, Page and Kerswell2023b), PDI can emerge at a length scale roughly of 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 regularization 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.
Funding
We are grateful for the support of EPSRC under grant EP/V027247/1.
Declaration of interests
The authors report no conflict of interest.
Open access
