Subcritical instabilities in plane Poiseuille flow of an Oldroyd-B fluid

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.

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 [2], 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 [2] 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 [3], 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 th 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 [4,5]. 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. 1 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' [6]: 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 Wi, 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 [7] 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 [8,9], 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 [10,11,12,13] 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 [13], and explicit but approximate nonlinear stability calculations for the case of Couette flow, led us to propose [8,12,13] 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 [14,15] 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 [12] 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 Amplitudeequation-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 [12]. 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.
Steady-state amplitude of the travelling-wave solution in plane Couette flow of UCM fluids found in [16] for kx = kz = 1; the ratio Re/Wi = 10 −3 was kept constant. Re-plotted from Ref. [16]. For definitions see Section 4.
2 Main result of perturbative expansion for viscoelastic plane Couette and plane Poiseuille flow Figure 1 summarizes our earlier result [12] for plane Couette flow and our extensions to the case of Poiseuille flow between plates (right panel) as a function of increasing Weissenberg number Wi 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 kx associated with the wavenumber of the modulation along the flow direction, and a wavenumber kz 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 stabilityor increase in time (meaning instability). The linear stability calculation, already done long ago [7], is based on assuming that the mode amplitude Φ (for a given kx and kz) 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 Ψ = |Φc| 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 th (red, label Ψ 2 ), 7 th (orange, label Ψ 3 ) , 9 th (light blue, label Ψ 4 ) and even 11 th order (dark blue, label Ψ 5 ) in our expansion for Φc for the case of plane Couette flow. Quite surprisingly -we had no reason to expect this a priori -the results for the critical value Φc 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.

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 where V is a d-dimensional vector of fields,L andÂ are linear operators, and N 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 (x, y, z), where x and z denote the unbounded directions, while y is the direction between two confining plates. The vector V is thought of to contain the velocity, stress components, and pressure in the fluid. We further assume thatL,Â, and N only contain spatial derivatives, and that N ( , where V 1 and V 2 are arbitrary functions of space, while f (t) 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 [1,17], for example). Our expansion is based on the eigenfunctions, e i(kxx+kzz) V (n) 0 (y), of the linear operator:L The index n 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 where e iξ V 0 is one of the eigenfunctions, and we have introduced ξ = kxx + kzz; " * " 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 [18] that we used in our previous work [12]. One can view Eq.(3) as a Fourier expansion in x and z with an extra assumption about the form and interrelation between the coefficients; Φ(t) is the basic amplitude of the mode whose non-linear behaviour we aim to study perturbatively. In the spirit of amplitude-equation techniques [1,17], the time-dependent amplitude Φ(t) follows the linear dynamics governed by Eq.
(2) at short times, which is replaced by nonlinear dynamics at large t. In the following, we derive the evolution equation for the amplitude Φ(t) 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 Φ(t), and that the higher harmonics Un's follow it adiabatically (they are 'slaved' to Φ(t)). This implies that the Fourier components Un's can only arise as a result of nonlinear interactions of the eigenmode with itself and other Applying this power-counting argument to all Fourier modes, we obtain where u (m) n (y) are unknown functions, which will be determined below; the subscript n and superscript m correspond to the Fourier harmonic and the order in Φ, respectively.
To derive the time-evolution equation for Φ(t), we substitute Eq.(3) into Eq.(1) and separate the terms proportional to e iξ . This yieldŝ where we used Eq.(2) and introducedN (A, Although the eigenmodes of the linear operators involved in flow stability problems often form full sets, they are typically non-normal [19] and their eigenmodes are not orthogonal. Therefore, to calculate the contribution of the r.h.s. of Eq.(5) to dΦ/dt, we employ the adjoint operatorL † , defined via where V 1 and V 2 are two arbitrary vectors satisfying proper boundary conditions. The scalar product · · · is given by where (A, B) = i A i B i is the Euclidean dot product, and 2d 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 e i(kxx+kzz) W (m) 0 (y) and the eigenvalues χm of the adjoint operator are given byL with m = 1, 2, . . . , 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 At the same time, Comparing Eqs. (9) and (10) we conclude that e iξ W (m) 0 (y)|Âe iξ V (n) 0 (y) = 0, unless λn = χ * m . The evolution equation for the amplitude Φ(t) is finally obtained by projecting Eq.(5) onto the eigenmode of the adjoint operator e iξ W 0 (y), selected such that its eigenvalue χ = λ * . According to Eq.(4), the r.h.s. of Eq.(5) is a polynomial in Φ(t) and its complex conjugate, and one obtains dΦ dt The coefficients C's are calculated by collecting terms of the corresponding order in Φ on the r.h.s. of Eq.(5) and projecting them on e iξ W 0 (y). Once these coefficients are known, we can then study whether there is a critical value Φc that separates decaying amplitudes, for small |Φ|, from the growing ones.
To illustrate the method, we show here how to calculate the coefficient C 3 , while the expressions for the higher-order coefficients are deferred to Appendix C.
To O Φ|Φ| 2 , we obtain from Eq.(5) by projecting it on e iξ W 0 (y) where Let us now determine the unknown functions u (2) 0 and u (2) 2 . Once again, we substitute Eq.(3) into Eq.(1), and separate the terms proportional to e 2iξ . To lowest order in Φ, it yieldŝ The time-derivative can be evaluated self-consistently with the help of the amplitude equation (11), and to O Φ 2 is equal to Therefore, u 0 (y) satisfies the following inhomogeneous ODÊ Similarly, u 2 (y) is given bŷ The equations for higher C's and u'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 C'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 Φ(t) = Ψ e i Ω t , 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 Ψ 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 Φ(t) is sufficiently small, where the scale is given by the coefficients Cn. To study convergence of series Eqs. (18) and (19), we employ a somewhat intuitive method based on the partial sums Sm, defined through We then solve a series of algebraic equations S 1 = 0, S 2 = 0, . . . , and obtain a corresponding series of solutions Ψ 1 , Ψ 2 , . . . , and Ω 1 , Ω 2 , . . . , 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 C 11 .

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 −∆P (plane Poiseuille flow). As in Section 3, we select a coordinate system with (x, y, z) being along the streamwise, vertical (gradient) and spanwise directions, respectively; the distance between the plates is 2d.
The velocity of the fluid v is governed by the Stokes equation and the incompressibility condition where p is the pressure, τ is the polymeric contribution to the total stress, and ηs 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 [20,14], for example). In both cases, the corresponding values of the Reynolds number (defined as Re = U d/ν, where U is the typical flow velocity and ν is the kinematic viscosity of the fluid) do not exceed 10 −2 − 10 −1 , and the inertial effects can be ignored.
To describe the dynamics of the polymer stress tensor τ , we employ the Oldroyd-B model given by where λ M is the Maxwell relaxation time of the fluid, ηp is the polymer viscosity, and T 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 [21,22,13]; for a detailed discussion of various viscoelastic equations of motion and their predictions see [23,24,25]. 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.
Equations (21)- (24) have as laminar solution the well-known parabolic profile The corresponding elastic stresses are given by 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 d as a unit of length, U = ∆P d 2 /(2 (ηs + ηp)) as a unit velocity, d/U as a unit of time, and ηp U/d 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 β = ηs/(ηs + ηp), and the Weissenberg number Wi = λ M U/d. 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, and introduce a perturbation vector 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 forL, N 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.

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. [26]. To calculate the eigenspectrum, we discretise Eq.(2) withL andÂ given in Appendix A using the fully spectral Chebyshev-tau method [27], and solve numerically the resulting generalised eigenvalue problem using Scientific Python [28].
In Fig.2 we present an example of the eigenvalue spectrum, plotted as Im(λ) vs Re(λ). The general structure of the spectrum for Wi = 1, β = 0.05, kx = 1, and kz = 2 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 Re(λ) < 0, indicating linear stability. This is confirmed numerically for all values of Wi, β, kx and kz; see also Wilson et al. [26]. The most prominent features in Fig.2, the balloon-like shapes at Re(λ) = 1/Wi and Re(λ) = 1/(βWi), are numerical approximations to the continuous spectrum of the linear operator, which corresponds to the singular points of the linear problem [26]. 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 λ 1 , λ 2 , and λ 3 . (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 [18] used in our previous work on plane Couette flow [12], while the other one is mostly present in the middle of the geometry. Next, we calculate the eigenmodes of the adjoint problem, Eq.(8), withL † 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 Im(χ) = −Im(λ). The adjoint eigenmodes, corresponding to the least stable eigenvalues λ 1 , λ 2 , and λ 3 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 C's as a function of the Weissenberg number Wi, 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.  Fig. 4 Spatial profiles of the adjoint eigenmode (arbitrarily normalised) corresponding to λ 1 in Fig.2b). For visualisation purposes, the velocities and sxx are scaled by the factors 10 4 and 30, correspondingly.  Table 2 Values of the non-linear coefficients C's corresponding to λ 3 . The values of β, kx and kz are the same as in Table 1.
and study their convergence. We found that using the eigenmode associated with λ 2 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 C's for selected values of Wi are given in Tables 1 and 2. For both states, the series of solutions Ψm share the same features. Here, m denotes the order 2m+1 to which the expansion in powers of Φ in Eq. (11) is taken. For any m > 1, the equation Sm = 0 has no real solutions at small Wi, two real solutions around the saddle-node bifurcation -the value of Wi at which two solution branches appears for the first time, and one real solution for larger values of Wi. 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 Wisn, and the implications of this behaviour are discussed below.
For both states, we observe that the upper-and lower-branch values of Ψm converge rapidly as m is increased from 1 to 5. In Fig.6 we assess this convergence more quantitatively, by plotting the saddle-node values Wisn for each Ψm as a    function of m (circles in Fig.6). In the range m = 1 . . . 5, corresponding to the amplitude equation expansion up to the eleventh order, the convergence is welldescribed by an exponential fit (dashed lines in Fig.6), approaching Wi  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 profilevx for State 1, calculated by averaging vx over x and z. In Fig.8a) we present the mean profile for Wi = 3.12, which is close enough to the saddle-node point for m = 5 so that Ψ 5 has two values, 0.15 and 0.22, corresponding to the lower-and upper-branches, respectively. While the lower-branch mean velocity profile is quite close to the laminar profile U 0 (y) = 1 − y 2 , the upper-branch profile looks very different, with sharp veloc- ity 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 x = 0 in the yz-plane, which is perpendicular to the streamwise direction: the in-plane components vy and vz are shown as arrows, while the streamwise velocity vx − U 0 (y) 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 z = 0 in the xy-plane, where arrows now trace the in-plane components vx − U 0 (y) and vy, while the colour gives the spatial profile of the spanwise velocity vz. Finally, in Fig.8d) we plot the largest component of the polymeric stress tensor, τxx, which is a deviation from the laminar stress, in the xy-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 kx and kz. 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 kx and kz are not too big (typically, smaller than kx ∼ kz ∼ 3. 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.

Discussion and Conclusion
Results presented in this work further corroborate our previous claim [10,8,12,13] 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. [9] and in particular the more recent detailed and systematic experiments by Arratia et al. [14,15] 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 [29,20,30].
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 [4,5]. 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 [12] might form the phase space scaffold of the purely elastic turbulence, and verifying their existence should be paramount to understanding its mechanism. The amplitudeequation 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, Wisn ∼ 3, that are somewhat lower than the onset of turbulence values reported in experiments in channel, Wi onset ∼ 5 [14,15], and pipe flows, Wi onset ∼ 4 [8,9]. This is not surprising since the fluids used in those works had significantly higher values of β than 0.05 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 Wisn than presented here. Also, the dynamical systems scenario of the transition, if applicable in the vis-coelastic case, would imply that Wisn should be smaller than the value of Wi onset at which sustained turbulence can be observed, providing yet another explanation for the difference. 2 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 [8,9,14,15], 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 [31], 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ć [32,33]. 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 z-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. [15].
Although we are confident that we obtained converged non-linear states, their existence needs to be verified numerically, by searching for steady-state, travellingwaves, 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 socalled High-Weissenberg Number Problem [34]). 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 [35,36], 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.
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 Obviously , N (A, B) = N (B, A).

C Higher-order coefficients for the amplitude equation
In Section 3 we demonstrated how the coefficient C 3 can be determined from Eq.(5). Here, we list the expressions for higher-order coefficients C's and the associated functions u (m) n (y) from Eq.(4), obtained in the manner discussed in Section 3.