eurm10 msam10

# Weakly nonlinear analysis of the viscoelastic instability in channel flow
for finite and vanishing Reynolds numbers

Gergely Buza1Jacob Page2andRich R. Kerswell1
gb643@cam.ac.ukjacob.page<EMAIL_ADDRESS>1Department of
Applied Mathematics and Theoretical Physics, University of Cambridge, CB3 0WA,
UK.
2School of Mathematics, University of Edinburgh, EH9 3FD, UK

(?; revised ?; accepted ?.)

###### Abstract

The recently-discovered centre-mode instability of rectilinear viscoelastic
shear flow (Garg et al. Phy. Rev. Lett. 121, 024502, 2018) has offered an
explanation for the origin of elasto-inertial turbulence (EIT) which occurs at
lower Weissenberg ($Wi$) numbers. In support of this, we show using weakly
nonlinear analysis that the subcriticality found in Page et al. (Phys. Rev.
Lett. 125, 154501, 2020) is generic across the neutral curve with the
instability only becoming supercritical at low Reynolds ($Re$) numbers and
high $Wi$. We demonstrate that the instability can be viewed as purely elastic
in origin even for $Re=O(10^{3})$, rather than ‘elasto-inertial’, as the
underlying shear does not energise the instability. It is also found that the
introduction of a realistic maximum polymer extension length, $L_{max}$, in
the FENE-P model moves the neutral curve closer to the inertialess $Re=0$
limit at a fixed ratio of solvent-to-solution viscosities, $\beta$. In the
dilute limit ($\beta\rightarrow 1$) with $L_{max}=O(100)$, the linear
instability can brought down to more physically-relevant $Wi\gtrsim 110$ at
$\beta=0.98$, compared with the threshold $Wi=O(10^{3})$ at $\beta=0.994$
reported recently by Khalid et al. (arXiv: 2103.06794) for an Oldroyd-B fluid.
Again the instability is subcritical implying that inertialess rectilinear
viscoelastic shear flow is nonlinearly unstable - i.e. unstable to finite
amplitude disturbances - for even lower $Wi$.

###### keywords:

viscoelasticity, shear flow, instability

##  1 Introduction

Viscoelastic flows have been of interest ever since the observation 70 years
ago that a substantial reduction in viscous drag on a wall of a pipe carrying
turbulent flow is possible after adding only a few parts per millon of long-
chain polymers (Toms, 1948). Just as curiously, adding further polymer quickly
saturates this effect when the so-called ‘maximum drag reduction’ regime (MDR)
is entered (Virk, 1970), with skin friction reduced by $\sim$$80$% relative to
its Newtonian value. Efforts to explain this phenomenon have naturally
focussed on understanding how low polymer concentrations moderate Newtonian
turbulence (NT) (e.g. Lumley, 1969; Tabor & de Gennes, 1986; Procaccia et al.,
2008; White & Mungal, 2008). However, the discovery of a new form of
viscoelastic turbulence - ‘elasto-inertial’ turbulence (EIT) - in 2013
(Samanta et al., 2013; Dubief et al., 2013; Sid et al., 2018) which exists at
large Reynolds number $Re=O(10^{3})$ and Weissenberg number $Wi=O(10)$ has
provided a competing and even less well understood possibility. Provided $Wi$
is large enough, EIT can exist at much lower $Re$ than NT explaining what has
been labelled in the past as ‘early turbulence’ (Jones & Maddock, 1966;
Goldstein et al., 1969; Hansen & Little, 1974; Draad et al., 1998; Samanta et
al., 2013; Choueiri et al., 2018; Chandra et al., 2018). At higher but fixed
$Re$, it is also possible, as the polymer concentration is steadily increased
from zero, to relaminarize NT before triggering EIT (Choueiri et al., 2018;
Chandra et al., 2018). In DNS, increasing $Wi$ from a state of EIT quenches
the flow down to a simple travelling wave solution and presumably laminar flow
if $Wi$ is large enough (e.g. see figure 2 in Page et al., 2020; Dubief et
al., 2020). At even higher $Re$, it is currently unclear whether the two types
of turbulence merge or co-exist, and how MDR fits into the situation remains
an outstanding issue (e.g. Xi & Graham, 2010, 2012; Graham, 2014; Samanta et
al., 2013; Choueiri et al., 2018, 2021; Lopez et al., 2019).

Further questions also exist as to how EIT relates to another form of
viscoelastic turbulence - ‘elastic’ turbulence (ET) - that was discovered a
decade earlier (Groisman & Steinberg, 2000). This is generated by the well-
known ‘elastic’ linear instability of curved streamlines (Larson et al., 1990;
Shaqfeh, 1996) and exists at vanishingly small Reynolds numbers so inertial
effects are unambiguously irrelevant for sustaining the turbulence. This
elastic instability is also possible in planar geometries, but requires
finite-amplitude disturbances to generate streamline curvature (Meulenbroek et
al., 2004; Morozov & Saarloos, 2007). In contrast to the inertialess ET, a
fairly large $Re$ is required for EIT, indicating that inertia is important
here. This suggests that EIT and ET are distinct phenomena (e.g. see figure 30
Chaudhary et al., 2021) yet they could still be two extremes of the same whole
(Samanta et al., 2013; Qin et al., 2019; Choueiri et al., 2021; Steinberg,
2021). Finally, the underlying mechanism which sustains EIT has yet to be
clarified (Dubief et al., 2013; Terrapon et al., 2015; Sid et al., 2018;
Shekar et al., 2018, 2020; Page et al., 2020; Chaudhary et al., 2021).

A major step forward in explaining the origin of EIT was made recently when a
linear instability was found at relatively high $Wi\gtrsim 20$ which could
reach down to a threshold $Re_{c}\approx 63$ in pipe flow (Garg et al., 2018;
Chaudhary et al., 2021). This finding overturned a long held view that no new
linear instability would appear by adding polymers to a Newtonian rectilinear
shear flow (see Chaudhary et al., 2019, 2021, for an extensive historical
discussion of this point). This instability was also confirmed in channel flow
(Khalid et al., 2021a) using an Oldroyd-B fluid but was found absent in an
upper-convected Maxwell (UCM) fluid (Chaudhary et al., 2019). The instability
is a centre-mode instability which has a phase speed close to the maximum
base-flow speed and appears to need inertia (finite $Re$) to exist: in a
channel with an experimentally-relevant $\beta$ (the ratio of solvent-to-
solution viscosities) of $0.9$ and elasticity number of $0.1$, the threshold
$Re_{c}\approx 200$. However, in the dilute limit ($\beta\rightarrow 1$) and
in contrast with pipe flow, Khalid et al. (2021a) found that $Re_{c}$ could be
pushed down to $\approx 5$ by the time $\beta$ reached $0.99$, albeit at very
large $Wi$. Further computations (Khalid et al., 2021b) have confirmed that
the elastic limit of $Re=0$ can indeed be reached at $\beta=0.9905$ and
$Wi\approx 2500$. Looking beyond the extreme value of $Wi$ –which is
apparently achievable experimentally (Vashney & Steinberg, 2018; Schnapp &
Steinberg, 2021)– this result has established a fascinating connection between
an instability which appears to need inertia, elasticity and solvent viscosity
(finite $(1-\beta)$) and a purely elastic instability when $(1-\beta)$ is
small enough (Khalid et al. (2021b) refer to this as an ‘ultra dilute’ polymer
solution).

However, EIT appears at lower $Wi$ (figure 2 in Page et al., 2020) and
sometimes lower $Re$ at a given $Wi$ (see figure 1b in Choueiri et al., 2021)
than the centre-mode instability. For example, in channel flow at $Re=1000$
and $\beta=0.9$ in a FENE-P fluid with $L_{max}=500$, EIT occurs around
$Wi=20$, whereas the centre-mode instability threshold is $Wi\approx 70$
(figure 2 (left) in Page et al., 2020). This means that if EIT is dynamically
connected to this instability, the hierarchy of nonlinear solutions which
emerge from the linear instability must be substantially subcritical, reaching
to $Wi$ far below those of the neutral curve (and similarly for $Re$ for high
enough $Wi$). This was confirmed in one specific case on the neutral curve -
$(Re,Wi,\beta)=(60,26.9,0.9)$ – where the bifurcation was shown to be strongly
subcritical with the branch of travelling waves solutions reaching down to
$Wi=8.77$ (Page et al., 2020). Moreover, the travelling wave solutions adopt a
distinctive ‘arrowhead’ form in the polymer stress when $Wi$ is small enough
which can be recognised as an intermittently-observed coherent structure in
the DNS of EIT (Dubief et al., 2020).

The primary purpose of this paper is to back this initial finding of
subcriticality up by carrying out a systematic survey of whether the centre
mode bifurcation is sub- or supercritical across the entire neutral curve for
a typical value of $\beta$ of $0.9$ using weakly nonlinear analysis (Stuart,
1960; Watson, 1960). In doing so, we also take the opportunity to confirm that
the instability is present for a FENE-P fluid with reasonable maximum polymer
extension $L_{max}$ (see 1d) and, spurred on by the recent results of Khalid
et al. (2021b), explore how the presence of finite $L_{max}$ affects the
dilute limit ($\beta\rightarrow 1$) where $Re=0$ can be reached. We also
examine the energetic source term, or terms, for the instability uncovering a
consistent picture even on the part of the neutral curve reaching to high
$Re$.

The plan of the paper is as follows. In §2, the FENE-P model is introduced and
the presence or not of polymer diffusion as indicated by a Schmidt number $Sc$
is discussed. The weakly nonlinear expansions are also introduced. While this
is now an established method in the fluid dynamicists’ toolbox, for
viscoelastic models where the (coarse-grained) local polymer configuration is
represented by a positive definite conformation tensor $\mathbf{C}$, there are
some technicalities which need some attention. We follow the framework
recently suggested by Hameduddin et al. (2018, 2019) to treat this issue which
requires a bit more formal development than is normal. Having set this up, §3
then presents the weakly nonlinear analysis which proceeds as usual albeit
with a proxy for $\mathbf{C}$ being expanded instead of $\mathbf{C}$ itself.
Results in §4 are arranged as follows: §4.1 and §4.2 consider
$(\beta,L_{max})=(0.9,500)$ with $Sc\to\infty$; §4.3 considers
$(\beta,L_{max},Sc)=(0.9,100,10^{6})$; §4.4 performs an energy analysis over
the neutral curves of §4.1 and 4.3; and finally §4.5 examines the $Re=0$
situation varying $\beta$ over the approximate range of $[0.97,0.99]$ for
$Wi\leq 200$ and $L_{max}\in[40,100]$ ($Sc\to\infty$). Lastly, §5 presents a
discussion of the paper’s results.

##  2 Formulation

We consider pressure-driven viscoelastic flow between two parallel,
stationary, rigid plates separated by a distance $2h$ and assume that the flow
is governed by the FENE-P model

$\displaystyle\partial_{t}\mathbf{u}+\left(\mathbf{u}\cdot\mbox{\boldmath$\nabla$}\right)\mathbf{u}+\mbox{\boldmath$\nabla$}p$
$\displaystyle=\frac{\beta}{Re}\Delta\mathbf{u}+\frac{(1-\beta)}{Re}\mbox{\boldmath$\nabla$}\cdot\mathbf{T}(\mathbf{C}),$
(1a) $\displaystyle\mbox{\boldmath$\nabla$}\cdot\mathbf{u}$ $\displaystyle=0,$
(1b)
$\displaystyle\partial_{t}\mathbf{C}+\left(\mathbf{u}\cdot\mbox{\boldmath$\nabla$}\right)\mathbf{C}+\mathbf{T}(\mathbf{C})$
$\displaystyle=\mathbf{C}\cdot\mbox{\boldmath$\nabla$}\mathbf{u}+\left(\mbox{\boldmath$\nabla$}\mathbf{u}\right)^{T}\cdot\mathbf{C}+\frac{1}{ReSc}\Delta\mathbf{C}.$
(1c) The constitutive relation for the polymer stress, $\mathbf{T}$, is given
by the Peterlin function
$\mathbf{T}(\mathbf{C}):=\frac{1}{Wi}\Big{(}f(\mathrm{tr}\,\mathbf{C})\mathbf{C}-{\mathbf{I}}\Big{)},\quad{\rm
where}\quad f(x):=\left(1-\frac{x-3}{L^{2}_{max}}\right)^{-1}$ (1d)

with $L_{max}$ denoting the maximum extensibility of polymer chains. Here
$\mathbf{C}\in\mathrm{Pos}(3)$ (the set of positive definite 3x3 matrices) is
the polymer conformation tensor and $\beta\in[0,1]$ denotes the viscosity
ratio, $\beta:=\nu_{s}/\nu$, where $\nu_{s}$ and $\nu_{p}=\nu-\nu_{s}$ are the
solvent and polymer contributions to the total kinematic viscosity $\nu$. The
equations are non-dimensionalised by $h$ and the bulk speed

$U_{b}:=\frac{1}{2h}\int^{h}_{-h}u_{x}\,dy$ (2)

which, through adjusting the pressure gradient appropriately, is kept constant
so that the Reynolds and Weissenberg numbers are defined as

$Re:=\frac{hU_{b}}{\nu},\quad Wi:=\frac{\tau U_{b}}{h}$ (3)

where $\tau$ is the polymer relaxation time. Polymer diffusion - the last term
in Eq. (1c) - is often omitted as the typical magnitude of the Schmidt number,
$Sc\sim O(10^{6})$. Here it is retained throughout the nonlinear analysis to:
1. allow a more realistic comparison with results from direct numerical
simulations (DNS), where a relatively low Schmidt number ($Sc\sim O(10^{3})$)
is required for the solver to converge (Page et al., 2020), and 2. assess its
importance more generally. Non-slip boundary conditions are imposed on the
velocity field. If an infinite Schmidt number $Sc$ is considered, no boundary
conditions for the conformation tensor $\mathbf{C}$ are needed. In the case of
finite Schmidt numbers, we apply $Sc\to\infty$ at the boundary to retain this
situation Sid et al. (2018).

In the course of this work, we compute neutral curves for the recently-
discovered centre mode instability in a pipe (Garg et al., 2018; Chaudhary et
al., 2021) and a channel flow (Khalid et al., 2021a, b). The marginally-stable
eigenfunctions form the basis of a weakly nonlinear expansion in the amplitude
of the bifurcating solution. The key objective here is to ascertain whether
the bifurcation is supercritical or subcritical. Subcriticality would indicate
that bifurcated solutions exist beyond the parameter domain of linear
instability, thereby implying that the flow is nonlinearly unstable - i.e.
unstable to sufficiently large amplitude disturbances - in new, potentially
more interesting parameter regimes. A case in point is the very recent
discovery that the centre mode instability still operates at $Re=0$ albeit at
very high $Wi=O(1000)$ and ultra-dilute polymer solutions of
$1-\beta=O(10^{-3})$ (Khalid et al., 2021b). While these extremes are on the
margins of physical relevance, a strongly-subcritical instability could still
see its consequences in the form of finite amplitude solutions at vastly
different $Wi$ and $\beta$.

###  2.1 Base state

The base state to (1a)-(1c) is the steady unidirectional solution and
satisfies the following reduced set of equations:

$\displaystyle\partial_{x}p=\frac{\beta}{Re}\partial_{yy}u_{x}+\frac{(1-\beta)}{ReWi}\left[\frac{\left(f(\mathrm{tr}\,\mathbf{C})\right)^{2}}{L^{2}_{max}}\mathrm{tr}(\partial_{y}C)C_{xy}+f(\mathrm{tr}\,\mathbf{C})\partial_{y}C_{xy}\right],$
(4a)
$\displaystyle\frac{1}{Wi}\left(f(\mathrm{tr}\,\mathbf{C})C_{xx}-1\right)=2C_{xy}\partial_{y}u_{x}+\frac{1}{ReSc}\partial_{yy}C_{xx},$
(4b)
$\displaystyle\frac{1}{Wi}\left(f(\mathrm{tr}\,\mathbf{C})C_{yy}-1\right)=\frac{1}{ReSc}\partial_{yy}C_{yy},$
(4c)
$\displaystyle\frac{1}{Wi}\left(f(\mathrm{tr}\,\mathbf{C})C_{zz}-1\right)=\frac{1}{ReSc}\partial_{yy}C_{zz},$
(4d)
$\displaystyle\frac{1}{Wi}\left(f(\mathrm{tr}\,\mathbf{C})C_{xy}\right)=C_{yy}\partial_{y}u_{x}+\frac{1}{ReSc}\partial_{yy}C_{xy}.$
(4e)

where $\mathbf{u}=u_{x}\hat{\bf x}+u_{y}\hat{\bf y}$. Since the $Re$ is based
on the bulk speed, the applied pressure gradient is adjusted until the bulk
speed is unity (after non-dimensionalisation) (e.g. Samanta et al., 2013;
Dubief et al., 2013; Sid et al., 2018; Dubief et al., 2020). Figure 1 displays
the base state $(\mathbf{u}_{b},p_{b},\mathbf{C}_{b})$ for a particular
parameter combination. It is worth remarking that $U_{max}$ is very nearly
$1.5$ in units of $U_{b}$ for the parameter settings considered so that a $Wi$
based upon the bulk velocity (as here) is very close to two thirds of a
Weissenberg number based on $U_{max}$ (Garg et al., 2018; Chaudhary et al.,
2021; Khalid et al., 2021a, b).

Figure 1: Laminar base state at $\beta=0.9$, $L_{max}=500$, $Sc\to\infty$,
$Wi=60$, $Re=68$. The components of the base flow conformation tensor,
$\mathbf{C}_{b}$, are normalized by their value at the bottom wall ($y=-1$).
($\left.C_{b,xx}\right|_{y=-1}=39235$,
$\left.C_{b,yy}\right|_{y=-1}=\left.C_{b,zz}\right|_{y=-1}=0.84$,
$\left.C_{b,xy}\right|_{y=-1}=129$.)

###  2.2 Perturbative Expansions

The weakly nonlinear expansions for the velocity and pressure components are
straightforwardly written in the form

$\mathbf{u}=\mathbf{u}_{b}+\sum_{k=1}^{N}\varepsilon^{k}\mathbf{u}_{(k)},\quad
p=p_{b}+\sum_{k=1}^{N}\varepsilon^{k}p_{(k)}.$ (5)

However, the conformation tensor, $\mathbf{C}$, calls for a more careful
treatment, since the set of positive definite $3\times 3$ matrices,
$\mathrm{Pos}(3)$, cannot be a vector space. Instead, it may be endowed with
the structure of a complete Riemannian manifold. Perturbations of order
$\varepsilon^{k}$ still make sense in this setting, but one has to interpret
the $\varepsilon^{k}$ distance in terms of the metric arising from the
Riemannian structure of the manifold $\mathrm{Pos}(3)$. In developing
perturbations for the conformation tensor, $\mathbf{C}$, we follow the
framework of Hameduddin et al. (2018, 2019) who focussed on precisely this
issue. We may view $\mathbf{C}$ as the left Cauchy-Green tensor associated to
the polymer deformation, i.e.,

$\mathbf{C}=\mathbf{F}\mathbf{F}^{T},$

where $\mathbf{F}$ denotes the deformation gradient with thermal equilibrium
taken as the reference configuration. A further decomposition of $\mathbf{F}$
into two successive deformations, which may be written as

$\mathbf{F}=\mathbf{F}_{b}\mathbf{L}$ (6)

separates the deformation corresponding to the perturbation, $\mathbf{L}$,
from the deformation associated with the base state, which may be expressed
as111This representation is not unique, any
$\mathbf{F}_{b}=\mathbf{C}_{b}^{\frac{1}{2}}{\bf R}$ works with ${\bf
R}\in\mathrm{SO}(3)$. The choice ${\bf R}={\bf I}$ is natural in the sense
that it allows for a geodesic between $\mathbf{C}_{b}$ and $\mathbf{C}$ to be
expressed solely in terms of $\mathbf{F}_{b}$ and $\mathbf{G}$.

$\mathbf{F}_{b}=\mathbf{C}_{b}^{\frac{1}{2}}.$

The fluctuating deformation gradient, $\mathbf{L}$, has an associated left
Cauchy-Green tensor $\mathbf{G}=\mathbf{L}\mathbf{L}^{T}$. Combining these
observations, we have that

$\mathbf{C}=\mathbf{F}_{b}\mathbf{G}\mathbf{F}_{b}^{T}.$ (7)

The tensor $\mathbf{G}$ is necessarily positive definite since $\mathbf{C}$
is, and by nature it acts as the conformation tensor representing the
fluctuations of $\mathbf{C}$ around $\mathbf{C}_{b}$.

The evolution equation (1c) for the conformation tensor can be rewritten in
terms of $\mathbf{G}$ as follows:

$\partial_{t}\mathbf{G}+(\mathbf{u}\cdot\mbox{\boldmath$\nabla$})\mathbf{G}=2\mathrm{sym}\left(\mathbf{G}h(\mathbf{u})\right)-\mathbf{F}_{b}^{-1}\mathbf{T}\mathbf{F}_{b}^{-T},$
(8)

with

$h(\mathbf{u})=\mathbf{F}_{b}^{T}\cdot\mbox{\boldmath$\nabla$}\mathbf{u}\cdot\mathbf{F}_{b}^{-T}-\left(\mathbf{F}_{b}^{-1}(\mathbf{u}\cdot\mbox{\boldmath$\nabla$})\mathbf{F}_{b}\right)^{T}.$

As described by Hameduddin et al. (2019), an additive expansion of the form
(5) no longer makes sense on $\mathrm{Pos}(3)$, since there is no a priori
guarantee that the resulting $\mathbf{C}$ remains positive definite. Instead,
Hameduddin et al. (2019) proposed a multiplicative expansion based on the
decomposition (6) that consists of a series of successively smaller
deformations, which may be written in the form

$\mathbf{L}_{wnl}=\mathbf{L}_{(1)}^{\varepsilon}\mathbf{L}_{(2)}^{\varepsilon^{2}}\cdots\mathbf{L}_{(N)}^{\varepsilon^{N}}.$

The matrix $\mathbf{L}_{wnl}$ may differ from $\mathbf{L}$ given in (6) by a
rotation only.

Under the additional assumption that the $\mathbf{L}_{(k)}^{\varepsilon^{k}}$
are rotation free with $\mathrm{det}(\mathbf{L}_{(k)})>0$, each
$\mathbf{L}_{k}$ is positive definite. The conformation tensors associated to
these deformations are then given by
$\mathbf{G}_{(k)}^{\varepsilon^{k}}=\mathbf{L}_{(k)}^{\varepsilon^{k}}\left(\mathbf{L}_{(k)}^{\varepsilon^{k}}\right)^{T}$.
To make sense of $\varepsilon$-magnitude perturbations, we make use of the
Riemannian manifold structure of $\mathrm{Pos}(3)$. In particular, the
$\mathbf{G}_{(k)}^{\varepsilon^{k}}$ may be thought of as length
$\sim|\varepsilon|^{k}$ geodesics emanating from $\mathbf{I}$ on the manifold
$\mathrm{Pos}(3)$. That is, we may take $\mathbfcal G_{(k)}\in
T_{\mathbf{I}}\mathrm{Pos}(3)=\mathrm{Sym}(3)$ such that

$\mathbf{G}_{(k)}^{\varepsilon^{k}}=\mathrm{exp}\left(\varepsilon^{k}\mathbfcal
G_{(k)}\right),$

with

$d(\mathbf{I},\mathbf{G}_{(k)}^{\varepsilon^{k}})=|\varepsilon|^{k}\|\mathbfcal
G_{(k)}\|_{F},$

where $d$ is the metric induced by the Riemannian structure of
$\mathrm{Pos}(3)$. Note that this is analogous to weakly nonlinear expansions
on vector spaces equipped with the Frobenius norm, only now we measure the
corresponding distance on $\mathrm{Pos}(3)$ with the Riemannian metric.

This approach eventually leads to an expansion of the form

$\displaystyle\mathbf{G}$
$\displaystyle=\mathrm{exp}\left(\varepsilon\frac{\mathbfcal
G_{(1)}}{2}\right)\cdots\mathrm{exp}\left(\varepsilon^{N-1}\frac{\mathbfcal
G_{(N-1)}}{2}\right)\mathrm{exp}\left(\varepsilon^{N}\mathbfcal
G_{(N)}\right)\mathrm{exp}\left(\varepsilon^{N-1}\frac{\mathbfcal
G_{(N-1)}}{2}\right)\cdots\mathrm{exp}\left(\varepsilon\frac{\mathbfcal
G_{(1)}}{2}\right)$ $\displaystyle=\mathbf{I}+\varepsilon\mathbfcal
G_{(1)}+\varepsilon^{2}\left(\mathbfcal G_{(2)}+\frac{\mathbfcal
G_{(1)}^{2}}{2}\right)+\varepsilon^{3}\left(\mathbfcal
G_{(3)}+\mathrm{sym}\left(\mathbfcal G_{(1)}\mathbfcal
G_{(2)}\right)+\frac{\mathbfcal G_{(3)}^{3}}{6}\right)+\ldots$ (9)

This representation of the weakly nonlinear terms is equivalent to a standard
expansion for $\mathbf{C}$ of the form (5), as the operation
$\mathbf{G}_{(j)}\mapsto\mathbf{F}_{b}\mathbf{G}_{(j)}\mathbf{F}_{b}^{T}$
serves as a bijection between the two solution sets, as long as
$\mathbf{F}_{b}\in\mathrm{GL}(3)$. Physically, this is always satisfied, as
$\mathrm{det}(\mathbf{F}_{b})=0$ would imply that material elements are
compressed to zero volume.

While the new formulation does not in practice modify the mechanics of
constructing a weakly nonlinear expansion, the mathematical consistency of the
approach yields a variety of tools for measuring perturbations on
$\mathrm{Pos}(3)$ in the only suitable manner, according to the corresponding
metric. One such measure, which we shall use frequently in the sections to
follow, is the geodesic distance from the mean, given by

$d(\mathbf{C}_{b},\mathbf{C})=d({\bf
I},\mathbf{G})=\sqrt{\mathrm{tr}\mathbfcal G^{2}}.$ (10)

##  3 Weakly nonlinear analysis

Let $\mbox{\boldmath$\varphi$}=(u_{x},u_{y},p,G_{xx},G_{yy},G_{zz},G_{xy})$
denote the vector composed of all state variables. This is further decomposed
into two parts: a contribution from the base state and a fluctuating part as
follows

$\mbox{\boldmath$\varphi$}=\mbox{\boldmath$\varphi$}_{b}+\hat{\mbox{\boldmath$\varphi$}},$

where the interest is now in solving the governing system (1) for the
perturbations $\hat{\mbox{\boldmath$\varphi$}}$. The Peterlin function (1d)
for $\mathbf{T}$ is first expanded around the base conformation state,
$\mathbf{C}_{b}$ as follows

$\mathbf{T}(\mathbf{C})=\mathbf{T}(\mathbf{C}_{b})+D\mathbf{T}(\mathbf{C}_{b})[\hat{\mathbf{C}}]+\frac{1}{2}D^{2}\mathbf{T}(\mathbf{C}_{b})[\hat{\mathbf{C}},\hat{\mathbf{C}}]+\frac{1}{6}D^{3}\mathbf{T}(\mathbf{C}_{b})[\hat{\mathbf{C}},\hat{\mathbf{C}},\hat{\mathbf{C}}]+\ldots.$
(11)

For the analysis which follows, it suffices to perform the above expansion
(11) up to third order and to compress the notation, we shall only consider
$Wi$ and $Re$ as varying parameters. The others, $\beta$ and $Sc$, are assumed
fixed but similar expansions for them may be obtained in an analogous fashion.
After a subtraction of the laminar solution, equation (1) can be written in an
operator form locally around the base state $(\mathbf{u}_{b},\mathbf{C}_{b})$
as

$\mathcal{L}\left(Re,Wi\right)\left[\hat{\mbox{\boldmath$\varphi$}}\right]+\mathcal{B}\left(Re,Wi\right)\left[\hat{\mbox{\boldmath$\varphi$}},\hat{\mbox{\boldmath$\varphi$}}\right]+\mathcal{T}\left(Re,Wi\right)\left[\hat{\mbox{\boldmath$\varphi$}},\hat{\mbox{\boldmath$\varphi$}},\hat{\mbox{\boldmath$\varphi$}}\right]={\bf
0},$ (12)

where $\mathcal{L}\left(Re,Wi\right)$ is linear,
$\mathcal{B}\left(Re,Wi\right)$ is bilinear and
$\mathcal{T}\left(Re,Wi\right)$ is symmetric trilinear. These are given
explicitly as

$\displaystyle\mathcal{L}\left(Re,Wi\right)\left[\hat{\mbox{\boldmath$\varphi$}}\right]=\begin{pmatrix}\partial_{t}\hat{\mathbf{u}}+(\mathbf{u}_{b}\cdot\mbox{\boldmath$\nabla$})\hat{\mathbf{u}}+(\hat{\mathbf{u}}\cdot\mbox{\boldmath$\nabla$})\mathbf{u}_{b}+\mbox{\boldmath$\nabla$}\hat{p}-\frac{\beta}{Re}\Delta\hat{\mathbf{u}}-\frac{1-\beta}{Re}\mbox{\boldmath$\nabla$}\cdot\left(D\mathbf{T}(\mathbf{C}_{b})\left[\mathbf{F}_{b}\hat{\mathbf{G}}\mathbf{F}_{b}^{T}\right]\right)\\\
\\\ \mbox{\boldmath$\nabla$}\cdot\hat{\mathbf{u}}\\\ \\\
\partial_{t}\hat{\mathbf{G}}+(\mathbf{u}_{b}\cdot\mbox{\boldmath$\nabla$})\hat{\mathbf{G}}-2\mathrm{sym}\left(h\left(\hat{\mathbf{u}}\right)+\hat{\mathbf{G}}h(\mathbf{u}_{b})\right)+\mathbf{F}_{b}^{-1}D\mathbf{T}(\mathbf{C}_{b})\left[\mathbf{F}_{b}\hat{\mathbf{G}}\mathbf{F}_{b}^{T}\right]\mathbf{F}_{b}^{-T}\\\
\hskip
170.71652pt-\frac{1}{ReSc}\mathbf{F}_{b}^{-1}\Delta\left(\mathbf{F}_{b}\hat{\mathbf{G}}\mathbf{F}_{b}^{T}\right)\mathbf{F}_{b}^{-T}\\\
\end{pmatrix},$
$\displaystyle\mathcal{B}\left(Re,Wi\right)\left[\hat{\mbox{\boldmath$\varphi$}}_{1},\hat{\mbox{\boldmath$\varphi$}}_{2}\right]=\begin{pmatrix}(\hat{\mathbf{u}}_{1}\cdot\mbox{\boldmath$\nabla$})\hat{\mathbf{u}}_{2}-\frac{1-\beta}{2Re}\mbox{\boldmath$\nabla$}\cdot\left(D^{2}\mathbf{T}(\mathbf{C}_{b})\left[\mathbf{F}_{b}\hat{\mathbf{G}}_{1}\mathbf{F}_{b}^{T},\mathbf{F}_{b}\hat{\mathbf{G}}_{2}\mathbf{F}_{b}^{T}\right]\right)\\\
\\\ 0\\\ \\\
(\hat{\mathbf{u}}_{1}\cdot\mbox{\boldmath$\nabla$})\hat{\mathbf{G}}_{2}-2\mathrm{sym}\left(\hat{\mathbf{G}}_{1}h(\hat{\mathbf{u}}_{2})\right)+\frac{1}{2}\mathbf{F}_{b}^{-1}D^{2}\mathbf{T}(\mathbf{C}_{b})\left[\mathbf{F}_{b}\hat{\mathbf{G}}_{1}\mathbf{F}_{b}^{T},\mathbf{F}_{b}\hat{\mathbf{G}}_{2}\mathbf{F}_{b}^{T}\right]\mathbf{F}_{b}^{-T}\end{pmatrix},$
$\displaystyle\mathcal{T}\left(Re,Wi\right)\left[\hat{\mbox{\boldmath$\varphi$}}_{1},\hat{\mbox{\boldmath$\varphi$}}_{2},\hat{\mbox{\boldmath$\varphi$}}_{3}\right]=\begin{pmatrix}-\frac{1-\beta}{6Re}\mbox{\boldmath$\nabla$}\cdot\left(D^{3}\mathbf{T}(\mathbf{C}_{b})\left[\mathbf{F}_{b}\hat{\mathbf{G}}_{1}\mathbf{F}_{b}^{T},\mathbf{F}_{b}\hat{\mathbf{G}}_{2}\mathbf{F}_{b}^{T},\mathbf{F}_{b}\hat{\mathbf{G}}_{3}\mathbf{F}_{b}^{T}\right]\right)\\\
\\\ 0\\\ \\\
\frac{1}{6}\mathbf{F}_{b}^{-1}D^{3}\mathbf{T}(\mathbf{C}_{b})\left[\mathbf{F}_{b}\hat{\mathbf{G}}_{1}\mathbf{F}_{b}^{T},\mathbf{F}_{b}\hat{\mathbf{G}}_{2}\mathbf{F}_{b}^{T},\mathbf{F}_{b}\hat{\mathbf{G}}_{3}\mathbf{F}_{b}^{T}\right]\mathbf{F}_{b}^{-T}\end{pmatrix}.$

It’s worth remarking that the base state $(\mathbf{u}_{b},\mathbf{C}_{b})$ in
the above operators depends on all parameter values $(Wi,Re,\beta,Sc)$ through
(4). Linear stability theory is concerned with the eigenvalue problem arising
from the linearized equations,
$\mathcal{L}\left(Re,Wi\right)\left[\hat{\mbox{\boldmath$\varphi$}}\right]={\bf
0}$. In practice, this is formally addressed by assuming a specific form of
the disturbance, and solving

$\mathcal{L}\left(Re,Wi\right)\left[\mbox{\boldmath$\varphi$}_{(1,1)}(y)\exp(ikx-i\omega
t)\right]=0,$ (13)

for pairs $(\omega,\mbox{\boldmath$\varphi$}_{(1,1)})$, where
$\omega=\omega_{r}+i\omega_{i}$ is the a priori unknown complex frequency,
$\mbox{\boldmath$\varphi$}_{(1,1)}$ is the associated eigenmode and $k$ is the
prespecified wave number.

Assume now that a bifurcation occurs at a certain triple $(Wi_{L},Re_{L},k)$,
i.e., there exists an eigenmode of (13) such that its associated
eigenfrequency is real (subsequently denoted by $\omega_{L}=\omega_{L,r}$),
which marks the state of marginal stability in the temporal sense. We wish to
uncover how the eigenfunction $\varphi_{(1,1)}$ evolves as we move slightly
away from the bifurcation point. For this, consider small perturbations to all
relevant parameters of the form

$(Wi,Re,\omega_{r})=(Wi_{L},Re_{L},\omega_{r,L})+\varepsilon^{2}(Wi_{1},Re_{1},\omega_{r,1})+\ldots,$

and formally expand the operator $\mathcal{L}$ around $(Re_{L},Wi_{L})$ as

$\mathcal{L}(Re_{L}+\varepsilon^{2}Re_{1},Wi_{L}+\varepsilon^{2}Wi_{1})=\mathcal{L}(Re_{L},Wi_{L})+\varepsilon^{2}Re_{1}\mathcal{L}^{\prime}_{Re}(Re_{L},Wi_{L})+\varepsilon^{2}Wi_{1}\mathcal{L}^{\prime}_{Wi}(Re_{L},Wi_{L}).$

The subtle difference here from standard weakly nonlinear expansions lies in
the fact that now the base state obtained from (4) depends on the parameters
$Wi$ and $Re$. To make this clear and explicit, we write

$\displaystyle\mathcal{L}^{\prime}_{Re}(Re_{L},Wi_{L})=\left.\frac{d}{dRe}\right|_{(Re_{L},Wi_{L})}\mathcal{L}=\left.\left(\frac{\partial}{\partial
Re}+\frac{\partial u_{b,i}}{\partial Re}\frac{\partial}{\partial
u_{b,i}}+\frac{\partial F_{b,ij}}{\partial Re}\frac{\partial}{\partial
F_{b,ij}}\right)\right|_{(Re_{L},Wi_{L})}\mathcal{L},$
$\displaystyle\mathcal{L}^{\prime}_{Wi}(Re_{L},Wi_{L})=\left.\frac{d}{dWi}\right|_{(Re_{L},Wi_{L})}\mathcal{L}=\left.\left(\frac{\partial}{\partial
Wi}+\frac{\partial u_{b,i}}{\partial Wi}\frac{\partial}{\partial
u_{b,i}}+\frac{\partial F_{b,ij}}{\partial Wi}\frac{\partial}{\partial
F_{b,ij}}\right)\right|_{(Re_{L},Wi_{L})}\mathcal{L},$

with

$\frac{\partial\mathcal{L}}{\partial
Re}(Re_{L},Wi_{L})[\hat{\mbox{\boldmath$\varphi$}}]=\begin{pmatrix}\frac{\beta}{Re_{L}^{2}}\Delta\hat{\mathbf{u}}+\frac{1-\beta}{Re_{L}^{2}}\mbox{\boldmath$\nabla$}\cdot\left(D\mathbf{T}(\mathbf{C}_{b})\left[\mathbf{F}_{b}\hat{\mathbf{G}}\mathbf{F}_{b}^{T}\right]\right)\\\
0\\\
\frac{1}{Re_{L}^{2}Sc}\mathbf{F}_{b}^{-1}\Delta\left(\mathbf{F}_{b}\hat{\mathbf{G}}\mathbf{F}_{b}^{T}\right)\mathbf{F}_{b}^{-T}\end{pmatrix},$

and

$\frac{\partial\mathcal{L}}{\partial
Wi}(Re_{L},Wi_{L})[\hat{\mbox{\boldmath$\varphi$}}]=\begin{pmatrix}\frac{1-\beta}{Re_{L}Wi_{L}}\mbox{\boldmath$\nabla$}\cdot\left(D\mathbf{T}(\mathbf{C}_{b})\left[\mathbf{F}_{b}\hat{\mathbf{G}}\mathbf{F}_{b}^{T}\right]\right)\\\
0\\\
-\frac{1}{Wi_{L}}\mathbf{F}_{b}^{-1}D\mathbf{T}(\mathbf{C}_{b})\left[\mathbf{F}_{b}\hat{\mathbf{G}}\mathbf{F}_{b}^{T}\right]\mathbf{F}_{b}^{-T}\end{pmatrix}.$

Due to the complexity of the laminar equations (4), the base flow’s dependence
on the parameters is sought numerically, i.e., the terms $\partial
u_{b,i}/\partial Re$ and $\partial F_{b,ij}/\partial Re$ \- and the
corresponding terms in the $Wi$ direction - are computed via a finite
difference scheme. We note here that alternatively one could also compute the
entirety of $\mathcal{L}^{\prime}_{Re}$ (and $\mathcal{L}^{\prime}_{Wi}$) with
a finite difference scheme.

To explore how the $\mbox{\boldmath$\varphi$}_{(1,1)}$ wave develops as these
parameters change, we seek solutions of (12) as a weakly nonlinear expansion
of the form

$\mbox{\boldmath$\varphi$}(t,x,y)=\mbox{\boldmath$\varphi$}_{b}(y)+\sum_{l=1}^{N}\sum_{q\in
J_{l}}\varepsilon^{l}\left(\mbox{\boldmath$\varphi$}_{(l,q)}+\tilde{\mbox{\boldmath$\varphi$}}_{(l,q)}\right)(y)\,\mathrm{exp}\big{(}iq(kx-\omega_{r}t)\big{)}+O(\varepsilon^{N+1}),$
(14)

where $J_{l}=\\{-l,-l+2,\ldots,l-2,l\\}$, and
$\tilde{\mbox{\boldmath$\varphi$}}_{(l,q)}$ is the term that represents the
dependence of $O(\varepsilon^{l})$ perturbations on the lower order
$\mathbfcal G_{(j)}$ terms in (9). For instance,
$\tilde{\mbox{\boldmath$\varphi$}}_{(1,q)}=0,$ $q\in\\{-1,1\\}$, and

$\tilde{\mbox{\boldmath$\varphi$}}_{(2,2)}=\frac{1}{2}\left(0,0,0,\left(\mathcal{G}_{(1,1)}^{2}\right)_{xx},\left(\mathcal{G}_{(1,1)}^{2}\right)_{yy},\left(\mathcal{G}_{(1,1)}^{2}\right)_{zz},\left(\mathcal{G}_{(1,1)}^{2}\right)_{xy}\right).$

To simplify the notation, let

$E_{q}:(t,x)\mapsto\mathrm{exp}\big{(}iq(kx-\omega_{r,L}t)\big{)},$

and

$\mathcal{L}_{q}[\mbox{\boldmath$\varphi$}]:=\mathcal{L}[\mbox{\boldmath$\varphi$}E_{q}].$

Now, upon substituting the specific form of $\hat{\mbox{\boldmath$\varphi$}}$
from Eq. (14) into (12), we obtain a hierarchy of problems as follows:

$\displaystyle O(\varepsilon):$
$\displaystyle\mathcal{L}_{1}[\mbox{\boldmath$\varphi$}_{(1,1)}]={\bf 0},$
(15a) $\displaystyle O(\varepsilon^{2}):$
$\displaystyle\mathcal{L}_{0}[\mbox{\boldmath$\varphi$}_{(2,0)}+\tilde{\mbox{\boldmath$\varphi$}}_{(2,0)}]+\mathcal{B}[\mbox{\boldmath$\varphi$}_{(1,1)}E_{1},\mbox{\boldmath$\varphi$}_{(1,-1)}E_{-1}]+\mathcal{B}[\mbox{\boldmath$\varphi$}_{(1,-1)}E_{-1},\mbox{\boldmath$\varphi$}_{(1,1)}E_{1}]={\bf
0},$ (15b)
$\displaystyle\mathcal{L}_{2}[\mbox{\boldmath$\varphi$}_{(2,2)}+\tilde{\mbox{\boldmath$\varphi$}}_{(2,2)}]+\mathcal{B}[\mbox{\boldmath$\varphi$}_{(1,1)}E_{1},\mbox{\boldmath$\varphi$}_{(1,1)}E_{1}]={\bf
0},$ (15c) $\displaystyle O(\varepsilon^{3}):$
$\displaystyle\mathcal{L}_{1}[\mbox{\boldmath$\varphi$}_{(3,1)}+\tilde{\mbox{\boldmath$\varphi$}}_{(3,1)}]+\mathcal{B}[\mbox{\boldmath$\varphi$}_{(1,-1)}E_{-1},(\mbox{\boldmath$\varphi$}_{(2,2)}+\tilde{\mbox{\boldmath$\varphi$}}_{(2,2)})E_{2}]$
$\displaystyle+\mathcal{B}[(\mbox{\boldmath$\varphi$}_{(2,2)}+\tilde{\mbox{\boldmath$\varphi$}}_{(2,2)})E_{2},\mbox{\boldmath$\varphi$}_{(1,-1)}E_{-1}]+\mathcal{B}[\mbox{\boldmath$\varphi$}_{(1,1)}E_{1},\mbox{\boldmath$\varphi$}_{(2,0)}+\tilde{\mbox{\boldmath$\varphi$}}_{(2,0)}]$
$\displaystyle+\mathcal{B}[\mbox{\boldmath$\varphi$}_{(2,0)}+\tilde{\mbox{\boldmath$\varphi$}}_{(2,0)},\mbox{\boldmath$\varphi$}_{(1,1)}E_{1}]+3\mathcal{T}[\mbox{\boldmath$\varphi$}_{(1,1)}E_{1},\mbox{\boldmath$\varphi$}_{(1,1)}E_{1},\mbox{\boldmath$\varphi$}_{(1,-1)}E_{-1}]$
$\displaystyle+Re_{1}\mathcal{L}^{\prime}_{Re}[\mbox{\boldmath$\varphi$}_{(1,1)}E_{1}]+Wi_{1}\mathcal{L}^{\prime}_{Wi}[\mbox{\boldmath$\varphi$}_{(1,1)}E_{1}]-i\omega_{r,1}\mbox{\boldmath$\varphi$}_{(1,1)}=:\mathcal{L}_{1}[\mbox{\boldmath$\varphi$}_{(3,1)}]+\mbox{\boldmath$\eta$}={\bf
0},$ (15d) $\displaystyle\qquad\vdots$

where $\eta$ is the known part of the last equation (15d). One subtlety in
solving the hierarchy of problems is maintaining the constancy of the
volumetric flux. This boils down to introducing a constant correction to the
pressure gradient, $\partial_{x}p_{(2,0)}$, to ensure
$\mbox{\boldmath$\varphi$}_{(2,0)}$ has zero flux. Provided that the
bifurcation is of codimension one, equation (15a) (equivalent to the linear
problem, (13)) has a non-unique solution of the form

$A\frac{\mbox{\boldmath$\varphi$}_{(1,1)}}{\|\mbox{\boldmath$\varphi$}_{(1,1)}\|_{L^{2}\left([-1,1];\mathbb{C}^{7}\right)}},\qquad
A\in\mathbb{C}.$ (16)

The aim is to map out the possible values of the steady-state amplitude $A$ in
the parameter space $(Wi,Re)$. Once an eigenmode of the form (16) is pushed
through equations (15a) to (15d), an explicit solvability condition can be
derived, as detailed in the following.

###  3.1 Solvability condition

Let us view the functions
$\mbox{\boldmath$\varphi$}_{(i,j)}:[-1,1]\rightarrow\mathbb{C}^{7}$ as
elements of $L^{2}\left([-1,1];\mathbb{C}^{7}\right)$. The inner product on
$L^{2}\left([-1,1];\mathbb{C}^{7}\right)$ is given by222 In the following, we
use an $L^{2}$ inner product on matrix valued functions as well. In this case,
we simply identify the matrices with vectors in the canonical way (i.e., we
replace the $\mathbb{C}^{7}$ inner product below the integral with a Frobenius
one).

$\langle\mbox{\boldmath$\varphi$},\mbox{\boldmath$\psi$}\rangle_{L^{2}\left([-1,1];\mathbb{C}^{7}\right)}=\int_{[-1,1]}\langle\mbox{\boldmath$\varphi$}(y),\mbox{\boldmath$\psi$}(y)\rangle_{\mathbb{C}^{7}}\,dy.$

The linear problem (15a) implies that $\mathcal{L}_{1}$ has a nontrivial
kernel. Therefore, the Fredholm alternative theorem (for elliptic PDEs)
implies the existence of a finite dimensional subspace of solutions to the
adjoint homogeneous problem

$\mathcal{L}^{*}_{1}[\mbox{\boldmath$\psi$}]={\bf 0},$

subject to the appropriate boundary conditions (matching those of the original
problem). Moreover, the original equation (15d) has a solution,
$\mbox{\boldmath$\varphi$}_{(3,1)}$, if and only if

$\langle\mbox{\boldmath$\eta$},\mbox{\boldmath$\psi$}\rangle_{L^{2}\left([-1,1];\mathbb{C}^{7}\right)}=0,\qquad\forall\mbox{\boldmath$\psi$}\in\mathrm{ker}\,\mathcal{L}^{*}_{1}\text{
satisfying the boundary conditions.}$ (17)

Assuming that the bifurcation is of codimension one, we know that
$\mathrm{dim}\left(\mathrm{ker}\,\mathcal{L}^{*}_{1}\right)=1$, so it suffices
to check (17) for any
$\mbox{\boldmath$\psi$}_{1}\in\mathrm{ker}\,\mathcal{L}^{*}_{1}$ that
satisfies the boundary conditions. With this procedure, we obtain the complex
solvability condition

$aRe_{1}+bWi_{1}+c|A|^{2}+d\omega_{r,1}=0,$ (18)

where

$\displaystyle a:=$
$\displaystyle\left\langle\mathcal{L}^{\prime}_{Re}[\mbox{\boldmath$\varphi$}_{(1,1)}E_{1}],\mbox{\boldmath$\psi$}_{1}\right\rangle_{L^{2}\left([-1,1];\mathbb{C}^{7}\right)}$
$\displaystyle b:=$
$\displaystyle\left\langle\mathcal{L}^{\prime}_{Wi}[\mbox{\boldmath$\varphi$}_{(1,1)}E_{1}],\mbox{\boldmath$\psi$}_{1}\right\rangle_{L^{2}\left([-1,1];\mathbb{C}^{7}\right)}$
$\displaystyle c:=$
$\displaystyle\Big{\langle}\mathcal{B}[\mbox{\boldmath$\varphi$}_{(1,-1)}E_{-1},(\mbox{\boldmath$\varphi$}_{(2,2)}+\tilde{\mbox{\boldmath$\varphi$}}_{(2,2)})E_{2}]+\mathcal{B}[(\mbox{\boldmath$\varphi$}_{(2,2)}+\tilde{\mbox{\boldmath$\varphi$}}_{(2,2)})E_{2},\mbox{\boldmath$\varphi$}_{(1,-1)}E_{-1}]$
$\displaystyle+\mathcal{B}[\mbox{\boldmath$\varphi$}_{(1,1)}E_{1},\mbox{\boldmath$\varphi$}_{(2,0)}+\tilde{\mbox{\boldmath$\varphi$}}_{(2,0)}]+\mathcal{B}[\mbox{\boldmath$\varphi$}_{(2,0)}+\tilde{\mbox{\boldmath$\varphi$}}_{(2,0)},\mbox{\boldmath$\varphi$}_{(1,1)}E_{1}]$
$\displaystyle+3\mathcal{T}[\mbox{\boldmath$\varphi$}_{(1,1)}E_{1},\mbox{\boldmath$\varphi$}_{(1,1)}E_{1},\mbox{\boldmath$\varphi$}_{(1,-1)}E_{-1}],\mbox{\boldmath$\psi$}_{1}\Big{\rangle}_{L^{2}\left([-1,1];\mathbb{C}^{7}\right)},$
$\displaystyle d:=$
$\displaystyle\left\langle-i\mbox{\boldmath$\varphi$}_{(1,1)},\mbox{\boldmath$\psi$}_{1}\right\rangle_{L^{2}\left([-1,1];\mathbb{C}^{7}\right)}.$

Equation (18) gives the desired relationship between the parameters
$(Wi_{1},Re_{1})$ and the steady state amplitude $A$, which allows us to track
how these finite amplitude states emerge from the bifurcation point.

##  4 Results

As indicated above, we are interested in uncovering the nature of the initial
bifurcation associated to the centre-mode instability first identified by Garg
et al. (2018) in pipe flow and, most relevantly for us, later by Khalid et al.
(2021a) in channel flow. This previous work assumed an Oldroyd-B fluid which
allows infinite polymer extension i.e. $L_{max}\to\infty$ for the FENE-P model
(1d). Given this, our objectives in what follows are two-fold. On the one
hand, we want to explore the effects of finite extensibility on the
aforementioned instability. And on the other, with the aid of the weakly
nonlinear analysis, we aim to identify parameter regions where the instability
persists beyond the neutral curve to lower $Wi$ in particular.

###  4.1 $\beta=0.9$, $L_{max}=500$ & $Sc\to\infty$

In order to test the weakly nonlinear analysis, we begin by examining the
parameter regime considered by Page et al. (2020) where $\beta=0.9$ and
$L_{max}=500$. Using $Sc=10^{3}$ to stabilise their time-stepping code, Page
et al. (2020) observed substantial subcriticality at $(Re,Wi,k)=(60,26.9,2)$
on the upper branch of the neutral curve since they were able to continue the
branch of solutions down to $Wi=8.77$. Figure 2 shows the neutral curve at
$\beta=0.9$, $L_{max}=500$ with $Sc\to\infty$: see appendix A for numerical
details. The neutral curve is insensitive to the choice of $Sc$ on the scale
of Figure 2 provided it is $\gg 10^{2}$. Alongside the neutral curve, we
display the results of the weakly nonlinear analysis by plotting a curve
corresponding to a finite (small) steady state amplitude $|A|$, as obtained
from the solvability condition (18). The linear instability is a Hopf
bifurcation and so the steady state solutions are travelling waves (in $x$)
with phase speed $\omega_{r}/k$ and a constant amplitude which decreases to
zero at the neutral curve. This finite amplitude curve in Figure 2 clearly
indicates subcriticality along the upper branch of the neutral curve.
Proceeding down to the lower branch of the curve, the Hopf bifurcation
switches to being supercritical for $Wi\gtrsim 40$ (the red dashed line
crosses the black neutral curve).

Figure 2: (Left) Neutral curve corresponding to marginal linear stability at
$\beta=0.9$, $L_{max}=500$, $Sc\to\infty$. Results of the weakly nonlinear
analysis are shown in the form of a curve at steady state amplitude $|A|=0.4$.
(Right) The development of the critical wave number, $k_{crit}$, along the
neutral curve. Since $k_{crit}$ varies monotonically along the neutral curve
its provides a convenient parametrization of it in subsequent figures. Figure
3: Bifurcation diagrams at $(Wi,Re,k)\approx(27,60,2)$ ($\Box$). Only the
unstable branch is displayed with a comparison of two methods for the
expansion for the conformation tensor.

Figure 2 confirms the subcritical behaviour observed by Page et al. (2020) at
the point $(Wi,Re,k)\approx(27,60,2)$, which is marked by a shaded square
$\Box$. The corresponding bifurcation diagrams with respect to model
parameters $Wi$ and $Re$ are shown in Figure 3. In the figure, the newly
developed approach for perturbative expansions of Hameduddin et al. (2019)
(described in §2.2) is compared with a standard expansion in the conformation
tensor, $\mathbf{C}$. As mentioned in §2.2, the two should be equivalent when
considering quantities that only depend on $\mathbf{C}$. In this context, the
real advantage of using the form of expansions established in §2.2 is that we
now have immediate access to quantities with tangible physical meaning (cf.
Figures 6, 7 and 8).

As a final check, results of the weakly nonlinear analysis are compared with a
full branch continuation computation (see appendix A for details of the
method) in Figure 4. A finite but large Schmidt number of $Sc=10^{3}$ had to
be selected for this comparison, as the latter method requires a diffusion
term to produce reliable results. The curves are in good agreement - on top of
each other near the bifurcation point but then diverging slightly (not visible
on the plots) as the amplitude increases (as they should). This divergence, of
course, is because the weakly nonlinear analysis is based upon a 3-Fourier-
mode expansion whereas the branch continuation curve is from a 40-mode Fourier
expansion.

Figure 4: Validation of the weakly nonlinear analysis at
$(Wi,Re,k)\approx(27,60,2)$ ($\Box$) with a full branch continuation
prediction. The $L^{2}$ norms are taken over the whole domain
$\Omega=[0,2\pi/k]\times[-1,1]$.

We now examine the bifurcation on the lower branch of the neutral curve for
$Wi>40$ to confirm the supercriticality predicted by the weakly nonlinear
analysis. In Figure 5, bifurcation diagrams resulting from the weakly
nonlinear analysis for the point $\triangle$ in Figure 2 are plotted with
result from the the Fourier-Chebyshev based branch continuation algorithm. The
clear agreement we observe in the vicinity of the critical point confirms the
existence of a stable supercritical state, and validates the weakly nonlinear
predictions along the lower branch of the neutral curve.

Figure 5: Bifurcation diagrams at point $\triangle$. The $L^{2}$ norms are
taken over the whole domain $\Omega=[0,2\pi/k]\times[-1,1]$. Figure 6: Real (
) and imaginary ( ) parts of the unstable eigenfunction
$\mbox{\boldmath$\varphi$}_{(1,1)}$ at the point $\triangle$. (Left) Axial
(streamwise) velocity $u_{(1,1),x}$ and vertical velocity $u_{(1,1),y}$.
(Middle) All four nonzero components of $\mathbfcal G_{(1,1)}\in
T_{\mathbf{I}}\mathrm{Pos}(3)$, the tangent form of the polymer strain
perturbation tensor. (Right) All four nonzero components of
$\mathbf{C}_{(1,1)}\in\mathrm{Pos}(3)$, the corresponding fluctuation tensor
from a standard expansion. Figure 7: The nonlinear mean correction
$\mbox{\boldmath$\varphi$}_{(2,0)}$ at the point $\triangle$. (Left) Axial
(streamwise) velocity $u_{(2,0),x}$ and vertical velocity $u_{(2,0),y}=0$.
(Middle) All four nonzero components of $\mathbfcal G_{(2,0)}\in
T_{\mathbf{I}}\mathrm{Pos}(3)$, the mean correction to the conformation tensor
in its tangent form. (Right) All four nonzero components of
$\mathbf{C}_{(2,0)}\in\mathrm{Pos}(3)$, the corresponding tensor from a
standard expansion. Figure 8: Real ( ) and imaginary ( ) parts of the
nonlinear correction $\mbox{\boldmath$\varphi$}_{(2,2)}$ at the point
$\triangle$. (Left) Axial (streamwise) velocity $u_{(2,2),x}$ and vertical
velocity $u_{(2,2),y}$. (Middle) All four nonzero components of $\mathbfcal
G_{(2,2)}\in T_{\mathbf{I}}\mathrm{Pos}(3)$. (Right) All four nonzero
components of $\mathbf{C}_{(2,2)}\in\mathrm{Pos}(3)$, the corresponding tensor
from a standard expansion.

###  4.2 Flow and polymer field prediction

The various flow and polymer fields generated as part of the weakly nonlinear
analysis can be used to generate an approximation to the solution near to a
bifurcation point. The structure of the critical eigenfunction at the
$\triangle$ in Figure 2 is shown in Figure 6. The flow and conformation tensor
structures are familiar from previous studies (Garg et al., 2018; Khalid et
al., 2021a) whereas the Cauchy-Green perturbation tensor $\mathbfcal
G_{(1,1)}$ hasn’t been shown before. Figure 6 shows that all components of
$\mathbfcal G_{(1,1)}$ are confined to the centerline of the channel. For
instance, ${\mathcal{G}}_{(1,1),xx}$ only develops a noticeable magnitude
above $y=-0.2$. On the other hand, $C_{(1,1),xx}$ indicates that the
streamwise normal stretch reaches its maximum towards the bottom of the
channel. This difference is explained by the shape of the laminar base state
(cf. Figure 1). $\mathbf{C}_{b}$ is smaller near the centerline, thus
computing $\mathbfcal
G_{(1,1)}=\mathbf{F}_{b}^{-1}\mathbf{C}_{(1,1)}\mathbf{F}_{b}^{-T}$ amplifies
changes in that region, i.e., $\mathbfcal G_{(1,1)}$ recognizes deformations
that are large relative to $\mathbf{C}_{b}$. Again, this is an immediate
consequence of the fact that the Riemannian metric on $\mathrm{Pos}(3)$
depends on the base point $\mathbf{C}_{b}$. Physically, the new formulation
highlights that the polymeric disturbance caused by the centre mode
instability is confined to a small layer around the centerline, which would
not be immediate from a standard expansion in $\mathbf{C}$ (cf. right panel of
Figure 6 or Figure 17 in Khalid et al. (2021a)).

Higher order disturbances are more difficult to interpret on
$\mathrm{Pos}(3)$, but up to $O(\varepsilon^{2})$ can still be thought of as
consecutive geodesic perturbations (Hameduddin et al., 2019). The
$O(\varepsilon^{2})$ terms from the weakly nonlinear expansion are given in
Figure 7, which displays the first nonlinear mean correction
$\mbox{\boldmath$\varphi$}_{(2,0)}$, and Figure 8, which shows
$\mbox{\boldmath$\varphi$}_{(2,2)}$. With these fields known, the full flow
state can be approximated by evaluating the weakly nonlinear expansion (14) up
to second order, including $|A|$ in the shape functions as necessary. This low
order approximation is compared with a full state from the continuation tool
in Figure 9 at the point $\Diamond$ on the supercritical bifurcation branch
(see Figure 5).

Figure 9: Comparison of the supercritical state at point $\Diamond$
(identified in Figure 5) as predicted by the weakly nonlinear analysis (left
panel) and branch continuation (right panel) techniques. Contours show the
geodesic distance between the base and full states
$d(C_{b},C)=\sqrt{\mathrm{tr}\,\mathcal{G}^{2}}$, the lines correspond to the
perturbation streamfunction.

###  4.3 $\beta=0.9$, $L_{max}=100$ & $Sc=10^{6}$

In this subsection, we reduce $L_{max}$ to $100$ to explore less extensible
(more realistic) polymers and reintroduce the conformation tensor diffusion
term into the governing equations (1c) by considering a finite Schmidt number,
$Sc=10^{6}$. Figure 10 shows the corresponding marginal stability curve
complemented with a finite amplitude curve from the weakly nonlinear analysis.
The $L_{max}=500$ neutral curve is also displayed for comparison in bright
grey. All visible changes are caused by the adjustment of $L_{max}$: the
introduction of finite $Sc$ alone has no visual effect. The key observation
from Figure 10 is that reducing $L_{max}$ shifts the neutral curve down in
$Re$, and reduces the slope of the lower branch. In particular, lowering
$L_{max}$ has a destabilizing effect in the elastic regime (low Reynolds
numbers). This counter-intuitive finding is the primary motivation for
examining the $Re=0$ instability recently found by Khalid et al. (2021b) at
finite $L_{max}$ in §4.5.

Figure 10: (Left) Neutral curve corresponding to marginal linear stability at
$\beta=0.9$, $L_{max}=100$, $Sc=10^{6}$. Results of the weakly nonlinear
analysis are shown in the form of a curve at steady state amplitude $|A|=0.4$.
The $L_{max}=500$ neutral curve from Figure 2 is also shown for comparison.
(Right) The development of the critical wave number, $k_{crit}$, along the
neutral curve (corresponding curve for $L_{max}=500$ again shown in grey).

###  4.4 Energy analysis

We now examine the energetic contributions of the different terms in the
equations (1) in order to examine the mechanisms driving the centre mode
instability. This approach has proved useful to diagnose the character of
instabilities - for example Joo & Shaqfeh (1991, 1992) identified purely
elastic instabilities in curved channel flows with this procedure (see also
Zhang et al. (2013) and Agarwal et al. (2014)). Taking an $L^{2}$ inner
product of the momentum equations at $O(\varepsilon)$ and the disturbance
velocity field
$\mathbf{u}_{(1,1)}=(\mbox{\boldmath$\varphi$}_{(1,1),1},\mbox{\boldmath$\varphi$}_{(1,1),2})$
gives (for more details see e.g. Zhang et al., 2013) the disturbance kinetic
energy equation

$\partial_{t}E:=\frac{1}{2}\partial_{t}\|\mathbf{u}_{(1,1)}\|_{L^{2}}^{2}=\mathscr{P}+\mathscr{E}+\mathscr{W},$
(19)

where

$\mathscr{P}:=-\frac{1}{2}\left\langle\mbox{\boldmath$\nabla$}\mathbf{u}_{b},\mathbf{u}_{(1,1)}\otimes\bar{\mathbf{u}}_{(1,1)}+\bar{\mathbf{u}}_{(1,1)}\otimes\mathbf{u}_{(1,1)}\right\rangle_{L^{2}}$

($\bar{\mathbf{u}}$ is the complex conjugate of $\mathbf{u}$) is the
disturbance kinetic energy production due to the underlying shear of
$\mathbf{u}_{b}$,

$\mathscr{E}:=-\frac{\beta}{Re}\|\mbox{\boldmath$\nabla$}\mathbf{u}_{(1,1)}\|_{L^{2}}^{2}$

represents the viscous dissipation and is strictly negative, and

$\mathscr{W}:=-\frac{(1-\beta)}{2Re}\left(\left\langle\mbox{\boldmath$\nabla$}\mathbf{u}_{(1,1)},\mathbf{T}_{(1,1)}\right\rangle_{L^{2}}+\left\langle\mathbf{T}_{(1,1)},\mbox{\boldmath$\nabla$}\mathbf{u}_{(1,1)}\right\rangle_{L^{2}}\right)$

indicates the rate of work done on the fluid by the polymeric stresses, with

$\mathbf{T}_{(1,1)}:=D\mathbf{T}(\mathbf{C}_{b})\left[\mathbf{F}_{b}\mathbfcal
G_{(1,1)}\mathbf{F}_{b}^{T}\right].$

Extending this procedure to identify the mechanisms behind the growth of
elastic energy stored in the polymer is well known to be problematic (Doering
et al., 2006). The underlying issue is that the elastic potential energy,
which is a function of $\mathrm{tr}\,\mathbf{C}$, does not correspond to a
norm in the obvious fashion that the kinetic energy does. Once again, this
essentially comes down to the fact that the set $\mathrm{Pos}(3)$ does not
constitute a linear vector space, and there is no notion of norm available.
This may be overcome by measuring disturbances in $\mathbf{C}$ along geodesics
in $\mathrm{Pos}(3)$, according to the metric induced by the Riemannian
structure. The work of Hameduddin et al. (2018) suggests that

$\left(d(\mathbf{C}_{b},\mathbf{C})\right)^{2}=\left(d({\bf
I},\mathbf{G})\right)^{2}=\mathrm{tr}\left(\mathbfcal G^{H}\mathbfcal
G\right),$

which immediately gives us a way of quantifying the evolution of polymer
disturbances as

$J:=\|d(\mathbf{C}_{b},\mathbf{C})\|_{L^{2}}^{2}=\int_{[-1,1]}\mathrm{tr}\left(\mathbfcal
G^{H}(y)\mathbfcal G(y)\right)\,dy,$ (20)

a formulation which was originally proposed in Hameduddin et al. (2019). This,
in fact, is the main advantage of relying on the alternative formulation of
the governing equations given in equation (8). This newly defined quantity $J$
in (20) is equal to $\|\mathbfcal G\|_{L^{2}}^{2}$ which is a natural
generalization of the kinetic energy from (19).

Adopting this polymer energy measure $J$, an energetic evolution equation for
the polymer disturbances can now be obtained by taking an $L^{2}$ inner
product of $\mathbfcal G_{(1,1)}$ with the linearized disturbance equation (in
a symmetric fashion) to obtain:

$\partial_{t}J=\mathscr{A}_{b}+\mathscr{A}_{1}+\mathscr{T}+\mathscr{E}_{p},$
(21)

where

$\mathscr{A}_{b}:=\left\langle\mathbfcal
G_{(1,1)},2\mathrm{sym}\left(\mathbfcal
G_{(1,1)}h(\mathbf{u}_{b})\right)\right\rangle_{L^{2}}+\left\langle
2\mathrm{sym}\left(\mathbfcal G_{(1,1)}h(\mathbf{u}_{b})\right),\mathbfcal
G_{(1,1)}\right\rangle_{L^{2}}$

represents the contribution due to the base velocity field,

$\mathscr{A}_{1}:=\left\langle\mathbfcal
G_{(1,1)},2\mathrm{sym}\left(h\left(\mathbf{u}_{(1,1)}\right)\right)\right\rangle_{L^{2}}+\left\langle
2\mathrm{sym}\left(h\left(\mathbf{u}_{(1,1)}\right)\right),\mathbfcal
G_{(1,1)}\right\rangle_{L^{2}}$

is the corresponding term capturing the effect of the disturbance velocity
field $\mathbf{u}_{(1,1)}$,

$\mathscr{T}:=-\left\langle\mathbfcal
G_{(1,1)},\mathbf{F}_{b}^{-1}\mathbf{T}_{(1,1)}\mathbf{F}_{b}^{-T}\right\rangle_{L^{2}}-\left\langle\mathbf{F}_{b}^{-1}\mathbf{T}_{(1,1)}\mathbf{F}_{b}^{-T},\mathbfcal
G_{(1,1)}\right\rangle_{L^{2}}$

is the polymeric relaxation term and

$\mathscr{E}_{p}:=\left\langle\mathbfcal
G_{(1,1)},\frac{1}{ReSc}\mathbf{F}_{b}^{-1}\Delta\left(\mathbf{F}_{b}\mathbfcal
G_{(1,1)}\mathbf{F}_{b}^{T}\right)\mathbf{F}_{b}^{-T}\right\rangle_{L^{2}}+\left\langle\frac{1}{ReSc}\mathbf{F}_{b}^{-1}\Delta\left(\mathbf{F}_{b}\mathbfcal
G_{(1,1)}\mathbf{F}_{b}^{T}\right)\mathbf{F}_{b}^{-T},\mathbfcal
G_{(1,1)}\right\rangle_{L^{2}}$

is the polymeric diffusion contribution.

The contribution of each individual term along the neutral curve of subsection
4.1 (parameterised by the wavenumber $k_{crit}$ which varies monotonically
along the curve) is shown in Figure 11 for both the kinetic energy equation
(19) (left) and polymer ‘energy’ equation (21) (right). Based on the recent
discovery of an inertialess linear instability that stems from the lower
branch of the neutral curve (Khalid et al., 2021b), it was anticipated that
the underlying destabilizing effects would be elastically driven along this
branch. This is exactly what is seen: the polymer stress term is the sole
energising term for the disturbance kinetic energy. Figure 11, however,
indicates that this holds over the upper branch as well so that the centre-
mode instability remains purely elastic - i.e., the rate of polymer work
$\mathscr{W}$ is the only positive contribution to $\partial_{t}E$ \-
throughout the entirety of the neutral curve shown. Not even at $Re=3000$ do
we have a positive contribution from the turbulence production term,
$\mathscr{P}$, which is the term that represents inertial effects and is
responsible for the onset of instability in Newtonian turbulence. In inertia-
dominated flows, $\mathscr{P}$ is the primary cause of turbulent kinetic
energy growth Zhang et al. (2013).

In the $J$ equation, the base flow ($\mathscr{A}_{b}$) barely contributes so
that the effect of polymeric relaxation processes, $\mathscr{T}$, is balanced
by the input of the perturbation velocity field through $\mathscr{A}_{1}$
(${\mathscr{E}}_{P}=0$ as $Sc\to\infty$ and so is not plotted). The dominance
of $\mathscr{A}_{1}$ which relies on the base polymer stretch rather than any
base flow shear confirms that the instability mechanism is entirely elastic.

Choosing large but finite $Sc$ does not change this conclusion. Figure 12
shows the energy analysis results for the neutral curve at $L_{max}=100$ in
Figure 10 of §4.3. Again, the polymeric viscous dissipation term,
$\mathscr{E}_{p}$, does not contribute to the growth of $J$
(${\mathscr{E}}_{p}$ only starts to become significant for $Sc\sim O(10^{2})$)
and the energy source for the instability is solely elastic.

Figure 11: Energy analysis results across the $L_{max}=500$ neutral curve
shown in Figure 2. (Left) Components contributing to the production of the
turbulent kinetic energy, $E$. (Right) Components contributing towards the
evolution of the polymeric disturbance, $J$. All values are normalized across
the neutral curve. Figure 12: Energy analysis results across the $L_{max}=100$
neutral curve shown in Figure 10. (Left) Components contributing to the
production of the turbulent kinetic energy, $E$. (Right) Components
contributing towards the evolution of the polymeric disturbance, $J$. All
values are normalized across the neutral curve.

###  4.5 Inertialess limit

In this section we explore the low-$Re$ elastic limit of the centre mode
instability motivated by the finding in §4.3 that decreasing $L_{max}$ makes
the instability move to lower $Re$. Recent work (Khalid et al., 2021b) has
found the centre mode instability for $Re=0$ in the Oldroyd-B model, albeit at
very high $Wi$ and very small $(1-\beta)$ i.e. the dilute limit. Our aim here
is to see if we can find this instability at a lower, more realistic $Wi$ by
varying $L_{max}$ in the FENE-P model. The effect of the viscosity ratio,
$\beta$, for $L_{max}\to\infty$ (an Oldroyd-B fluid) is already known (see
inset (B) of Figure 2 in Khalid et al., 2021b). The instability first appears
at $\beta=0.9905$ with the critical $Wi$ decreasing as $\beta$ increases to
0.994, reaching a minimum of $Wi\approx 649$ (note their value
$Wi^{\prime}=973.8$ is defined using the base centreline speed) and then
increases again as $\beta$ continues to increase beyond 0.994 towards 1. Thus,
the lowest $\beta$ for which the $Re=0$ instability still exists (limited by
the slope of the lower branch on the marginal curve) could also be decreased
if the threshold $Wi$ for instability is decreased through adjusting
$L_{max}$. This is what we find: see figure 13, which shows that instability
at $Re=0$ is possible at just over $Wi=100$ for $\beta=0.98$ and
$L_{max}=100$. The finite amplitude curves generated by weakly nonlinear
analysis and shown in figure 13 further imply the existence of an unstable
subcritical state in this inertialess regime. That is, the flow continues to
be nonlinearly unstable when lowering $Wi$ below the threshold for linear
instability.

Figure 13: Neutrally stable curves ( ) around the inertialess ($Re=0$) limit
for ultra-dilute polymer solutions at $\beta=0.99$ (left) and $\beta=0.98$
(right). The dashed lines ( ) are finite amplitude curves that show the
nonlinear behaviour indicated by the weakly nonlinear analysis.

Figure 13 suggests further reduction in the threshold $Wi$ for instability may
be possible by making $L_{max}$ even smaller. Neutral curves in the $Wi-\beta$
plane at $Re=0$ for $L_{max}=40,70$ and $100$ are shown in Figure 14 along
with the concomitant finite amplitude curves. Two important features are
evident from this figure. Firstly, the destabilizing effect of $L_{max}$ has a
limit, which appears to be in the interval $L_{max}\in[40,100]$ for $Re=0$.
Secondly, the weakly nonlinear analysis indicates that the bifurcation is
subcritical with respect to $\beta$, except for high $Wi$ along the lower
branch (in the $Wi-\beta$ plane) of the neutral curve where it becomes
supercritical.

Figure 14: (Left) Neutrally stable curves ( ) at the inertialess limit $Re=0$
for ultra-dilute polymer solutions. The dashed lines ( ) are finite amplitude
curves that show the nonlinear behaviour indicated by the weakly nonlinear
analysis. (Right) Changes in the critical wave number, $k_{crit}$, as the
neutral curves are traversed.

The results of an energy budget analysis are shown in Figure 15 for this
$Re=0$ instability at $Wi=115$ and $\beta=0.98$ \- the $\bigtriangledown$ in
Figure 14 \- as a function of $L_{max}$. The kinetic energy evolution equation
(19) is unable to handle the vanishing $Re$ situation and so we exclusively
focus on the budget in $J$, the measure introduced for polymeric
perturbations. Figure 15 tracks how the disturbance growth rate,
$\partial_{t}J$, and each term contributing to it changes as $L_{max}$ is
varied at point $\bigtriangledown$ (cf. Figure 14) ($\partial_{t}J=0$ indicate
points on the neutral curve e.g. there is no instability at $L_{max}=100$ at
$\bigtriangledown$). The contribution stemming from the base flow,
$\mathscr{A}_{b}$, is still negligible, which indicates that stability is
determined by the balance between (destabilizing) $\mathscr{A}_{1}$ and
(stabilizing) $\mathscr{T}$. The dissipation rate associated with polymeric
relaxation processes, $\mathscr{T}$, becomes increasing negative as $L_{max}$
is decreased, ultimately causing stabilisation. As expected from Figure 14, an
optimal $L_{max}$ exists ($\approx 60$) for this particular pairing of $Wi$
and $\beta$. That it exists at all - i.e. the FENE-P model is more unstable
than the Oldroyd-B model to this inertialess centre mode instability - is a
surprise.

Figure 15: The energy budget for the polymeric disturbance, $J$, at the
inertialess limit, at point $\bigtriangledown$ ($Wi=115$, $\beta=0.98$) in
Figure 14. Note that the scale for $\partial_{t}J$ (left axis) is enlarged to
improve visibility.

##  5 Discussion

In this paper, we have considered the character of the bifurcation of a
recently-discovered centre-mode (Garg et al., 2018; Khalid et al., 2021a) in
rectilinear viscoelastic channel flow for large $Re=O(10^{3})$ down to the
inertialess limit of $Re=0$. Using weakly nonlinear analysis within a formal
framework which respects the positive definiteness of the conformation tensor
$\mathbf{C}$ (Hameduddin et al., 2018, 2019), we find that the subcriticality
found by Page et al. (2020) for one point of the neutral curve at
$L_{max}=500$ is generic across the neutral curve and for different $L_{max}$.
Supercriticality is only found at large $Wi$ on the ‘lower’ (low-$Re$) branch
of the neutral curve in the $(Wi,Re)$ plane otherwise the branch of travelling
waves arising from the neutral curve reach down to lower $Wi$ and the region
where EIT is found. In this extended region of parameter space, the base flow
is nonlinearly unstable to disturbances of sufficient amplitude. The threshold
amplitude to trigger this instability is determined by the minimal amplitude
of approach of the stable manifold of the lower branch of travelling waves to
the base flow. This is bounded from above by the amplitude of the lower branch
itself and the one branch-tracking calculation done so far (see Figure 3 in
Page et al., 2020) indicates that this is small: the volume-averaged
tr$\mathbf{C}$ of the travelling wave solutions stays within 5% of the base
flow value even when the $Wi$ is reduced to 50% of its value at the
bifurcation. Hence, for practical purposes, the base flow may well appear
linearly unstable below the neutral curve in $Wi$ (recent experiments suggest
a similar situation in $Re$ (Choueiri et al., 2021)). Assessing how far this
situation continues as $Wi$ is decreased requires, of course, a full branch
continuation procedure to map out the surface of travelling wave solutions.

By using a FENE-P fluid we have also confirmed that the centre-mode
instability persists for maximum polymer extension down to $L_{max}=40$ at
least. Somewhat counterintuitively, the introduction of finite $L_{max}$ is
found to move the neutral curve closer to the inertialess $Re=0$ limit at
fixed $\beta$. Pursuing this further by entering the dilute ($\beta\rightarrow
1$) limit, we also find that finite $L_{max}$ can bring the linear instability
recently found by Khalid et al. (2021b) down to more physically-relevant
$Wi\gtrsim 110$ at $\beta=0.98$ compared with their threshold of $Wi\approx
649$ (based on the bulk velocity) at $\beta=0.994$ for $L_{max}\to\infty$.
Again the instability is subcritical implying that inertialess rectilinear
viscoelastic shear flow is nonlinearly unstable for even lower $Wi$. Assessing
exactly how low again requires locating the saddle node (turning point) of the
travelling waves as $Wi$ decreases which requires a branch continuation code.

Finally, by considering the various energy terms in the disturbance kinetic
energy equation, we have found that the centre-mode instability is purely
elastic in origin even for $Re=O(10^{3})$, rather than ‘elasto-inertial’, as
the underlying shear does not energise the instability. This finding is
consistent with the recent smooth connection found by Khalid et al. (2021b) to
an entirely elastic instability at $Re=0$ and suggests that EIT and ET may
indeed be two different extremes of the same whole. Given that this
instability is being suggested as the origin of EIT (Garg et al., 2018; Page
et al., 2020; Chaudhary et al., 2021; Khalid et al., 2021a), the importance of
inertia must emerge at finite amplitude and is perhaps already there in the
travelling wave solutions especially when they establish their ‘arrowhead’
form familiar from DNS at higher amplitudes and lower $Wi$ (Dubief et al.,
2020).

In terms of experiments, the centre mode instability has recently been
investigated in both channel (Schnapp & Steinberg, 2021) and pipe flow
(Choueiri et al., 2021). In a pipe, (Choueiri et al., 2021) observed evidence
of the centre mode instability at high $Wi=O(100)$ and low (subcritical) $Re$.
More relevant to the current results are the essentially inertialess
($Re\lesssim 0.3$) channel flow experiments of Schnapp & Steinberg (2021),
which were conducted at very high $Wi\in(100,1700]$. Finite amplitude
traveling waves (or ‘elastic’ waves in their terminology) were triggered by
‘small’ disturbances – in contrast to the ‘large’ disturbances used in Pan et
al. (2013) for $Re\lesssim 0.01$ and $Wi\lesssim 10$. Interestingly for the
calculations performed here, they estimate the presence of a linear
instability at $Wi=125\pm 25$. However, both studies were performed at
considerably lower values of $\beta$ than those studied in the bulk of this
paper ($\beta=0.74$ in Schnapp & Steinberg (2021) and $\beta=0.56$ in Choueiri
et al. (2021)). We have examined both of these solvent viscosities in Appendix
B and find that the significant reduction in $\beta$ leads to both (i) a
smaller unstable region in the $Wi$-$Re$ plane and (ii) almost uniformly
supercritical behaviour around the neutral curve. This does not preclude the
possibility that the branch may bend back down towards lower $Re$ and $Wi$,
which cannot be captured in our third order weakly nonlinear analysis but
which can be studied by branch continuation of the travelling waves.

The obvious next steps after the analysis described here – and particularly
important in the context of the experimental observations at low $\beta$ – is
to employ a branch continuation procedure to track the travelling waves
produced by the centre-mode instability to finite amplitudes and then to
explore where they exist in parameter space. The inertialess limit is perhaps
the most interesting but hardest to access numerically. These travelling
waves, of course, provide their own launchpad for further (secondary)
bifurcations from which subsequent solutions then suffer tertiary bifurcations
and so forth. Establishing that this bifurcation cascade occurs precisely
where EIT is observed in parameter space would provide convincing evidence of
the importance of the centre-mode instability. We hope to report on further
progress in this direction in the near future.

Acknowledgements: GB gratefully acknowledges the support of the Harding
Foundation through a PhD scholarship
(https://www.hardingscholars.fund.cam.ac.uk).

##  Appendix A Numerical methods

The eigenvalue problem (15a) and the subsequent nonlinear equations in (15)
were solved using a Chebyshev discretization. Exploiting the symmetries of the
centre eigenmode, the expansions were performed over half the channel width,
$y\in[-1,0]$, with appropriate boundary conditions to enforce the symmetry of
$u_{x}$, antisymmetry of $u_{y}$ and appropriate symmetries for the various
components of $\mathbf{C}$. This approach crucially concentrates the
collocation points near both the channel boundary and the centreline where the
eigenmode is localised so that manageable truncations prove adequate. For the
$\beta=0.9$ neutral curves, $200$ Chebyshev were sufficient while higher
$\beta$ values needed $300$-$400$ Chebyshev modes due to the increasing
localization of the unstable eigenmode (see Khalid et al. (2021b)). The
neutral curves were obtained using a continuation scheme that relies on the
tangent that the weakly nonlinear analysis yields. Specifically, in the $Wi-
Re$ plane, this is given by substituting $|A|=0$ into (18):

$\frac{Re_{1}}{Wi_{1}}=-\frac{\mathrm{Im}(\bar{d}b)}{\mathrm{Im}(\bar{d}a)}.$

In solving the eigenvalue problem, a shift-inverse spectral transformation
Meerbergen et al. (1994) was employed, using the eigenvalue at the previous
continuation step, to isolate the critical eigenmode. The unstable mode was
then obtained via standard power iteration. All results were cross-checked
using two grid resolutions.

Results of the weakly nonlinear analysis were validated by an independently-
developed branch continuation routine. In this, the flow solution is assumed
to be steady in an appropriately-chosen Galilean frame (i.e. a travelling
wave) which allows the time derivatives to be replaced by a spatial derivative
in $x$ premultiplied by an a priori unknown phase speed $c=\omega_{r}/k$. The
governing equations are then discretized in space using Fourier modes in $x$
and Chebyshev modes in $y$ across the domain $(x,y)\in[0,2\pi/k]\times[-1,1]$
to leave a high-dimensional - typically $O(10^{5})$ degrees of freedom -
nonlinear system of equations for the expansion coefficients. A good starting
guess for the solution and $c$ can be generated near the neutral curve and
then the solver propagates along the solution surface via a pseudo-arclength
continuation algorithm based on a Newton-Raphson iterative scheme (e.g.
Dijkstra et al. (2014)). Simulations for the curves appearing in Figures 4 and
5 were run at $80$ Chebyshev and $40$ Fourier modes (so $7\times 40\times
80\times 2=44,800$ real degrees of freedom). Resolution independence was
carefully checked at the terminal point of each branch shown (using up to
$80,000$ degrees of freedom). In this paper, the branch continuation code was
only used to confirm the weakly nonlinear analysis. A future report will
describe it in detail when the results of using it to explore solution
morphology a finite distance from the neutral curve will be presented.

##  Appendix B Results at moderate $\beta$

Figure 16:  Neutrally stable curves ( ) for low solvent viscosities
$\beta=0.74$ (left) and $\beta=0.56$ (right). The dashed lines ( ) are finite
amplitude curves that show the nonlinear behaviour indicated by the weakly
nonlinear analysis.

Motivated by recent experimental results (Choueiri et al., 2021; Schnapp &
Steinberg, 2021) at higher polymer concentrations, we briefly discuss the
impact of reducing $\beta$ on both the linear instability and the predictions
of our weakly nonlinear analysis. We consider two solvent viscosities,
$\beta=0.74$ and $\beta=0.56$, which match the values obtained in Schnapp &
Steinberg (2021) and Choueiri et al. (2021) respectively (note the latter
study was done in a pipe precluding any direct comparison here). Neutral
curves and the weakly nonlinear results are reported in figure 16 for both
Oldroyd-B fluids and FENE-P fluid with relatively high $L_{max}$. The
reduction in $\beta$ noticably shrinks the region of instability in the
$Wi$-$Re$ plane, notably bending the lower part of the curve - which connects
to $Re=0$ at high $\beta$ \- upwards. Moreover, in contrast to the dilute
($\beta\geq 0.9$) results in the bulk of this paper, the introduction of
finite extensibility has a uniformly stabilizing effect – this behaviour is
perhaps more typical of the more realistic polymer model; in many cases the
reduction in the base-state normal stress tends to suppress more ‘interesting’
Oldroyd-B results (e.g. see the linear analyses in Ray & Zaki, 2014; Page &
Zaki, 2015).

In addition, the weakly nonlinear results (dashed lines in figure 16) indicate
almost uniformly supercritical behaviour around the neutral curve (note the
small exception at high $Wi$ for $\beta=0.74$ and $L_{max}=250$). This finding
should be contrasted to the recent experimental results at extreme $Wi\geq
100$ of Schnapp & Steinberg (2021), who have observed finite amplitude
traveling waves at very low $Re$ at $\beta=0.74$, and motivates further study
via branch continuation of exactly where nonlinear traveling waves are
predicted to exist in the parameter space.

## References

  * Agarwal et al. (2014) Agarwal, A., Brandt, L. & Zaki, T. A. 2014 Linear and nonlinear evolution of a localized disturbance in polymeric channel flow. Journal of Fluid Mechanics 760, 278–303.
  * Chandra et al. (2018) Chandra, B., Shankar, V. & Das, D. 2018 Onset of transition in the flow of polymer solutions in microtubes. Journal of Fluid Mechanics 844, 1052–1083.
  * Chaudhary et al. (2019) Chaudhary, Indersh, Garg, Piyush, Shankar, V. & Subramanian, Ganesh 2019 Elasto-inertial wall mode instabilities in viscoelastic plane poiseuille flow. Journal of Fluid Mechanics 881, 119–163.
  * Chaudhary et al. (2021) Chaudhary, I., Garg, P., Subramanian, G. & Shankar, V. 2021 Linear instability of viscoelastic pipe flow. Journal of Fluid Mechanics 908, A11.
  * Choueiri et al. (2018) Choueiri, G. H., Lopez, J. M. & Hof, B. 2018 Exceeding the Asymptotic Limit of Polymer Drag Reduction§. Physical Review Letters 120, 124501.
  * Choueiri et al. (2021) Choueiri, G. H., Lopez, J. M., Varshey, A., Sankar, S. & Hof, B. 2021 Experimental observation of the origin and structure of elasto-inertial turbulence. arXiv:2103.00023 .
  * Dijkstra et al. (2014) Dijkstra, Henk A., Wubs, Fred W., Cliffe, Andrew K., Doedel, Eusebius, Dragomirescu, Ioana F., Eckhardt, Bruno, Gelfgat, Alexander Yu., Hazel, Andrew L., Lucarini, Valerio, Salinger, Andy G. & et al. 2014 Numerical bifurcation methods and their application to fluid dynamics: Analysis beyond simulation. Communications in Computational Physics 15 (1), 1–45.
  * Doering et al. (2006) Doering, C. R., Eckhardt, B. & Schumacher, J. 2006 Failure of energy stability in Oldroyd-B fluids at arbitrarily low Reynolds numbers. Journal of Non-Newtonian Fluid Mechanics 135, 92–96.
  * Draad et al. (1998) Draad, A. A., Kuiken, G. D. C. & Nieuwstadt, F. T. M. 1998 Laminar-turbulent transition in pipe flow for newtonian and non-newtonian fluids. Journal of Fluid Mechanics 377, 267–312.
  * Dubief et al. (2020) Dubief, Y., Page, J., Kerswell, R. R., Terrapon, V. E. & Steinberg, V. 2020 A first coherent structure in elasto-inertial turbulence. arXiv 2006.06770 .
  * Dubief et al. (2013) Dubief, Y., Terrapon, V. E. & Soria, J. 2013 On the mechanism of elasto-inertial turbulence. Physics of Fluids 25 (11), 110817.
  * Garg et al. (2018) Garg, P., Chaudhary, I., Khalid, M., Shankar, V. & Subramanian, G. 2018 Viscoelastic Pipe Flow is Linearly Unstable. Physical Review Letters 121, 024502.
  * Goldstein et al. (1969) Goldstein, R. J., Adrian, R. J. & Kreid, D. K. 1969 Turbulent and transition pipe flow of dliye aqueous polymer solutions. Industrial and Engineering Chemistry Fundamentals 8, 498.
  * Graham (2014) Graham, M. D. 2014 Drag reduction and the dynamics of turbulence in simple and complex fluids. Physics of Fluids 26, 101301.
  * Groisman & Steinberg (2000) Groisman, A. & Steinberg, V. 2000 Elastic turbulence in a polymer solution flow. Nature 405, 53–55.
  * Hameduddin et al. (2019) Hameduddin, Ismail, Gayme, Dennice F. & Zaki, Tamer A. 2019 Perturbative expansions of the conformation tensor in viscoelastic flows. Journal of Fluid Mechanics 858, 377–406.
  * Hameduddin et al. (2018) Hameduddin, Ismail, Meneveau, Charles, Zaki, Tamer A. & Gayme, Dennice F. 2018 Geometric decomposition of the conformation tensor in viscoelastic turbulence. Journal of Fluid Mechanics 842, 395–427.
  * Hansen & Little (1974) Hansen, R. J. & Little, R. C. 1974 Early turbulence and drag reduction phenomena in larger pipes. Nature 252, 690.
  * Jones & Maddock (1966) Jones, W. & Maddock, J. L. 1966 Onset of instabilities and reduction of drag in flow of relaxing liquids through tubes and porous beds. Nature 212, 388.
  * Joo & Shaqfeh (1991) Joo, Yong Lak & Shaqfeh, Eric S. G. 1991 Viscoelastic poiseuille flow through a curved channel: A new elastic instability. Physics of Fluids A: Fluid Dynamics 3 (7), 1691–1694.
  * Joo & Shaqfeh (1992) Joo, Y. L. & Shaqfeh, E. S. G. 1992 A purely elastic instability in Dean and Taylor-Dean flow. Physics of Fluids 4, 524.
  * 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, A43.
  * Khalid et al. (2021b) Khalid, Mohammad, Shankar, V. & Subramanian, Ganesh 2021b A continuous pathway between the elasto-inertial and elastic turbulent states in viscoelastic channel flow. arXiv:2103.06794 .
  * Larson et al. (1990) Larson, R. G., Shaqfeh, E. S. G. & Muller, S. J. 1990 A purely elastic instability in Taylor-Couette flow. Journal of Fluid Mechanics 218, 573–600.
  * Lopez et al. (2019) Lopez, J. M., Choueiri, George H. & Hof, Bjoern 2019 Dynamics of viscoelastic pipe flow at low reynolds numbers in the maximum drag reduction limit. Journal of Fluid Mechanics 874, 699–719.
  * Lumley (1969) Lumley, J. L. 1969 Drag reduction by additives. Annual Review of Fluid Mechanics 656 (33), 367–384.
  * Meerbergen et al. (1994) Meerbergen, Karl, Spence, A & Roose, Dirk 1994 Shift-invert and cayley transforms for detection of rightmost eigenvalues of nonsymmetric matrices. Bit Numerical Mathematics 34 (3), 409–423.
  * Meulenbroek et al. (2004) Meulenbroek, B., Storm, C., Morozov, A. N. & van Saarloos, W. 2004 Weakly nonlinear subcritical instability of visco-elastic Poiseuille flow. Journal of Non-Newtonian Fluid Mechanics 116 (2-3), 235–268.
  * Morozov & Saarloos (2007) Morozov, A. N. & Saarloos, W. V. 2007 An introductory essay on subcritical instabilities and the transition to turbulence in visco-elastic parallel shear flows. Physics Reports 447, 112–143.
  * Page et al. (2020) Page, Jacob, Dubief, Yves & Kerswell, Rich R. 2020 Exact traveling wave solutions in viscoelastic channel flow. Phys. Rev. Lett. 125, 154501\.
  * Page & Zaki (2015) Page, J. & Zaki, T. A. 2015 The dynamics of spanwise vorticity perturbations in homogeneous viscoelastic shear flow. Journal of Fluid Mechanics 777, 327–363.
  * Pan et al. (2013) Pan, L., Morozov, A., Wagner, C. & Arratia, P. E. 2013 Nonlinear elastic instability in channel flows at low Reynolds numbers. Physical Review Letters 110, 174502.
  * Procaccia et al. (2008) Procaccia, I., Lvov, Victor & Benzi, Roberto 2008 Colloquium: Theory of drag reduction by polymers in wall-bounded turbulence. Reviews of Modern Physics 1, 225–247.
  * Qin et al. (2019) Qin, B. Y., Salipante, P. F., Hudson, S. D. & Arratia, P. E. 2019 Flow resistance and structure in viscoelastic channel flows at low ‘re’. Physical Review Letters 123, 194501.
  * Ray & Zaki (2014) Ray, P. K. & Zaki, T. A. 2014 Absolute instability in viscoelastic mixing layers. Physics of Fluids 26 (1), 014103.
  * Samanta et al. (2013) Samanta, D. S., Dubief, Y., Holzner, H., Schäfer, C., Morozov, A. N., Wagner, C. & Hof, B. 2013 Elasto-inertial turbulence. Proceedings of the National Academy of Sciences of the United States of America 110, 10557–10562.
  * Schnapp & Steinberg (2021) Schnapp, R. & Steinberg, V. 2021 Elastic waves above elastically driven instabiity in weakly perturbed channel flow. arXiv:2106.01817 .
  * Shaqfeh (1996) Shaqfeh, E. S. G. 1996 Purely Elastic Instabilities in Viscometric Flows. Annual Review of Fluid Mechanics 28, 129–185.
  * Shekar et al. (2020) Shekar, Ashwin, MucMullen, Ryan M., McKeon, B. J. & Graham, M. D. 2020 Self-sustained elastoinertial tommien-schlichting waves. Journal of Fluid Mechanics 897, A3.
  * Shekar et al. (2018) Shekar, Ashwin, MucMullen, Ryan M., Wang, S. N., McKeon, B. J. & Graham, M. D. 2018 Critical-layer structures and mechanisms in elastoinertial turbulence. Physical Review Letters 122, 124503.
  * Sid et al. (2018) Sid, S., Terrapon, V. E. & Dubief, Y. 2018 Two-dimensional dynamics of elasto-inertial turbulence and its role in polymer drag reduction. Physical Review Fluids 3, 01130(R).
  * Steinberg (2021) Steinberg, Victor 2021 Elastic turbulence: an experimental view on inertialess random flow. Annual Review of Fluid Mechanics 53, 27\.
  * Stuart (1960) Stuart, J. T. 1960 On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. part 1. the basic behaviour in plane poiseuille flow. Journal of Fluid Mechanics 9, 353–370.
  * Tabor & de Gennes (1986) Tabor, M. & de Gennes, P. G. 1986 A cascade theory of drag reduction. Europhysics Letters 2 (7), 519–522.
  * Terrapon et al. (2015) Terrapon, V., Dubief, Yves & Soria, Julio 2015 On the role of pressure in elasto-inertial turbulence. Journal of Turbulence 16, 26–43.
  * Toms (1948) Toms, B. A. 1948 Observation on the flow of linear polymer solutions through straight tubes at large reynolds numbers. First Intern. Congr. on Rheology 11, 135–141.
  * Vashney & Steinberg (2018) Vashney, A. & Steinberg, V. 2018 Drag enhancement and drag reduction in viscoelastic flow. Physical Review Fluids 3, 103302.
  * Virk (1970) Virk, P. S. 1970 Drag reduction fundamentals. AIChE Journal 21, 625–656.
  * Watson (1960) Watson, J. 1960 On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. part 2. the development of a solution for plane poiseuille flow and for plane couette flow. Journal of Fluid Mechanics 9, 372–389.
  * White & Mungal (2008) White, C. M. & Mungal, M. G. 2008 Mechanics and prediction of turbulent drag reduction with polymer additives. Annual Review of Fluid Mechanics 40, 235–256.
  * Xi & Graham (2010) Xi, L. & Graham, M. D. 2010 Turbulent drag reduction and multistage transitions in viscoelastic minimal flow units. Journal of Fluid Mechanics 647, 421–452.
  * Xi & Graham (2012) Xi, L. & Graham, M. D. 2012 Intermittent dynamics of turbulence hibernation in Newtonian and viscoelastic minimal channel flows. Journal of Fluid Mechanics 693, 433–472.
  * Zhang et al. (2013) Zhang, M., Lashgari, I., Zaki, T. A. & Brandt, L. 2013 Linear stability analysis of channel flow of viscoelastic Oldroyd-B and FENE-P fluids. Journal of Fluid Mechanics 737, 249–279.