Finite-amplitude elastic waves in viscoelastic channel flow from large to zero Reynolds number

Abstract

Using branch continuation in the FENE-P model, we show that finite-amplitude travelling waves borne out of the recently discovered linear instability of viscoelastic channel flow (Khalid et al., J. Fluid Mech., vol. 915, 2021, A43) are substantially subcritical reaching much lower Weissenberg ($Wi$) numbers than on the neutral curve at a given Reynolds ($Re$) number over $Re \in [0,3000]$. The travelling waves on the lower branch are surprisingly weak indicating that viscoelastic channel flow is susceptible to (nonlinear) instability triggered by small finite-amplitude disturbances for $Wi$ and $Re$ well below the neutral curve. The critical $Wi$ for these waves to appear in a saddle node bifurcation decreases monotonically from, for example, $\approx 37$ at $Re=3000$ down to $\approx 7.5$ at $Re=0$ at the solvent-to-total-viscosity ratio $\beta =0.9$. In this latter creeping flow limit, we also show that these waves exist at $Wi \lesssim 50$ for higher polymer concentrations, $\beta \in [0.5,0.97)$, where there is no known linear instability. Our results therefore indicate that these travelling waves, found in simulations and named ‘arrowheads’ by Dubief et al. (Phys. Rev. Fluids, vol. 7, 2022, 073301), exist much more generally in $(Wi,Re, \beta )$ parameter space than their spawning neutral curve and, hence, can either directly, or indirectly through their instabilities, influence the dynamics seen far away from where the flow is linearly unstable. Possible connections to elastic and elasto-inertial turbulence are discussed.

1. Introduction

It is now well known that even small concentrations of long-chain polymers in a Newtonian solvent can give rise to interesting new behaviour (e.g. Larson 1988). Perhaps the most extreme demonstration of this is the existence of ‘elastic’ turbulence (ET) at vanishingly small Reynolds numbers ($Re$) where inertia is minimal (Groisman & Steinberg 2000, 2001; Steinberg 2021). In 2013, a further multiscale, time-dependent state, ‘elasto-inertial’ turbulence (EIT), was found which differs from Newtonian turbulence (NT) in being predominantly two-dimensional and seems to require finite Reynolds number ($Re=O(10^{3}$)) and Weissenberg number $Wi=O(10)$ to exist (Dubief, Terrapon & Soria 2013; Samanta et al. 2013; Sid, Terrapon & Dubief 2018). Understanding exactly how these different types of turbulence relate to each other remains an outstanding challenge. Work at the NT–EIT interface has so far focussed on the possible sustenance of elastically modified Tollmien–Schlichting waves at least for very dilute solutions and weak elasticity (Shekar et al. 2019, 2020). Our focus here is the possible relationship between EIT and ET: are they two extremes of one whole (Samanta et al. 2013; Qin et al. 2019; Choueiri et al. 2021; Steinberg 2021) or distinct flow responses (see, e.g., figure 30 of Chaudhary et al. (2021) and figures 21 and 22 of Datta et al. (2021)). Finding the dynamical origin for either could help in resolving this question.

The very recent discovery of a new linear instability in dilute viscoelastic rectilinear flows at high $Wi=O(20)$ (in pipes by Garg et al. (2018) and channels by Khalid et al. (2021a)) seems highly relevant. Such ‘straight’ flows had always been believed linearly stable due to the absence of curved streamlines (see, e.g., Chaudhary et al. (2019, 2021), Datta et al. (2021) and Castillo-Sanchez et al. (2022) for extensive discussions of this) although there had been some evidence of instability to finite-amplitude disturbances at low $Re$ (Bertola et al. 2003; Pan et al. 2013; Choueiri et al. 2021; Jha & Steinberg 2021). Significantly, the neutral curve for this instability lies in a region of the $(Wi,Re)$ parameter space between where EIT and ET are believed to exist. The instability was initially only found above $Re \approx 63$ in pipe flow in the Oldroyd-B model (Garg et al. 2018; Chaudhary et al. 2021), suggesting that it needs some inertia to function. However, the corresponding instability in channel flow was found to have no such finite-$Re$ threshold, although for Oldroyd-B fluids, the instability is restricted to ultra dilute solutions with $\beta \gtrsim 0.99$, and very large $Wi=O(10^{3})$ (Khalid et al. 2021a; Khalid, Shankar & Subramanian 2021b). Subsequently, these conditions have been relaxed to a more physically relevant critical $Wi\gtrsim 110$ at $\beta \approx 0.98$, by limiting the maximum extension of the polymers ($L_{max}=70$) in a FENE-P model (Buza, Page & Kerswell 2022). This suggests that a purely elastic instability can smoothly morph into an elasto-inertial instability, where inertia plays a role but the instability is found to only derive its energy through elastic terms. This remains the case even as high as $Re=O(1000)$ (Buza et al. 2022). While this new instability is active for a wide range of parameter values, it does not appear to overlap with areas where either EIT or ET have been found, consistently appearing at much higher $Wi$ at a given $Re$. Therefore, the question of its relevance to these nonlinear states remains open.

A key issue is whether the branch of travelling wave (TW) solutions which emerge from the neutral curve is subcritical and so exist down to some saddle node at Weissenberg number $Wi_{sn}$ below the critical value $Wi_c$, thereby potentially connecting the instability to either EIT and/or ET in parameter space. Page, Dubief & Kerswell (2020) demonstrated the existence of substantial subcriticality albeit at $Re=60$ (and $\beta =0.9$) where $Wi_{sn}=8.8$ is much lower than $Wi_c=26.7$. Despite EIT not existing at this low $Re$, the upper branch TWs found there clearly resembled the ‘arrowhead’ states found in the simulations of Dubief et al. (2022) at $Re=1000$ when EIT was annealed by increasing the elasticity. Weakly nonlinear analysis (in the channel by Buza et al. (2022) and pipe flow by Wan, Sun & Zhang (2021)) has confirmed the general subcritical nature of the instability but cannot give global information about how far $Wi_{sn}(Re,\beta )$ is below $Wi_c(Re,\beta )$. Our purpose here is to answer this by performing an investigation using branch continuation to track where the TWs exist in $(Wi,Re,\beta )$-parameter space. This turns out to be feasible, as a branch continuation procedure based on solving an algebraic set of equations derived directly from the governing equations is much more efficient than branch continuing using a direct numerical simulation (DNS) code as done in Page et al. (2020). There are two reasons for this. First, the TW is highly symmetric: it is two-dimensional and has a symmetry around the channel's midplane. Second, far fewer degrees of freedom are needed to resolve the flow algebraically compared with the number needed to keep a time-stepping code stable. For example, the algebraic formulation needs only ${\approx }50$ Chebyshev modes in the cross-stream direction for convergence at the parameters considered while the DNS code needs ${\approx }128$ modes to remain time stable. There have been previous theoretical attempts to generate nonlinear solutions to viscoelastic flow in channels and pipes but without an anchoring bifurcation point. These have centred on constructing a high-order expansion assuming the solution is dominantly streamwise and temporally monochromatic and taking the leading state to be one of the least-damped linear modes of the base state (Meulenbroek et al. 2003; Morozov & Saarloos 2005; Morozov & van Saarloos 2019).

This approach has produced some interesting signs of convergence with an increasing number of terms included in the expansion. In particular, by taking expansions up to 11th order in the amplitude, Morozov & Saarloos (2005) and Morozov & van Saarloos (2019) (plane Couette and channel flow, respectively) see apparent convergence to non-trivial TW solutions in creeping ($Re \ll 1$) flows of upper-convected Maxwell (UCM) fluids ($\beta = 0$) as well as Oldroyd-B fluids at low $\beta$. The branch continuation used here is similar in spirit but closer to classical weakly nonlinear theory, and differs in two significant ways: (1) it is firmly rooted in the neutral curve found by Khalid et al. (2021a), i.e. the zero-amplitude limit smoothly leads to the neutral curve (unknown in Morozov & van Saarloos 2019); and (2) the order of the expansion is taken as high as necessary (typically 50–80 Fourier modes) to obtain convergence.

The rest of the paper is organised as follows. Section 2 briefly recaps the formulation of viscoelastic channel flow described in our earlier work (Buza et al. 2022). Section 3 then outlines the branch continuation approach, with the technical details relegated to a series of Appendices. The results are presented in two sections: § 4 considers finite inertia $Re >0$ and § 5 deals with inertialess flows at $Re=0$. Section 4 exclusively concentrates on $\beta =0.9$ and considers how the subcritical TW branches behave as: (1) $Re$ varies over the range $Re \in [0,3000]$; and (2) as the domain size varies at $Re=30$. For the analysis of creeping flow in § 5 at $Re=0$, we explore the existence of the TWs over the $(Wi,\beta )$ plane for $\beta \in [0.5,1)$ and $Wi < 50$ and $L_{max} \in \{70,100,500\}$, the maximum polymer extensibility in the FENE-P model. Morozov (2022) has concurrently found TWs in viscoelastic channel flow at $Re=0.01$ by the complementary approach of time stepping in the Phan–Thien–Tanner model. These waves correspond to the attracting upper branch of the curves shown here. Finally, a discussion follows in § 6.

2. Formulation

As in Buza et al. (2022), we consider pressure-driven flow of an incompressible viscoelastic fluid in a channel bounded by two parallel, stationary, rigid plates separated by a distance of $2h$. We model viscoelasticity using the FENE-P model so that the governing equations are

(2.1a)$$\begin{gather} Re [ \partial_t \boldsymbol u + ( \boldsymbol u \boldsymbol{\cdot} \boldsymbol{\nabla} ) \boldsymbol u ]+ \boldsymbol{\nabla} p = \beta {\rm \Delta} \boldsymbol u + (1-\beta) \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}(\boldsymbol{\mathsf{C}}) + \begin{pmatrix} F_x \\ 0 \end{pmatrix}, \end{gather}$$

(2.1b)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol u =0, \end{gather}$$

(2.1c)$$\begin{gather}\partial_t \boldsymbol{\mathsf{C}} + ( \boldsymbol u \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{\mathsf{C}} + \boldsymbol{\mathsf{T}}(\boldsymbol{\mathsf{C}}) = \boldsymbol{\mathsf{C}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol u + ( \boldsymbol{\nabla} \boldsymbol u )^{\rm T} \boldsymbol{\cdot} \boldsymbol{\mathsf{C}} + \frac{1}{Re Sc} {\rm \Delta} \boldsymbol{\mathsf{C}}. \end{gather}$$

The constitutive relation for the polymer stress, $\boldsymbol{\mathsf{T}}$, is given by the Peterlin function

(2.1d)\begin{equation} \boldsymbol{\mathsf{T}}(\boldsymbol{\mathsf{C}}) := \frac{1}{Wi} ( f (\mathrm{tr} \, \boldsymbol{\mathsf{C}}) \boldsymbol{\mathsf{C}} - \boldsymbol{\mathsf{I}}), \quad {\rm where} \ f(s) := \left(1-\frac{s-3}{L^{2}_{max}}\right)^{{-}1}, \end{equation}

with $L_{max}$ denoting the maximum extensibility of polymer chains. Here $\boldsymbol{\mathsf{C}} \in {Pos}(3)$ is the positive-definite polymer conformation tensor and $\beta := \nu _s/\nu \in [0,1]$ denotes the viscosity ratio 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-dimensionalized by $h$ and the bulk speed

(2.2)\begin{equation} U_b:= \frac{1}{2h} \int^{h}_{{-}h} u_x\, {{\rm d} y}, \end{equation}

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

(2.3a,b)\begin{equation} Re:= \frac{hU_b}{\nu}, \quad Wi:= \frac{\tau U_b}{h}, \end{equation}

where $\tau$ is the polymer relaxation time.

The Schmidt number $Sc$, appearing solely in the polymer diffusion term and defined as the ratio between the solvent kinematic viscosity and polymer diffusivity (Sid et al. 2018), is typically of the order of $O (10^{6})$ in physical applications. In this work, enhanced diffusion (i.e. lower $Sc$) had to be employed to regularize the hyperbolic equation (2.1c), as is customarily done in other works involving viscoelastic DNS (Sid et al. 2018; Dubief et al. 2022). This is also necessitated by the spectral methods embedded in our branch continuation scheme, which behave slightly worse than finite-difference methods in this respect, pushing the maximum admissible Schmidt number down to $Sc = 250$ from typically $1000$ (cf. Dubief et al. (2022) and Appendix C).

Equation (2.1) is supplemented with non-slip boundary conditions on the velocity field. For the conformation tensor $\boldsymbol{\mathsf{C}}$, we impose

(2.4)\begin{equation} \partial_t \boldsymbol{\mathsf{C}} + ( \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{\mathsf{C}} + \boldsymbol{\mathsf{T}}(\boldsymbol{\mathsf{C}}) = \boldsymbol{\mathsf{C}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol u + ( \boldsymbol{\nabla} \boldsymbol u )^{\rm T} \boldsymbol{\cdot} \boldsymbol{\mathsf{C}} + \frac{1}{Re Sc} \partial_{xx} \boldsymbol{\mathsf{C}} , \end{equation}

at the wall, i.e. we minimize the deviation from the $Sc \to \infty$ limit, where no boundary conditions are necessary. In the streamwise ($x$) direction, periodic boundary conditions are imposed on both $\boldsymbol u$ and $\boldsymbol{\mathsf{C}}$. Solutions to (2.1) of the form

(2.5)\begin{equation} \boldsymbol{\varphi}(x,y,t;Wi,Re,\beta) = \boldsymbol{\varphi}_b(y;Wi,Re,\beta) + \hat{\boldsymbol{\varphi}}(X:=x-ct,y;Wi,Re,\beta), \end{equation}

(where $\boldsymbol {\varphi } = (u_X,u_y,p,C_{XX},C_{yy},C_{zz},C_{Xy})$ is the vector composed of all variables) are sought in two consecutive steps. First, the steady, one-dimensional base state $\boldsymbol {\varphi }_b(y;Wi,Re,\beta )$ is solved for numerically at a given $Wi$, $Re$ and $\beta$ (with other model parameters such as $L_{max}$ suppressed for clarity). Then a possibly large, two-dimensional perturbation $\hat {\boldsymbol {\varphi }}(X,y;Wi,Re,\beta )$ is sought which is steady in a frame travelling at some a priori unknown phase speed $c$ in the $\hat {\boldsymbol x}$ direction.

3. Numerical methods

For TWs, time derivatives can be replaced by $-c\partial _X$ in the governing equations and the problem then becomes elliptic with a ‘nonlinear’ eigenvalue $c$. This approach circumvents the need for time integration but at the price of specialising to steady solutions viewed from some Galilean frame. Writing the various terms of the governing equations (2.1) for $\hat {\boldsymbol {\varphi }}$ according to their degree of nonlinearity gives

(3.1)\begin{equation} \mathcal{L}[\hat{\boldsymbol{\varphi}}] + \mathcal{B}[\hat{\boldsymbol{\varphi}}, \hat{\boldsymbol{\varphi}}] + \mathcal{N}[\hat{\boldsymbol{\varphi}}] + \boldsymbol F = \boldsymbol{0},\end{equation}

where

(3.2)\begin{equation} \mathcal{L}[\hat{\boldsymbol{\varphi}}] := \begin{pmatrix} Re ({-}c \partial_X \hat{\boldsymbol u} + ( \boldsymbol u_b \boldsymbol{\cdot} \boldsymbol{\nabla} ) \hat{\boldsymbol u} + ( \hat{\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol u_b ) + \boldsymbol{\nabla} \hat{p} - \beta {\rm \Delta} \hat{\boldsymbol u}\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \hat{\boldsymbol u}\\ -c \partial_X \hat{\boldsymbol{\mathsf{C}}} + ( \boldsymbol u_b \boldsymbol{\cdot} \boldsymbol{\nabla} ) \hat{\boldsymbol{\mathsf{C}}} + ( \hat{\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{\mathsf{C}}_b - 2 \mathrm{sym} ( \boldsymbol{\mathsf{C}}_b \boldsymbol{\cdot} \boldsymbol{\nabla} \hat{\boldsymbol u} + \hat{\boldsymbol{\mathsf{C}}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol u_b) - \dfrac{1}{Re Sc} {\rm \Delta} \hat{\boldsymbol{\mathsf{C}}}\\ \end{pmatrix}, \end{equation}

collects the linear contributions,

(3.3)\begin{equation} \mathcal{B}[\hat{\boldsymbol{\varphi}}_1,\hat{\boldsymbol{\varphi}}_2] := \begin{pmatrix} Re \, ( \hat{\boldsymbol u}_1 \boldsymbol{\cdot} \boldsymbol{\nabla} ) \hat{\boldsymbol u}_2\\ 0\\ (\hat{\boldsymbol u}_1 \boldsymbol{\cdot} \boldsymbol{\nabla}) \hat{\boldsymbol{\mathsf{C}}}_2 - 2 \mathrm{sym} ( \hat{\boldsymbol{\mathsf{C}}}_1 \boldsymbol{\cdot} \boldsymbol{\nabla} \hat{\boldsymbol u}_2 ) \end{pmatrix}, \end{equation}

forms the bilinear part of the nonlinearity and

(3.4a,b)\begin{align} \mathcal{N} [ \boldsymbol{\varphi} ] := \begin{pmatrix} -(1- \beta) \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}} (\hat{\boldsymbol{\mathsf{C}}})\\ 0\\ \boldsymbol{\mathsf{T}} (\hat{\boldsymbol{\mathsf{C}}}) \end{pmatrix} \quad \text{and} \quad \boldsymbol F := \begin{pmatrix} (1- \beta) \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}(\boldsymbol{\mathsf{C}}_b) +\begin{pmatrix} F_X \\ 0 \end{pmatrix}\\ 0\\ - \boldsymbol{\mathsf{T}}(\boldsymbol{\mathsf{C}}_b) \end{pmatrix}, \end{align}

contain the remainder of the terms, with $\mathcal {N}$ representing the general nonlinearity that originates from the constitutive relation ($\boldsymbol u_b$ is the base flow and $\boldsymbol{\mathsf{C}}_b$ is the base conformation tensor which, along with the base pressure, make up $\boldsymbol {\varphi }_b$). The channel is the two-dimensional domain $\varOmega = S^{1} \times [-1,1]$, with $S^{1} = \mathbb {R} / (2 {\rm \pi}/ k)\mathbb {Z}$ denoting a $2 {\rm \pi}/ k$-periodic domain that represents the streamwise ($X$) direction.

The bifurcating eigenfunction has a symmetry about the midplane, $(u_X,C_{XX},C_{yy},C_{zz})$ are symmetric in $y$ and $(u_y,C_{Xy})$ are antisymmetric, which is preserved at finite amplitude in the subsequent ‘arrowhead’-type TWs. This is exploited in what follows by only solving the flow in the lower half of the channel $y \in [-1,0]$ and assuming appropriate symmetry conditions at the midplane $y=0$.

3.1. Branch continuation

All dependent variables are approximated using a Fourier–Chebyshev basis $\{ \phi _n(X) \psi _m(y) \}_{n,m \in \mathbb {N}}$, where

(3.5a,b)\begin{equation} \phi_n(X):= \sqrt{k/(2 {\rm \pi})} \, {\rm e}^{{\rm i}nkX} \quad {\rm and} \quad \psi_m(y) := \cos [m \cos^{{-}1}(2y+1)]. \end{equation}

Corresponding to this basis, a TW truncated at order $N_X \times N_y$ may be written as

(3.6)\begin{equation} \hat{\boldsymbol{\varphi}}(X,y) = \sum_{n ={-} N_X}^{N_X} \sum_{m=0}^{N_y} \boldsymbol a_{nm} \phi_n(X) \psi_m(y), \end{equation}

where $\boldsymbol a_{nm} \in \mathbb {C}^{7}$ is the vector of coefficients satisfying

(3.7)\begin{equation} \boldsymbol a_{({-}n)m} = \bar{\boldsymbol a}_{nm},\end{equation}

($\bar {\boldsymbol a}_{nm}$ is the complex conjugate of $\boldsymbol a_{nm}$). Substituting (3.6) into (3.1), gives

(3.8)\begin{align} \sum_{n,m} \mathcal{L}[\boldsymbol a_{nm} \phi_n \psi_m] + \sum_{n,m} \sum_{p,q} \mathcal{B} [ \boldsymbol a_{nm} \phi_n \psi_m, \boldsymbol a_{pq} \phi_p \psi_q ] + \mathcal{N} \left[\sum_{n,m} \boldsymbol a_{nm} \phi_n \psi_m\right] + \boldsymbol F = \boldsymbol{0}. \end{align}

A projection onto the $j$th Fourier mode now yields ($\mathcal {L}_\ell [\varphi ]$ is a slight abuse of notation that stands for $(\mathrm {vec}(\mathcal {L}[\varphi ]))_\ell$)

(3.9)\begin{align} &\sum_{m} \mathcal{L}^{j}_\ell[\boldsymbol a_{jm} \psi_m] + \sum_{m} \sum_{q} \sum_r \mathcal{B}^{j-r}_\ell [ \boldsymbol a_{rm} \psi_m, \boldsymbol a_{(j-r)q} \psi_q ]\nonumber\\ &\quad +\left\langle \mathcal{N}_\ell \left[\sum_{n,m} \boldsymbol a_{nm} \phi_n \psi_m\right], \phi_j \right\rangle_{L^{2}(S^{1};\mathbb{C})} + F_\ell \delta_{0j} = 0. \end{align}

where $\mathcal {L}^{j}$ (and, similarly, $\mathcal {B}^{j}$) is the operator $\mathcal {L}$ (and $\mathcal {B}$) modified such that derivatives in the streamwise direction $\partial _X$ are replaced by multiplications with $ikj$. Thus, the $X:=x-ct$ dependence is now fully eliminated from the equations. To treat the $y$ direction, a collocation method is employed over the Gauss–Lobatto points given by

(3.10)\begin{equation} y_s = \frac{1}{2}\left[\cos \left(\frac{s {\rm \pi}}{N_y} \right)-1\right] \in [{-}1,0], \quad s = 0, \ldots, N_y, \end{equation}

Crucially, these are concentrated near both the channel boundary and the centreline where the resolution is generally most needed. The exception is near the saddle node where the ‘arrowhead’ polymer structure significantly extends into the region between the midplane and boundary of the channel and is therefore the most challenging to resolve (see, e.g., figure 7 later). The resulting system of complex algebraic equations are

(3.11)\begin{align} &\sum_{m} \mathcal{L}^{j}_\ell[\boldsymbol a_{jm} \psi_m](y_s) + \sum_{m} \sum_{q} \sum_r \mathcal{B}^{j-r}_\ell [ \boldsymbol a_{rm} \psi_m, \boldsymbol a_{(j-r)q} \psi_q ](y_s) \nonumber\\ &\quad + \left\langle \mathcal{N}_\ell \left[\sum_{n,m} \boldsymbol a_{nm} \phi_n \psi_m \right](y_s), \phi_j \right\rangle_{L^{2}(S^{1};\mathbb{C})} \nonumber\\ &\quad + F_\ell (y_s) \delta_{0j}= 0, \quad \text{for } j = 0,\ldots,N_X, \quad s = 0,\ldots,N_y, \quad \ell = 1,\ldots,7, \end{align}

for the coefficients $\boldsymbol a_{nm} \in \mathbb {C}^{7}$ with $n,m \geq 0$. The remainder of the coefficients in (3.6) are computed via (3.7).

Two further equations are needed to determine the wave speed $c$ and the applied pressure gradient $F_X$. As indicated previously, $F_X$ is determined by ensuring the perturbation volume flux vanishes,

(3.12)\begin{equation} \int^{1}_{{-}1}\hat{u}_X \, {{\rm d} y}=0,\end{equation}

and

(3.13)\begin{equation} {Im} \int^{2{\rm \pi}/k}_0 \, {\rm e}^{-{\rm i}kX} \hat{u}_X(X,y_{15}) \, {{\rm d}\kern0.06em x} =0,\end{equation}

is imposed to eliminate the phase degeneracy of the TW and thereby determine the wave speed (the exact collocation point $y_{15}$ is chosen arbitrarily; see, e.g., Wedin & Kerswell (2004)). The resulting nonlinear, complex, algebraic system of equations comprising (3.11), (3.12) and (3.13) reads

(3.14)\begin{equation} \mathcal{F} (\boldsymbol a,c,F_X;k,Wi,Re,\beta,Sc) = \boldsymbol{0},\end{equation}

where $\boldsymbol a = \mathrm {vec}((a_{nm})_\ell )$. System (3.14) gives rise to $Q:=2+2 \times 7 \times N_X \times (N_y+1) + 7 \times (N_y+1) \sim 14N_X N_y$ real nonlinear equations to be solved simultaneously (by slight abuse of notation we shall denote the real parts of $\mathcal {F}$ and $\boldsymbol a$ in (3.14) by the same letters in what follows). Steady states of interest may now be extracted from (3.14) using a Newton–Raphson root finding scheme given a good enough initial guess. The neutral curve found by Khalid et al. (2021a) and the weakly nonlinear analysis in Buza et al. (2022) are used to generate this initially. Then pseudo arclength continuation (see Appendix A) is used to proceed along the solution branch to higher amplitudes away from the neutral curve.

Simulations were typically run at $(N_X,N_y) = (50,60)$ where $Q \approx 43\ 000$ real degrees of freedom or $(N_X,N_y) = (40,50)$ ($Q \approx 29\ 000$), depending on the complexity of the tracked states, with occasional grid-convergence checks at much higher resolutions up to $(N_X,N_y) = (80,80)$ ($Q \approx 91\ 000$); see Appendix B. In general, lower branch solutions were less resolution dependent, and required about half the Fourier modes of their upper branch counterparts. The minimum requirement for the number of Chebyshev modes, $N_y$, was around $40$ across all parameter regimes, with a slight increase to $50$ around saddle-node points due to the suboptimal placement of collocation points for this region. Reducing the polymer diffusion increases the requirements both in $N_X$ and $N_y$, and adjustments in $k$, the domain size, necessitate equivalent adjustments in $N_X$.

3.2. Direct numerical simulations

The Dedalus codebase (Burns et al. 2020) was used to time step (2.1) in order to examine the stability of the TWs found. To allow this DNS code to interface seamlessly with the branch continuation code, the simulations were also performed on the half-channel using exactly the same symmetry boundary conditions described previously and using the same spectral expansions. This allowed an unstable lower branch solution of the branch continuation procedure to be used directly as an initial condition for the DNS and the fact that this remained steady under time-stepping provided a valuable cross-check of the two approaches.

In the DNS, the full state $\boldsymbol {\varphi }$ was minimally expanded into $N_x=128$ Fourier modes in the periodic $x$-direction and into $N_y=128$ Chebyshev modes in the wall-normal direction with higher resolutions of 256 and 512 available in either or both dimensions to check truncation robustness. The equations were advanced in time using a third-order semi-implicit backward differentiation formula (BDF) scheme (Wang & Ruuth 2008) and a constant timestep ${\rm \Delta} t=5\times 10^{-3}$.

4. Results: TWs at finite $Re$ for $(\beta, {L}_{max},Sc) = (0.9,500,250)$

As in Page et al. (2020) and Buza et al. (2022), we fix $\beta = 0.9$ and $L_{max} = 500$ for the initial set of computations. The Schmidt number had to be chosen slightly smaller than that of Page et al. (2020) and Dubief et al. (2022) at $Sc = 250$ due to the considerations given in § 2 and Appendix C.

Upon supplying the weakly nonlinear predictions as initial conditions to the continuation routine, any branch of solutions emanating from the neutral curve can be tracked starting directly from its bifurcation point. Three branches were launched downwards in $Wi$ at fixed $Re = 1000,200$ and $60$, starting from their respective bifurcation points at $k_{opt} = 4.7, 2.7$ and $1.8$. These wave numbers are optimal in the sense of marginal stability and so do not necessarily minimise $Wi_{sn}(Re,\beta )$, but do provide a good upper estimate of it. An additional, fourth branch was initialized from the lowest point on the neutral curve at $Wi = 30$ and $k_{opt} = 1.6$, continued down to $Re = 30$ at fixed $Wi$, then, after a switch in direction, towards decreasing $Wi$ at fixed $Re$. A schematic depiction of these branches is given in the left panel of figure 1.

Figure 1. (a) Linearly and nonlinearly unstable regions in the $Wi\unicode{x2013}Re$ plane for $\beta =0.9$, $L_{max}=500$ and $Sc=250$. The saddle-node Weissenberg numbers $Wi_{sn}(Re,k)$ shown are $Wi_{sn}(30,1.6)=7.6$, $Wi_{sn}(60,1.8)=8.7$, $Wi_{sn}(200,2.7)=14.3$, $Wi_{sn}(1000,4.7)=24.9$ and $Wi_{sn}(3000,4.7)=36.7$. Coloured horizontal lines correspond to branches on the right panel and symbols indicate the saddle-node points. (b) Solution branches tracking TWs as $Wi$ varies at constant $Re \in \{30,60,200,1000\}$ (note horizontal axis is $Re=20$).

The right panel of figure 1 shows the amplitude evolution of these four branches. As a measure of amplitude, we chose the volume-averaged trace of the polymer conformation relative to the laminar value, i.e.

(4.1)\begin{equation} \mathcal{A}: = \frac{\langle \mathrm{tr} \, \boldsymbol{\mathsf{C}} \rangle_\varOmega}{\langle \mathrm{tr} \, \boldsymbol{\mathsf{C}}_b \rangle_\varOmega}, \end{equation}

again, to remain consistent with Page et al. (2020). The lower $Re$ branches of figure 1(b) are reminiscent of the branch shown in Page et al. (2020) (in fact, the green branch is at the same Reynolds number ($Re = 60$), albeit with different $k$ and $Sc$), and the higher $Re$ ones are lower-amplitude variants of these. This shrink in relative amplitude can be attributed to the increase in both $Re$ and $k$, with the latter playing a non-negligible role through the accompanying change in domain size (see § 4.2).

We explore one of these states (the saddle node from the $Re=200$ branch) further in figure 3, where we report the perturbation velocities as a fraction of the local base streamwise velocity, $u_{b,X}$. The arrowhead of polymer stretch close to the centreline is associated with a backwards ‘jet’ in the perturbation streamwise velocity where the horizontal flow is locally slower than the base parallel solution. Both this and the contours of vertical velocity are consistent with the counter-rotating vortex pair observed by Morozov (2022) in the co-moving frame, and include a change in sign of $\partial _y \hat {u}_y$ across a stagnation point leading to a locally extensional flow.

4.1. General interpretation of solution branches

Qualitatively, all solution branches behave in a similar way. A sample case is depicted in figure 4, showcasing the main features. The lower branches emanating from the neutral curve are all unstable until reaching their respective saddle-node points, labeled by a variety of symbols in figure 1 (circle for our sample branch in figure 4), with the corresponding states shown in figure 2. Points on the (unstable) lower branches are found to be edge states which are attracting states on a codimension-one manifold separating two different basins of attraction (Skufca, Yorke & Eckhardt 2006; Duguet, Willis & Kerswell 2008; Schneider et al. 2008). This is illustrated by figure 4(b) at $Wi=20$ and $Re=60$ where an edge-tracking procedure, applied between the upper branch and laminar states, converges on the lower branch state. The lower branch state is a saddle but with only one unstable direction either pointing to the laminar or upper branch state. Upper branch states, at least $Re = 60$ (Page et al. 2020), start as stable nodes as $Wi$ increases away from $Wi_{sn}$ but quickly experience Hopf bifurcations to tertiary states (if the base state is the ‘primary’). These bifurcations and where these tertiary states lead are interesting questions beyond the scope of this paper.

Figure 2. Snapshots of full states, in terms of $\mathrm {tr} \, \boldsymbol{\mathsf{C}}$ (contours), at the saddle-node bifurcation points from figure 1. Lines correspond to level sets of the perturbation stream function. The noticeable thinning of the arrowhead structure with increasing $Re$ (from bottom right to top left) is partially due to the corresponding increase in $k_{opt}$ (and decrease in domain length; see § 4.2 and compare with figure 5). The domain size is $2 {\rm \pi}/ k_{opt}$ in the $X$ direction, with $k_{opt} \in \{4.7,2.7,1.8,1.6\}$ in panel order.

Figure 3. (a) Contours of $\hat {u}_X/u_{b,X}$ at the saddle-node point of the $Re = 200$ branch (blue triangle in figures 1 and 2). The domain length is $2 {\rm \pi}/ 2.7$ as in figure 2. (b) Same as (a) but with $\hat {u}_y/u_{b,X}$ displayed. (c) Mean velocity profiles $u_X/u_{b,X}$ at the saddle-node points marked on figures 1 and 2, with colour codes matching that of figure 1.

Figure 4. (a) The $Re = 60$ branch from figure 1, with edge states indicated by the dashed line. (b) Results of the independent edge-tracking algorithm at $Wi=20$. The blue/red lines indicate evolutions which start close to the lower branch state and converge to the laminar/upper branch state, respectively.

Based on these observations, we have the following picture: if the laminar state is disturbed with a perturbation large enough to reach a certain threshold, determined by the minimal amplitude of approach of the stable manifold of a lower branch state, the flow will evolve towards the upper branch, forming a stable TW. The threshold amplitude to trigger growth is bounded above by the amplitude of the lower branch state itself, which remains $\mathcal {A} < 1.05$ across the domain of existence of TWs. In other words, this domain (shaded bright grey on figure 1) is nonlinearly unstable when subjected to finite but small amplitude disturbances.

4.2. Influence of domain length

This section is dedicated to studying the effect of $k$, and thus the influence of domain size $[0,2 {\rm \pi}/k]$, on the TWs. Figure 5 shows how a single branch of TWs at $Re = 30$ (purple in figure 1) changes with $k$. It has already been established that the steady arrowhead structure is highly sensitive to domain length in the EIT regime (Dubief et al. 2022). There, through capturing larger-scale motions, an increase in domain length was found to unveil structures of increasing complexity, with the possibility of inducing chaotic dynamics at certain parameter combinations. Similar tendencies can be observed in our case (cf. figure 5): An increase in $k$ (and, thus, decrease in domain length) has a considerable weakening effect on the arrowhead structure, eventually resulting in a complete eradication of TWs and a subsequent relaminarization. Despite this observation, the location of the saddle-node points seems largely unaffected by $k$ (cf. figure 5), making the marked ‘nonlinearly unstable’ region on figure 1 robust to changes in the assumed periodicity and domain size.

Figure 5. (a) Branches of TWs at different wave numbers, and thus domain sizes, at $k = 1$ (full line), $k$ $(= k_{opt})= 1.6$ (dashed line) and $k = 2$ (dotted line) with fixed $Re=30$ ($k=1.6$ is also shown coloured purple in figure 1). (b) Snapshots of $k = 1$ and $k = 2$ upper branch TWs at fixed $Wi= 19.814$ (note the difference in domain size). Contours correspond to $\mathrm {tr} \, \boldsymbol{\mathsf{C}}$ and lines correspond to level sets of the perturbation stream function. All visible differences are contained in the superimposed wave solutions: the laminar state does not depend on $k$.

4.3. High-elasticity regime: $Re \rightarrow 0$

The high-elasticity regime is difficult to access using time-stepping as it becomes increasing stiff as $Re \rightarrow 0$. The algebraic approach taken here suffers no such problems and we can approach and even consider $Re=0$ (see the next section) without difficulty.

The existence of the centre-mode linear instability at $Re = 0$ is already known in the limit of very dilute polymer solutions ($\beta \to 1$) for $Wi = O (10^{3})$ in Oldroyd-B fluids (Khalid et al. 2021b) and for $Wi = O (10^{2})$ in FENE-P fluids at finite extensibility ($L_{max}$) (Buza et al. 2022). To substantiate its connection to ET, the time evolution of these growing modes has to be tracked to see whether they are able to produce turbulent behaviour, presumably after transitioning through a cascade of intermediate states. Our goal here is to see where the first level of intermediate state, the TWs, exist at low and vanishing $Re$.

Weakly nonlinear theory predicts supercritical behaviour in the high elasticity ($Wi/Re$) regime, i.e. along the lower boundary of the linearly unstable region. To probe this, a fifth branch was initiated at fixed $Wi = 60$, starting upwards in $Re$ as indicated by the weakly nonlinear analysis (Buza et al. 2022) from a marginally stable point in this region (indicated by orange on figure 1). The resulting branch of solutions is shown in figure 6. Given the supercriticality, this branch starts off as a stable node, moving up in $Re$. Almost immediately after leaving the initial bifurcation point (of linear stability), it reaches a saddle-node bifurcation point, turns around and proceeds to advance towards decreasing $Re$, maintaining a relatively low amplitude until reaching a second saddle node and transitioning to the upper branch. This is an example of how local information provided by weakly nonlinear analysis can be misleading. In fact, the neutral curve gives rise to TWs which reach to lower $Wi$ at fixed $Re$ and lower $Re$ at fixed $Wi$ as shown by figure 1. Upon further inspection, it turns out that the lower (secondary) fold shown in figure 6 at $Re \approx 30$ is purely a feature of the polymer diffusion term $1/(Re\, Sc) {\rm \Delta} \boldsymbol{\mathsf{C}}$ growing artificially large (as $Re$ is decreased), the effect of which is already known to destroy small-scale dynamics (Dubief et al. 2022). It turns out that the point at which the saddle-node bifurcation occurs can be delayed arbitrarily by adjusting $Sc$ in accordance with the variations in $Re$ to keep the polymer diffusion finite. Numerical experimentation suggested a revised polymer diffusion term of the form

(4.2)\begin{equation} \frac{\lambda}{Wi} {\rm \Delta} \boldsymbol{\mathsf{C}}, \end{equation}

for some fixed number $\lambda$. The choice of an inverse scaling with $Wi$ is motivated by observations at $Re = 0$ shown in Appendix C. If $\lambda = 0.005$ is enforced for the branch in question, which amounts to fixing the coefficient $1/(Re\, Sc)$ at the point marked by ‘$+$’ in figure 6, the branch of solutions can be followed down to $Re = 0$ along the lower branch (cf. the dashed line in figure 6).

Figure 6. Solution branch at fixed $Wi = 60$, indicated by orange on figure 1 ($k = 1$). The solid black line shows the weakly nonlinear prediction of supercriticality. In orange are the results from branch continuation with the $1/(Re\, Sc)$ formulation (full line) and with the $\lambda$ formulation (dashed line) ($\lambda = 0.005$ and $Sc = 250$) which shows that the branch of TWs quickly turns around and heads to lower $Re$, i.e. the TWs are substantially subcritical.

5. Results: TWs in the creeping flow limit $Re = 0$

Once $Re = 0$ is reached, we redirect the continuation tool towards decreasing $Wi$. The resulting branch takes the familiar shape (from the $Re>0$ cases) shown in figure 7, attaining its saddle-node bifurcation point at $Wi \approx 7.5$, which serves as a lower bound for the region where TWs exist (note the waves found by Morozov (2022), at $Re=0$ are all above $Wi=20$, albeit with a different model). Figure 7 gives a detailed description of this branch, containing snapshots of full states that illustrate how these waves evolve as $Wi$ is varied. Arrowhead-shaped structures are still clearly visible at low $Wi$ (cf. panels on the left-hand side of figure 7), establishing their prevalence even in the high-elasticity regime. We also report the phase speed, $c$, of the inertialess TWs as a function of $Wi$ in figure 7. Similar to the elasto-inertial case reported in Page et al. (2020), the phase speed is always faster than the bulk velocity. On the lower branch, $c$ tends to a constant value which matches that of the linear eigenfunction. The phase speed of the upper branch solution is notably slower, and is presumably set by the self-sustaining mechanism proposed by Morozov (2022).

Figure 7. (Middle right) Branch of TWs at $Re = 0$, $\beta = 0.9$, $L_{max} = 500$, $\lambda = 0.005$, $k=1$ in terms of the amplitude $\mathcal {A}$ and wavespeed $c$ (here the unstable lower branch is distinguished via dashes). All other panels correspond to states at the locations marked via different symbols. In these plots, contours correspond to $\mathrm {tr} \, \boldsymbol{\mathsf{C}}$ and lines correspond to level sets of the perturbation stream function. Domain length is $2 {\rm \pi}$ in all state plots.

For the particular case considered in figure 7 (specifically $k=1$, $L_{max}=500$, $\lambda =0.005$) the stability of steady states was examined along the upper branch using DNS. At four points, $Wi = 10,20,30$ and $50$, solutions of the branch continuation tool were transferred into the Dedalus-based DNS code, and were subsequently subjected to disturbances of finite amplitude. The perturbations were constructed from snapshots extracted from separate simulations of EIT at high-$Re$, which were pre-multiplied by $10^{-6}$ and added to the TWs. All perturbed states returned to their respective stable upper-branch solutions after a period of transient growth, suggesting that two-dimensional ET cannot be initiated from these TWs in a direct manner.

5.1. $\beta \in [0.5,1)$: relation to recent experiments

The first experiments claiming to see nonlinear instability in viscoelastic channel flow were performed by Arratia and colleagues (Pan et al. 2013; Qin & Arratia 2017; Qin et al. 2019). Finite-amplitude perturbations were induced by an array of obstacles placed upstream, with the number of obstacles serving as a measure of amplitude. Based on measurements taken further downstream, far away from the initial disturbances, they conclude the existence of a subcritical nonlinear instability that persists down to $Wi \approx 5.4$. With the caveats that their channel had a square cross-section and FENE-P is an approximation, their results are encouragingly comparable to the two-dimensional channel prediction made here of $Wi_{sn}=7.5$. In Figure S1 of the supplementary material to Pan et al. (2013), the authors indicated the boundary to the observed instability in a $Wi$ versus perturbation amplitude plane, essentially matching our predictions for the threshold of nonlinear instability given by the $Re = 0$ lower branch (shown in the middle-right panel of figure 7) (that the unstable region in our case is bounded by the branch is an immediate consequence of the discussion in § 4.1: once the threshold amplitude of a lower branch edge state is reached, solutions continue to grow). However, in later proceedings, the authors claim that the unstable flow remains time dependent with features reminiscent of ET (Qin & Arratia 2017; Qin et al. 2019), as opposed to the upper branch TW scenario described here.

More recently, (Schnapp & Steinberg 2022) have experimentally examined the stability of viscoelastic channel flow to small disturbances using a channel of width-to-height ratio of 7. The viscosity ratio was $\beta = 0.74$ and so significantly smaller than the above presented $\beta = 0.9$. They describe seeing a complex transition which they call ‘an elastic non-normal mode instability’ characterised by the direct transition to many modes at once. This is different from a modal instability of a TW solution discussed here but may also be the effect of a comparatively large initial disturbance so the transition is immediately strongly nonlinear. Other recent experiments by Choueiri et al. (2021) in a pipe have also considered smaller $\beta =0.57$ albeit at larger $Re$. Motivated by this, we also performed a few TW branch continuations (at $Wi = 10, 20, 30$ and $50$) with decreasing $\beta$ in an attempt to map out the nonlinearly unstable domain in the $Wi\unicode{x2013}\beta$ plane (all with zero inertia). Results from these computations are shown in figure 8, with the left panel indicating the unstable region and the right panel containing the solution branches found with varying $\beta$. It transpires that at lower $\beta$, the solutions are a little more sensitive to the artificial diffusion (see Appendix C for further details), necessitating multiple simulations at different $\lambda$ values, all of which are also shown in figure 8.

Figure 8. (a) Nonlinearly unstable regions in the $Wi$–$\beta$ plane at $Re = 0$ for $\lambda = 0.005$ (full line), $\lambda = 0.003$ (dashed line) and $\lambda = 0.002$ (dotted line). Horizontal and vertical lines correspond to solution branches computed via the continuation routine. (b) Branches with respect to $\beta$, obeying the same colour code as on the left panel (at $\lambda = 0.005$).

Figure 9 tracks changes in the lower branch TW amplitude, which is an upper bound on the threshold for growth, as $Wi$ is increased. For the parameter combination in question ($\beta = 0.74$, $Re = 0$), both our standard amplitude measure $\mathcal {A}$ and a separate measure for the magnitude of velocity perturbations, $\Vert \hat {\boldsymbol u} \Vert _{L^{2} (\varOmega )}$, are shown. For $Wi > 30$, $\mathcal {A}$ remains below $1.005$, implying that perturbations amounting to just $0.5\,\%$ of the laminar conformation field $\mathrm {tr} \, \boldsymbol{\mathsf{C}}_b$ are sufficient to trigger growth, making this scenario practically indistinguishable from a linearly unstable one. Here $\mathcal {A}$ reaches its minimum roughly around $Wi = 100$, then, despite setting off in an increasing trend, remains negligibly small up to $Wi=500$.

Figure 9. (a) The $\beta = 0.74$ lower branch from figure 8, at $\lambda = 0.001$, with respect to the standard amplitude measure $\mathcal {A}$ (in black) and the perturbation velocity magnitude $\Vert \hat {\boldsymbol u} \Vert _{L^{2} (\varOmega )}$ (in blue). (b) State snapshot at $Wi = 20$ from the branch above. Contours show $\hat {u}_X/u_{b,X}$. (c) State snapshot at $Wi = 407$ from the branch above. Contours show $\hat {u}_X/u_{b,X}$.

In addition to the solution branch, figure 9 also displays two state plots at $Wi = 20$ and $Wi = 407$, now in terms of $\hat {u}_X/u_{b,X}$ to aid comparison with experimental results. The former, at $Wi = 20$, still resembles the structural composition of the linearly unstable centre eigenmode at higher $Re$: chevron-shaped streaks remain visible, but are now disconnected at the centreline. However, by $Wi = 407$ (and similarly for all $Wi>100$ lower branch states), these structures straighten out and form streamwise counter-propagating streaks placed symmetrically around the centreline.

In their pipe set-up, Choueiri et al. (2021) detected fluctuations at $Re \approx 5$ for $\beta = 0.57$, $Wi = 104$ while the flow remained laminar at $Re \approx 3$. Computationally, the $Re=0$ branch may be continued to $\beta = 0.57$ (for $\lambda$ sufficiently small) in our channel flow setting (cf. figure 8 again), albeit with a slightly larger threshold amplitude than that of figure 9. This may explain the need for a finite $Re$ for ‘unperturbed’ instabilities (Choueiri et al. 2021). Interestingly, the flow states shown by Choueiri et al. (2021) have the connected, chevron-shaped streaks reminiscent of the centre-mode eigenfunction. In our setting, ‘connectedness’ of the chevrons is lost shortly after leaving the initial bifurcation (the general shape is still retained for moderate $Wi$, see the middle panel of figure 9), but the scenario could be quite different in pipes.

5.2. $L_{max} \in \{70,100,500\}$

For completeness, the effect of varying the last outstanding parameter, $L_{max}$, on the region of nonlinear instability is depicted in figure 10. In contrast to the observations of Buza et al. (2022), which indicated that lowering $L_{max}$ has a destabilizing role in the elastically dominated regime (cf. their figure 13), here we see that the nonlinearly unstable region shrinks with decreasing $L_{max}$. This suppressing effect is in line with past observations of the effect of finite extensibility; see, e.g., the linear analyses in Ray & Zaki (2014), Page & Zaki (2015) or even figure 16 of Buza et al. (2022). It should be noted that, in our experience, solution branches became difficult to extract at lower $L_{max}$, with convergence issues appearing along upper branches. In fact, we can only reliably obtain upper branches for $L_{max}>150$, but lower branches remain accessible due to their lower resolution requirements (see Appendix B).

Figure 10. The effect of varying $L_{max}$ on the nonlinearly unstable region at fixed $\lambda = 0.003$. Horizontal and vertical lines correspond to solution branches computed via the continuation routine with $L_{max} = 70$ (black solid lines), $L_{max} = 100$ (red dashed lines) and $L_{max} = 500$ (blue dotted lines).

6. Discussion

In this paper, we have used branch continuation to track two-dimensional, finite-amplitude TWs in a viscoelastic channel flow, using the FENE-P model. The TWs are borne out of the centre mode instability (Khalid et al. 2021a) which is known to be subcritical over large areas of the $Re-Wi$ parameter space (Buza et al. 2022). Here, we have shown that the TW solution branches extend to significantly lower $Re$ and $Wi$ than the curve of marginal stability. For instance, we showed that the saddle node at $Re=1000$ drops as low as $Wi\approx 25$, while the associated linear instability at this $Re$ does not occur until $Wi\approx 80$. Most significantly, we demonstrated the persistence of the nonlinear TWs at $Re=0$ for a large range of $Wi$, with the saddle node sitting at $Wi \approx 7.5\unicode{x2013}10$ in dilute ($\beta =0.9$) solutions for a range of $L_{max}$. Across a broad range of the parameter space, including at $Re=0$, the upper branch TWs resemble the arrowhead structures observed in DNS at higher $Re$ (Page et al. 2020; Dubief et al. 2022).

A key feature of the solution branch at $Re=0$ is the low amplitude of the lower branch TW across a very large range of $Wi$. This suggests that only a very weak disturbance would be required to cause the flow to transition to a potentially complicated time-dependent state, a scenario which would potentially be indistinguishable from a linear instability in an experiment. The small amplitude of the lower branch also poses the question whether the amplitude expansion pursued in Morozov & Saarloos (2005) and Morozov & van Saarloos (2019) was attempting to resolve it; see, e.g., figure 8(c) of Morozov & van Saarloos (2019).

While our results demonstrate the existence of finite-amplitude arrowhead TWs over a very large range of the parameter space, a direct connection to either EIT or ET has yet to be established even within the same FENE-P model. In inertia-dominated EIT at Reynolds numbers $Re=O(1000)$, instantaneous arrowhead-shaped flow structures resembling the stable upper branch states have been observed numerically (Dubief et al. 2022). While the region of nonlinear instability found here overlaps with these observations, further work is required to assess if there is a direct connection between the exact coherent states and EIT (e.g. through a sequence of successive bifurcations). The picture is complicated by observations from both numerics and experiments at higher $Re$, where near-wall activity has been observed to be dominant in both channels and pipes (Shekar et al. 2019, 2020; Choueiri et al. 2021). At present, no exact coherent states connected to other linear instabilities (e.g. Tollmien–Schlichting waves in channels, see Shekar et al. 2019) have been discovered, and the possible role of such states in EIT remains an interesting open problem. At low $Re$, our numerical experiments indicate that the upper branch TWs are linearly stable, and we have been unable to trigger ET. Our computations have been restricted to two-dimensional flows, and so this does not rule out the existence of ET in three dimensions, or a direct route from the upper branch TWs to such a state. We hope to report results on this soon.