Subcritical instabilities in plane Poiseuille flow of an Oldroyd-B fluid A.M. acknowledges support from the UK Engineering and Physical Sciences Research Council (grant number EP/I004262/1).

# Subcritical instabilities in plane Poiseuille flow of an Oldroyd-B fluid ††thanks: A.M. acknowledges support from the UK Engineering and Physical Sciences Research Council (grant number EP/I004262/1).

Alexander Morozov Alexander Morozov SUPA, School of Physics and Astronomy, The University of Edinburgh, James Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, United Kingdom
Tel.: +44 131 650 5882
22email: Alexander.Morozov@ed.ac.ukWim van Saarloos Instituut–Lorentz, Universiteit Leiden, Postbus 9506, 2300 RA Leiden, The Netherlands
44email: saarloos@lorentz.leidenuniv.nl
Wim van Saarloos Alexander Morozov SUPA, School of Physics and Astronomy, The University of Edinburgh, James Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, United Kingdom
Tel.: +44 131 650 5882
22email: Alexander.Morozov@ed.ac.ukWim van Saarloos Instituut–Lorentz, Universiteit Leiden, Postbus 9506, 2300 RA Leiden, The Netherlands
44email: saarloos@lorentz.leidenuniv.nl
###### Abstract

Recently, detailed experiments on visco-elastic channel flow have provided convincing evidence for a nonlinear instability scenario which we had argued for based on calculations for visco-elastic Couette flow. Motivated by these experiments we extend the previous calculations to the case of visco-elastic Poiseuille flow, using the Oldroyd-B constitutive model. Our results confirm that the subcritical instability scenario is similar for both types of flow, and that the nonlinear transition occurs for Weissenberg numbers somewhat larger than one. We provide detailed results for the convergence of our expansion and for the spatial structure of the mode that drives the instability. This also gives insight into possible similarities with the mechanism of the transition to turbulence in Newtonian pipe flow.

###### Keywords:
Subrcritical instability visoelastic Poiseuille flow transition to turbulence
journal: Journal of Statistical Physics

## 1 Introduction

Pierre Hohenberg played an important role in the scientific life and career of one of us, Wim van Saarloos, who was a junior colleague of Pierre at Bell Labs from 1982-1990. During these years Pierre worked steadily on his magnum opus on Pattern Formation, the review Pattern formation outside of equilibrium with Mike Cross, published in 1993 cross1993 (). Although the two actually only started to work together during the last two of these years, indirectly and through informal discussions and by Pierre acting as a soundboard, Wim profited a lot from Pierre’s wisdom and insight in pattern formation. Moreover, Pierre’s attention for understanding real experiments, for translating one’s theoretical analysis into experimentally testable predictions, and his emphasis on the importance of writing longer papers and reviews, has had a lasting influence on him. It is therefore a privilege to contribute to this special issue honoring Pierre Hohenberg, with a topic that has many elements of Pierre’s interests and ways of doing physics.

In his personal reflections which will appear in a memorial book for Pierre Hohenberg vansaarloos2018 (), Wim has extensively described Pierre’s influence on him. Moreover, it is explained there in detail how the topic which the two of us, authors of this paper, explored together some ten to fifteen years ago, exhibits traces of the earlier collaborations of Wim and Pierre. We refer to these personal reminiscences vansaarloos2018 () for this background. Connections relevant to the present paper are: the relation with their unpublished explorations of the transition to turbulence in Newtonian pipe flow, how this culminated in their long paper on the quintic Complex Ginzburg-Landau equation vansaarloos1992 (), how this in the end relates to the topic of interest here, and how much Pierre enjoyed it when we presented this whole storyline at the Rutgers meeting celebrating his 80 birthday in December 2014.

The issue at stake is the question whether viscoelastic channel flow exhibits a nonlinear flow-instability in the small Reynolds number limit.

It is well known that Newtonian fluids flowing through a pipe, so-called Poiseuille flow, exhibit a nonlinear instability to turbulence Eckhardt2008 (); barkley2016 (). The transition must be nonlinear, since the laminar Poiseuille flow state is linearly stable. The same holds for Couette flow, flow induced by having two plates moving in opposite direction.111Planar Poiseuille flow actually does become linearly unstable at sufficiently high Reynolds numbers, but this linear instability is actually irrelevant for the transition to turbulence, which happens at much lower Reynolds numbers.

Adding polymers to a fluid can have drastic effects. With polymers in it, the fluid behaves as a so-called viscoelastic fluid. Because shear causes stretching of the polymers, a viscoelastic polymer fluid exhibits elasticity, relaxation and anisotropy: each stretched polymer acts like a little elastic rubber band which is oriented by the shear direction and which takes time to respond (relax) when the local shear rate varies. These effects are stronger, of course, the longer the polymer molecules are, and the higher their concentration. A well-know dramatic manifestation of practical interest of these viscoelastic effects is ’turbulent drag reduction’ White2008 (): when sufficiently long polymers are added in small amounts to a fluid, the drag experienced by turbulent flow through the pipe is reduced. The precise origin of this effect is still a matter of active research.

We are here interested in another limit, the limit in which both the polymer concentration and their length are large enough that the effective fluid viscosity is large, so much so that the flow can be considered as small Reynolds number flow. The convective nonlinearity in the Navier-Stokes equation can then be neglected, and the only nonlinearities come from the so-called constitutive equation that relates the polymer stress tensor and the flow field. In this limit, the only dimensionless quantity for simple shear flows is the so-called Weissenberg number , which is a measure of the normal stress effects in the fluid, resulting from the orientation and stretching of the polymers.

It has been known since 1977 Ho1977 () that viscoelastic Poiseuille flow, as modeled by the so-called Oldroyd-B constitutive equation, is linearly stable in the small Reynolds number limit. Motivated by experimental work by Bonn and co-workers Bertola2003 (); Bonn2011 (), and by the similarity of the questions concerning the nonlinear transition to turbulence of Newtonian fluids, we therefore explored theoretically in a series of papers Meulenbroek2003 (); Meulenbroek2004 (); Morozov2005prl (); Morozov2007 () the question whether parallel viscoelastic shear flows would show a nonlinear flow instability as well.

A combination of analytical arguments, based on insights into the Newtonian case and the theory of dynamical systems and pattern formation Morozov2007 (), and explicit but approximate nonlinear stability calculations for the case of Couette flow, led us to propose Bertola2003 (); Morozov2005prl (); Morozov2007 () that viscoelastic shear flows should indeed show this putative nonlinear transition for Weissenberg numbers somewhat larger than one. We speculated that the transition would lead to turbulent flow.

Several initial attempts to test our proposed scenario experimentally gave negative results. Some numerical investigations did show behavior reminiscent of the proposed scenario, but it is well known to experts that simulations of these bilinear type of constitutive equations are very prone to numerical instabilities, resulting from flow regions which lead to exponential divergence of stress components.Thus, also numerical investigations were considered to give inconclusive results.

The issue therefore remained unsettled till Arratia and co-workers Pan2013 (); Qin2017 () convincingly showed in 2013, using microfluidics experiments, that indeed small Reynolds number viscoelastic channel flow does indeed exhibit a nonlinear transition to turbulence. Qualitatively, their finding is very much in line with what we had proposed: by perturbing the viscoelastic flow in very long micro-channels behind the inlet in a controlled way, by putting a variable number of cylindrical obstacles in the channel, they found that sufficiently strong perturbations lead to a well-developed asymptotic turbulent state far down the channel. Moreover, the experiments showed that the transition occurs for Weissenberg numbers somewhat larger than one, very much like what our own approximate analysis had suggested.

The unequivocal experimental evidence for the nonlinear viscoelastic instability of channel flow leads us to revisit and extend, in this paper, our earlier theoretical analysis. The earlier calculations Morozov2005prl () were focussed on Couette flow (two plates moving in opposite directions); we here extend the results to the experimentally relevant case of Poiseuille flow. We also take the opportunity to present in more detail the nature of our analysis, which is based on an asymptotic Amplitude-equation-like expansion taken to unusually large order (which, incidentally, was an aspect that Pierre liked very much when we presented this at the Rutgers meeting celebrating his 80th birthday). This also allows us to put our results into perspective and to substantiate our previous claim that our results converge to well-defined values and flow profiles. Indeed, we also show for the first time the spatial structure of the nonlinear modulated waves which according to our analysis determine the instability threshold.

In the next section, we simply summarise the main result of our analysis, which shows within the context of the same approximate analysis that viscoelastic Poiseuille flow between plates exhibits a nonlinear instability which is very similar to the one for plane Couette flow published previously Morozov2005prl (). Also the transition amplitude and the values of the Weissenberg number where the instability is found, is similar. In section 3 we then present the technical details of the expansion underlying our results, which we then apply to the case of plane Poiseuille flow in section 5. In our concluding section 6 we put our results into perspective, and speculate on the transition to turbulence mechanism of viscoelastic flow.

## 2 Main result of perturbative expansion for viscoelastic plane Couette and plane Poiseuille flow

Figure 1 summarizes our earlier result Morozov2005prl () for plane Couette flow and our extensions to the case of Poiseuille flow between plates (right panel) as a function of increasing Weissenberg number which, as stated above, is the only dimensionless number for these flow configurations in the zero Reynolds number limit. What do these results mean?

As detailed in later sections, these calculations are based on perturbing the laminar viscoelastic flow with a finite-wavelength mode of amplitude . For flow between plates, one does this by picking a wavenumber associated with the wavenumber of the modulation along the flow direction, and a wavenumber in the transverse direction. We imagine these to be fixed, and then ask ourselves whether, if we think of the flow to be perturbed by a modulation of this type and of a given initial amplitude , this amplitude would decay — meaning stability — or increase in time (meaning instability). The linear stability calculation, already done long ago Ho1977 (), is based on assuming that the mode amplitude (for a given and ) is infinitesimally small, and the finding from these calculations is that when the initial amplitude is infinitesimally small, it will indeed decay in time: the flow is linearly stable.

Our calculations are based on probing whether will grow or decay in time perturbatively, by performing an expansion in powers of , for fixed wavenumbers. In other words, we only probe growth or decay of the amplitude, without allowing for spatial instabilities of the modes themselves.

In Figure 1, we plot along the vertical axis the critical dimensionless value of the amplitude — which we can think of the the maximum relative change in shear rate at the walls — separating decay (for smaller amplitudes) from growth (for larger amplitudes) for particular values of the wavenumbers indicated.

Clearly, these results indicate that within the limitations of our perturbative calculations, we expect that nonlinear instability to be very similar for plane Couette and Poiseuille flow. This is of course what one would expect physically as the instability is driven by the self-amplifying nature of viscoelastic flow along curved streamlines. This should depend little on the details of the way the flow is driven.

Of course, as explained above our analysis is based on an amplitude expansion, which itself is an asymptotic expansion. It is legitimate to ask why we were confident to draw conclusions from such an expansion. To illustrate why we did dare to do so, we show in Figure 1 the results up to 5 (red, label ), 7 (orange, label ) , 9 (light blue, label ) and even 11 order (dark blue, label ) in our expansion for for the case of plane Couette flow. Quite surprisingly — we had no reason to expect this a priori — the results for the critical value are quite robust. Moreover, even the ’nose’ of the curves on the left seems to converge: whenever the expansion is extended by including another order, the curve shifts by about half the amount that it did in the previous step.

## 3 Derivation of the amplitude equation

In this Section we present a method of approximating non-linear solutions for a class of non-linear partial differential equations in the following form

 ^LV+^A∂V∂t=N(V,V), (1)

where is a -dimensional vector of fields, and are linear operators, and is a quadratic non-linear operator. Our method is primarily designed for problems of hydrodynamic stability in parallel shear flows, and, from the onset, we incorporate some of the main features of the corresponding equations of motion directly into Eq.(1). First, we introduce a Cartesian coordinate system , where and denote the unbounded directions, while is the direction between two confining plates. The vector is thought of to contain the velocity, stress components, and pressure in the fluid. We further assume that , , and only contain spatial derivatives, and that , where and are arbitrary functions of space, while is an arbitrary function of time.

Our goal is to demonstrate that Eq.(1) can have finite-amplitude solutions besides the trivial (laminar) one. In the absence of a linear instability, there is no systematic way of finding such a solution and we attempt to construct it perturbatively, in a way motivated by amplitude expansion methods developed for studying near-threshold behaviour of pattern-forming instabilities (see cross1993 (); vanHecke1994 (), for example). Our expansion is based on the eigenfunctions, , of the linear operator:

 ^Lei(kxx+kzz)V(n)0(y)=−λn^Aei(kxx+kzz)V(n)0(y). (2)

The index enumerates the eigenfunctions; their number and form depends on the problem at hand. The desired finite-amplitude solution for the real-valued fields is then assumed to be represented by

 V(x,y,z,t)=Φ(t)eiξV0(y)+Φ∗(t)e−iξV∗0(y) +U0(y,t)+∞∑n=2[Un(y,t)einξ+U∗n(y,t)e−inξ], (3)

where is one of the eigenfunctions, and we have introduced ; ”” denotes complex conjugation; the choice of a particular eigenfunction will be discussed below, but we note here that for plane Couette flow these are the Gorodtsov-Leonov modes Gorodtsov1967 () that we used in our previous work Morozov2005prl (). One can view Eq.(3) as a Fourier expansion in and with an extra assumption about the form and interrelation between the coefficients; is the basic amplitude of the mode whose non-linear behaviour we aim to study perturbatively. In the spirit of amplitude-equation techniques cross1993 (); vanHecke1994 (), the time-dependent amplitude follows the linear dynamics governed by Eq.(2) at short times, which is replaced by nonlinear dynamics at large . In the following, we derive the evolution equation for the amplitude assuming that it is small in some sense. This assumption will be checked for self-consistency, once the solution is obtained.

To proceed, we assume that the dominant dynamics in Eq.(3) are given by the time-evolution of , and that the higher harmonics ’s follow it adiabatically (they are ’slaved’ to ). This implies that the Fourier components ’s can only arise as a result of nonlinear interactions of the eigenmode with itself and other ’s. For instance, if we denote , then has contributions from the interactions of a) with , b) , , , and , etc. Applying this power-counting argument to all Fourier modes, we obtain

 U0(y,t)=|Φ(t)|2u(2)0(y)+|Φ(t)|4u(4)0(y)+⋯, U2(y,t)=Φ2(t)u(2)2(y)+Φ2(t)|Φ(t)|2u(4)2(y)+⋯, (4) U3(y,t)=Φ3(t)u(3)3(y)+⋯, ⋯,

where are unknown functions, which will be determined below; the subscript and superscript correspond to the Fourier harmonic and the order in , respectively.

To derive the time-evolution equation for , we substitute Eq.(3) into Eq.(1) and separate the terms proportional to . This yields

 ^L(Φ(t)eiξV0(y))+dΦ(t)dt^A(eiξV0(y))=(dΦ(t)dt−λΦ(t))^A(eiξV0(y)) =¯N(Φ(t)eiξV0(y),U0(y,t))+¯N(Φ∗(t)e−iξV∗0(y),U2(y,t)e2iξ) +∞∑n=2¯N(Un+1(y,t)ei(n+1)ξ,U∗n(y,t)e−inξ), (5)

where we used Eq.(2) and introduced . Although the eigenmodes of the linear operators involved in flow stability problems often form full sets, they are typically non-normal SchmidHenningson2000 () and their eigenmodes are not orthogonal. Therefore, to calculate the contribution of the r.h.s. of Eq.(5) to , we employ the adjoint operator , defined via

 ⟨V1|^LV2⟩=⟨^L†V1|V2⟩, (6)

where and are two arbitrary vectors satisfying proper boundary conditions. The scalar product is given by

 ⟨V1|V2⟩=limLx,Lz→∞12Lx∫Lx−Lxdx12d∫d−ddy12Lz∫Lz−Lzdz(V∗1,V2), (7)

where is the Euclidean dot product, and is the distance between the confining plates. The actual form of the adjoint operator is obtained by performing integration-by-parts in Eq.(6), as will be demonstrated later. The eigenmodes and the eigenvalues of the adjoint operator are given by

 ^L†ei(kxx+kzz)W(m)0(y)=−χm^A†ei(kxx+kzz)W(m)0(y), (8)

with , and being the operator adjoint to ; for the problem we consider in Section 4, operator is self-adjoint. Since Eq.(2) is a generalised eigenvalue problem, the orthogonality condition between the eigenmodes of the linear and adjoint operators is somewhat unusual, and we state it here explicitly. We consider

 ⟨eiξW(m)0(y)|^LeiξV(n)0(y)⟩=−λn⟨eiξW(m)0(y)|^AeiξV(n)0(y)⟩. (9)

At the same time,

 ⟨eiξW(m)0(y)|^LeiξV(n)0(y)⟩=⟨^L†eiξW(m)0(y)|eiξV(n)0(y)⟩ =−χ∗m⟨^A†eiξW(m)0(y)|eiξV(n)0(y)⟩=−χ∗m⟨eiξW(m)0(y)|^AeiξV(n)0(y)⟩. (10)

Comparing Eqs.(9) and (10) we conclude that , unless .

The evolution equation for the amplitude is finally obtained by projecting Eq.(5) onto the eigenmode of the adjoint operator , selected such that its eigenvalue . According to Eq.(4), the r.h.s. of Eq.(5) is a polynomial in and its complex conjugate, and one obtains

 dΦdt=λΦ+C3Φ|Φ|2+C5Φ|Φ|4+C7Φ|Φ|6+C9Φ|Φ|8+C11Φ|Φ|10⋯. (11)

The coefficients ’s are calculated by collecting terms of the corresponding order in on the r.h.s. of Eq.(5) and projecting them on . Once these coefficients are known, we can then study whether there is a critical value that separates decaying amplitudes, for small , from the growing ones.

To illustrate the method, we show here how to calculate the coefficient , while the expressions for the higher-order coefficients are deferred to Appendix C. To , we obtain from Eq.(5) by projecting it on

 C3=1Δ⟨eiξW0(y)∣∣∣¯N(eiξV0(y),u(2)0(y))+¯N(e−iξV∗0(y),e2iξu(2)2(y))⟩, (12)

where

 Δ=⟨eiξW0(y)|eiξ^AV0(y)⟩. (13)

Let us now determine the unknown functions and . Once again, we substitute Eq.(3) into Eq.(1), and separate the terms proportional to . To lowest order in , it yields

 ^L(|Φ(t)|2u(2)0(y))+d|Φ(t)|2dt^Au(2)0(y)=¯N(Φ(t)eiξV0(y),Φ∗(t)e−iξV∗0(y)). (14)

The time-derivative can be evaluated self-consistently with the help of the amplitude equation (11), and to is equal to

 d|Φ|2dt=dΦdtΦ∗+ΦdΦ∗dt=(λ+λ∗)|Φ|2. (15)

Therefore, satisfies the following inhomogeneous ODE

 ^Lu(2)0(y)+(λ+λ∗)^Au(2)0(y)=¯N(eiξV0(y),e−iξV∗0(y)). (16)

Similarly, is given by

 ^L(e2iξu(2)2)+2λ^Ae2iξu(2)2=N(eiξV0,eiξV0). (17)

The equations for higher ’s and ’s are obtained by a straightforward, though lengthy, generalisation to higher orders in of the procedure outlined above. The corresponding expressions can be found in Appendix C.

Equation (11), together with the procedure to systematically determine the coefficients ’s, is the central result underlying our non-linear analysis of the flow stability. It allows one to calculate the amplitude of a non-linear solution to Eq.(1) in the form given by Eq.(3). In what follows, we will be particularly interested in simple solutions to this equation, either in the form of stationary points or travelling waves. As we will see below, for the problem discussed in Section 4, the least stable eigenvalues are complex, and, therefore, the relevant solutions are of the latter type, given by , where and are real numbers. Substituting this ansatz into the amplitude equation (11), and separating the real and imaginary parts, we obtain for a non-trivial solution with a stationary amplitude

 0=Re(λ)+Re(C3)Ψ2+Re(C5)Ψ4+⋯, (18) Ω=Im(λ)+Im(C3)Ψ2+Im(C5)Ψ4+⋯. (19)

The asymptotic nature of Eqs.(18) and (19) imply that only converging series can represent a physical solution. In turn, this translates into the requirement that the solution amplitude is sufficiently small, where the scale is given by the coefficients . To study convergence of series Eqs.(18) and (19), we employ a somewhat intuitive method based on the partial sums , defined through

 Sm≡Re(λ)+Re(C3)Ψ2+⋯+Re(C2m+1)Ψ2m. (20)

We then solve a series of algebraic equations , , , and obtain a corresponding series of solutions , , , and , , , using Eq.(19). If both series approach limiting values, the latter are recognised as representing a physical solution; see Section 2 for discussion.

In the next Sections we adopt this method to parallel shear flows of model viscoelastic fluids and calculate the coefficients in Eq.(11) up to .

## 4 Channel flow of a viscoelastic fluid

Here we adopt the method presented in the previous Section to the study of non-linear stability of parallel shear flows of model polymer fluids. We consider a flow between two parallel plates forced by either a constant pressure gradient (plane Poiseuille flow). As in Section 3, we select a coordinate system with being along the streamwise, vertical (gradient) and spanwise directions, respectively; the distance between the plates is .

The velocity of the fluid is governed by the Stokes equation

 −∇p+ηs∇2v+∇⋅τ=0, (21)

and the incompressibility condition

 ∇⋅v=0, (22)

where is the pressure, is the polymeric contribution to the total stress, and is the solvent viscosity. In Eq.(21), we neglected the inertial terms as typical experiments on elastic instabilities and turbulence are usually performed either in microfluidic devices or with high-viscosity fluids (see Groisman2004 (); Pan2013 (), for example). In both cases, the corresponding values of the Reynolds number (defined as , where is the typical flow velocity and is the kinematic viscosity of the fluid) do not exceed , and the inertial effects can be ignored.

To describe the dynamics of the polymer stress tensor , we employ the Oldroyd-B model given by

 τ+λ[∂τ∂t+v⋅∇τ−(∇v)T⋅τ−τ⋅(∇v)]=ηp[(∇v)+(∇v)T], (23)

where is the Maxwell relaxation time of the fluid, is the polymer viscosity, and denotes the matrix transpose. The Oldroyd-B model is the simplest equation incorporating normal-stress effects that are the driving force for many viscoelastic flow instabilities Larson1992 (); Shaqfeh1996 (); Morozov2007 (); for a detailed discussion of various viscoelastic equations of motion and their predictions see Bird1987 (); Larson1988book (); Morozov2015 (). Finally, the velocity field is assumed to satisfy the no-slip boundary condition

 vx(x,y=±d,z)=vy(x,y=±d,z)=vz(x,y=±d,z)=0. (24)

Equations (21)-(24) have as laminar solution the well-known parabolic profile with , where

 U0(y)=ΔPd22(ηs+ηp)[1−(yd)2]. (25)

The corresponding elastic stresses are given by

 τ(0)xx=2ηpλMU′20(y), τ(0)xy=ηpU′0(y), (26) τ(0)xz=τ(0)yy=τ(0)yz=τ(0)zz=0.

In what follows we render equations (21)-(24) dimensionless by re-scaling the variables by appropriately chosen units, suggested by the laminar solution above. We use as a unit of length, as a unit velocity, as a unit of time, and as a unit of stress. The material properties of the fluid are controlled by two dimensionless parameters: the ratio of the solvent to total viscosities , and the Weissenberg number . The latter controls the strength of non-Newtonian effects in Eq.(23) and plays in elastic instabilities and turbulence the same role as the Reynolds number does in Newtonian fluid mechanics.

Next, we split the dimensionless velocity and stress fields into the laminar part and a deviation from the laminar solution,

 τij=2WiU′20(y)δixδjx+U′0(y)(δixδjy+δiyδjx)+τ(1)ij, (27) v=(U0(y),0,0)+v(1), (28)

and introduce a perturbation vector , where is the pressure perturbation. Substituting these expressions into the equations of motion, we arrive at the compact form, Eq.(1), introduced in Section 3. The explicit expressions for , and are given in Appendix A. In the next Section we present the finite-amplitude solutions of Eqs.(21)-(24) found by the method presented in Section 3.

## 5 Results

As the eigenmodes of the linear operator in Eq.(2) are the starting point of our analysis, here we briefly present the linear stability analysis of plane Poiseuille flow of an Oldroyd-B fluid; for a detailed discussion see Wilson et al. Wilson1999 (). To calculate the eigenspectrum, we discretise Eq.(2) with and given in Appendix A using the fully spectral Chebyshev-tau method Canuto:book (), and solve numerically the resulting generalised eigenvalue problem using Scientific Python scipy ().

In Fig.2 we present an example of the eigenvalue spectrum, plotted as vs . The general structure of the spectrum for , , , and is shown in Fig.2a), while Fig.2b) presents a zoom-in on the least stable part of the spectrum. As can be seen, all eigenvalues have , indicating linear stability. This is confirmed numerically for all values of , , and ; see also Wilson et al. Wilson1999 (). The most prominent features in Fig.2, the balloon-like shapes at and , are numerical approximations to the continuous spectrum of the linear operator, which corresponds to the singular points of the linear problem Wilson1999 (). The least unstable eigenvalues, which are used in the non-linear analysis below, are the right-most modes in Fig.2b), which we denote by , , and . (There is another discreet eigenvalue which is very close to the continuous spectrum, and we do not attempt to perform calculations for the associated eigenmode as it would require a very high spectral resolution.) The corresponding eigenmodes are presented in Fig.3. As can be seen there, two of the eigenmodes are mostly pronounced close to the walls, and are related to the Gorodtsov-Leonov modes Gorodtsov1967 () used in our previous work on plane Couette flow Morozov2005prl (), while the other one is mostly present in the middle of the geometry.

Next, we calculate the eigenmodes of the adjoint problem, Eq.(8), with defined in Appendix B, by using the same numerical method, as above. As demonstrated in Section 3, every eigenvalue of the linear problem has its adjoint counterpart , such that . Therefore, for the same parameters as above, the adjoint spectrum looks like Fig.2, with . The adjoint eigenmodes, corresponding to the least stable eigenvalues , , and share the same spatial features as their linear counterparts and are either confined to the walls or the bulk of the system; see Fig.4, for example.

In what follows, we use the three least stable modes as the starting point of the non-linear analysis presented in Section 3, and assess whether they result in converged non-linear states. For each mode, we calculate the coefficients ’s as a function of the Weissenberg number , and use them to solve Eq.(18) for the amplitude of the travelling-wave state. As explained in Section 3, we construct a series of solutions at progressively higher orders in the amplitude, see Eq.(20), and study their convergence. We found that using the eigenmode associated with does not lead to a converging series of amplitudes , while the other two eigenmodes lead to converging non-linear solutions, which we refer to as ’State 1’ and ’State 3’, and below we focus on these two modes.

The results of this procedure are presented in Fig.5, while the values of the coefficients ’s for selected values of are given in Tables 1 and 2. For both states, the series of solutions share the same features. Here, denotes the order to which the expansion in powers of in Eq. (11) is taken. For any , the equation has no real solutions at small , two real solutions around the saddle-node bifurcation – the value of at which two solution branches appears for the first time, and one real solution for larger values of . The lower branches define the threshold amplitude required to destabilise the flow, while the upper branches are supposed to set the saturated value of the amplitude at the instability. As can be seen from Fig.5, the upper branches of all solutions seem to diverge rapidly close to the saddle-node value , and the implications of this behaviour are discussed below.

For both states, we observe that the upper- and lower-branch values of converge rapidly as is increased from to . In Fig.6 we assess this convergence more quantitatively, by plotting the saddle-node values for each as a function of (circles in Fig.6). In the range , corresponding to the amplitude equation expansion up to the eleventh order, the convergence is well-described by an exponential fit (dashed lines in Fig.6), approaching and for States 1 and 3, correspondingly. The phase speed of the travelling-wave solutions, , calculated using the values of presented above are shown in Fig.7. Again, we observe that the position of the saddle-node and the lower-branch values converge rapidly. Using the criterion presented in Section 3, we conclude that our consecutive approximations to States 1 and 3 converge towards physical solutions that we now examine in more detail.

To gain insight into the spatial structure of the travelling-wave solutions reported above, we now plot the mean velocity profile for State 1, calculated by averaging over and . In Fig.8a) we present the mean profile for , which is close enough to the saddle-node point for so that has two values, and , corresponding to the lower- and upper-branches, respectively. While the lower-branch mean velocity profile is quite close to the laminar profile , the upper-branch profile looks very different, with sharp velocity gradients close the walls and a shifted parabolic-like profile in the bulk of the system. Surprisingly, the centre-line velocity appears to be larger than its laminar counterpart. In Fig.8b) we present the velocity profile at in the -plane, which is perpendicular to the streamwise direction: the in-plane components and are shown as arrows, while the streamwise velocity is given by the colour. One can clearly see two arrays of streamwise vortices next to each wall superimposed onto the corresponding arrays of high- and low-velocity streamwise streaks. Fig.8c) shows a similar velocity profile at in the -plane, where arrows now trace the in-plane components and , while the colour gives the spatial profile of the spanwise velocity . Finally, in Fig.8d) we plot the largest component of the polymeric stress tensor, , which is a deviation from the laminar stress, in the -plane. Most of the stress is concentrated close to the boundaries consistent with the presence of sharp velocity gradients there.

The profiles presented in Fig.8b)-d) correspond to the upper branch of State 1. The lower branch profiles have a very similar structure, albeit with a significantly smaller amplitude, and are not presented here. In a similar fashion, the mean profiles and the spatial structure of State 3 bear strong similarities with State 1, and are, therefore, also omitted here.

Finally, the two states presented here are shown for the somewhat randomly selected values of the wavenumbers and . Preliminary studies show that for small values of , converged solutions similar to States 1 and 3 persist for a wide range of wavenumbers, provided both and are not too big (typically, smaller than . For larger values of , this region shrinks, which is either related to the convergence properties of our technique, or is connected to the actual shrinking of the region of existence of such travelling-wave solutions. This point is deferred to future studies.

## 6 Discussion and Conclusion

Results presented in this work further corroborate our previous claim Meulenbroek2003 (); Bertola2003 (); Morozov2005prl (); Morozov2007 () that while parallel shear flows of viscoelastic fluids are linear stable, they exhibit sub-critical instabilities that lead directly to a chaotic state. Early experiments by Bonn et al. Bonn2011 () and in particular the more recent detailed and systematic experiments by Arratia et al. Pan2013 (); Qin2017 () confirm the existence of such sub-critical instabilities in channel flows of dilute polymer solutions and demonstrates that the chaotic state observed there is related to the phenomenon of purely elastic turbulence previously only reported in shear flows with curved streamlines Groisman2000 (); Groisman2004 (); Burghelea2007 ().

The emergent scenario of the transition in parallel shear flows of viscoelastic fluids parallels that for their Newtonian counterparts. Recently, significant progress was made in understanding the transition to turbulence in pipe, plane Couette and channel flows of Newtonian fluids by studying it from the dynamical systems’ point of view Eckhardt2008 (); barkley2016 (). The key ingredients there are the so-called coherent structures, the exact solutions of the Navier-Stokes equations, either travelling waves or periodic orbits, that are three-dimensional but relatively simple. These solutions appear through a sub-critical bifurcation and correspond to saddles in the phase space of the flow: while their upper and lower branches are linearly unstable, their vicinity is attractive; also, their number increases with the Reynolds number. When the phase space is sufficiently populated with such solutions, a turbulent trajectory performing a random walk between a large number of saddle-like states gets trapped for a very long time. Coupled to the phenomenon of splitting of localised exact coherent states, this scenario firmly places the transition to Newtonian turbulence within the directed percolation universality class.

Our results indicate that a similar scenario might also be at work in the viscoelastic case. The solutions presented here and in our previous work Morozov2005prl () might form the phase space scaffold of the purely elastic turbulence, and verifying their existence should be paramount to understanding its mechanism. The amplitude-equation type technique employed here attempts to construct a non-linear solution as an asymptotic series, and, as such, it has a limited radius of convergence. While the lower branches of our solutions, indicative of the amplitude threshold required to trigger the turbulent state, are well-converged, the upper branches disappear in a close vicinity of the saddle-node bifurcation. This is either related to the radius of convergence of the asymptotic series, Eq.(11), as mentioned above, or can be the direct consequence of the upper-branch non-linear state being turbulent, and the failure of the technique to capture it.

The solutions that we found in this work appear at the values of the Weissenberg number, , that are somewhat lower than the onset of turbulence values reported in experiments in channel, Pan2013 (); Qin2017 (), and pipe flows, Bertola2003 (); Bonn2011 (). This is not surprising since the fluids used in those works had significantly higher values of than used in this work, and, moreover, exhibited various degrees of shear thinning, which is not included in the present analysis. Both factors would result in higher values of than presented here. Also, the dynamical systems scenario of the transition, if applicable in the viscoelastic case, would imply that should be smaller than the value of at which sustained turbulence can be observed, providing yet another explanation for the difference.222One should also keep in mind that we study a mode which is fully periodic in the streamwise direction, whereas in the experiments the flow is perturbed by a small number of cylindrical obstacles in the flow channel. We do not argue in favour of either of the explanations and, instead, draw a more conservative conclusion that our work provides further evidence for a subcritical (nonlinear) instability scenario at moderate values of the Weissenberg number, and that it demonstrates that exact coherent solutions do exist in this type of viscoelastic flows and proper numerical investigation of such states is required.

While experiments provide very limited information about the spatial structure of purely elastic turbulence in parallel shear flows Bertola2003 (); Bonn2011 (); Pan2013 (); Qin2017 (), our work sheds light on what profiles might be expected in such geometries. First, we observe that, according to our calculations, the turbulent mean velocity profiles should exhibit a larger centre-line velocity than their laminar counterparts, see Fig.8a). This is in stark contrast with the Newtonian turbulent mean profiles, which are plug-like due to the momentum re-distribution between the walls and the bulk, and are always slower than the corresponding laminar velocity field in the middle of the gap. Intriguingly, a similar profile was reported under certain conditions in recent experiments on pipe and channel flows of rather concentrated polymer solutions Poole2016 (), where it was attributed to shear-thinning. The profiles reported here suggest that it might be a more generic feature related to the elasticity of the fluid. Our second observation comprises the existence of streamwise vortices and streaks in a plane perpendicular to the flow direction, Fig.8b). These structures are a hallmark of Newtonian coherent structures, and it appears that they also play a role in viscoelastic non-linear solutions. This is, perhaps, not surprising since they feature prominently in non-normal growth analysis by Kumar and Jovanović Jovanovic2010 (); Jovanovic2011 (). The associated stresses (not shown) are in line with the positions of large velocity gradients in-between the streamwise vortices, although both structures are tilted due to the travelling-wave nature of whole the state. Finally, the plane perpendicular to the spanwise -direction contains widely-spaced co-rotating vortices interlaid with expanding-contracting streamlines and large, wall-localised stresses. Although such structures have not been reported in the literature, they are consistent with the near-wall mixing patterns observed by Qin et al. Qin2017 ().

Although we are confident that we obtained converged non-linear states, their existence needs to be verified numerically, by searching for steady-state, travelling-waves, and periodic orbits, using a Newton-Raphson-type algorithm. Until recently, such calculations were not feasible as a Newton-Raphson step is akin to a time-iteration step for the same equations, and viscoelastic constitutive models are notoriously difficult to time-step at sufficiently high Weissenberg numbers (the so-called High-Weissenberg Number Problem OwensPhillips ()). In the past years, there emerged a class of numerical techniques to ensure positive-definiteness of the conformation tensor (absence thereof was implicated as a cause of the High- Weissenberg Number Problem), led by the log-conformation algorithm Fattal2005 (); Pimenta2017 (), and a combination of such techniques with a Newton-Raphson-type algorithm should be able to overcome this problem. The states predicted in this work can then serve as a good initial guess for the search algorithm. Once found, upper branches of these solutions should be studied for their linear stability that will assess whether the transition mechanism in viscoelastic fluids bears similarities with its Newtonian counterpart. This work will be a subject of our future studies.

## Appendix A Explicit form of the matrix equation

As discussed in the main text, the equations of motion Eqs.(21)-(24) can be written in the matrix form

 ^LV+^A∂V∂t=N(V,V), (29)

where is the dimensionless deviation of the hydrodynamic fields from their laminar values (we have dropped the superscript ”1” in Eqs.(27) and (28) to simplify notation). The operator is given by a constant diagonal matrix

 Widiag(0,0,0,0,1,1,1,1,1,1), (30)

while the linear operator is define by it action on

 ^LV= (31)

Here, and , and the dimensionless laminar velocity profile is given by .

The r.h.s. of Eq.(29) represents the non-linear terms in the original equations Eqs.(21)-(24), and is given by a bilinear form

 N(V(A),V(B))=−(v(A)⋅∇)^AV(B) +Wi⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝00002[τ(A)xx∂xv(B)x+τ(A)xy∂yv(B)x+τ(A)xz∂zv(B)x]τ(A)xx∂xv(B)y−τ(A)xy∂zv(B)z+τ(A)xz∂zv(B)y+τ(A)yy∂yv(B)x+τ(A)yz∂zv(B)xτ(A)xx∂xv(B)z+τ(A)xy∂yv(B)z−τ(A)xz∂yv(B)y+τ(A)yz∂yv(B)x+τ(A)zz∂zv(B)x2[τ(A)xy∂xv(B)y+τ(A)yy∂yv(B)y+τ(A)yz∂zv(B)y]τ(A)xy∂xv(B)z+τ(A)xz∂xv(B)y+τ(A)yy∂yv(B)z−τ(A)yz∂xv(B)x+τ(A)zz∂zv(B)y2[τ(A)xz∂xv(B)z+τ(A)yz∂yv(B)z+τ(A)zz∂zv(B)z]\par⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

Obviously, .

The actual form of the adjoint operator is obtained by using in Eq.(6) the expression for the linear operator , defined above, and performing integration-by-parts. Again, we define by its action on a mode , which is given by

 ^L†W= (32)

where . The adjoint velocities are subject to the boundary conditions

 νx=νy=νz=0,y=±1, (33)

which follow from the requirement that the boundary terms from the integration by parts in Eq.(6) vanish for any satisfying Eq.(24).

## Appendix C Higher-order coefficients for the amplitude equation

In Section 3 we demonstrated how the coefficient can be determined from Eq.(5). Here, we list the expressions for higher-order coefficients ’s and the associated functions from Eq.(4), obtained in the manner discussed in Section 3.

To , we have

 ^L(e3iξu(3)3)+3λ^Ae3iξu(3)3=¯N(eiξV0,e2iξu(2)2), (34)
 ^L(e2iξu(4)2) +(3λ+λ∗)^Ae2iξu(4)2 =¯N(e−iξV∗0,e3iξu(3)3)+¯N(u(2)0,e2iξu(2)2)−2C3e2iξu(2)2, (35)
 ^Lu(4)0+2(λ+λ∗)^Au(4)0=N(u(2)0,u(2)0)+¯N(e2iξu(2)2,e−2iξu(2)∗2)−(C3+C∗3)u(2)0, (36)

which yields

 C5=1Δ⟨eiξW0∣∣∣¯N(eiξV0,u(4)0)+¯N(e−iξV∗0,e2iξu(4)2)+¯N(e−2iξu(2)∗2,e3iξu(3)3)⟩, (37)

where, as in Section 3,

 Δ=⟨eiξW0(y)|eiξ^AV0(y)⟩. (38)

To , we have

 ^L(e4iξu(4)4)+4λ^Ae4iξu(4)4=¯N(eiξV0,e3iξu(3)3)+N(e2iξu(2)2,e2iξu(2)2), (39)
 ^L(e3iξu(5)3) +(4λ+λ∗)^Ae3iξu(5)3 =¯N(eiξV0,e2iξu(4)2)+¯N(e−iξV∗0,e4iξu(4)4)+¯N(u(2)0,e3iξu(3)3) −3C3e3iξ<