eurm10 msam10

# Asymptotics of the centre mode instability in viscoelastic channel flow:
with and without inertia

Rich R. Kerswell1and Jacob Page2
r.r.kerswell<EMAIL_ADDRESS>1DAMTP, Centre for
Mathematical Sciences, University of Cambridge, CB3 0WA, UK.
2School of Mathematics, University of Edinburgh, EH9 3FD, UK

(?; revised ?; accepted ?.)

###### Abstract

Motivated by the recent numerical results of Khalid et al., Phys. Rev. Lett.,
127 134502 (2021), we consider the large-Weissenberg-number ($W$) asymptotics
of the centre mode instability in inertialess viscoelastic channel flow. The
instability is of the critical layer type in the distinguished ultra-dilute
limit where $W(1-\beta)=O(1)$ as $W\rightarrow\infty$ ($\beta$ is the ratio of
solvent-to-total viscosity). In contrast to centre modes in the Orr-Sommerfeld
equation, $1-c=O(1)$ as $W\rightarrow\infty$ where $c$ is the phase speed
normalised by the centreline speed as a central ‘outer’ region is always
needed to adjust the non-zero cross-stream velocity at the critical layer down
to zero at the centreline. The critical layer acts as a pair of intense
‘bellows’ which blows the flow streamlines apart locally and then sucks them
back together again. This compression/rarefaction amplifies the streamwise-
normal polymer stress which in turn drives the streamwise flow through local
polymer stresses at the critical layer. The streamwise flow energises the
cross-stream flow via continuity which in turn intensifies the critical layer
to close the cycle. We also treat the large-Reynolds-number ($Re$) asymptotic
structure of the upper (where $1-c=O(Re^{-2/3})$) and lower branches of the
$Re$-$W$ neutral curve confirming the inferred scalings from previous
numerical computations. Finally, we argue that the viscoelastic centre mode
instability was actually first found in viscoelastic Kolmogorov flow by
Boffetta et al., J. Fluid Mech. 523, 161-170 (2005).

###### keywords:

##  1 Introduction

It is well known that the addition of long-chain polymers to a Newtonian fluid
introduces elasticity which can give rise to fascinating new ‘viscoelastic’
flow phenomena. Prime examples of this are a new form of spatiotemporal chaos
- dubbed ‘elastic’ turbulence (ET) (Groisman & Steinberg, 2000) \- which
exists in inertialess curvilinear flows and ‘elasto-inertial’ turbulence (EIT)
(Samanta et al., 2013) which can occur in 2D rectilinear flows (Sid et al.,
2018) where inertia and elasticity balance each other. While ET is assumed
triggered by a linear ‘hoop stress’ instability of curved streamlines Larson
et al. (1990); Shaqfeh (1996), the origin of EIT remains unclear (Datta et
al., 2022; Dubief et al., 2023) as does any possible relationship to ET.

The breakdown of viscoelastically-modified Tollmien-Schlichting modes has been
suggested as a cause of EIT (Shekar et al., 2019, 2021) at least at high
Reynolds number, $Re$ and low Weissenberg number, $W$. At low $Re$ and high
$W$, however, the recent discovery of a new linear instability of rectilinear
viscoelastic shear flow seems more viable (Garg et al., 2018; Chaudhary et
al., 2021; Khalid et al., 2021a). This instability occurs at higher $W$ than
generally associated with EIT but has been shown to be subcritical (Page et
al., 2020; Wan et al., 2021; Buza et al., 2022b, a). In particular, travelling
wave solutions, which have a distinctive ‘arrowhead’ structure, originating
from the neutral curve reach down in $W$ to where EIT exists in parameter
space (Page et al., 2020; Buza et al., 2022a; Dubief et al., 2022). This
instability is of centre-mode type being localised either at the centre of a
pipe (Garg et al., 2018; Chaudhary et al., 2021) or midplane of a channel
Khalid et al. (2021a, b) but is notably absent in plane Couette flow (Garg et
al., 2018). Perhaps most intriguingly, the instability can be traced down to
$Re=0$ in channel flow (Khalid et al., 2021b) in the ultra-dilute limit of the
solvent-to-total viscosity ratio approaching 1 while a minimum $Re\approx 63$
exists in pipe flow (Chaudhary et al., 2021). Subsequently, travelling wave
solutions have been numerically computed in 2 dimensions and at $Re=0$ (Buza
et al., 2022a; Morozov, 2022) and their instability examined (Lellep et al.,
2023a, b).

Apart from numerically-inferred scaling relationships (Garg et al., 2018;
Chaudhary et al., 2021; Khalid et al., 2021a, b), the only work to unpick the
asymptotic structure of the centre-mode instability is that of Dong & Zhang
(2022) in the pipe flow. They identify the asymptotic structure on the upper
branch of the neutral curve characterised by $W\sim Re^{1/3}$ as
$Re\rightarrow\infty$ and consider the long wavelength limit but stop short of
treating the lower branch of the neutral curve. Here we do both for the
channel and go further to examine the inertialess regime in channel flow which
is absent in pipe flow. Unravelling the $Re=0$ situation asymptotically is
actually our main motivation here as it differs fundamentally from all the
classical Orr-Sommerfeld work performed for Newtonian shear flows (Drazin &
Reid, 1981). In particular, the regularizing feature of the critical layer
formed (e.g. figure 3 of Khalid et al. (2021b) and figure 4 below) is the
presence of elastic relaxation rather than viscosity and the ‘outer’
relaxation-free solutions satisfy a fourth order differential equation rather
than the classical, inviscid, second order Rayleigh equation in Newtonian
flows. This means that matching conditions across the critical layer need to
be sought down to the third order derivative in the cross-stream velocity (or
streamfunction) and, due to a logarithmic singularity in the first order
derivative, computations need to go beyond double precision accuracy to
achieve a convincing correspondence between numerical results and the
asymptotic predictions; see Table 4. A particularly interesting feature of
this viscoelastic centre mode instability is the critical layer does not
approach the midplane as $W\rightarrow\infty$, that is, the phase speed of the
instability approaches a non-trivial value very close to but distinct from 1,
the maximum speed of the base flow. Ultimately, though, the point of the
asymptotic analysis is to identify the mechanism of the instability and to
understand, if possible, why it does not manifest in plane Couette flow.

The plan of this paper is first to introduce the channel flow problem in §2
and the viscoelastic model (Oldroyd-B) used by Khalid et al. (2021a, b). The
first results section, §3, then examines the large $Re$-asymptotics of the
upper (§3.1) and lower branches (§3.2) of the neutral curve in the $Re$-$W$
plane for fixed $\beta$: see figure 1. Reduced eigenvalue problems based only
on O(1) quantities (relative to $Re$) can be straightforwardly derived for
both upper and lower branches. Interestingly, if $\beta\gtrsim 0.9905$, Khalid
et al. (2021b) showed that the lower branch crosses the $Re=0$ axis and the
appropriate (mathematical) limit is then $Re\rightarrow-\infty$. Figure 3
indicates that nothing mathematically untoward happens as the neutral curve
swings around from pointing at $Re\rightarrow\infty$ to $Re\rightarrow-\infty$
although, of course, negative $Re$ makes little physical sense. The special
case of $Re=0$ or vanishing inertia, however, does and the asymptotics as
$W\rightarrow\infty$ is studied in §4. The work of Khalid et al. (2021b) has
already indicated that the appropriate distinquished limit is that in which
$\beta$ simultaneously approaches 1 such that $W(1-\beta)$ stays finite.
Section 5 goes on to use the asymptotic solution to discuss the mechanics of
the inertialess instability and §6 describes some numerical experiments to
understand how the instability responds to the problem becoming a bit more
plane-Couette like. A brief §7 presents evidence that the centre mode
instability was actually found first in viscoelastic Kolmogorov flow (Boffetta
et al., 2005) before a final discussion follows in §8.

##  2 Formulation

We consider pressure-driven, incompressible channel flow between two walls
$y=\pm h$ in the $x$-direction. Using the half-channel height, $h$, and the
base centreline speed $U_{max}$ to non-dimensionalize the problem, the
governing equations become

$Re\left(\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\mbox{\boldmath$\nabla$}\mathbf{u}\right)=-\nabla{\cal
P}+\beta\nabla^{2}\mathbf{u}+\mbox{\boldmath$\nabla$}\cdot\mbox{\boldmath{${\cal
T}$}},$ (1) $\nabla\cdot\mathbf{u}=0,$ (2)
$\frac{\partial\mbox{\boldmath{${\cal T}$}}}{\partial
t}+\mathbf{u}\cdot\mbox{\boldmath$\nabla$}\mbox{\boldmath{${\cal T}$}}-2\,{\rm
sym}(\mbox{\boldmath{${\cal
T}$}}.\mbox{\boldmath$\nabla$}\mathbf{u})=-\frac{1}{W}\mbox{\boldmath{${\cal
T}$}}+\frac{1-\beta}{W}(\,\mbox{\boldmath$\nabla$}\mathbf{u}+\mbox{\boldmath$\nabla$}\mathbf{u}^{T}\,)$
(3)

where $\mathbf{u}$ is the velocity field, ${\cal P}$ the pressure and ${\cal
T}$ the polymer stress following Khalid et al. (2021a). Here an Oldroyd-B
fluid has been assumed so

$\mbox{\boldmath{${\cal T}$}}=\frac{1-\beta}{W}({\bf C}-{\bf I})$ (4)

where ${\bf C}$ is the conformation tensor. The parameters of the problem are
the Reynolds number, Weissenberg number and the solvent-to-total viscosity
ratio,

$Re:=\frac{U_{max}h}{\nu},\quad W:=\frac{\lambda
U_{max}}{h}\quad\&\quad\beta:=\frac{\nu_{s}}{\nu}$ (5)

respectively, where $\lambda$ is the microstructural relaxation time,
$\nu_{s}$ is the solvent kinematic viscosity and $\nu$ is the total kinematic
viscosity (following Khalid et al. (2021a, b)). The scaling of the pressure
has been done in anticipation of setting $Re=0$ in §4.

The 1-dimensional base state is

$\mathbf{U}=U(y)\mathbf{\hat{x}}:=(1-y^{2})\mathbf{\hat{x}},\quad\mbox{\boldmath$\nabla$}P=-2\mathbf{\hat{x}},$
(6) ${\bf T}=\left[\begin{array}[]{ll}T_{11}&T_{12}\\\
T_{12}&T_{22}\end{array}\right]=(1-\beta)\left[\begin{array}[]{ll}2WU^{\prime
2}&U^{\prime}\\\ U^{\prime}&0\end{array}\right]$ (7)

where $U^{\prime}:=dU/dy$ and henceforth, the analysis is entirely
2-dimensional. The linearized equations for small perturbations

$\mathbf{v}=u\mathbf{\hat{x}}+v\mathbf{\hat{y}}:=\mathbf{u}-\mathbf{U},\quad
p:={\cal
P}-P\quad\&\quad\mbox{\boldmath{$\tau$}}=\left[\begin{array}[]{ll}\tau_{11}&\tau_{12}\\\
\tau_{12}&\tau_{22}\end{array}\right]:=\mbox{\boldmath{${\cal T}$}}-{\bf T}$
(8)

which are all assumed proportional to $e^{ik(x-ct)}$ where $k\in\mathbb{R}$ is
a real wavenumber but the frequency $c=c_{r}+ic_{i}\in\mathbb{C}$ can be
complex ($c_{i}>0$ indicates instability; $c_{r},c_{i}\in\mathbb{R}$), are the
momentum and incompressibility equations,

$\displaystyle Re\left[\,ik(U-c)u+U^{\prime}v\,\right]$
$\displaystyle=-ikp+\beta(D^{2}-k^{2})u+ik\tau_{11}+D\tau_{12}$ (9)
$\displaystyle Re\left[\,ik(U-c)v\,\right]$
$\displaystyle=-Dp+\beta(D^{2}-k^{2})v+ik\tau_{12}+D\tau_{22},$ (10)
$\displaystyle iku+Dv$ $\displaystyle=0$ (11)

for the velocity field and

$\displaystyle\left[\frac{1}{W}+ik(U-c)\right]\tau_{11}$
$\displaystyle=-vDT_{11}+2ikT_{11}u+2T_{12}Du+2U^{\prime}\tau_{12}$
$\displaystyle+\frac{2ik(1-\beta)}{W}u,$ (12)
$\displaystyle\left[\frac{1}{W}+ik(U-c)\right]\tau_{12}$
$\displaystyle=-vDT_{12}+ikT_{11}v+U^{\prime}\tau_{22}$
$\displaystyle+\frac{1-\beta}{W}(Du+ikv),$ (13)
$\displaystyle\left[\frac{1}{W}+ik(U-c)\right]\tau_{22}$
$\displaystyle=\,\,\,\,\,2ikT_{12}v$ $\displaystyle+\frac{2(1-\beta)}{W}Dv$
(14)

for the polymer field where $D:=d/dy$, $T_{11}=2\Lambda U^{\prime 2}$,
$T_{12}=\Lambda U^{\prime}/W$ and $T_{22}=0$ and $\Lambda:=W(1-\beta)$ (see
(7)) in preparation for §4. The pressure $p$ can be eliminated between (9) and
(10) to produce the vorticity equation

$\displaystyle\beta(D^{2}-k^{2})^{2}v=-k^{2}D(\tau_{11}-\tau_{22})$
$\displaystyle+ik(D^{2}+k^{2})\tau_{12}$ $\displaystyle\hskip
28.45274pt+ikRe\left[(U-c)(D^{2}-k^{2})v-U^{\prime\prime}v\right]$ (15)

This equation is good for (asymptotic) analysis but not for a numerical
solution where discretizing two 2nd order equations rather than one 4th order
equation is far better conditioned process.

Figure 1:  The centre mode neutral curve (black dashed line) for $\beta=0.9$
and the asymptotic predictions (upper/lower branch - red/blue solid lines) for
$W$ (top left), $k$ (top right) and $1-c$ (bottom left). The bottom right plot
shows how the lower branch of the neutral curve rotates around and crosses the
inertialess limit of $Re=0$ as $\beta$ increases from $\beta=0.9$ (black
again) through $\beta=0.98$ (green) to $\beta=0.994$ (dark red). The solid
blue lines indicate the asymptotic prediction for the various lower branches
(the $\beta=0.994$ prediction is reached for $Re\rightarrow-\infty$; note only
$Re=O(10)$ is shown here). The dark red dot at $(W,Re)=(973.8,0)$ is the
lowest point reached by the neutral curve for $Re=0$ for an Oldroyd-B fluid
(Khalid et al. (2021b)). The $W$-$Re$ neutral curve in the top left of this
figure is the channel flow equivalent of the pipe flow curve shown on the left
in figure 1 of Dong & Zhang (2022).

##  3 $Re\rightarrow\infty$ Asymptotics for Channel Flow

A natural starting point for examining the centre mode instability is to
consider the neutral curve in the $Re$-$W$ plane for fixed $\beta$ (e.g.
figure 2 of Page et al. (2020), figure 1 here for
$\beta\in\\{0.9,0.98,0.994\\}$ and figure 1 in Dong & Zhang (2022) for pipe
flow). The upper and lower branches of this neutral curve have
$|Re|\rightarrow\infty$ limits which are now explored.

###  3.1 Upper branch in $Re$ vs $W$ plane at fixed $\beta$

Numerical calculations by Khalid et al. (2021a) on the upper branch neutral
curve suggest the scaling behaviour

$(W,k,y,1-c)=\left(\frac{\hat{W}}{\delta},\frac{\hat{k}}{\delta},\delta\hat{Y},\hat{a}\delta^{2}\right)$
(16)

where all hatted variables are $O(\delta^{0})$ for some
$\delta=\delta(Re)\rightarrow 0$ as $Re\rightarrow\infty$. This is the channel
flow equivalent of the short-wavelength scalings for pipe flow studied by Dong
& Zhang (2022) in their §4.1. Rescaling the variables (18)-(23) as follows

$(u,v,p,\tau_{11},\tau_{12},\tau_{22})=\left(\hat{u},\hat{v},\frac{\hat{p}}{\delta},\frac{\hat{\tau}_{11}}{\delta},\frac{\hat{\tau}_{12}}{\delta},\frac{\hat{\tau}_{22}}{\delta}\right)$
(17)

leaves the problem

$\displaystyle
Re\delta^{3}\left[\,i\hat{k}(\hat{a}-\hat{Y}^{2})\hat{u}-2\hat{Y}\hat{v}\,\right]$
$\displaystyle=-i\hat{k}\hat{p}+\beta(\hat{D}^{2}-\hat{k}^{2})\hat{u}+i\hat{k}\hat{\tau}_{11}+\hat{D}\hat{\tau}_{12}$
(18) $\displaystyle
Re\delta^{3}\left[\,i\hat{k}(\hat{a}-\hat{Y}^{2})\hat{v}\,\right]$
$\displaystyle=-\hat{D}\hat{p}+\beta(\hat{D}^{2}-\hat{k}^{2})\hat{v}+i\hat{k}\hat{\tau}_{12}+\hat{D}\hat{\tau}_{22},$
(19) $\displaystyle i\hat{k}\hat{u}+\hat{D}\hat{v}$ $\displaystyle=0,$ (20)
$\displaystyle\left[\frac{1}{\hat{W}}+i\hat{k}(\hat{a}-\hat{Y}^{2})\right]\hat{\tau}_{11}$
$\displaystyle=-16(1-\beta)\hat{W}\hat{Y}\hat{v}+16i\hat{k}(1-\beta)\hat{W}\hat{Y}^{2}\hat{u}$
$\displaystyle\hskip
28.45274pt-4(1-\beta)\hat{Y}\hat{D}\hat{u}-4\hat{Y}\hat{\tau}_{12}+\frac{2ik(1-\beta)}{\hat{W}}\hat{u},$
(21)
$\displaystyle\left[\frac{1}{\hat{W}}+i\hat{k}(\hat{a}-\hat{Y}^{2})\right]\hat{\tau}_{12}$
$\displaystyle=2(1-\beta)\hat{v}+8i\hat{k}(1-\beta)\hat{W}\hat{Y}^{2}\hat{v}-2\hat{Y}\hat{\tau}_{22}$
$\displaystyle\hskip
113.81102pt+\frac{1-\beta}{\hat{W}}(\hat{D}\hat{u}+i\hat{k}\hat{v}),$ (22)
$\displaystyle\left[\frac{1}{\hat{W}}+i\hat{k}(\hat{a}-\hat{Y}^{2})\right]\hat{\tau}_{22}$
$\displaystyle=\,\,\,\,\,-4i\hat{k}(1-\beta)\hat{Y}\hat{v}+\frac{2(1-\beta)}{\hat{W}}\hat{D}\hat{v},$
(23)

where $\hat{D}:=\partial/\partial\hat{Y}=\delta D$ and no terms have been
dropped. The polymer equations are invariant under this scaling (no terms are
dropped) regardless of $\delta$ but $\delta:=Re^{-1/3}$ is forced by the
momentum equation if inertia and viscous effects are to be balanced in the
usual Newtonian way near a critical layer (where $\Re e(c)=U(y)$). With this
choice, no terms are also dropped in the momentum equation so the scaling
transformation is exact here for a parabolic base profile. The one change
going from the original eigenvalue problem to this scaled version is the
position of the boundary which is transformed to $\hat{Y}=\pm\infty$. Solving
the asymptotic eigenvalue problem on the neutral curve is then one of finding
a neutral eigenfunction which decays away at infinity. Given the symmetry of
the centre mode (Khalid et al., 2021a),

$(u,v,p,\tau_{11},\tau_{12},\tau_{22})(-y)=(u,-v,p,\tau_{11},-\tau_{12},\tau_{22})(y),$
(24)

it is sufficient to just impose $\hat{u}=\hat{v}=0$ at some large distance
$\hat{Y}=-L$ ($L\gg 1$) solving over the lower half of the channel and
imposing appropriate symmetry across $y=0$ ($L=15$ to $50$ were used to
explore convergence at $\beta=0.9$): see eigenfunctions in figure 2.

The asymptotic properties $(\hat{W},\hat{k},\,\Re e(\hat{a}))$ of the upper
branch neutral curve are given by seeking

$\max_{\hat{k}}\,\\{\,\hat{W}\,|\,\min_{\hat{k}}\Im
m[\hat{a}(\hat{k},\hat{W})]=0\,\\}$ (25)

in the eigenvalue problem (18)-(23), that is, by finding the largest value of
$\hat{W}$ for which there are no unstable eigenfunctions (the growth rate
$kc_{i}=-\Im m(\hat{a})/Re$) ). The required $\hat{k}$ and $\hat{a}$ are
defined by the neutral eigenfunction at this maximum. The results of this
procedure for $\beta=0.9$ are that

$W\sim 9.1725Re^{1/3},\quad k\sim 0.5023Re^{1/3},\quad c_{r}\sim
1-0.0908Re^{-2/3}$ (26)

on the upper branch neutral curve as $Re\rightarrow\infty$. Table 1 and figure
1 show that this asymptotic result is good down to at least $Re=100$. Using
the elasticity number $E:=W/Re$ as in Khalid et al. (2021a), these scalings
are equivalent to $Re\sim O(E^{-3/2}$), $k\sim O(E^{-1/2})$ and $1-c\sim O(E)$
as $E\rightarrow 0$ which is consistent with figure 11 in Khalid et al.
(2021a).

$Re$ |  | $\hat{W}$ |  | $\hat{k}$ |  | $\hat{a}$
---|---|---|---|---|---|---
100 |  | 8.9409 |  | 0.4912 |  | 0.1018
200 |  | 9.1809 |  | 0.5010 |  | 0.0908
1,000 |  | 9.1730 |  | 0.5020 |  | 0.0909
10,000 |  | 9.1727 |  | 0.5022 |  | 0.0908
$\infty$ |  | 9.1725 |  | 0.5023 |  | 0.0908
Table 1:  Upper branch neutral curve characteristics for $\beta=0.9$ as
$Re\rightarrow\infty$. At a given finite $Re$, $\hat{W}=W/Re^{1/3}$,
$\hat{k}=k/Re^{1/3}$ and $\hat{a}=(1-c)Re^{2/3}$. The bottom line shows the
values given in (26) found by solving the asymptotic problem in (25).

Figure 2:  The (asymptotic) eigenfunction on the upper branch neutral curve
for $\beta=0.9$: $\hat{u}$ (red) and $\hat{v}$ (blue) on the left;
$\hat{\tau}_{11}$ (blue), $\hat{\tau}_{12}$ (red) and $\hat{\tau}_{22}$
(black) on the right (in both real/imaginary parts are solid/dashed). The
eigenfunction has been normalised so that $\hat{u}(0)=1$. The calculation has
been done over the domain $\hat{Y}\in[-15,0]$ for clarity but larger domains
(e.g. $[-50,0]$) were used for convergence purposes.

###  3.2 Lower branch in $Re$ vs $W$ plane at fixed $\beta$

Numerical calculations on the lower branch neutral curve (Khalid et al.
(2021a)) suggest a long-wavelength limit scaling of the following form

$(u,v,p,\tau_{11},\tau_{12},\tau_{22},W,k)=\left(\hat{u},\frac{\hat{v}}{Re},Re\,\hat{p}_{0}+\hat{p}_{1}+O(Re^{-1}),Re\,\hat{\tau}_{11},\hat{\tau}_{12},\frac{\hat{\tau}_{22}}{Re},Re\,\hat{W},\frac{\hat{k}}{Re}\right)$
(27)

where all hatted variables are $O(1)$ as $Re\rightarrow\infty$, $\hat{p}_{0}$
is a constant ($D\hat{p}_{0}=0$) and $c$ stays $O(1)$ and bounded away from
$0$ (the wall advection speed) and $1$ (the centreline advection speed): see
figure 1 (bottom left). In their long-wavelength analysis for pipe flow (their
§3), Dong & Zhang (2022) only consider $1/Re\ll k\ll 1$ and so do not treat
the lower branch neutral curve where again $k=O(1/Re)$ is found numerically
(Garg et al. (2018)). With these rescalings, (18)-(23) become

$\displaystyle i\hat{k}(1-y^{2}-c)\hat{u}-2y\hat{v}$
$\displaystyle=-i\hat{k}\hat{p}_{0}+\beta
D^{2}\hat{u}+i\hat{k}\hat{\tau}_{11}+D\hat{\tau}_{12}$ (28) $\displaystyle
i\hat{k}(1-y^{2}-c)\hat{v}\hskip 28.45274pt$
$\displaystyle=-D\hat{p}_{1}+\beta
D^{2}\hat{v}+i\hat{k}\hat{\tau}_{12}+D\hat{\tau}_{22},$ (29) $\displaystyle
i\hat{k}\hat{u}+D\hat{v}$ $\displaystyle=0,$ (30)

for the velocity field and, for the polymer stress,

$\displaystyle\left[\frac{1}{\hat{W}}+i\hat{k}(1-y^{2}-c)\right]\hat{\tau}_{11}$
$\displaystyle=-16(1-\beta)\hat{W}y\,\hat{v}+16i\hat{k}(1-\beta)\hat{W}y^{2}\,\hat{u}-4y\,\hat{\tau}_{12}$
$\displaystyle\hskip 156.49014pt-4(1-\beta)yD\hat{u},$ (31)
$\displaystyle\left[\frac{1}{\hat{W}}+i\hat{k}(1-y^{2}-c)\right]\hat{\tau}_{12}$
$\displaystyle=\,2(1-\beta)\hat{v}+8i\hat{k}(1-\beta)\hat{W}y^{2}\,\hat{v}-2y\,\hat{\tau}_{22}+\frac{1-\beta}{\hat{W}}D\hat{u},$
(32)
$\displaystyle\left[\frac{1}{\hat{W}}+i\hat{k}(1-y^{2}-c)\right]\hat{\tau}_{22}$
$\displaystyle=\,\,\,\,\,-4i\hat{k}(1-\beta)y\,\hat{v}+\frac{2(1-\beta)}{\hat{W}}D\hat{v}.$
(33)

Since $D\hat{p}_{0}=0$, differentiating (28) leads directly to the vorticity
equation

$\beta
D^{4}\hat{v}=-\hat{k}^{2}D\hat{\tau}_{11}+i\hat{k}D^{2}\hat{\tau}_{12}+i\hat{k}\left[(1-y^{2}-c)D^{2}\hat{v}+2\hat{v}\right]$
(34)

and (29), which just defines $\hat{p}_{1}$, can be ignored. The problem
defined by (30)-(34) is then an eigenvalue problem for $c$. Since there is no
rescaling of the spatial dimension, the neutral eigenfunction is global and
easily resolved. The asymptotic properties $(\hat{W},\hat{k},\,c$ of the lower
branch neutral curve are given by seeking

$\max_{\hat{k}}\,\\{\,\hat{W}\,|\,\max_{\hat{k}}\Im
m[c(\hat{k},\hat{W})]=0\,\\}$ (35)

and the results are shown in Table 2. The asymptotic scalings for $\beta=0.9$
where $Re\rightarrow\infty$ are the same as for $\beta=0.994$ where
$Re\rightarrow-\infty$ but the eigenfunctions looks distinctly different in
the polymer stress field - see figure 3.

The lower branch scalings are apparent in figure 13 of Khalid et al. (2021a)
(see also their figure 18). The upper branch asymptote is reached for
$Re_{c}\rightarrow\infty$ and $E\rightarrow 0$ whereas the vertical asymptote
$E\rightarrow E_{\infty}$ (a finite value) as $Re_{c}\rightarrow\infty$
corresponds to the lower branch asymptote and the results in Table 2 can be
used to predict $E_{c}$. For example $(1-\beta)E_{c}=(1-\beta)\hat{W}=0.1058$
at $\beta=0.9$ and $(1-\beta)E_{c}=(1-\beta)\hat{W}=0.352$ at $\beta=0.98$
(note $E_{c}\lesssim\max_{Re_{c}}{E}$ for a given $\beta$ in figure 13 in
Khalid et al. (2021a)). Khalid et al. (2021b) (their figure 4) show that there
is no asymptote for $\beta>0.990552$. Instead the asymptote has to flip to
$Re\rightarrow-\infty$ and $E_{c}<0$ as shown for example with $\beta=0.994$
in figure 1 (bottom right).

$\beta$ |  | $\hat{W}$ |  | $\hat{k}$ |  | $c$
---|---|---|---|---|---|---
0.9 |  | 1.058 |  | 62.17 |  | 0.999553
0.98 |  | 17.60 |  | 13.32 |  | 0.998788
0.994 |  | -128.1 |  | -5.97 |  | 0.999198
Table 2: Lower branch neutral curve characteristics for $\beta=0.9$ and $0.98$
as $Re\rightarrow\infty$ and for $\beta=0.994$ as $Re\rightarrow-\infty$
(hence $\hat{W}$ and $\hat{k}$ are both negative). A comparison with results
from finite $Re$ calculations are shown in figure 1 (bottom right).

Figure 3:  The (asymptotic) eigenfunctions on the lower branch neutral curve
for $\beta=0.9$ (top row) and $\beta=0.994$ (bottom row): $\hat{u}$ (red) and
$\hat{v}$ (blue) on the left; $\hat{\tau}_{11}$ (blue), $\hat{\tau}_{12}$
(red) and $\hat{\tau}_{22}$ (black) on the right (in both real/imaginary parts
are solid/dashed). Both eigenfunctions have been normalised so that
$\hat{u}(0)=1$.

##  4 $W\rightarrow\infty$ Asymptotics for Inertialess ($Re=0$) Channel Flow

Without inertia ($Re=0$), the relevant asymptotic limit is
$W\rightarrow\infty$ and $\beta\rightarrow 1$ such that $W(1-\beta)=\Lambda$
is an $O(1)$ constant, e.g. see insets A and B of figure 2 in Khalid et al.
(2021b) (or figure 8 in their supplementary material) which suggest
$3.5\lesssim\Lambda\lesssim 10$ for instability. Physically, of course, this
means $W$ is large but finite whereas $Re$ can be considered separately as
small as desired but is strictly not zero as there is flow. The point is
mathematically, the $Re\rightarrow 0$ limit is regular so it is convenient to
set $Re=0$ to get the true limiting values of key dependencies (e.g. how $1-c$
scales with $W$ on the neutral curve).

As already mentioned in §3 and $Re>0$, the centre mode instability has a
certain symmetry about the midplane $y=0$: $u$ is symmetric and $v$
antisymmetric: see (24). Henceforth we only consider $y\in[-1,0]$ and impose
non-slip boundary conditions $v(-1)=Dv(-1)=0$ at the solid lower plate and
symmetry conditions $v(0)=D^{2}v(0)=0$ at the midplane. Numerically (see
Appendix A for details), we find on the neutral curve that the eigenfunction
has a critical layer near the midplane across which $v$ is continuous but
there are jumps in $\Re e(Dv)$ and $\Re e(D^{2}v)$ and singular-looking
behaviour for $\Im m(Dv)$ where the phase of the eigenfunction is set by
making $\Re e(Dv)=1$ at the midplane: see figure 4.

A key issue is whether the critical layer at $y=y^{*}$, defined by
$U(y^{*})=c_{r}$ so $y^{*}:=-\sqrt{1-c_{r}}\in[-1,0]$, approaches the midplane
as it does in classical Orr-Sommerfeld analysis for Newtonian shear flows
(Drazin & Reid, 1981) or not. Certainly figure 4 suggests ‘not’ and earlier
(pre-shooting code) attempts to develop the asymptotic structure could not
reconcile $c_{r}\rightarrow 1$ as $W\rightarrow\infty$ with an O(1) jump in
$D^{2}v$ across the critical layer. This means a novel aspect of the
asymptotics here is that despite $1-c_{r}$ being very small, it does in fact
remain O(1) as $W\rightarrow\infty$: see Table 3.

| $W$ |  | $c_{r}$ | $k\Lambda_{-}$ |  | $c_{r}$ | $k\Lambda_{+}$
---|---|---|---|---|---|---|---
$k=0.1$ | 32,000 |  | 0.998390 | 3.41331 |  | 0.999893 | 7.62705
| 128,000 |  | 0.998408 | 3.39349 |  | 0.999943 | 9.25320
| 512,000 |  | 0.998412 | 3.38849 |  | 0.999955 | 9.76344
| *2,048,000 |  | 0.998413 | 3.38724 |  | 0.999957 | 9.80010
| *8,192,000 |  | 0.998413 | 3.38693 |  | 0.999957 | 9.81447
| *32,768,000 |  | 0.998413 | 3.38686 |  | 0.999957 | 9.81826
| $\infty$ |  | 0.998413 | 3.38683 |  | 0.999957 | 9.81954
$k=1.1$ | 64,000 |  | 0.9992186 | 4.526249 |  | 0.999853 | 7.63580
| 128,000 |  | 0.9992183 | 4.523676 |  | 0.999855 | 7.65132
| 256,000 |  | 0.9992181 | 4.522393 |  | 0.999855 | 7.65932
| *1,024,000 |  | 0.9992180 | 4.521432 |  | 0.999856 | 7.66538
| *8,192,000 |  | 0.9992180 | 4.521154 |  | 0.999856 | 7.66715
| $\infty$ |  | 0.9992180 | 4.521114 |  | 0.999856 | 7.66740
Table 3:  $c_{r}$ and $k\Lambda:=k(1-\beta)W$ as plotted in inset A of figure
2 in Khalid et al. (2021b) on the neutral curve for $k=0.1$ and $1.1$ at
various large $W$ ($\pm$ indicates upper and lower parts of the neutral curve
and * results computed using quadruple rather than double precision). The
$\infty$ entry comes from Richardson extrapolation eliminating the leading
$O(1/W)$ error. In all cases $\Lambda$ is slower to converge than $c_{r}$ with
the upper curve calculations on the right highlighting this.

We introduce a small parameter

$\varepsilon:=\frac{1}{W}$ (36)

and take the distinquished limit $\beta=1-\varepsilon\Lambda$ where $\Lambda$
is an O(1) number to be determined. The eigenfunction plotted in figure 4
shows abrupt changes in the solution as it crosses a critical layer at
$y=c_{r}$. The solution either side of the critical layer - the ‘outer’
solution - must satisfy the governing equations with $\beta=1$,
$\varepsilon=0$ and $\Lambda=O(1)$ with the critical layer supplying
appropriate ‘matching’ conditions between the two parts. The purpose of the
asymptotic analysis developed below is to identify analytic expressions for
these matching conditions so that the two parts of the outer solution can be
fitted together in the limit of $W\rightarrow\infty$ without solving for the
critical layer. This defines the leading solution to the problem which
includes the leading $O(W^{0})$ value of $c$.

Figure 4:  Neutral eigenfunction at W=32,000 (left full plot), 128,000 (middle
partial plot) and 512,000 (right partial plot), $k=1.1$ and $\Lambda\approx
4.11=W(1-\beta)$ and $c_{r}\approx 0.99218$. Vertical dotted line near $y=0$
is the critical layer where $U=1-y^{2}=c_{r}$. $v$ is blue (real/imaginary
parts solid/dashed respectively), real part of $Dv$ is red and imaginary part
is black. Notice the O(1) jump in $\Re e(Dv)$ (and in $\Re e(D^{2}v)$) across
the critical layer and the increasingly singular behaviour of $\Im m(Dv)$
(black solid line) as $W$ increases.

###  4.1 Outer solution

We refer to the ‘outer’ solution as the solution in the regions
$y-y^{*}=O(1)$. Assuming $v=O(1)$, then $\tau_{11}$ and $\tau_{12}$ are $O(1)$
whereas $\tau_{22}=O(1/W)$. The leading order outer problem is then

$\displaystyle(D^{2}-k^{2})^{2}v$
$\displaystyle=-k^{2}D\tau_{11}+ik(D^{2}+k^{2})\tau_{12},$ (37) $\displaystyle
iku+Dv$ $\displaystyle=0,$ (38) $\displaystyle ik(U-c)\tau_{11}$
$\displaystyle=-vDT_{11}+2ikT_{11}u+2U^{\prime}\tau_{12},$ (39) $\displaystyle
ik(U-c)\tau_{12}$ $\displaystyle=ikT_{11}v,$ (40)

which can be simplified to the 4th order problem

$(D^{2}-k^{2})^{2}v=-ikD\left[T_{11}D\left(\frac{v}{U-c}\right)\right]+\frac{ik^{3}T_{11}v}{U-c}.$
(41)

For $y<y^{*}$, outer boundary conditions are $v(-1)=0,Dv(-1)=0$ and for
$y>y^{*}$ $v(0)=0=D^{2}v(0)$ with the critical layer supplying 4 matching
conditions (for $v$, $Dv$, $D^{2}v$ and $D^{3}v$ respectively). This will
produce a well posed eigenvalue problem for $c(\Lambda,k)$. Ultimately, the
problem is to find $\min\Lambda$ over all pairs $(\Lambda,k)$ where
$c_{i}(\Lambda,k)=0$.

###  4.2 Inner solution

The inner solution is the solution in the critical layer
$y-y^{*}=O(\varepsilon)$ where the thickness comes from balancing polymer
advection and relaxation processes: see the left hand sides of the perturbed
polymer stress equations (12)-(14). We therefore define a critical layer
variable $Y$ and the corresponding derivative $\hat{D}$,

$Y:=\frac{y-y^{*}}{\varepsilon},\qquad
D=\frac{1}{\varepsilon}\hat{D}:=\frac{1}{\varepsilon}\frac{\partial}{\partial
Y},$ (42)

so that, for example,
$1/W+ik(U-c)=\varepsilon(1+ikU_{*}^{{}^{\prime}}Y+ikU_{*}^{{}^{\prime\prime}}\varepsilon
Y^{2})$ where $U_{*}^{{}^{\prime}}:=U^{{}^{\prime}}(y^{*})=-2y^{*}>0.$ Then
the appropriate expansions turn out to be

$\displaystyle v$ $\displaystyle=\hskip
14.22636pt\hat{v}_{0}+\varepsilon\log\varepsilon\,\Pi
Y+\varepsilon\,\hat{v}_{1}(Y)+\varepsilon^{2}\log\varepsilon\biggl{(}\hat{v}_{\ell}(Y)+\tfrac{1}{2}\Omega
Y^{2}\biggr{)}+\varepsilon^{2}\,\hat{v}_{2}(Y)$ $\displaystyle\hskip
142.26378pt+\varepsilon^{3}\log\varepsilon\biggl{(}\hat{v}_{2\ell}+\tfrac{1}{6}\Upsilon
Y^{3}\biggr{)}+\varepsilon^{3}\,\hat{v}_{3}(Y)+\ldots$ (43) $\displaystyle u$
$\displaystyle=\frac{i}{k}\left[\hskip
31.2982pt\log\varepsilon\,\Pi+\hat{D}\hat{v}_{1}(Y)+\varepsilon\log\varepsilon\biggl{(}\hat{D}\hat{v}_{\ell}(Y)+\Omega
Y\biggr{)}+\varepsilon\hat{D}\hat{v}_{2}(Y)\right.$ $\displaystyle\hskip
128.0374pt+\varepsilon^{2}\log\varepsilon\biggl{(}\hat{D}\hat{v}_{2\ell}+{\textstyle{1\over
2}}\Upsilon Y^{2}\biggr{)}+\varepsilon^{2}\hat{D}\hat{v}_{3}(Y)+\ldots\ $ (44)
$\displaystyle Du$ $\displaystyle=\frac{i}{k}\left[\hskip
65.44142pt\frac{1}{\varepsilon}\hat{D}^{2}\hat{v}_{1}(Y)+\log\varepsilon\biggl{(}\hat{D}^{2}\hat{v}_{\ell}(Y)+\Omega\biggr{)}+\hat{D}^{2}\hat{v}_{2}(Y)\right.$
$\displaystyle\hskip
142.26378pt+\varepsilon\log\varepsilon\biggl{(}\hat{D}^{2}\hat{v}_{2\ell}+\Upsilon
Y\biggr{)}+\varepsilon\hat{D}^{2}\hat{v}_{3}(Y)+\ldots$ (45) $\displaystyle
D^{2}u$ $\displaystyle=\frac{i}{k}\left[\hskip
65.44142pt\frac{1}{\varepsilon^{2}}\hat{D}^{3}\hat{v}_{1}(Y)+\frac{\log\varepsilon}{\varepsilon}\hat{D}^{3}\hat{v}_{\ell}(Y)+\frac{1}{\varepsilon}\hat{D}^{3}\hat{v}_{2}(Y)\right.$
$\displaystyle\hskip
156.49014pt+\log\varepsilon\biggl{(}\hat{D}^{3}\hat{v}_{2\ell}+\Upsilon\biggr{)}+\hat{D}^{3}\hat{v}_{3}(Y)+\ldots$
(46)

for the velocity fields and

$\displaystyle\tau_{11}$
$\displaystyle=\frac{1}{\varepsilon^{2}}\hat{\tau}_{11}^{0}(Y)+\frac{\log\varepsilon}{\varepsilon}\hat{\tau}_{11}^{\ell}(Y)+\frac{1}{\varepsilon}\hat{\tau}_{11}^{1}(Y)+\log\varepsilon\,\hat{\tau}_{11}^{2\ell}+\hat{\tau}_{11}^{2}(Y)+\ldots,$
(47) $\displaystyle\tau_{12}$
$\displaystyle=\frac{1}{\varepsilon}\hat{\tau}_{12}^{0}(Y)+\log\varepsilon\,\hat{\tau}_{12}^{\ell}(Y)+\hat{\tau}_{12}^{1}(Y)+\varepsilon\log\varepsilon\,\hat{\tau}_{12}^{2\ell}+\varepsilon\hat{\tau}_{12}^{2}(Y)+\ldots,$
(48) $\displaystyle\tau_{22}$
$\displaystyle=\hat{\tau}_{22}^{0}(Y)+\varepsilon\log\varepsilon\,\hat{\tau}_{22}^{\ell}(Y)+\varepsilon\hat{\tau}_{22}^{1}(Y)+\varepsilon^{2}\log\varepsilon\,\hat{\tau}_{22}^{2\ell}+\varepsilon^{2}\hat{\tau}_{22}^{2}(Y)+\ldots$
(49)

for the polymer stresses where $\Pi$, $\Omega$ and $\Upsilon$ are complex
constants which will emerge below. The reason we need to go so deep into these
expansions is the outer problem is 4th order and therefore requires jump
conditions down to $D^{3}v$ plus there is a singularity at the critical layer.
Together these require considering the equation for
$\hat{D}^{3}\hat{v}_{3}(Y)$ which is $O(\varepsilon^{2})$ down in the
expansion of $D^{3}v$ i.e. we need to go to third order in the expansion. The
intermediate $O(\log\varepsilon)$ terms are formally needed to fix up the
logarithmic terms which arise but turn out to be unimportant for deriving the
matching conditions.

To keep track of the influence of $U^{\prime}$ and $U^{{}^{\prime\prime}}$
when we probe later why plane Couette flow doesn’t have a neutral curve, it is
useful to expand the base state around $y=y_{*}$ in the critical layer as
follows,

$\displaystyle U$ $\displaystyle=U_{*}+\varepsilon
YU_{*}^{{}^{\prime}}+{\textstyle{1\over
2}}\varepsilon^{2}Y^{2}U_{*}^{{}^{\prime\prime}}+\ldots,$ (50) $\displaystyle
DU$ $\displaystyle=U_{*}^{{}^{\prime}}+\varepsilon
YU_{*}^{{}^{\prime\prime}}+\ldots,$ (51) $\displaystyle T_{11}$
$\displaystyle=T_{11}^{*(0)}+\varepsilon
YT_{11}^{*(1)}+\varepsilon^{2}Y^{2}T_{11}^{*(2)}+\ldots,$ (52) $\displaystyle
DT_{11}$ $\displaystyle=T_{11}^{*(1)}+2\varepsilon YT_{11}^{*(2)}+\dots,$ (53)
$\displaystyle T_{12}$ $\displaystyle=\hskip 28.45274pt+\varepsilon
T_{12}^{*(1)}\hskip 7.11317pt+\varepsilon^{2}YT_{12}^{*(2)}+\ldots,$ (54)
$\displaystyle DT_{12}$ $\displaystyle=\hskip 28.45274pt+\varepsilon
T_{12}^{*(2)}+\ldots$ (55)

where

$T_{11}^{*(0)}:=2\Lambda U_{*}^{{}^{\prime}2},\quad T_{11}^{*(1)}:=4\Lambda
U_{*}^{{}^{\prime}}U_{*}^{{}^{\prime\prime}},\quad T_{11}^{*(2)}:=2\Lambda
U_{*}^{{}^{\prime\prime}2}$ (57)

and

$T_{12}^{*(1)}:=\Lambda U_{*}^{{}^{\prime}},\quad T_{12}^{*(2)}:=\Lambda
U_{*}^{{}^{\prime\prime}}$ (58)

assuming that $U_{*}^{{}^{\prime\prime\prime}}=0$ for simplicity (true for
both channel and Couette flow).

Now we substitute expansions (43)-(46) for the perturbation velocity field,
(47)-(49) for the perturbation polymer stresses and (50)-(54) for the base
state into the equations (11)-(15) and collect similar order terms to create a
hierarchy of problems in the usual way.

###  4.3 Leading order: $O(1/\varepsilon^{3})$ in the Stokes equation

At leading order,

$\displaystyle X\hat{\tau}_{22}^{0}$
$\displaystyle=2ikT_{12}^{*(1)}\hat{v}_{0}$ (59) $\displaystyle
X\hat{\tau}_{12}^{0}$
$\displaystyle=ikT_{11}^{*(0)}\hat{v}_{0}+U_{*}^{{}^{\prime}}\hat{\tau}_{22}^{0}$
(60) $\displaystyle X\hat{\tau}_{11}^{0}$
$\displaystyle=2U_{*}^{{}^{\prime}}\hat{\tau}_{12}^{0}$ (61)
$\displaystyle\hat{D}^{4}\hat{v}_{1}$
$\displaystyle=-k^{2}\hat{D}\hat{\tau}_{11}^{0}+ik\hat{D}^{2}\hat{\tau}_{12}^{0}$
(62)

where $X:=(1+ikU_{*}^{{}^{\prime}}Y)$. Integrating the last equation twice
with respect to $Y$ gives

$D^{2}v=\frac{1}{\varepsilon}\hat{D}^{2}\hat{v}_{1}=-\frac{k^{2}}{\varepsilon}\int\hat{\tau}_{11}^{0}dY+\frac{ik}{\varepsilon}\hat{\tau}_{12}^{0}=\frac{k^{2}T_{11}^{*(0)}\hat{v}_{0}}{\varepsilon
X}\sim\frac{-ikT_{11}^{*(0)}\hat{v}_{0}}{U_{*}^{{}^{\prime}}(y-y^{*})}\quad{\rm
as}\quad Y\rightarrow\pm\infty$ (63)

as the $O(1/\varepsilon)$ integration constants must be zero otherwise the
solution cannot be matched with the central region. This leading solution for
$D^{2}v$ immediately suggests that the asymptotic matching will be a
challenge. There is a simple pole singularity in $D^{2}v$ at the critical
layer and as a consequence, a double pole singularity in $D^{3}v$. Since a
double pole is symmetric across the critical layer, it does not enter into the
matching conditions for $D^{3}v$ but will certainly obscure any matching
criterion present involving higher order, less singular behaviour.

Forewarned, we press on and integrate once to give

$Dv=\hat{D}\hat{v}_{1}=\frac{-ikT_{11}^{*(0)}\hat{v}_{0}}{U_{*}^{{}^{\prime}}}\biggl{[}\log
X+\alpha_{1}\biggr{]}$ (64)

where $\alpha_{1}$ is another complex constant). Matching to the exterior
requires that a
$-ikT_{11}^{*(0)}\hat{v}_{0}/U_{*}^{{}^{\prime}}\log\varepsilon$ term is
present in the expression for $Dv$ so then

$Dv\sim\frac{-ikT_{11}^{*(0)}\hat{v}_{0}}{U_{*}^{{}^{\prime}}}\biggl{[}\log[ikU_{*}^{{}^{\prime}}(y-y^{*})]+\alpha_{1}\biggr{]}\quad{\rm
for}\quad\varepsilon\ll|y-y^{*}|$ (65)

This means that the first complex coefficient $\Pi$ in the expansions
(43)-(46) has to be
$\Pi=-ikT_{11}^{*(0)}\hat{v}_{0}/U_{*}^{{}^{\prime}}=-2ik\Lambda
U_{*}^{{}^{\prime}}\hat{v}_{0}$ and so

$Dv\rightarrow\left\\{\begin{array}[]{cc}\Pi\biggl{(}\log|kU_{*}^{{}^{\prime}}(y-y^{*})|+\alpha_{1}-\tfrac{1}{2}i\pi\,\biggr{)}&y\rightarrow
y^{*-}\\\ &\\\
\Pi\biggl{(}\log|ikU_{*}^{{}^{\prime}}(y-y^{*})|+\alpha_{1}+\tfrac{1}{2}i\pi\,\biggr{)}&y\rightarrow
y^{*+}\end{array}\right.$ (66)

This gives a jump in $Dv$ across the critical layer of $i\pi\Pi$. Integrating
(64) gives

$\hat{v}_{1}=\Pi\left[\frac{X\log
X}{ikU_{*}^{{}^{\prime}}}+(\alpha_{1}-1)Y\right]$ (67)

since $\hat{v}_{1}(y_{*})=0$ as $v(y_{*})=\hat{v}_{0}$ by definition.

###  4.4 Next order: $O(\log\varepsilon/\varepsilon^{2})$ in the Stokes
equation

Working to next order

$\displaystyle\hat{\tau}_{22}^{\ell}$ $\displaystyle=\frac{2ikT_{12}^{*(1)}\Pi
Y+2\Lambda\Pi}{X}=2\Lambda\Pi$ (68) $\displaystyle\hat{\tau}_{12}^{\ell}$
$\displaystyle=\frac{ikT_{11}^{*(0)}\Pi
Y+U_{*}^{{}^{\prime}}\hat{\tau}_{22}^{\ell}}{X}=2\Lambda\Pi
U_{*}^{{}^{\prime}}$ (69) $\displaystyle\hat{\tau}_{11}^{\ell}$
$\displaystyle=\frac{2\Pi(-T_{11}^{*(0)}X+2ikU_{*}^{{}^{\prime}2}T_{12}^{*(1)}+2U_{*}^{{}^{\prime}2}\Lambda)}{X^{3}}=0$
(70)

so $\hat{D}^{4}\hat{v}_{\ell}=0$ and therefore $\hat{v}_{\ell}=0$.

###  4.5 $O(1/\varepsilon^{2})$ in the Stokes equation

This is the order which will complete the jump condition in $D^{2}v$. At
$O(1/\varepsilon^{2})$, we get

$\displaystyle X\hat{\tau}_{22}^{1}$ $\displaystyle=-{\textstyle{1\over
2}}ikU_{*}^{{}^{\prime\prime}}Y^{2}\hat{\tau}_{22}^{0}+2ikT_{12}^{*(2)}\hat{v}_{0}+2ikT_{12}^{*(1)}\hat{v}_{1}+2\Lambda\hat{D}\hat{v}_{1}$
(71) $\displaystyle X\hat{\tau}_{12}^{1}$ $\displaystyle=-{\textstyle{1\over
2}}ikU_{*}^{{}^{\prime\prime}}Y^{2}\hat{\tau}_{12}^{0}+(ikT_{11}^{*(1)}Y-T_{12}^{*(2)})\hat{v}_{0}+ikT_{11}^{*(0)}\hat{v}_{1}+U_{*}^{{}^{\prime}}\hat{\tau}_{22}^{1}+U_{*}^{{}^{\prime\prime}}Y\hat{\tau}_{22}^{0}+\tfrac{i\Lambda}{k}\hat{D}^{2}\hat{v}_{1}$
(72) $\displaystyle X\hat{\tau}_{11}^{1}$ $\displaystyle=-{\textstyle{1\over
2}}ikU_{*}^{{}^{\prime\prime}}Y^{2}\hat{\tau}_{11}^{0}-T_{11}^{*(1)}\hat{v}_{0}-2T_{11}^{*(0)}\hat{D}\hat{v}_{1}+\tfrac{2i}{k}T_{12}^{*(1)}\hat{D}^{2}\hat{v}_{1}+2U_{*}^{{}^{\prime}}\hat{\tau}_{12}^{1}+2U_{*}^{{}^{\prime\prime}}Y\hat{\tau}_{12}^{0}$
(73)

with

$\hat{D}^{4}\hat{v}_{2}=\Lambda\hat{D}^{4}\hat{v}_{1}-k^{2}\hat{D}\hat{\tau}_{11}^{1}+ik\hat{D}^{2}\hat{\tau}_{12}^{1}.$
(74)

Integrating twice

$\hat{D}^{2}\hat{v}_{2}=\Lambda\hat{D}^{2}\hat{v}_{1}-k^{2}\int\hat{\tau}_{11}^{1}\,dY+ik\hat{\tau}_{12}^{1}+\Theta
Y+\Phi.$ (75)

where $\Theta$ and $\Phi$ are complex constants. $\Theta Y$ is unmatchable in
the interior so $\Theta$ must be $0$. In terms of deriving jump conditions
across the critical layer, the presence of the constant $\Phi$ means we are
only interested in the asymptotic behaviour (as $Y\rightarrow\pm\infty$) of
the RHS of (75) which gives rise to jumps across the layer. With this in mind,
it is straightforward to show from (71) and (72) that

$ik\hat{\tau}_{12}^{1}=\frac{k^{2}[T_{11}^{*(0)}]^{2}}{U_{*}^{{}^{\prime}2}}\hat{v}_{0}\log
X+const\qquad{\rm as}\quad Y\rightarrow\pm\infty$ (76)

The other term on the RHS,

$\displaystyle-k^{2}\int\hat{\tau}_{11}^{1}\,dY=\int$
$\displaystyle\overbrace{\frac{{\textstyle{1\over
2}}k^{2}U_{*}^{{}^{\prime\prime}}Y^{2}\hat{\tau}_{11}^{0}}{U_{*}^{{}^{\prime}}X}}^{(i)}\overbrace{-\frac{ikT_{11}^{*(0)}\hat{v}_{0}}{U_{*}^{{}^{\prime}}X}}^{(ii)}\overbrace{-\frac{2ikT_{11}^{*(0)}\hat{D}\hat{v}_{1}}{U_{*}^{{}^{\prime}}X}}^{(iii)}$
$\displaystyle-\frac{2T_{12}^{*(1)}\hat{D}^{2}\hat{v}_{1}}{U_{*}^{{}^{\prime}}X}+\underbrace{\frac{2ikU_{*}^{{}^{\prime\prime}}Y\hat{\tau}_{12}^{0}}{U_{*}^{{}^{\prime}}X}}_{(iv)}+\underbrace{\frac{2ik\hat{\tau}_{12}^{1}}{X}}_{(v)},$
(77)

is more involved with each labelled term contributing. Respectively, as
$Y\rightarrow\pm\infty$,

$\displaystyle(i)$
$\displaystyle\rightarrow-\frac{ikU_{*}^{{}^{\prime\prime}}T_{11}^{*(0)}}{U_{*}^{{}^{\prime}2}}\hat{v}_{0}\log
X$ $\displaystyle(ii)$
$\displaystyle\rightarrow-\frac{kT_{11}^{*(1)}}{U_{*}^{{}^{\prime}}}\hat{v}_{0}\log
X$ $\displaystyle(iii)$
$\displaystyle\rightarrow-\frac{2k^{2}[T_{11}^{*(0)}]^{2}\,\hat{v}_{0}}{U_{*}^{{}^{\prime}2}}\biggl{[}{\textstyle{1\over
2}}(\log X)^{2}+\alpha_{1}\log X\biggr{]}$ $\displaystyle(iv)$
$\displaystyle\rightarrow\frac{2ikU_{*}^{{}^{\prime\prime}}T_{11}^{*(0)}\hat{v}_{0}}{U_{*}^{{}^{\prime}2}}\log
X$

The $(v)^{th}$ term has to be further subdivided as follows

$\displaystyle\frac{ik}{U_{*}^{{}^{\prime}}}\int\frac{2U_{*}^{{}^{\prime}}\hat{\tau}_{12}^{1}}{X}\,dY=\int$
$\displaystyle\overbrace{\frac{k^{2}U_{*}^{{}^{\prime\prime}}Y^{2}\hat{\tau}_{12}^{0}}{X^{2}}}^{a}-\frac{2ik\Lambda
U_{*}^{{}^{\prime\prime}}\hat{v}_{0}}{X^{2}}+\overbrace{\frac{-2k^{2}T_{11}^{*(1)}Y\hat{v}_{0}}{X^{2}}}^{b}$
$\displaystyle+\underbrace{\frac{-2k^{2}T_{11}^{*(0)}\hat{v}_{1}}{X^{2}}}_{c}+\frac{2ikU_{*}^{{}^{\prime}}\hat{\tau}_{22}^{1}}{X^{2}}+\frac{2ikU_{*}^{{}^{\prime\prime}}Y\hat{\tau}_{22}^{0}}{X^{2}}-\frac{2\Lambda\hat{D}^{2}\hat{v}_{1}}{X^{2}},$
(78)

with the respective asymptotic behaviour s $Y\rightarrow\pm\infty$,

$\displaystyle(v)_{a}$
$\displaystyle\rightarrow-\frac{ikU_{*}^{{}^{\prime\prime}}T_{11}^{*(0)}\hat{v}_{0}}{U_{*}^{{}^{\prime}2}}\log
X$ $\displaystyle(v)_{b}$
$\displaystyle\rightarrow\frac{2ikT_{11}^{*(1)}\hat{v}_{0}}{U_{*}^{{}^{\prime}}}\log
X$ $\displaystyle(v)_{c}$
$\displaystyle\rightarrow\frac{2k^{2}[T_{11}^{*(0)}]^{2}\,\hat{v}_{0}}{U_{*}^{{}^{\prime}2}}\biggl{[}{\textstyle{1\over
2}}(\log X)^{2}+(\alpha_{1}-1)\log X\biggr{]}$

Adding all the contributions

$\displaystyle\hat{D}^{2}\hat{v}_{2}$
$\displaystyle\sim\biggl{[}\frac{-k^{2}[T_{11}^{*(0)}]^{2}}{U_{*}^{{}^{\prime}2}}+\frac{ikT_{11}^{*(1)}}{U_{*}^{{}^{\prime}}}\biggr{]}\log
X+const$ $\displaystyle\hskip 85.35826pt=(\,4ik\Lambda
U_{*}^{{}^{\prime\prime}}-4k^{2}\Lambda^{2}U_{*}^{{}^{\prime}2}\,)\hat{v}_{0}\log[ikU_{*}^{{}^{\prime}}(y-y_{*})]+const,$
$\displaystyle\hskip
85.35826pt=\Omega\log[ikU_{*}^{{}^{\prime}}(y-y_{*})]+const.$ (79)

Therefore the second complex coefficient $\Omega$ in the expansions (43)-(46)
has to be $\Omega=(\,4ik\Lambda
U_{*}^{{}^{\prime\prime}}-4k^{2}\Lambda^{2}U_{*}^{{}^{\prime}2}\,)\hat{v}_{0}$
and this indicates the finite jump in $D^{2}v$ across the critical layer.
Integrating (79) twice gives

$\displaystyle\hat{D}\hat{v}_{2}$
$\displaystyle=-\frac{i\Omega}{kU_{*}^{{}^{\prime}}}\biggl{(}X\log
X-X\biggr{)}+\alpha_{2},$ $\displaystyle\hat{v}_{2}$
$\displaystyle=-\frac{\Omega}{k^{2}U_{*}^{{}^{\prime}2}}\biggl{(}\frac{1}{2}X^{2}\log
X-\frac{3}{4}X^{2}\biggr{)}+\alpha_{2}Y+\beta_{2},$ (80)

where $\alpha_{2}$ and $\beta_{2}$ are constants. Since $\hat{v}_{2}=0$ at
$Y=0\,(X=1)$ by definition of $\hat{v}_{0}$,
$\beta_{2}=-3\Omega/(4k^{2}U_{*}^{{}^{\prime}2})$.

###  4.6 $O(\log\varepsilon/\varepsilon)$ in the Stokes equation

The $O(\log\varepsilon/\varepsilon)$ Stokes equation balance integrated once
gives

$\hat{D}^{3}\hat{v}_{2\ell}=-k^{2}\hat{\tau}_{11}^{2\ell}+ik\hat{D}\hat{\tau}_{12}^{2\ell}+const$
(81)

as $\hat{v}_{\ell}=0$. Now, since

$\displaystyle\hat{\tau}_{22}^{2\ell}$
$\displaystyle=\frac{1}{X}\biggl{(}-{\textstyle{1\over
2}}ikU_{*}^{{}^{\prime\prime}}Y^{2}\hat{\tau}_{22}^{\ell}+ik\Omega
T_{12}^{*(1)}Y^{2}+2ikT_{12}^{*(2)}Y^{2}\Pi+2\Lambda\Omega
Y\biggr{)}\sim\frac{Y^{2}}{X},$ $\displaystyle\hat{\tau}_{12}^{2\ell}$
$\displaystyle=\frac{1}{X}\biggl{(}-{\textstyle{1\over
2}}ikU_{*}^{{}^{\prime\prime}}Y^{2}\hat{\tau}_{12}^{\ell}+(ikT_{11}^{*(1)}Y-T_{12}^{*(2)})\Pi
Y+{\textstyle{1\over 2}}ik\Omega T_{11}^{*(0)}Y^{2}$ $\displaystyle\hskip
170.71652pt+U_{*}^{{}^{\prime\prime}}Y\hat{\tau}_{22}^{\ell}+U_{*}^{{}^{\prime}}\hat{\tau}_{22}^{2\ell}+\frac{i\Lambda\Omega}{k}\biggr{)}\sim\frac{Y^{2}}{X},$
$\displaystyle\hat{\tau}_{11}^{2\ell}$
$\displaystyle=\frac{1}{X}\biggl{(}-3\Pi YT_{11}^{*(1)}-2T_{11}^{*(0)}\Omega
Y+\frac{2iT_{12}^{*(1)}\Omega}{k}+2U_{*}^{{}^{\prime\prime}}Y\hat{\tau}_{12}^{\ell}+2U_{*}^{{}^{\prime}}\hat{\tau}_{12}^{2\ell}\biggr{)}\sim
O\biggl{(}\frac{Y^{2}}{X^{2}},\frac{Y}{X}\biggr{)}$

as $|Y|\rightarrow\infty$, the RHS of (89) generates at best constant terms as
$Y\rightarrow\pm\infty$ which can be removed by the ‘const’. Hence there is no
consequence outside the critical layer as expected (the $\log\varepsilon$
ordered terms are there just to fix up the logarithmic terms). Note, however,
$\hat{v}_{2\ell}$ is non trivial in the critical layer but it is just doesn’t
drive a jump across it.

###  4.7 $O(1/\varepsilon)$ in the Stokes equation

This is the order at which a finite jump in $D^{3}v$ across the critical layer
is determined. The $O(1/\varepsilon)$ Stokes equation balance integrated once
gives

$\hat{D}^{3}\hat{v}_{3}=\underbrace{2k^{2}\hat{D}\hat{v}_{1}}_{(a)}+\Lambda\hat{D}^{3}\hat{v}_{2}-\underbrace{k^{2}\hat{\tau}_{11}^{2}}_{(d)}+k^{2}\hat{\tau}_{22}^{0}+\overbrace{ik\hat{D}\hat{\tau}_{12}^{2}}^{(c)}+\overbrace{ik^{3}\int\hat{\tau}_{12}^{0}dY}^{(b)}$
(82)

Since we are only interested in jumps across the critical layer, we focus on
the logarithmic terms appearing on the RHS of (82) (only the labelled terms
contribute). Terms (a) and (b) follow immediately

$\displaystyle(a)$
$\displaystyle\sim-\frac{2ik^{3}T_{11}^{*(0)}\hat{v}_{0}}{U_{*}^{{}^{\prime}}}\log
X\quad({\rm from}\,\,(\ref{Dv_lead}))$ $\displaystyle(b)$
$\displaystyle\sim\frac{ik^{3}T_{11}^{*(0)}\hat{v}_{0}}{U_{*}^{{}^{\prime}}}\log
X\quad({\rm from}\,\,(\ref{tau_12^0}))$ (83)

as $|Y|\rightarrow\infty$. Terms (c) and (d) require more work going to yet
higher order in the polymer stress equations (12)-(14). Starting with (c),

$\displaystyle ik\hat{D}\hat{\tau}_{12}^{2}$
$\displaystyle=\overbrace{\hat{D}\biggl{[}\frac{{\textstyle{1\over
2}}k^{2}U_{*}^{{}^{\prime\prime}}Y^{2}\hat{\tau}_{12}^{*(1)}}{X}\biggr{]}}^{(i)}+\overbrace{\hat{D}\biggl{[}\frac{(-k^{2}T_{11}^{*(1)}Y-ikT_{12}^{*(2)})\hat{v}_{1}}{X}\biggr{]}}^{(ii)}+\overbrace{\hat{D}\biggl{[}-\frac{k^{2}T_{11}^{*(0)}\hat{v}_{2}}{X}\biggr{]}}^{(iii)}\biggr{.}$
$\displaystyle\hskip
42.67912pt+\hat{D}\biggl{[}\frac{ikU_{*}^{{}^{\prime}}\hat{\tau}_{22}^{2}}{X}+\frac{-k^{2}(\Lambda+Y^{2}T_{11}^{*(2)})\hat{v}_{0}}{X}+\frac{ikU_{*}^{{}^{\prime\prime}}Y\hat{\tau}_{22}^{1}}{X}-\frac{\Lambda\hat{D}^{2}\hat{v}_{2}}{X}\biggr{]}$
(84)

Then the labelled terms have the following logarithmic behaviour

$\displaystyle(i)$
$\displaystyle\sim-\frac{k^{2}U_{*}^{{}^{\prime\prime}}[T_{11}^{*(0)}]^{2}\hat{v}_{0}}{2U_{*}^{{}^{\prime}3}}\log
X$ $\displaystyle(ii)$
$\displaystyle\sim\frac{k^{2}T_{11}^{*(0)}T_{11}^{*(1)}\hat{v}_{0}}{U_{*}^{{}^{\prime}2}}\log
X$ $\displaystyle(iii)$
$\displaystyle\sim\frac{ikT_{11}^{*(0)}\Omega}{2U_{*}^{{}^{\prime}}}\log X$
(85)

as $|Y|\rightarrow\infty$. Now considering (d),

$\displaystyle-k^{2}\hat{\tau}_{11}^{2}$
$\displaystyle=\frac{{\textstyle{1\over
2}}ik^{3}U_{*}^{{}^{\prime\prime}}Y^{2}\hat{\tau}_{11}^{1}}{X}+\frac{2k^{2}T_{11}^{*(2)}Y\hat{v}_{0}}{X}+\overbrace{\frac{k^{2}T_{11}^{*(1)}\hat{v}_{1}}{X}}^{(a)}+\overbrace{\frac{k^{2}T_{11}^{*(0)}\hat{D}\hat{v}_{2}}{X}}^{(b)}+\overbrace{\frac{k^{2}T_{11}^{*(1)}Y\hat{D}\hat{v}_{1}}{X}}^{(c)}$
$\displaystyle\hskip
14.22636pt-2\frac{ik^{2}T_{12}^{*(1)}\hat{D}^{2}\hat{v}_{2}}{X}-2\frac{ik^{2}T_{12}^{*(2)}Y\hat{D}^{2}\hat{v}_{1}}{X}+\underbrace{\frac{-2k^{2}U_{*}^{{}^{\prime}}\hat{\tau}_{12}^{2}}{X}}_{(d)}+\underbrace{\frac{-2k^{2}U_{*}^{{}^{\prime\prime}}Y\hat{\tau}_{12}^{1}}{X}}_{(e)}$
(86)

with labelled terms having the following logarithmic behaviour

$\displaystyle(a)$
$\displaystyle\sim-\frac{k^{2}T_{11}^{*(0)}T_{11}^{*(1)}\hat{v}_{0}}{U_{*}^{`2}}\log
X$ $\displaystyle(b)$
$\displaystyle\sim-\frac{2ikT_{11}^{*(0)}\Omega}{U_{*}^{{}^{\prime}}}\log X$
$\displaystyle(c)$
$\displaystyle\sim-\frac{2k^{2}T_{11}^{*(0)}T_{11}^{*(1)}\hat{v}_{0}}{U_{*}^{{}^{\prime}2}}\log
X$ $\displaystyle(d)$
$\displaystyle\sim\biggl{[}-\frac{k^{2}U_{*}^{{}^{\prime\prime}}[T_{11}^{*(0)}]^{2}\hat{v}_{0}}{U_{*}^{{}^{\prime}3}}+\frac{2k^{2}T_{11}^{*(0)}T_{11}^{*(1)}\hat{v}_{0}}{U_{*}^{{}^{\prime}2}}+\frac{ikT_{11}^{*(0)}\Omega}{U_{*}^{{}^{\prime}}}\biggr{]}\log
X$ $\displaystyle(e)$
$\displaystyle\sim\frac{2k^{2}U_{*}^{{}^{\prime\prime}}[T_{11}^{*(0)}]^{2}\hat{v}_{0}}{U_{*}^{{}^{\prime}3}}\log
X$ (87)

as $|Y|\rightarrow\infty$. Bringing the expressions from (83), (85) & (87)
together,

$\hat{D}^{3}\hat{v}_{3}\sim\Upsilon\log
X:=\biggl{[}k\Pi(k-2i\Lambda)-ik\Lambda U_{*}^{{}^{\prime}}\Omega\biggr{]}\log
X\quad{\rm as}\,\,Y\rightarrow\pm\infty$ (88)

which gives a finite jump in $D^{3}v$ across the critical layer. Integrating
repeatedly gives

$\displaystyle\hat{D}^{2}\hat{v}_{3}$
$\displaystyle=-\frac{i\Upsilon}{kU_{*}^{{}^{\prime}}}\biggl{(}X\log
X-X\biggr{)}+\alpha_{3},$ $\displaystyle\hat{D}\hat{v}_{3}$
$\displaystyle=-\frac{\Upsilon}{k^{2}U_{*}^{{}^{\prime}2}}\biggl{(}\frac{1}{2}X^{2}\log
X-\frac{3}{4}X^{2}\biggr{)}+\alpha_{3}Y+\beta_{3},$ $\displaystyle\hat{v}_{3}$
$\displaystyle=\frac{i\Upsilon}{k^{3}U_{*}^{{}^{\prime}3}}\biggl{(}\frac{1}{6}X^{3}\log
X-\frac{11}{36}X^{3}\biggr{)}+\frac{1}{2}\alpha_{3}Y^{2}+\beta_{2}Y+\gamma_{3}.$
(89)

where $\alpha_{3},\beta_{3}$ and $\gamma_{3}$ are constants with $\gamma_{3}$
set by ensuring that $\hat{v}_{3}=0$ at $Y=0\,(X=1)$).

Figure 5:  The ‘error’ in $c$ and $\Lambda$ as a function of finite $\delta$
compared to their $\delta=0$ limiting values for the asymptotically matched
solution and for the numerical solution as a function of $1/W$ compared to the
$W\rightarrow\infty$ limit for $k=1.1$ and $(c,\Lambda)=(0.999218,4.1101)$.
Circles are used for the asymptotically matched solutions (filled red circles
for $c$ and open blue circles for $\Lambda$); for the numerical
approximations, red $+$ is used for $c$ and blue filled squares for $\Lambda$.
The ‘truth’ $(c(0),\Lambda(0))$ is estimated by assuming
$c(0)=c(\delta)+a\delta+b\delta^{2}+\ldots$ and eliminating the leading error
between the two most accurate predictions - i.e.
$c(0):=(\,10c(10^{-7})-c(10^{-6})\,)/9$ (Note $1/W\leq 1/974$ for an Oldroyd-B
fluid (Khalid et al. (2021b))).

Figure 6:  The outer solution at $k=1.1$ and $\Lambda\approx 4.11=W(1-\beta)$
found by applying the matching conditions (96)-(101) with $\delta=2\times
10^{-4}$ (full plot) and $\delta=5\times 10^{-6}$ (partial right plot). As in
figure 1, $v$ is blue (real/imaginary parts solid/dashed respectively), real
part of $Dv$ is red and imaginary part is black. Notice the increasingly
singular behaviour in the imaginary part of $Dv$ at the critical point as
$\delta$ decreases.

Figure 7:  A comparison of the full stress field at $k=1.1$, $\Lambda\approx
4.11$ computed numerically at $W=32,000$ (in blue both left and right columns)
with the outer solutions (left column) and leading inner asymptotic solutions
given by (59)-(61) (right column - note the rescaling of both axes), both
marked in red. Solid lines indicate real parts and dashed lines imaginary
parts.) Note there is no $\tau_{22}$ component in the leading order outer
solution hence the absence of a red line in the top left plot.

###  4.8 Matching conditions across the critical layer

All the work above has built up a high order (in $\varepsilon$) representation
for $\hat{v}_{0}$ as follows

$\displaystyle v$
$\displaystyle=\hat{v}_{0}+\varepsilon\,\biggl{[}\hat{v}_{1}(Y)+\Pi
Y\log\varepsilon\biggr{]}+\varepsilon^{2}\,\biggl{[}\hat{v}_{2}(Y)+\tfrac{1}{2}\Omega
Y^{2}\
\log\varepsilon\biggr{]}+\varepsilon^{3}\log\varepsilon\,\hat{v}_{2\ell}$
$\displaystyle\hskip
199.16928pt+\varepsilon^{3}\biggl{[}\,\hat{v}_{3}(Y)+\tfrac{1}{6}\Upsilon
Y^{3}\,\biggr{]}+\ldots.$ (90)

What’s important is the behaviour as $y\rightarrow y^{*}$ for matching to the
outer solutions; specifically

$\displaystyle v$
$\displaystyle=\hat{v}_{0}+\biggl{[}\Pi(y-y^{*})+\frac{1}{2}\Omega(y-y^{*})^{2}+\frac{1}{6}\Upsilon(y-y^{*})^{3}+\ldots\biggr{]}\log[ikU_{*}^{{}^{\prime}}(y-y^{*})]$
$\displaystyle\hskip
56.9055pt+\Pi(y-y^{*})(\alpha_{1}-1)-\frac{3}{4}\Omega(y-y^{*})^{2}-\frac{11}{36}\Upsilon(y-y^{*})^{3}+O(\varepsilon,(y-y^{*})^{4})$
(91)

where all the constants disappear at leading order in $\varepsilon$ with the
exception of $\alpha_{1}$. From this, the various jump conditions

$[A]^{+}_{-}:=A(y^{*}+\delta)-A(y^{*}-\delta),$ (92)

where $\delta\rightarrow 0$, can be deduced as

$\displaystyle[v]^{+}_{-}$
$\displaystyle=2\Pi\delta[\,\log|kU_{*}^{{}^{\prime}}\delta|+\alpha_{1}-1\,]+\tfrac{1}{2}i\pi\Omega\delta^{2}+\Upsilon\delta^{3}[\,\tfrac{1}{3}\log|kU_{*}^{{}^{\prime}}\delta|-\tfrac{11}{18}\,],$
(93) $\displaystyle[Dv]^{+}_{-}$
$\displaystyle=i\pi\Pi+2\delta\Omega[\,\log|kU_{*}^{{}^{\prime}}\delta|-1\,]+{\textstyle{1\over
2}}i\pi\Upsilon\delta^{2}$ (94) $\displaystyle[D^{2}v]^{+}_{-}$
$\displaystyle=\frac{2\Pi}{\delta}+i\pi\Omega+2\delta\Upsilon[\,\log|kU_{*}^{{}^{\prime}}\delta|-1\,],$
(95) $\displaystyle[D^{3}v]^{+}_{-}$
$\displaystyle=\frac{2\Omega}{\delta}+i\pi\Upsilon.$ (96)

with

$\Pi=-2ik\Lambda U_{*}^{{}^{\prime}}\hat{v}_{0},\quad\Omega=(\,4ik\Lambda
U_{*}^{{}^{\prime\prime}}-4k^{2}\Lambda^{2}U_{*}^{{}^{\prime}2}\,)\hat{v}_{0},\quad\Upsilon=k\Pi(k-2i\Lambda)-ik\Lambda
U_{*}^{{}^{\prime}}\Omega.$ (97)

The conditions (93)-(96) are correct to 3 orders in $\delta$. This is clear in
all but the last jump condition. Here $D^{3}v\sim-\Pi/\delta^{2}$ which, since
it is even in $\delta$, does not appear in the required jump.

The interior problem (41) is 4th order with 2 boundary conditions at each
boundary (non-slip at $y=-1$ and symmetry conditions at $y=0$). So matching
the interior solutions across the critical layer requires 4 (complex)
conditions to determine 4 complex constants. Since the problem is linear and
so the amplitude and phase are indeterminate, one of these conditions instead
sets the (complex) frequency $c$ for given (real) $\Lambda$ (in the case of
the neutral curve the determination for $\Im m(c)=c_{i}$ is replaced by that
for $\Lambda$ to make $c_{i}=0$). However, there is an extra (complex) unknown
in the jump conditions (93)-(96), $\alpha_{1}$, which means a 5th (complex)
condition is needed.

In practice, it is convenient to choose $\hat{v}_{0}=-1$ to set the amplitude
and phase of the eigenfunction thereby upping the matching requirement to 6
complex equations (for the 4 complex constants specifying the two outer
solutions, $c$ and $\alpha_{1}$). An extra equation comes from now being able
to impose the behaviour of $v$

$\displaystyle v^{-}+\delta\Pi\alpha_{1}$
$\displaystyle=-1+(-\delta\Pi+\tfrac{1}{2}\delta^{2}\Omega-\tfrac{1}{6}\delta^{3}\Upsilon)[\,\log|kU_{*}^{{}^{\prime}}\delta|-{\textstyle{1\over
2}}i\pi\,]+\delta\Pi-\tfrac{3}{4}\delta^{2}\Omega+\tfrac{11}{36}\delta^{3}\Upsilon,$
(98) $\displaystyle v^{+}-\delta\Pi\alpha_{1}$
$\displaystyle=-1+(\delta\Pi+\tfrac{1}{2}\delta^{2}\Omega+\tfrac{1}{6}\delta^{3}\Upsilon)[\,\log|kU_{*}^{{}^{\prime}}\delta|+{\textstyle{1\over
2}}i\pi\,]-\delta\Pi-\tfrac{3}{4}\delta^{2}\Omega-\tfrac{11}{36}\delta^{3}\Upsilon$
(99)

as the critical layer is approached from either side
($v^{\pm}=v(y^{*}\pm\delta)$). The final extra equation comes from also doing
the same for $Dv$

$\displaystyle Dv^{-}-\Pi\alpha_{1}$
$\displaystyle=(\Pi-\delta\Omega+{\textstyle{1\over
2}}\delta^{2}\Upsilon)[\,\log|kU_{*}^{{}^{\prime}}\delta|-{\textstyle{1\over
2}}i\pi\,]+\delta\Omega-\tfrac{3}{4}\delta^{2}\Upsilon,$ (100) $\displaystyle
Dv^{+}-\Pi\alpha_{1}$ $\displaystyle=(\Pi+\delta\Omega+{\textstyle{1\over
2}}\delta^{2}\Upsilon)[\,\log|kU_{*}^{{}^{\prime}}\delta|+{\textstyle{1\over
2}}i\pi\,]-\delta\Omega-\tfrac{3}{4}\delta^{2}\Upsilon.$ (101)

which crucially does not introduce any new unknown constants. These 4
conditions together with the jump conditions (95) and (96) give the required 6
conditions for the 4 unknowns specifying the interior solution below and above
the critical layer, the two real numbers, $c_{r}$ and $\Lambda$, and the
complex constant $\alpha_{1}$.

| $\delta$ |  | $c_{r}$ | $k\Lambda_{-}$ |  | $c_{r}$ | $k\Lambda_{+}$
---|---|---|---|---|---|---|---
$k=0.1$ | $32\times 10^{-4}$ |  | 0.99855267713 | 3.446546 |  |  |
| $8\times 10^{-4}$ |  | 0.99844946362 | 3.400230 |  |  |
| $2\times 10^{-4}$ |  | 0.99842217083 | 3.389936 |  |  |
| $*1\times 10^{-4}$ |  | 0.99841758893 | 3.388332 |  | 0.99995720674 | 9.856305
| $*1\times 10^{-5}$ |  | 0.99841347771 | 3.386965 |  | 0.99995700930 | 9.823193
| $*1\times 10^{-6}$ |  | 0.99841306978 | 3.386842 |  | 0.99995698930 | 9.819909
| $*3\times 10^{-7}$ |  | 0.99841304032 | 3.386833 |  | 0.99995698769 | 9.819643
| 0 |  | 0.99841302769 | 3.386829 |  | 0.99995698700 | 9.819529
$k=1.1$ | $32\times 10^{-4}$ |  | 0.99923583992 | 4.525599 |  | 0.99988318362 | 8.328840
| $8\times 10^{-4}$ |  | 0.99922208240 | 4.519946 |  | 0.99986369555 | 7.823761
| $2\times 10^{-4}$ |  | 0.99921892090 | 4.520468 |  | 0.99985801319 | 7.706180
| $*1\times 10^{-4}$ |  | 0.99921842343 | 4.520720 |  | 0.99985702092 | 7.686784
| $*1\times 10^{-5}$ |  | 0.99921799797 | 4.521053 |  | 0.99985611212 | 7.669353
| $*1\times 10^{-6}$ |  | 0.99921795926 | 4.521106 |  | 0.99985602002 | 7.667615
| $*1\times 10^{-7}$ |  | 0.99921795576 | 4.521113 |  | 0.99985600515 | 7.667345
| 0 |  | 0.99921795537 | 4.521114 |  | 0.99985600350 | 7.667315
Table 4:  $c_{r}$ and $k\Lambda$ found by asymptotic matching on the neutral
curve for $k=0.1$ and $1.1$ with various choices of $\delta$ ($\pm$ indicates
upper and lower parts of the neutral curve and * results computed using
quadruple rather than double precision). The quadruple precision calculations
were done using the multi-precision extension package ‘Advanpix’ for Matlab.
The ‘$0$’ predictions are found by applying Richardson extrapolation to the
last two entries assuming a leading $O(\delta)$ error. All four cases compare
well with the numerical predictions marked by ‘$\infty$’ in Table 3.

In practice, we impose the conditions (96) - (101) to specify the velocity
fields and $\alpha_{1}$, and then find (real) $c$ and $\Lambda$ using Newton-
Raphson on the remaining $D^{2}v$ jump condition (95). The matching is not
surprisingly quite delicate because up to the third derivative in $v$ has to
be matched, there is singular behaviour and the critical layer can get very
close to the centreline which imposes the condition $\delta\ll\sqrt{1-c}$.
Invariably $\delta$ needs to be smaller than $10^{-4}$ to see convergence
which, since terms up to $O(\delta^{3})$ need to be resolved, requires
quadruple precision arithmetic. Table 4 show the results of matching at
$k=0.1$, which is away from the neutral curve nose (see inset A of figure 2 of
Khalid et al. (2021b)), and $k=1.1$ which is close to it. For the upper
neutral curve at $k=0.1$ where $c=0.99957$, the critical layer is so close to
the centreline that only matching with quadruple precision is possible.

Figure 5 plots the error in estimating the (lower neutral curve) limiting
values of $c$ and $\Lambda$ at $k=1.1$ via the numerical solution taking
$W\rightarrow\infty$ and the asymptotic matching approach taking
$\delta\rightarrow 0$. The asymptotically-matched outer velocity fields over
$y\in[-1,y^{*}-\delta]\cup[y^{*}+\delta,0]$ compare excellently with the full
numerical solution shown in figure 4 \- see figure 6. Reducing $\delta$ from
$2\times 10^{-4}$ to $5\times 10^{-6}$ shows the singularity in $\Im m(Dv)$ at
the critical layer when the phase of the eigenfunction is set by $Dv=1$ at the
midplane. The numerically-computed stress field at $W=32,000$ is compared in
figure 7 with the matched outer solution (using $\delta=2\times 10^{-4}$) and
the leading inner asymptotic solution, again with excellent agreement.

The key realisation from the asymptotic analysis is that the structure of the
critical layer is built upon a non-vanishing cross-stream velocity there. This
is reflected in the fact that everything scales with it - specifically
$\hat{v}_{0}$ in the expansion (43). For example, the 3 matching constants in
(97) are proportional to $\hat{v}_{0}$: if these are zero, the critical layer
has no effect on the outer solutions. The cross-stream velocity, however, must
vanish at the midplane by the symmetry conditions and so there has to be an
$O(1)$ outer layer between the critical layer and the centreline to bring
about this adjustment (derivatives in the outer regions are $O(1)$ and
$\hat{v}_{0}=O(1)$ to set the normalisation of the eigenfunction). This
explains why the critical layer cannot approach the centreline as
$W\rightarrow\infty$ or equivalently why $c$ converges to a finite value which
is not 1. It remains unclear why this finite value is numerically so close to
1 (e.g. $1-c=4.3\times 10^{-5}$ for $k=0.1$ in Table 3) but the plausible
hypothesis is that the critical layer can only manifest in a low shear region
compared to the rest of the domain. This is probed a little in §6 below but
first we give a discussion on the instability mechanism.

Figure 8:  Left: $P(y)$ plotted against $y$ for a solution ($k=1.1$,
$\Lambda=4.520468$ and $c=0.99921892$) matched across the critical layer using
$\delta=2\times 10^{-4}$ (the position of the critical layer is shown as the
vertical dashed black line). Right: $C(y)$ (green) and $D(y)$ (black dashed)
plotted against $y$ for the same matched solution. These figures show that the
velocity field is driven locally by the critical layer whereas the polymer
perturbation is driven globally by the gradual relaxation of the streamline
distortion caused by the critical layer.

##  5 Instability Mechanism

The asymptotic analysis above separates the ‘inner’ solution in the critical
layer from the ‘outer’ solution, allowing scrutinization of how the
instability manifests in the latter. The role of the critical layer is then
viewed as generating ‘energising’ internal conditions for the outer solution
(or boundary conditions for the two parts of the outer solution). Before
proceeding in this manner, we simplify the outer equations by rewriting some
terms using the streamline displacement $\phi:=v/ik(U-c)$ instead of $v$
following Rallison & Hinch (1995) to get

$\displaystyle 0$ $\displaystyle=-ikp+(D^{2}-k^{2})u-ikT_{11}D\phi,$ (102)
$\displaystyle 0$ $\displaystyle=-Dp+(D^{2}-k^{2})v-k^{2}T_{11}\phi,$ (103)
$\displaystyle 0$ $\displaystyle=iku+Dv,$ (104) $\displaystyle\tau_{11}$
$\displaystyle=-2T_{11}D\phi-\phi DT_{11},$ (105) $\displaystyle\tau_{12}$
$\displaystyle=ikT_{11}\phi.$ (106)

While the above analysis shows that $v$ is continuous across the critical
layer, $\phi\sim 1/(U-c)$ as the critical layer is approached. Hence the
streamline displacement is maximal there and the question is how this drives
the polymer field which must in turn offset the viscous dissipation in the
inertialess momentum equation.

Approaching the neutral curve, $c_{i}\rightarrow 0$, means that $\phi$ is
exactly out of phase with $v$ so that $-k^{2}T_{11}\phi$ can do no work in
energising $v$ in (103), that is,

$\frac{k}{2\pi}\int^{2\pi/k}_{0}\Re e[ve^{ik(x-ct)}]\Re e[-k^{2}T_{11}\phi
e^{ik(x-ct)}]\,dx=0.$ (107)

This means the viscous dissipation in the cross-stream variable $v$ can only
be offset by the pressure term and the polymer stress driving of the velocity
field must occur in (102). To confirm this, the power input

$P(y):=\frac{k}{2\pi}\int^{2\pi/k}_{0}\Re e[ue^{ik(x-ct)}]\Re e[-ikT_{11}D\phi
e^{ik(x-ct)}]\,dx$ (108)

is plotted in figure 8 which clearly shows that the polymer stress energising
of $u$ is localised at the critical layer. The specific solution used here is
shown in figure 9.

Figure 9:  Contour plots of the outer solution corresponding to $k=1.1$,
$\Lambda=4.520468$ and $c=0.99921892$ matched across the critical layer using
$\delta=2\times 10^{-4}$. Top left $u$, top right $v$, bottom left $\tau_{11}$
and bottom right $\phi$. The last plot focusses on the area around the
critical layer (the white gap is just the excluded critical layer region
$[y^{*}-\delta,y^{*}+\delta]$). For $u$ and $v$, 20 equally spaced contour
levels between the maximum and minimum values are used. For $\tau_{11}$, the
contour levels are $\pm[0,10,20,30,40,50,60,100,200,1000,2000,5000]$ and 50
equally spaced levels are used for $\phi$ with all the contours essentially
contained in the shown area around the critical layer..

The polymer stress equations are particularly simple when written in terms of
the streamline displacement and interpretable. The first term on the right of
(105) represents the increase in the streamwise-normal stress due to
streamline compression (in $y$), the second term represents the maintenance of
the initial basic stress on displaced streamlines, and the third term (the the
RHS of (106) ) is simply the generation of tangential polymer stress due to
tilting of the base polymer stress lines (Rallison & Hinch, 1995). Multiplying
(105) by $ik(U-c)$ recovers the time-derivative and advection terms on the LHS
of (105) and thereby recovers the proper driving terms on the RHS,

$ik(U-c)\tau_{11}=-2ik(U-c)T_{11}D\phi-ik(U-c)\phi DT_{11}.$ (109)

Their energising effects,

$\displaystyle C(y)$ $\displaystyle:=\frac{k}{2\pi}\int^{2\pi/k}_{0}\Re
e[\tau_{11}e^{ik(x-ct)}]\Re e[-2ik(U-c)T_{11}D\phi e^{ik(x-ct)}]\,dx$ (110)
$\displaystyle\&\quad D(y)$ $\displaystyle:=\frac{k}{2\pi}\int^{2\pi/k}_{0}\Re
e[\tau_{11}e^{ik(x-ct)}]\Re e[-ik(U-c)\phi DT_{11}e^{ik(x-ct)}]\,dx$ (111)

are shown in figure (8). Both terms are global and barely register the
critical layer with streamline $C$ompression causing the streamwise-normal
stress to increase ($C(y)$) while $D$isplacement across the base shear field
($D(y)$) works negatively to balance it on the neutral curve. It is worth
remarking that the reintroduction of the factor $(U-c)$ into (109) is
responsible for desensitizing $C(y)$ to the critical layer as otherwise a
significant component - $T_{11}D\phi$ \- is the same as that in $P$.

The mechanism of the instability can therefore be seen as one in which the
critical layer acts like a pair of ‘bellows’ periodically sucking the flow
streamlines together - see $\phi(\pi/2,y)$ in figure 9 \- and then blowing
them apart - see $\phi(3\pi/2,y)$ same figure. This amplifies the base
streamwise-normal stress field which in turn exerts a streamwise stress on the
flow locally at the critical layer. The streamwise flow drives the cross-
stream flow through continuity which then intensifies the critical layer
closing the loop.

The one outstanding question is why the critical layer has to be so close to
the centreline as $W\rightarrow\infty$. The asymptotic analysis above
indicates that the shear at the critical layer needs to be $O(W^{0})$ as
$W\rightarrow\infty$ but can tell us nothing about the size of this $O(W^{0})$
number relative to 1. In particular, given the complicated analysis, it is
still not clear why the instability does not manifest in plane Couette flow.
To help answer this, we conduct some simple experiments in the next section.

##  6 Moving towards plane Couette flow

The complexity of the matching analysis means it is difficult to discern the
importance of $U_{*}^{{}^{\prime\prime}}$ or the size of $U_{*}^{{}^{\prime}}$
(e.g. just look at the 3 constants that emerge in (97) ) despite our best
efforts to keep them separated. So here, we perform a series of numerical
experiments exploring the effect of small changes which would bring the
channel flow closer to plane Couette flow (pCf). To keep things manageable,
these experiments concentrate on studying how the lowest $W$ point on the
neutral curve, $W_{min}$, varies as the problem is changed slightly. Four
experiments are undertaken in which this minimum point is tracked as a
homotopy parameter $\lambda$ is reduced from 1 which is the channel flow
problem studied above where $W_{min}=974$. These are as follows.

  * •

Expt. 1 explores the importance of $U^{{}^{\prime\prime}}$ by (artificially)
changing it while keeping $U^{{}^{\prime}}$ fixed. Specifically
$U^{{}^{\prime\prime}}:=-2\lambda$ everywhere so $U^{{}^{\prime\prime}}$ can
be reduced without changing $U=1-y^{2}$ or $U^{{}^{\prime}}=-2y$ in the code.

  * •

Expt 2 explores the effect of changing the boundary condition at the midplane
towards a solid boundary to mimick the monotonic increase in $U$ across the
domain of pCf. The boundary condition at $y=0$ is set to $\lambda
D^{2}v+(1-\lambda)Dv=0$ so $\lambda=1$ corresponds to the stress-free/symmetry
conditions considered above and $\lambda=0$ to a non-slip solid wall.

  * •

Expt. 3 explores the effect of increasing the minimum shear across the domain
$y\in[-1,0]$. The base flow is set to $U:=\lambda(1-y^{2})+(1-\lambda)y$ which
mixes in pCf in a way to gradually generate a (minimum) non-zero shear
$=1-\lambda$ at the midplane.

  * •

Expt 4 explores the effect of moving the $U^{{}^{\prime\prime}}=0$ point into
the interior. The base flow is set to $U:=\lambda(1-y^{2})-(1-\lambda)y$ which
mixes in pCf in a way to move the zero-shear point at
$y=-(1-\lambda)/2\lambda$ away from the midplane and into the interior
$y\in[-1,0]$.

The results are shown in figure 10. The effect of changing the boundary
conditions (blue dashed line) at the midplane is minimal and reducing
$U^{{}^{\prime\prime}}$ (black dash-dot line) is destabilizing. The presence
of vanishing shear, however, seems crucial: removing it (yellow line) quickly
stabilizes the instability whereas moving it from the midplane (red line) is
destabilizing. The plausible conclusion is that the instability needs small
shear relative to the rest of the domain to nucleate. This is certainly absent
in pCf.

Figure 10: Left: the effect on $W_{min}$ of decreasing $\lambda$ from 1 to
0.96 (the common point on $\lambda=1$ is the undisturbed value of
$W_{min}=974$): Expt. 1 black dash-dot line; Expt. 2 blue dashed line; Expt. 3
yellow solid line; Expt. 4 red solid line. Right: the base flow profile
$U:=\lambda(1-y^{2})+\lambda y$ in Expt. 3 (yellow line) and
$U:=\lambda(1-y^{2})-\lambda y$ in Expt. 4 (red line) for $\lambda=0.96$
showing how the former eliminates the zero-shear point whereas the latter
moves it into the interior (notice only a small part of the flow domain is
being shown close to the midplane. In all cases the problem is solved over
$y\in[-1,0]$)

Figure 11:  Top left: Viscoelastic Kolmogorov base flow as given in (112) in
green and channel flow $(1-y^{2})$ in dashed black. Neutral vKf eigenfunctions
for $k=1$ and $(W,\beta,c)=(2000,0.99954,0.98029)$ (top right);
$(W,\beta,c)=(200,0.99530,0.97897)$(bottom left) and
$(W,\beta,c)=(20,0.93972,0.96726)$(bottom right). $v$ is blue (real/imaginary
parts solid/dashed respectively), real part of $Dv$ is red and imaginary part
is black.

##  7 Viscoelastic Kolmogorov flow

Finally, the similarity of the base flow shape in viscoelastic Kolmogorov flow
(vKf) to that in channel flow suggests that the viscoelastic linear
instability found there (Boffetta et al., 2005) could be the centre mode
instability of Garg et al. (2018); Khalid et al. (2021a, b). It is
straightforward to confirm this by renormalising the base flow in vKf to the
form

$U=\tfrac{1}{2}(1+\cos\pi y)$ (112)

(to most closely match $(1-y^{2})$ \- see Figure 11) and considering
disturbances which are periodic over $y\in[-1,1]$ and have the same symmetry
(24) as the centre mode instability about $y=0$. These properties actually
imply that the disturbance satisfies stress-free boundary conditions at $y=\pm
1$ and so all the numerical codes developed for channel flow can trivially be
reapplied to vKf by just i) changing $U(y)$ and ii) imposing stress-free
boundary conditions on the perturbation at $y=-1$.

Using a shooting code which identifies the neutral curve at a given $k$ and
$W$ after a guess for $\beta$ and $c$ (based on the channel flow solution)
quickly identifies a neutral vKf eigenfunction for $(W,k)=(2000,1)$ at
$(\beta,c)=(0.99953538,0.98029137)$. Reducing $W$ to $200$, keeping $k=1$ and
using the values of $\beta$ and $c$ found at $W=2000$ as initial guesses, the
neutral eigenfunction at $W=200$ converges easily to
$(\beta,c)=(0.99530288,0.97896974)$. Repeating this procedure, reducing $W$
from $200$ to $20$, again converges smoothly to
$(\beta,c)=(0.93972375,0.96725546)$: see figure 11. The similarity of the
neutral eigenfunctions in figure 11 to those in channel flow (modulo the
different boundary conditions at $y=-1$) is striking and suggests that
Boffetta et al. (2005) were actually the first to find the centre mode
instability in rectilinear viscoelastic flow. To further support this
conclusion, Berti et al. (2008); Berti & Boffetta (2010) see ‘arrowhead’
solutions when tracking their vKf instability to finite amplitude (e.g. see
figures 7 and 8 in Berti & Boffetta (2010)) just like those found in channel
flow originating from the centre mode instability Page et al. (2020); Buza et
al. (2022a); Beneitez et al. (2023a).

Lastly, it is also worth remarking that $W$ can be reduced down to at least
$2.1$ at $k=1$ where $(\beta,c)=0.035543817,0.82144631)$ in vKf. These $W$ and
$\beta$ are an order of magnitude smaller than those required for centre mode
instability in channel flow presumably because of the rigid boundaries
present.

##  8 Discussion

We first summarise the findings of the paper. The first part of these concern
the $Re\rightarrow\infty$ asymptotics of the upper (§3.1) and lower branches
(§3.2) of the centre-mode neutral curve in the $Re$-$W$ plane for viscoelastic
channel flow.

Along the upper branch

$W\sim Re^{1/3},\quad k\sim Re^{1/3},\quad 1-c\sim Re^{-2/3}$ (113)

with numerical coefficients given in (26) for $\beta=0.9$. These scalings are
equivalent to $Re\sim O(E^{-3/2}$), $k\sim O(E^{-1/2})$ and $1-c\sim O(E)$ as
the elasticity number $E:=W/Re\rightarrow 0$ consistent with figure 11 in
Khalid et al. (2021a).

Along the lower branch,

$W\sim Re,\quad k\sim\frac{1}{Re},\quad c=O(1)$ (114)

with numerical coefficients computed for $\beta=0.9,0.98$ and $0.994$ given in
Table 2. These lower branch scalings are apparent in figure 13 of Khalid et
al. (2021a) (see also their figure 18).

The second part of the findings described in §4 concern the inertialess limit
of viscoelastic channel flow. By $\beta=0.994$ as $\beta$ increases, the lower
branch has swung sufficiently clockwise in the $Re$-$W$ plane to cross the
$Re=0$ axis (see figure 1). This reveals the existence of an inertialess
($Re=0$) centre mode instability and the $W\rightarrow\infty$ asymptotic
problem in the ultra-dilute limit where $W(1-\beta)=O(1)$ is then treated. A
matched asymptotic analysis is performed in which a critical layer region is
resolved sufficiently to extract matching conditions, (93)-(96), to connect up
outer regions either side. Interestingly, the outer problem is 4th order as
opposed to the usual 2nd order problem for the Orr-Sommerfeld problem and so
requires matching conditions all the way down to the 3rd order derivative in
the cross-stream velocity. This leads to a particularly delicate matching
procedure (§4.8) where the matching conditions need to be resolved to third
order in the small matching parameter and quadruple precision is needed to
make contact with numerical solutions. The completed analysis is successful in
revealing, again unlike the Orr-Sommerfeld problem, that

$1-c=O(1)\qquad{\rm as}\quad W\rightarrow\infty$ (115)

This $O(1)$ number can be deceivingly small compared to 1 (e.g. $4.3\times
10^{-5}$ for $k=0.1$ in Table 3) but nevertheless remains finite as
$W\rightarrow\infty$. That this has to be so is clear from the structure of
the critical layer that is built around an $0(1)$ cross-stream velocity which
has to be brought to zero at the midplane by an $O(1)$ outer region separating
the critical layer from it. Quite why $c$ has to be so close to 1 or
equivalently why the critical layer sets up in a region of small shear is
unclear (and unknowable from the asymptotic analysis). Some simple numerical
experiments (§6) suggest that the lack of this small shear region is the
likely reason the instability does not manifest in plane Couette flow.

The asymptotic analysis also clarifies that the instability mechanism (§5) is
one in which the critical layer acts like a pair of ‘bellows’ periodically
sucking the flow streamlines together and then blowing them apart (see figure
9). This amplifies the streamwise-normal polymer stress field which in turn
exerts a streamwise stress on the flow locally at the critical layer. The
streamwise flow drives the cross-stream flow by continuity which then
intensifies the critical layer closing the loop.

Finally in §7, a connection is made between the centre mode instability of
channel flow and an earlier linear instability found in viscoelastic
Kolmogorov flow by Boffetta et al. (2005). The fact that the instability in
Kolmogorov flow was discovered at much lower $Re$ and $W$ to the extent it was
viewed as a purely ‘elastic’ instability disguised its connection to the work
of Garg et al. (2018). They worked at $Re=O(100)$-$O(1000)$ and $W\gtrsim 20$
in a very different geometry and so viewed their instability as ‘elasto-
inertial’ in origin. It is clear in hindsight that the apparent difference in
the regimes is more a function of the boundary conditions - a solid wall in
the pipe verses periodicity in Kolmogorov flow - than any deeper dynamical
difference as evident in the Newtonian versions of the respective problems
($Re_{crit}=O(10)$ for Kolmogorov flow while $Re_{crit}=5772$ for channel
flow).

The importance of the centre mode instability for elasto-inertial turbulence
(EIT) or indeed elastic turbulence (ET) is still an area of much current
speculation (e.g. Datta et al., 2022; Dubief et al., 2023). While computations
have confirmed that the instability leads to travelling wave solutions dubbed
‘arrowheads (Page et al., 2020; Dubief et al., 2022; Buza et al., 2022a), it
remains unclear what these lead to via their own bifurcations. Recently
Beneitez et al. (2023a) have found that the arrowhead solutions coexist with
EIT rather breaking down to it. The situation, however, is slightly clearer at
$Re=0$ (perhaps because the parameter space is one dimension less) where the
(2-dimensional) arrowhead solution can become unstable to 3-dimensional
disturbances (Lellep et al., 2023a). Very recent calculations using a large
domain indicate that this instability can lead to a 3-dimensional chaotic
state (Lellep et al., 2023b).

Further complicating the picture is the very recent emergence of another
viscoelastic instability - dubbed ‘PDI’ for polymer diffusive instability -
when polymer stress diffusion is present (Beneitez et al., 2023b; Couchman et
al., 2023; Lewy & Kerswell, 2023). This is a ‘wall’ mode which also exists for
all $Re$ including $0$ and any shear flow with a solid wall is susceptible.
Beneitez et al. (2023b) have already found that PDI can lead to a chaotic 3D
state in inertialess plane Couette flow using the FENE-P model.

Going forward, the challenge is to try to unpick which process of the current
contenders - viscoelastic Tollmien Schlichting instability, the centre mode
instability or the PDI - triggers EIT and ET in what part of parameter space.
This will assist in simulating EIT and ET and ultimately in manipulating those
states as required for industrial applications.

Acknowledgements: The authors acknowledge EPSRC support for this work under
grant EP/V027247/1.

Declaration of interests: The authors report no conflict of interest.

##  Appendix A Numerical methods for eigenvalue problem

Two complementary approaches were developed: a matrix formulation and a
shooting technique.

###  A.1 Matrix

A generalised eigenvalue code was written to solve for all 6 variables
$(u,v,p,\tau_{11},\tau_{12},\tau_{22})$ building in the symmetry of the
unstable eigenfunction around the midline to improve efficiency. This was done
by mapping the lower half of the channel $y\in[-1,0]$ to the full Chebyshev
domain $[-1,1]$ so collocation points are concentrated at the wall and the
centreline, and imposing symmetry boundary conditions at the centreline $y=0$
which are $\partial u/\partial y=v=0$ (no b.c.s are imposed on $p$ or the
polymer stress $\tau$ anywhere). The fields are expanded using individual
functions which incorporate these conditions: specifically

$\left[\begin{array}[]{c}u\\\ v\\\ p\\\ \tau_{11}\\\ \tau_{12}\\\
\tau_{22}\end{array}\right]=e^{ik(x-ct)}\sum_{n=1}^{N}\left[\begin{array}[]{c}u_{n}\\{\,T_{n+1}(2y+1)+\alpha_{n}T_{n}(2y+1)+(\alpha_{n}-1)T_{n-1}(2y+1)\,\\}\\\
v_{n}\\{\,T_{n+1}(2y+1)-T_{n-1}(2y+1)\,\\}\\\ p_{n}T_{n-1}(2y+1)\\\
t_{n}T_{n-1}(2y+1)\\\ r_{n}T_{n-1}(2y+1)\\\
s_{n}T_{n-1}(2y+1)\end{array}\right]$

where $T_{n}(z):=\cos(n\cos^{-1}z)$ is the $n$th Chebyshev polynomial,
$(u_{n},v_{n},p_{n},t_{n},r_{n},s_{m})\in{\mathbb{C}}^{6}$ and
$\alpha_{n}:=-4n/(2n^{2}-2n+1)$ so that $u(x,-1,t)=0=\partial u/\partial
y(x,0,t)$. A complementary inverse iteration code was also written which could
take an eigenvalue from the generalised eigenvalue code and converge it at
much higher resolution (e.g. Table 1). This was important as the generalised
eigenvalue problem is not well-conditioned with increasing resolution: see the
drift in the eigenvalue for $N\gtrsim 300$ in Table 1 using eig in Matlab).
This lack of conditioning gets worse near the neutral curve where an interior
critical layer is present. Inverse iteration treats exactly the same matrices
but is far better conditioned - there is no drift in Table 1 even increasing
$N$ to 2000.

###  A.2 Shooting

Two shooting codes were also written based on different integrators. The first
used RK4 over a uniformly spaced grid (e.g. 50,000 points across [-1,0] in
Table 1) with inbuilt re-orthogonalisation of shooting solutions across the
domain and a second used Matlab’s ODE15s with relative and absolute tolerances
set at $3\times 10^{-14}$ which did not. For the solutions sought, re-
orthogonalisation was not needed and so the latter, which was more efficient
as it has locally adaptive stepping, was used for all subsequent calculations.
The eigenvalue problem is 4th order so the usual shooting code approach takes
a guess for the complex phase speed $c$ and searches for the 2 unknown
velocity boundary conditions at one wall which mean that the required boundary
conditions at the other are obeyed. This can be readily adapted to search for
the neutral curve directly by setting $c_{i}=0$ and instead adjusting one
(real) parameter of the problem. Here we chose to vary $\Lambda=(1-\beta)W$
keeping $W$ fixed. This can be used to recreate the upper inset in figure 2 of
Khalid et al. (2021b): see Tables 2 and 3 for sample computations at $k=0.1$
and $k=1.1$.

| $N$ |  | $c_{r}$ |  | $c_{i}$
---|---|---|---|---|---
Eigenvalue code | 100 |  | 0.999608051011 |  | 8.2367630455$\times 10^{-5}$
| 200 |  | 0.999608009040 |  | 8.2751793055$\times 10^{-5}$
| 300 |  | 0.999607999400 |  | 8.2756092074$\times 10^{-5}$
| 500 |  | 0.999608029589 |  | 8.2771153285$\times 10^{-5}$
| 1000 |  | 0.999608103554 |  | 8.2644316360$\times 10^{-5}$
| 2000 |  | 0.999607409613 |  | 8.3760120903$\times 10^{-5}$
Inverse iteration | 100 |  | 0.999608051169 |  | 8.2367612043$\times 10^{-5}$
| 200 |  | 0.999608007042 |  | 8.2751523639$\times 10^{-5}$
| 300 |  | 0.999608007115 |  | 8.2751700997$\times 10^{-5}$
| 500 |  | 0.999608007115 |  | 8.2751701045$\times 10^{-5}$
| 1000 |  | 0.999608007115 |  | 8.2751701039$\times 10^{-5}$
| 2000 |  | 0.999608007115 |  | 8.2751701041$\times 10^{-5}$
Shooting | ODE15s |  | 0.999608007115 |  | 8.2751701043$\times 10^{-5}$
| 50,000 |  | 0.999608007115 |  | 8.2751701040$\times 10^{-5}$
Table 5: Check of codes with the eigenvalue shown in figure 1(inset) of Khalid
et al. (2021b). The unstable centre mode at $\beta=0.997$, $k=0.75$ and
$W=2500$ is shown there with $c_{r}\approx 0.9996$ and $c_{i}\approx 8\times
10^{-5}$. Two shooting codes were written based on different integrators. The
first used RK4 over a uniformly spaced grid (here 50,000 points across [-1,0])
and the second used Matlab’s ODE15s with relative and absolute tolerances set
at $3\times 10^{-14}$.

## References

  * Beneitez et al. (2023a) Beneitez, M., Page, J., Dubief, Y. & Kerswell, R. R. 2023a Multistability of elasto-inertial two-dimensional channel flow. https://arxiv.org/abs/2308.11554 .
  * Beneitez et al. (2023b) Beneitez, M., Page, J. & Kerswell, R. R. 2023b Polymer diffusive instability leading to elastic turbulence in plane Couette flow. Phys. Rev. Fluids 8, L101901.
  * Berti et al. (2008) Berti, S., Bistagnino, A., Boffetta, G., Celani, A. & Musacchio, S. 2008 Two-dimensional elastic turbulence. Phys. Rev. E 77, 055306.
  * Berti & Boffetta (2010) Berti, S. & Boffetta, G. 2010 Elastic waves and transition to elastic turbuelnce in a two-dimensional viscoelastic Kolmogorov flow. Phys. Rev. E 82, 036314.
  * Boffetta et al. (2005) Boffetta, G., Celani, A., Mazzino, A., Puliafito, A. & Vergassola, M. 2005 The viscoelastic Kolmogorov flow: eddy viscosity and linear instability. J. Fluid Mech. 523, 161–170.
  * Buza et al. (2022a) Buza, Gergely, Beneitez, Miguel, Page, Jacob & Kerswell, Rich R 2022a Finite-amplitude elastic waves in viscoelastic channel flow from large to zero Reynolds number. Journal of Fluid Mechanics 951, A3.
  * Buza et al. (2022b) Buza, Gergely, Page, Jacob & Kerswell, Rich R 2022b Weakly nonlinear analysis of the viscoelastic instability in channel flow for finite and vanishing Reynolds numbers. Journal of Fluid Mechanics 940, A11.
  * Chaudhary et al. (2021) Chaudhary, Indresh, Garg, P., Subramanian, Ganesh & Shankar, V. 2021 Linear instability of viscoelastic pipe flow. Journal of Fluid Mechanics 908, A11.
  * Couchman et al. (2023) Couchman, M. M. P., Beneitez, M., Page, J. & Kerswell, R. R. 2023 Inertial enhancement of the polymer diffusive instability. http://arxiv.org/abs/2308.14879 .
  * Datta et al. (2022) Datta, Sujit S, Ardekani, Arezoo M, Arratia, Paulo E, Beris, Antony N, Bischofberger, Irmgard, McKinley, Gareth H, Eggers, Jens G, López-Aguilar, J Esteban, Fielding, Suzanne M, Frishman, Anna et al. 2022 Perspectives on viscoelastic flow instabilities and elastic turbulence. Physical Review Fluids 7 (8), 080701.
  * Dong & Zhang (2022) Dong, M. & Zhang, M. 2022 Asymptotic study of linear instability in a viscoelastic pipe flow. Journal of Fluid Mechanics 935, A28.
  * Drazin & Reid (1981) Drazin, P. G. & Reid, W. H. 1981 Hydrodynamic Stability .
  * Dubief et al. (2022) Dubief, Yves, Page, Jacob, Kerswell, Richard R, Terrapon, Vincent E & Steinberg, Victor 2022 First coherent structure in elasto-inertial turbulence. Physical Review Fluids 7 (7), 073301.
  * Dubief et al. (2023) Dubief, Y., Terrapon, V. E. & Hof, B. 2023 Elasto-inertial turbulence. Annual Review of Fluid Mechanics 55, 675–705.
  * Garg et al. (2018) Garg, Piyush, Chaudhary, Indresh, Khalid, Mohammad, Shankar, V & Subramanian, Ganesh 2018 Viscoelastic pipe flow is linearly unstable. Physical Review Letters 121 (2), 024502.
  * Groisman & Steinberg (2000) Groisman, Alexander & Steinberg, Victor 2000 Elastic turbulence in a polymer solution flow. Nature 405 (6782), 53–55.
  * Khalid et al. (2021a) Khalid, Mohammad, Chaudhary, Indresh, Garg, Piyush, Shankar, V & Subramanian, Ganesh 2021a The centre-mode instability of viscoelastic plane Poiseuille flow. Journal of Fluid Mechanics 915.
  * Khalid et al. (2021b) Khalid, Mohammad, Shankar, V & Subramanian, Ganesh 2021b Continuous pathway between the elasto-inertial and elastic turbulent states in viscoelastic channel flow. Physical Review Letters 127 (13), 134502.
  * Larson et al. (1990) Larson, Ronald G, Shaqfeh, Eric SG & Muller, Susan J 1990 A purely elastic instability in Taylor–Couette flow. Journal of Fluid Mechanics 218, 573–600.
  * Lellep et al. (2023a) Lellep, M., Linkmann, M. & Morozov, A. 2023a Linear instability analysis of purely elastic travelling waves in pressure-driven channel flows. J. Fluid Mech. 959, R1.
  * Lellep et al. (2023b) Lellep, M., Linkmann, M. & Morozov, A. 2023b Purely elastic turbulence in pressure-driven channel flows. http://arxiv.org/arXiv.2312.0891 .
  * Lewy & Kerswell (2023) Lewy, T. & Kerswell, R. R. 2023 The polymer diffusive instability in highly concentrated polymer fluids. https://arxiv.org/abs/2311.05251 .
  * Morozov (2022) Morozov, Alexander 2022 Coherent structures in plane channel flow of dilute polymer solutions with vanishing inertia. Physical Review Letters 129 (1), 017801.
  * Page et al. (2020) Page, Jacob, Dubief, Yves & Kerswell, Rich R 2020 Exact traveling wave solutions in viscoelastic channel flow. Physical Review Letters 125 (15), 154501.
  * Rallison & Hinch (1995) Rallison, J. M. & Hinch, E. J. 1995 Instability of a high-speed submerged elastic jet. J. Fluid Mech. 228, 311–324.
  * Samanta et al. (2013) Samanta, Devranjan, Dubief, Yves, Holzner, Markus, Schäfer, Christof, Morozov, Alexander N, Wagner, Christian & Hof, Björn 2013 Elasto-inertial turbulence. Proceedings of the National Academy of Sciences 110 (26), 10557–10562.
  * Shaqfeh (1996) Shaqfeh, ES G 1996 Purely elastic instabilities in viscometric flows. Annual Review of Fluid Mechanics (28), 129–185.
  * Shekar et al. (2021) Shekar, A., McMullen, R. M., McKeon, B. J. & Graham, M. D. 2021 Tollmien-Schlichting route to elastoinertial turbulence in channel flow. Physical Review Fluids 6, 093301.
  * Shekar et al. (2019) Shekar, A., McMullen, R. M., Wang, S-N., McKeon, B. J. & Graham, M. D. 2019 Critical-layer structures and mechanisms in elastoinertial turbulence. Phys. Rev. Lett. 122, 124503.
  * Sid et al. (2018) Sid, Samir, Terrapon, VE & Dubief, Y 2018 Two-dimensional dynamics of elasto-inertial turbulence and its role in polymer drag reduction. Physical Review Fluids 3 (1), 011301.
  * Wan et al. (2021) Wan, D., Sun, G. & Zhang, M. 2021 Subcritical and suprcritical bifurcations in axisymmetric viscoelastic pipe flows. Journal of Fluid Mechanics 929, A16.