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

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. {\em J. Fluid Mech.} {\bf 915}, A43, 2021) 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 viscolastic 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. {\em arXiv}.2006.06770 (2020) - 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.

(Received xx; revised xx; accepted xx) 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. 915, A43, 2021) are substantially subcritical reaching much lower Weissenberg ( ) numbers than on the neutral curve at a given Reynolds ( ) number over ∈ [0, 3000]. The travelling waves on the lower branch are surprisingly weak indicating that viscolastic channel flow is susceptible to (nonlinear) instability triggered by small finite amplitude disturbances for and well below the neutral curve. The critical for these waves to appear in a saddle node bifurcation decreases monotonically from, for example, ≈ 37 at = 3000 down to ≈ 7.5 at = 0 at the solvent-to-total-viscosity ratio = 0.9. In this latter creeping flow limit, we also show that these waves exist at 50 for higher polymer concentrations -∈ [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. arXiv.2006.06770 (2020 -exist much more generally in ( , , ) 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.

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 ( ) where inertia is minimal (Groisman & Steinberg 2000, 2001Steinberg 2021). In 2013, a further multiscale, time-dependent state -'elasto-inertial' turbulence or EIT -was found which differs from Newtonian turbulence (NT) in being predominantly 2D and seems to require finite Reynolds number ( = (10 3 )) and Weissenberg number = (10) to exist (Samanta et al. 2013;Dubief et al. 2013;Sid et al. 2018). Understanding exactly how these different types of turbulence relate to each 2 other remains an outstanding challenge. Work at the NT-EIT interface has so far focussed on the possible sustenance of elastically-modified Tollmein-Schlicting waves at least for very dilute solutions and weak elasticity (Shekar et al. 2018(Shekar et al. , 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)  The very recent discovery of a new linear instability in dilute viscoelastic rectilinear flows at high = (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 (e.g. see Chaudhary et al. 2019Chaudhary et al. , 2021Datta et al. 2021;Castillo-Sanchez et al. 2022, for extensive discussion of this) although there had been some evidence of instability to finite-amplitude disturbances at low (Bertola et al. 2003;Pan et al. 2013;Choueiri et al. 2021;Jha & Steinberg 2020). Significantly, the neutral curve for this instability lies in a region of the ( , ) parameter space between where EIT and ET are believed to exist. The instability was initially only found above ≈ 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-threshold, although for Oldroyd-B fluids, the instability is restricted to ultra dilute solutions with 0.99, and very large = (10 3 ) (Khalid et al. 2021a,b). Subsequently, these conditions have been relaxed to a more physicallyrelevant critical 110 at ≈ 0.98, by limiting the maximum extension of the polymers ( = 70) in a FENE-P model (Buza et al. 2021). This suggests that a purely elastic instability can smoothly morph into an elasto-inertial one, 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 = (1000) (Buza et al. 2021). 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 at a given . Therefore, the question of its relevance to these nonlinear states remains open.
A key issue is whether the branch of travelling wave solutions which emerge from the neutral curve is subcritical and so exist down to some saddle node at Weissenberg number below the critical value , thereby potentially connecting the instability to either EIT and/or ET in parameter space. Page et al. (2020) demonstrated the existence of substantial subcriticality albeit at = 60 (and = 0.9) where = 8.8 is much lower than = 26.7. Despite EIT not existing at this low , the upper branch travelling waves found there clearly resembled the 'arrowhead' states found in the simulations of Dubief et al.
(2020) at = 1000 when EIT was annealed by increasing the elasticity. Weakly nonlinear analysis (in the channel by Buza et al. (2021) and pipe flow by Wan et al. (2021)) has confirmed the general subcritical nature of the instability but can not give global information about how far ( , ) is below ( , ). Our purpose here is to answer this by performing an investigation using branch continuation to track where the travelling waves exist in ( , , )-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 DNS code as done in Page et al. (2020). There are two reasons for this. Firstly, the travelling wave is highly symmetric: it is 2-dimensional and has a symmetry around the channel's midplane. Secondly, far fewer degrees of freedom are needed to resolve the flow algebraically compared to the number needed to keep a time-stepping code stable. For example, the algebraic formulation needs only ≈ 50 Chebyshev modes in the cross-stream direction for convergence at the parameters considered while the DNS code needs ≈ 128 modes to stay 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 nontrivial TW solutions in creeping ( 1) flows of upper-convected Maxwell (UCM) fluids ( = 0) as well as Oldroyd-B fluids at low . 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 get 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. 2021). 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: section 4 considers finite inertia > 0 and section 5 deals with inertialess flows at = 0. Section 4 exclusively concentrates on = 0.9 and considers how the subcritical travelling wave branches behave as: 1) varies over the range ∈ [0, 3000]; and 2) as the domain size varies at = 30. For the analysis of creeping flow in section 5 at = 0, we explore the existence of the travelling waves over the ( , ) plane for ∈ [0.5, 1) and < 50 and ∈ {70, 100, 500}, the maximum polymer extensibility in the FENE-P model. Morozov (2022) has concurrently found travelling waves in viscoelastic channel flow at = 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 section 6.

Formulation
As in Buza et al. (2021), 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 2ℎ. We model viscoelasticity using the FENE-P model so that the governing equations are which, through adjusting the imposed pressure gradient appropriately, is kept constant so that the Reynolds and Weissenberg numbers are defined as where is the polymer relaxation time. The Schmidt number , 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 order (10 6 ) in physical applications. In this work, enhanced diffusion (i.e. lower ) had to be employed to regularize the hyperbolic equation (2.1c), as is customarily done in other works involving viscoelastic direct numerical simulations (Dubief et al. 2020;Sid et al. 2018). 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 = 250 from typically 1000 (cf. Dubief et al. 2020, and Appendix C).
Equation (2.1) is supplemented with non-slip boundary conditions on the velocity field. For the conformation tensor C, we impose C + (u · ∇) C + T(C) = C · ∇u + (∇u) · C + 1 C at the wall, i.e. we minimize the deviation from the → ∞ limit, where no boundary conditions are necessary. In the streamwise ( ) direction, periodic boundary conditions are imposed on both u and C. Solutions to (2.1) of the form ) is the vector composed of all variables) are sought in two consecutive steps. First, the steady, 1-dimensional base state ( ; , , ) is solved for numerically at a given , and (with other model parameters such as suppressed for clarity). Then a possibly-large, 2-dimensional perturbationˆ( , ; , , ) is sought which is steady in a frame travelling at some a priori unknown phase speed in thex direction.

Numerical Methods
For travelling waves (TW), time derivatives can be replaced by − in the governing equations and the problem then becomes elliptic with a 'nonlinear' eigenvalue . 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ˆaccording to their degree of nonlinearity gives 1 · ∇û 2 forms the bilinear part of the nonlinearity, and and contain the remainder of the terms, with N representing the general nonlinearity that originates from the constitutive relation (u is the base flow and C is the base conformation tensor which, along with the base pressure, make up ). The channel is the 2-dimensional domain Ω = 1 × [−1, 1], with 1 = R/(2 / )Z denoting a 2 / -periodic domain that represents the streamwise ( ) direction.
The bifurcating eigenfunction has a symmetry about the midplane -( , , , ) are symmetric in and ( , ) are antisymmetric -which is preserved at finite amplitude in the subsequent 'arrowhead'-type travelling waves. This is exploited in what follows by only solving the flow in the lower half of the channel ∈ [−1, 0] and assuming appropriate symmetry conditions at the midplane = 0. (3.2) Corresponding to this basis, a TW truncated at order × may be written aŝ

Branch continuation
where a ∈ C 7 is the vector of coefficients satisfying A projection onto the -th Fourier mode now yields † where L (and similarly, B ) is the operator L (and B) modified such that derivatives in the streamwise direction are replaced by multiplications with . Thus, the := − † L ℓ [ ] is a slight abuse of notation that stands for (vec(L [ ])) ℓ . 6 dependence is now fully eliminated from the equations. To treat the direction, a collocation method is employed over the Gauss-Lobatto points given by

5)
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 (e.g. see Figure 7 later). The resulting system of complex algebraic equations are for the coefficients a ∈ C 7 with , 0. The remainder of the coefficients in (3.3) are computed via (3.4).
Two further equations are needed to determine the wave speed and the applied pressure gradient . As indicated above, is determined by ensuring the perturbation volume flux vanishes, and Im is imposed to eliminate the phase degeneracy of the travelling wave and thereby determine the wave speed (the exact collocation point 15 is chosen arbitrarily, e.g. see Wedin & Kerswell (2004)). The resulting nonlinear, complex, algebraic system of equations comprising (3.6), (3.7) and (3.8) reads F (a, , ; , , , , real nonlinear equations to be solved simultaneously (by slight abuse of notation we shall denote the real parts of F and a in (3.9) by the same letters in what follows). Steady states of interest may now be extracted from (3.9) using a Newton-Raphson root finding scheme given a good enough initial guess. The neutral curve found by Khalid et al. (2021a) and the the weakly nonlinear analysis in Buza et al. (2021) 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 ( , ) = (50, 60) where ≈ 43, 000 real degrees of freedom or ( , ) = (40, 50) ( ≈ 29, 000), depending on the complexity of the tracked states, with occasional grid-convergence checks at much higher resolutions up to ( , ) = (80, 80) ( ≈ 91, 000); see Appendix B. Generally, 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, , 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 and , and adjustments in , the domain size, necessitate equivalent adjustments in .

Direct numerical simulations
The Dedalus codebase (Burns et al. 2020) was used to time step eqs. (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 above 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 was minimally expanded into = 128 Fourier modes in the periodic -direction and into = 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 3rd-order semi-implicit BDF scheme (Wang & Ruuth 2008) and a constant timestep Δ = 5 × 10 −3 . 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 at fixed = 1000, 200, 60, starting from their respective bifurcation points at = 4.7, 2.7, 1.8. These wave numbers are optimal in the sense of marginal stability and so do not necessarily minimise ( , ), but do provide a good upper estimate of it. An additional, fourth branch was initialized from the lowest point on the neutral curve at = 30 and = 1.6, continued down to = 30 at fixed , then -after a switch in direction -towards decreasing at fixed . A schematic depiction of these branches is given in the left panel of Figure 1.

Results
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., We explore one of these states (the saddle node from the = 200 branch) further in figure 3, where we report the perturbation velocities as a fraction of the local base streamwise velocity, , . The arrowhead of polymer stretch close to the centreline is associated with a backwards 'jet' in the perturbation streamwise velocity. The contours of vertical velocity indicate a change in sign ofˆacross a stagnation point, which is consistent with the physical interpretation of the self-sustaining mechanism proposed by Morozov (2022).

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-1 manifold separating two different basins of attraction (Skufca et al. 2006;Schneider et al. 2008;Duguet et al. 2008). This is illustrated by Figure 4 (right) at = 20 and = 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 = 60 (Page et al. 2020), start as stable nodes as increases away from 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 manuscript. Based on the above 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 travelling wave. The threshold amplitude to trigger growth is bounded above by the amplitude of the lower branch state itself, which remains A < 1.05 across the domain of existence of travelling waves. In other words, this domain (shaded bright grey on Figure 1) is nonlinearly unstable when subjected to finite but small amplitude disturbances.

Influence of domain length
This section is dedicated to studying the effect of , and thus the influence of domain size [0, 2 / ] on the TWs. Figure 5 shows how a single branch of travelling waves at = 30 (purple in Figure 1) changes with . It has already been established that the steady arrowhead structure is highly sensitive to domain length in the EIT regime (Dubief et al. 2020). 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 (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 (cf. Figure 5), making the marked 'nonlinearly unstable' region on Figure 1 robust to changes in the assumed periodicity and domain size.

High elasticity regime:
→ 0 The high elasticity regime is difficult to access using time-stepping as it becomes increasing stiff as → 0. The algebraic approach taken here suffers no such problems and we can approach and even consider = 0 (see the next section) without difficulty. The existence of the centre-mode linear instability at = 0 is already known in the limit 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 .
Weakly nonlinear theory predicts supercritical behaviour in the high elasticity ( / ) regime, i.e., along the lower boundary of the linearly unstable region. To probe this, a fifth branch was initiated at fixed = 60, starting upwards in as indicated by the weakly nonlinear analysis (Buza et al. 2021) 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 . 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 , 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 at fixed and lower at fixed as shown by Figure 1. Upon further inspection, it turns out that the lower (secondary) fold shown in Figure 6 at ≈ 30 is purely a feature of the polymer diffusion term 1/( )Δ ⊗ growing artificially large (as is decreased), the effect of which is already known to destroy small scale dynamics (Dubief et al. 2020). It turns out that the point at which the saddle-node bifurcation occurs can be delayed arbitrarily by adjusting in accordance with the variations in to keep the polymer diffusion finite. Numerical experimentation suggested a revised polymer diffusion term of the form for some fixed number . The choice of an inverse scaling with is motivated by observations at = 0 shown in Appendix C. If = 0.005 is enforced for the branch in question, which amounts to fixing the coefficient 1/( ) at the point marked by '+' in  Figure 6, the branch of solutions can be followed down to = 0 along the lower branch (cf. the dashed line in Figure 6).

Results: Travelling waves in the creeping flow limit
= 0 Once = 0 is reached, we redirect the continuation tool towards decreasing . The resulting branch takes the familiar shape (from the > 0 cases) shown in Figure 7, attaining its saddle-node bifurcation point at ≈ 7.5, which serves as a lower bound for the region where travelling waves exist (note the waves found by Morozov 2022, at = 0 are all above = 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 is varied. Arrowhead-shaped structures are still clearly visible at low (cf. panels on the left side of Figure 7), establishing their prevalence even in the high elasticity regime.
For this particular case, specify = 1, = 500, = 0.005), the stability of steady states was examined along the upper branch using DNS. At four points, = 10, 20, 30, 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-, which were pre-multiplied by 10 −6 and added to the travelling waves. All perturbed states returned 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 ≈ 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 2D channel prediction made here of = 7.5. In Figure  S1 of the supplementary material to Pan et al. (2013), the authors indicate the boundary to the observed instability in a vs. perturbation amplitude plane, essentially matching our predictions for the threshold of nonlinear instability given by the More recently, Steinberg and coworkers (Jha & Steinberg 2020; Schnapp & Steinberg 2021) have obtained results in an experimental setup using a channel with a width/height ratio of 7 and so more approximately 2D. However, the viscosity ratio was = 0.74, significantly smaller than the above presented = 0.9. Other recent experiments in a pipe have also considered smaller (e.g. Choueiri et al. 2021, at = 0.56). Motivated by this, we also performed a few TW branch continuations (at = 10, 20, 30 and 50) with decreasing in an attempt to map out the nonlinearly-unstable domain in the − 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 . It transpires that at lower , the solutions are a little more sensitive to the  (2021)), arbitrarily small perturbations are sufficient to trigger these growing, elastic waves.
In an attempt to recover these observations in the present setting, Figure 9 tracks changes in the lower branch TW amplitude, which is an upper bound on the threshold for growth, as is increased. For the parameter combination in question ( = 0.74, = 0), both our standard amplitude measure A and a separate measure for the magnitude of velocity perturbations, û 2 (Ω) , are shown. For > 30, A remains below 1.005, implying that perturbations amounting to 0.5% of the laminar conformation field tr C are sufficient to trigger growth, making this scenario practically indistinguishable from a linearly unstable one. A reaches its minimum roughly around = 100, then -despite setting off in an increasing trend -remains negligibly small for the remainder. Perhaps more in line with experimental results, the minimal disturbance momentumû decreases steadily with respect to .
In addition to the solution branch, Figure 9 displays two state plots at = 20 and = 407, now in terms ofˆ/ , to aid comparison with experimental results. The former, at = 20, still resembles the structural composition of the linearly unstable center eigenmode at higher -chevron shaped streaks remain visible, but are now disconnected at the centerline. By  (2020) in particular). These similarities indicate that the observed stream-wise travelling waves are indeed dynamically connected to the center-mode instability, via a standard subcritical route. However, given that our analysis is 2D, we are unable to identify the span-wise travelling waves of Schnapp & Steinberg (2021), which might serve to explain the discrepancies in threshold amplitude A on Figure 9 -for instance, there might exist a branch of spanwise travelling solutions with a lower threshold amplitude for large enough. This is under investigation. Another recent experimental work observing instabilities in the elastic regime, Choueiri et al. (2021), was conducted in a pipe -precluding direct comparisons with channel computations performed here -at even lower viscosity ratios, ∈ [0.5, 0.6]. Without externally perturbing the system, Choueiri et al. (2021) detected fluctuations at ≈ 5 for = 0.57, = 104 -however, the flow remains laminar at ≈ 3. Computationally, the = 0 branch may be continued up to = 0.57 (for sufficiently small) in our channel flow setting (cf. Figure 8 again), albeit with a slightly larger threshold amplitude than that of Figure 9, which might serve as an explanation for the finite required for 'unperturbed' instabilities (Choueiri et al. 2021). Surprisingly, flow states in Choueiri et al. (2021) still retain their connected, chevron shaped streaks from the center-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 , see the middle panel of Figure 9), but the scenario could be quite different in pipes.

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

Discussion
In this paper, we have used branch continuation to track two-dimensional, finite-amplitude travelling waves in a viscoelastic channel flow, using the FENE-P model. The travelling waves are borne out of the centre mode instability (Khalid et al. 2021a) which is known to be subcritical over large areas of the − parameter space (Buza et al. 2021). Here, we have shown that the TW solution branches extend to significantly lower and than the curve of marginal stability. For instance, we showed that the saddle node at = 1000 drops as low as ≈ 25, while the associated linear instability at this does not occur until ≈ 80. Most significantly, we demonstrated the persistence of the nonlinear TWs at = 0 for a large range of , with the saddle node sitting at ≈ 7.5 − 10 in dilute ( = 0.9) solutions for a range of . Across a broad range of the parameter space, including at = 0, the upper branch TWs resemble the arrowhead structures observed in DNS at higher (Dubief et al. 2020;Page et al. 2020).
A key feature of the solution branch at = 0 is the low amplitude of the lower branch TW across a very large range of . This suggests that only a very weak disturbance would be required to cause the flow to transition to the higher-amplitude upper branch solution, a scenario which would potentially be indistinguishable from a linear instability in an experiment. This observation is consistent 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. In inertia-dominated EIT, instantaneous arrowhead-shaped flow structures resembling the stable upper branch states have been observed numerically (Dubief et al. 2020). 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). At low-, 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 2D 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. Appendix A. Pseudo-arclength continuation The branch continuation routine is launched from a point on the neutral curve (curve of marginal linear stability) with the aid of the weakly nonlinear theory, which provides the very first initial condition near the solution branch of interest. Beyond this point, as parameters are varied in larger increments, supplying sufficiently accurate initial conditions for (3.9) amounts to predicting the shape of the bifurcation branch, precisely the objective of pseudoarclength continuation (see e.g. Dĳkstra et al. 2014). If ∈ { , , } is the parameter allowed to vary in the continuation (all others are fixed), F restricts to a map F : R +1 → R and (3.9) translates to F (a, , , ) = 0.
(A 1) Arclength based techniques interpret the branch of solutions as a curve : R → R +1 embedded in configuration space, which is spanned by (a, , , ) in this particular case. Assume now that a steady state has been computed at point on the branch , which we write as ( ) = (a , , , , ) . A prediction for the solution at the next step, ( +1 ), is given by moving a distance of tangentially along , where denotes the a priori specified step size. Making use of the fact that a solution branch must satisfy F • ≡ 0, the tangent at point is computed according to If no tangent is available at the previous step, the last row is replaced by (0, . . . , 1). The initial prediction for the solution at step + 1 is given bỹ where the upper indices correspond to the number of completed Newton-Raphson iterates (see below), and the twiddle serves to distinguish the converged branch from its approximate counterpart˜. The next step is to employ the Newton-Raphson method using (A 2) as initial condition to obtain an exact solution of (A 1). Knowledge of the tangent may be used to aid this procedure, by means of constraining subsequent iterates to remain on the hyperplane orthogonal to ( ). Incorporating this condition into a standard Newton-Raphson scheme, we obtain the system F (˜( +1 )) ( ) Δ˜+ 1 ( +1 ) = −F (˜( +1 )) − ( ),˜( +1 ) − ( ) , supplemented with the update rulẽ +1 ( +1 ) =˜( +1 ) + Δ˜+ 1 ( +1 ).
Once the solutions have converged to a sufficient degree, i.e., ) = (20, 60, 0.9, 250, 500) ) with the base flow subtracted. The various truncations of branch continuation show very good agreement with the spectra from the DNS. It is worth remarking that to obtain converged spectra in the DNS it was necessary to reduce the timestep to Δ = 1 × 10 −3 . Further discrepancies for the higher order modes may be caused by the finite accuracy of the edge tracking algorithm in locating the lower branch.

Appendix C. Remarks on the polymer diffusion term
Memory limitations associated with increasing and provide an upper bound on how large a value of can be considered whereas taking too small is known to eradicate small scale dynamics. Thus, was selected between the two of these bounds: just above the smallest number where solution branches become independent of , but resolution requirements are still moderate enough for the physical memory to handle.
In the creeping flow limit ( → 0), a distinguished limit with → ∞ such that = 1/( ) with a constant clearly needs to be taken to retain finite polymer diffusion as follows C + (u · ∇) C + T(C) = C · ∇u + (∇u) · C + Δ ⊗ C.
The new polymer diffusion coefficient, , is selected analogously to , whereby independence of solution branches is sought while keeping the grid size manageable. Numerical observations indicate that there cannot be a universal that is optimal in this sense across all parameter regimes within the = 0 limit. However, the slight rescaling = (cf. Section 4.3) seems to help. Figure 12 displays three pairs of solution branches computed at different values of fixed (in green) and fixed (in blue). While the -based branches can be extended up to arbitrarily large (the = 0.005 one was continued up to > 1000 -not shown), the based ones form isolas at low Weissenberg numbers < 30, resulting in a loss of robustness with respect to . Using , robustness is recovered after a certain threshold is hit ( ≈ 0.007 for this branch), which served as the primary motivation behind choosing as the parameter to be fixed throughout the main text.