{'fraction_non_alphanumeric': 0.03828265498067263, 'fraction_numerical': 0.016454100088885214, 'mean_word_length': 4.340324539132355, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Linear Stability Analysis for Plane-Poiseuille Flow of an Elastoviscoplastic fluid with internal microstructure

We study the linear stability of Plane Poiseuille flow of an elastoviscoplastic fluid using a revised version of the model proposed by Putz and Burghelea (Rheol. Acta (2009)48:673-689). The evolution of the microstructure upon a gradual increase of the external forcing is governed by a structural variable (the concentration of solid material elements) which decays smoothly from unity to zero as the stresses are gradually increased beyond the yield point. Stability results are in close conformity with the ones of a pseudo-plastic fluid. Destabilizing effects are related to the presence of an intermediate transition zone where elastic solid elements coexist with fluid elements. This region brings an elastic contribution which does modify the stability of the flow.


Introduction
During the past several decades, yield stress fluids found an increasing number of practical applications for several major industries (which include cosmetics, foods, oil field etc.) and they are encountered in the daily life in various forms such as hair gels, food pastes, cement, mud and more. Such materials can be loosely and simplistically defined as materials that do not flow unless a minimal stress (referred as 'yield stress') is applied onto them. Understanding and controlling the hydrodynamic stability during flows of yield stress fluids is important for practical applications which involve flows of such materials through conduits. From a fundamental standpoint, yield stress materials continue triggering intensive debates and posing difficult challenges, of both theoretical and experimental nature. Undoubtedly, the best known and most cited debate is related to the very definition of such materials and the existence of a 'true yield stress', [1,2]. It was argued in Refs. [1,2] that the yield stress emerges as an artefact related to the inability of the rheometric equipment to properly identify a viscous flow regime in a range of negligibly small rates of deformation. Recent studies claim to have solved the "yield stress debate" by arguing that prior to yielding the viscosity of the material is infinite which demonstrates the existence of a true solid state, [4,23].

Yielding of a Carbopol R gel and its relevance to hydrodynamic studies
The existence of a "true" yield stress and the nature of the transition from a solid to fluid behaviour may look at a first glance of little or no relevance to the hydrodynamic stability of the flow of a yield stress material. Indeed, whereas the yielding transition occurs at low values of the Reynolds number (typically Re < 1) a loss of the hydrodynamic stability is typically observed at significantly larger Re (typically Re > 1000). On the other hand, the base flows usually considered in the linear analysis of the hydrodynamic stability of channel flows of yield stress fluids are characterized by a significant stratification of the velocity gradients: large values near the channel boundaries which are consistent with a yielded flow region and vanishing values near the center-line, which are consistent with a plug region. A recent experimental investigation of the laminar-turbulent transition in the pipe flow of a yield stress fluids demonstrates that the transition to turbulence occurs when the Reynolds stresses balance the yield stress of the fluid, that is when the solid plug is broken, [14]. These findings corroborate well with the idea that, in fact, the nature of the solid-fluid transition and the yielding scenario may play an important role in the stability problem. A rather new and certainly unexpected insight into the yielding transition has been brought by an experimental study on sedimentation of spherical particles in a physical gel (Carbopol R 940), [35]. Quite unexpectedly, the fore-aft symmetry of the flow patterns around the falling object was broken (in spite of the smallness of the Reynolds number during these experiments) and even more surprisingly a negative wake has been observed (see Fig. 8 in Ref. [35]). A similar non symmetric flow pattern has been observed experimentally during flows of Carbopol R gel past a cylinder, [41]. A first message of the sedimentation study presented in [35] was that the solid-fluid transition in the Carbopol R gel might not be reversible upon increasing/decreasing stresses, unlike it has been traditionally assessed. A second important message was that the elastic effects responsible for the emergence of the experimentally observed negative wake might be consistent with an intermediate transition regime where solid material elements (unyielded) coexist with fluid (yielded) ones. A more systematic investigation of these features of the solid-fluid transition in a Carbopol R gel has been presented in Ref. [34]. Indeed, as suggested by the sedimentation experiments presented in [35], a hysteresis has been found in the solid-fluid transition of a Carbopol R gel upon increasing/decreasing stresses (see Fig. 3 in Ref. [34]). A clear signature of the elasticity involved in the transition region in the form of a recoil effect (negative shear rates) on the decreasing stress branch of the flow curve has been observed as well. The hysteresis effect observed in the increasing/decreasing stress flow curves and reported in [34] has also been observed experimentally by Divoux et al., [9], and attributed to critical slowing down of the system dynamics close to yielding (and not to any memory effect).

Linear stability analysis of a yield stress fluid: a brief survey
The first study on linear stability analysis of a Bingham viscoplatic fluid is due to Frigaard et al. [10]. These authors concentrated on the asymptotic stability of two-dimensional traveling wave perturbations and found that the critical Reynolds number increases linearly with Bingham number. They also showed that the plug region remains unaffected by the disturbance field. In this work the authors imposed even symmetry at the yield surface and the results are incomplete. Some years later Nouar et. al. studied again this problem using the exact boundary conditions at the yield surface and they found that the plane-Bingham-Poiseuille flow is linearly stable. This being a consequence of the vanishing perturbation at the yield surface. They also conjectured that this result could be extended to Poiseuille flow of a viscoplastic fluid in a circular or annular pipe and to other rheological models such as the Herschel-Bulkley and Casson models. As mentioned before, recent studies have found that the solid-fluid transition in some yield stress materials, such as Carbopol R , is not reversible under increasing/decreasing stresses, thus the objective of this work is the linear stability analysis of an elastoviscoplastic fluid. In this work we consider a modified version of the structural model developed by Putz and Burghelea in [34]. This models consists of a non-linear Maxwell-type viscoelastic constitutive equation and a kinematic equation that governs the behaviour of the microstructure. The relaxation time of the fluid depends on the structural variable which is equal to one if the fluid is unyielded and zero if the fluid is fully yielded. We should note that the transition from viscoelastic solid to fluid is continuous and smooth. For the non-linear viscosity we consider a regularized viscoplastic model, its behaviour is the one of a pseudo-plastic fluid, that is, if the second invariant of the stress tensor is below the yield stress the fluid presents a very high viscosity.
Several authors have studied the linear stability problem of a flow with variations in viscosity across the channel [6,25,32]. Govindarajan et. al. [13] showed that by having a viscosity which is decreasing function of space it is sufficient to considerably delay the onset of two-dimensional instabilities. This result was latter confirmed by Nouar and Bottaro for the case of shear-thinning fluids [25]. Linear stability analysis for a regularized Bingham model were carried out by Frigaard and Nouar in [11]. They showed that the spectrum for the regularized problem had some physical spurious eigenvalues, not present in the Bingham problem. The stability of these eigenvalues depended on the size of the regularization parameter and unstable modes can be found even for moderate Reynolds numbers if this parameter is small enough.
For the case of Maxwell liquids, Gorodstov and Leonov [12] showed the existence of a stable continuous spectrum and two stable eigenvalues for each value of the wave number in the streamwise direction, α. In the same work they predicted instabilities at finite Reynolds numbers. Later, this result was proven wrong by several authors, e.g. [16,20,36]. Through careful numerical simulation of the generalised eigenvalue problem no unstable modes were found, even for high Reynolds numbers. Denn and coworkers [16,33,38] also studied the stability of plane Poiseuille flow and found that at high Reynolds numbers UCM fluids are significantly less stable than Newtonian fluids. Sureshkumar & Beris [40] confirmed this and also showed that destabilisation is reduced for the Oldroyd-B and Chilcott-Rallison models. Renardy [37] has proven the linear stability of Couette flow of a UCM fluid, but there is no proof of linear stability of Couette flows for more general fluids of this type. However, it is generally believed that plane Couette flow is linearly stable. Indeed much of the literature is focused at the study of interesting features of the eigenspectra, rather than marginal stability, e.g. [19,42].

Elastoviscoplastic model with internal microstructure
The elastovisplastic fluid is described by the following set of equations with unknowns (p, u, τ , Φ): where µ s denotes the solvent viscosity. The rate of strain tensor isγ =γ(u). The constitutive equation for the elastic stress tensor τ is: is the Upper-Convected-Derivative and D · /Dt is the usual material derivative. Finally, the concentration of the solid state, Φ, satisfies the following kinematic equation: with R d being the rate of destruction of solid units and R r is the rate of fluid recombination of fluid elements into a gelled structure. Note that when Φ ≡ 1 the model describes a Maxwell type viscoelastic fluid with a nonlinear viscosity µ(γ) and if Φ ≡ 0 the model represents a Generalized Newtonian fluid.
In both (3) and (4),γ and τ are the second invariant of the rate of strain and elastic stress tensors respectively. These are defined as follow: 2.1. Choice of functions for nonlinear viscosity µ (γ), destruction rate R d (Φ, τ (u)) and recombination rate R r (Φ, τ (u)).
In this section we present our choices for the nonlinear viscosity µ(γ) and the destruction and recombination rate functions. The choices are clearly non unique. Constitutive models for fluids with a yield stress trace back to the early 1900's . The most widely used are the ones due to Bingham [3], Herschel-Bulkley [15] and Casson [5]. Description of the rates of destruction and recombination is rather empirical and for this we follow closely [34].

Non-linear viscosity µ(γ)
The fluids we are interested in present a yield stress. It is well known that an "ideal" viscoplastic fluid behaves as a Newtonian (or Generalized Newtonian) fluid if τ is above the yield stress, if the contrary happens then the material behaves as a plastic solid and one says that the fluid is unyielded. In [34] it has been shown that the transition from an elastic solid behavior (unyielded) to a fluid behavior (yielded) is not direct but mediated by a viscoelastic like regime where solid and fluid behavior coexist. This effect, though previously unnoticed, is particularly pronounced in the case of "fast" yielding, i.e when the externally applied stress changes unsteadily and becomes vanishingly small in the limit of steady yielding, in other words, when the stresses vary infinitesimally slow. Although it has been argued that if one waits "long enough" the viscosity in the unyielded region actually tends to infinity, see [23], in an experimental (to investigate hydrodynamic stability for example) or industrial setting it is unlikely one will reach those time frames. Therefore for the rest of this study we will work with a regularized version of a viscoplastic constitutive model. Because our main interest is on linear stability analysis and this involves high Reynolds numbers it is better if we have a model with non-zero infinite shear viscosity, for this reason we choose the shear-thinning regularized Casson's model for the viscosity: where m ≥ 1 is the power law index and τ y is the yield stress and ǫ ≪ 1 is the regularization parameter . We recover Casson's model when m = 2, Bingham model when m = 1 and the model reduces to a Newtonian viscosity when m = 1 and τ y = 0. In Figure 1a show an example of (5), note that if τ < τ y the viscosity becomes very large but never infinity as long as ǫ is fixed and different from zero, thus we will be working with a pseudo-plastic viscosity (shear-thinnig).
2.1.2. R d (Φ, τ (u)) and R r (Φ, τ (u)) As discussed by Putz and Burghelea in [34] R d and R r have to satisfy the following assumptions: 1. R d (Φ, τ (u)) is proportional to the relative speed of neighboring solid units and the existing amount of solid, thus and where g(τ ) is an increasing function of the second invariant of the stress tensor, K d is the rate of destruction with units s −1 and w determines the "width" of the solid-fluid region and has units Pa.
2. R r (Φ, τ (u)) decreases with the relative speed of neighboring elements, begin practically zero in a fast enough flow. Therefore, and where f (τ ) is a decreasing function of the second invariant of the stress tensor. In Figure  1b we show an example of (7) and (9). Note that when τ is far below the yield stress, the material is fully "unyielded" which means it behaves as a viscoelastic solid. This because when τ < τ y , µ(γ) is of the order O(1/ǫ) and Φ ≡ 1, thus by a simple dominant balance analysis equation (3) reduces to the one for a Kelvin-Voight viscoelastic solid. If the opposite happens (τ far above the yield stress) Φ ≡ 0 and the fluid behaves as a shear-thinning fluid. It is important to note that in a neighborhood around τ y the material is neither a viscoelastic solid nor a viscous fluid, both solid and fluid structures coexist and this region will play a crucial role int the stability of the flow. The existence of this region is the main difference with respect to other proposed models for this type of elastoviscoplastic fluids. The main difference of the approach used here and previous work, see [22], is that we do not use the steady states of (4) but solve it numerically and thus describe a non-steady yielding behavior which we believe might be of some practical relevance (a)  (b) f (τ ) and g(τ ) with K d = K r = 0.3 and w = 0.1.

Validation of the model against experimental results
To validate the model, we use the rheological measurements presented in [34] conducted for various Carbopol R solutions with concentrations (by weight) ranging in between 0.05% and 0.2%. The measurements were conducted on a Bohlin (now Malvern) C−V OR rotational rheometer equipped with a digitally controlled temperature (within 0.1 • C ) bath (RT E − 111 from Neslab). The rheometer was operated in a controlled stress mode. To insure a good reproducibility of the measurements, prior to each experiment the sample has been pre-sheared at a constant stress (usually the largest stress applied during the test, so that yielding occurs in the entire volume of the material) for 300 s. After the pre-shear step, the samples were allowed to relax for 300 s prior to the rheological measurements. A first major concern during these experiments was the occurrence of the wall slip phenomenon at the contact with the measuring geometry, which is a well known and documented effect for Carbopol R gels, [31]. To prevent this, a serrated parallel plates geometry with an approximate roughness of 0.8 mm has been used. A second concern was related to the possible artefacts introduced by fluid evaporation during long experimental runs. To prevent this, a solvent trap has been placed around the free fluid meniscus. After each experimental run we have carefully checked (by visual inspection) that no significant changes in the shape of the meniscus occurred, and thus concluded that evaporation effects were either minimal or absent. The radius of the plates was R = 2 cm and the distance between them was d = 1 mm (measured between the plates protuberances). Two types of rheological tests have been performed: controlled stress linear ramps and controlled stress oscillatory sweeps at a fixed frequency. For each sample we have conducted controlled stress experiments for 19 different values of the total ramp time. Each constant stress rheological experiment started with a increasing stress ramp and ended with a decreasing stress ramp within the same range of stresses and the same stress step. The data averaging time per stress value, t 0 (referred in [34] as the characteristic time of forcing), has been varied between 0.2 s and 2 s. For each up-down stress ramps 1000 stress values have been explored ranging in between 0.1 P a and 20 P a. We emphasize here that a true steady state of deformation can only be inferred in the asymptotic limit t 0 → ∞, unlike in previous studies concerning Carbopol R gels where a steady state was a priori set by deliberately choosing large values of t 0 (an accurate description of such procedure is presented, for example, in Ref. [7]). Alternatively, the time dependent response of the samples has been tested via stress controlled oscillatory experiments at several frequencies and in the same range of stresses. We now turn our attention to the comparison of the predictions of the model described in Sec. 2 with the experimental results presented in Ref. [34].
Bearing in mind that the experiments provide a scalar data set, the equations (1)-(4) reduce to: Because the rheometer operates in a controlled stress mode, the above equations decouple. First we solve for (10) using the function ode15s provided with the MATLAB R distribution. The function ode15s performs well with stiff problems and its accuracy is recognized as medium to high. Equation (11) reduces to an algebraic equation forγ which we solve using the function fzero provided with the MATLAB R distribution as well. The function fzero uses a combination of bisection, secant, and inverse quadratic interpolation methods.
In Figure 2a we present the steady state solution of equation (10) which reads In Figure 2b we present the solution of (10) using the same stress ramp as in the experiments presented in Ref. [34]. Clearly, for intermediate values of τ , the decreasing stress branch of the solid concentration lags behind the increasing branch and a hysteresis is clearly visible. Having now both τ and Φ we can calculateγ from (11) using (5) for µ(γ).
The values for m, µ ∞ and ǫ are obtained from the fit performed in the rheometer. We obtain the elastic modulus G from a linear fit of the experimentally measured dependence τ = Gγ using the experimental data. The parameters K d , K r and w were chosen such that the numerical flow curve is in close agreement with the experimental flow curve.
In Figure 3a we compare the model discussed above with the experimental data. An excellent agreement is obtained, particularly at the point whereγ = 0. Finally, we test the model against to oscillatory flow measurements where a harmonic forcing τ = τ 0 sin(2πκt) and the strain response γ = γ(t) is monitored. This step is needed for a proper validation of the model, as the model involves three adjustable parameters (K d , K r and w). Note thatγ = dγ/dt. As before, having τ we solve for Φ using (10), thenγ is found using (11), and finally we can use standard numerical integration to find γ. The results are presented in Figure 3b.
One can conclude that using the values of the adjustable parameters that have validated the flow curves presented in Fig. 3a a good agreement with the oscillatory measurements is obtained without any additional parametric adjustments, Fig. 3b. We consider this result as a validation of the model initially proposed in [34] and revisited in Sec. 2 which further justifies employing it in the linear stability analysis presented through the rest of the paper.

Non-dimensionalization
To conclude this section we present the non-dimensional set of equations. Let where µ r is the ratio of the solvent viscosity to the total viscosity, i.e. µ s /µ. Finally: where The Reynolds (Re), Weissenberg (W e) and Bingham (B) numbers are defined: Through the rest of the paper only non-dimensional variables will be used and the "hat" shall be dropped for simplicity.

Linear stability analysis for shear-thinning regularized Casson's Model
First we consider the linear stability of the regularized Casson model. We note that our main goal is to investigate the effect the solid and solid-fluid regions have on the stability of the flow. It is known that a viscoplastic fluid is stable to all infinitesimal perturbations (at least no instabilities have been found numerically) see [27], therefore a natural second step towards our goal is to ask ourselves how the presence of a highly stratified viscosity affects the stability of the flow, thus we consider the case Φ ≡ 0. Linear stability analysis for a regularized Bingham model was carried out by Frigaard and Nouar in [11]. There they showed that as ǫ → 0 the velocity field using a regularized viscosity converges to the velocity field using the non-regularized viscoplastic viscosity. While considering the problem of linear stability they showed that the spectrum for the regularized problem had some physical spurious eigenvalues, not present in the Bingham problem. The stability of the eigenvalues depended on ǫ, we will discuss further about this in Sec. 3.4.3. Here we are not interested in the distinguished limit ǫ → 0. Based on experimental data we are aware of that typically O(10 −3 ) ≤ ǫ ≤ O(10 −2 ), therefore for the rest of the paper we fix ǫ = 0.002 unless otherwise stated.
For the rest of this section we set Φ = 0, µ s = 0 and µ = µ ∞ , then the equation of motion reads along with the incompressibility condition ∇ · u = 0.

Baseflow profile
We consider shear flows of the form U = (U (y), 0, 0) with y ∈ [−1, 1], then (13) reduces to: dτ xy dy = Re dp dx (21) with the usual non-slip boundary conditions U (−1) = U (1) = 0 at the walls. It is clear that in this situationγ = |dU/dy| and because U (y) reaches its maximum at the center of the channel we have:γ = −dU/dy if y ∈ [0, 1] andγ = dU/dy otherwise. Having this in mind and integrating once (21) we have the following algebraic equation forγ: We solve this non-linear equation forγ using fzero in MATLAB for each point of our discrete domain. Once we have the approximation forγ we integrate it with respect to y and use our boundary conditions to get U (y). We feed this result to an user-built function in MATLAB and use fzero again to find Re(dp/dx) such that max y (U (y)) = 1.
We present different examples of the base flow for different Bingham numbers in Figure  4, clearly as we increase B we can see the existence of a pseudo-plug region around the center line of the channel. Recall that this is not a true plug, i.e. the velocity there is not constant.

Linearised momentum equation and tangent viscosity
As is common in linear stability analysis we consider an infinitesimal perturbation (εu ′ , εp ′ ) superimposed upon the base flow and linearize the momentum equation (20) around (U , p) to get: where τ ′ is the stress perturbation given by: whereμ is the viscosity perturbation and it is given by: Given the fact that the flow we consider is unidirectional, it can be shown that (see [25]): where is the tangent viscosity. For a one-dimensional shear flow, the tangent viscosity is defined by µ t = dτ xy /dγ xy whereas the effective viscosity is defined as µ = τ xy /γ xy . In Figure 5 we show some examples of the non-linear viscosity µ(y) and µ t (y) for increasing Bingham numbers, note that the functions are smooth, µ ≥ µ t and that max y µ(y) ≡ max y µ t (y) ∼ B/ǫ. For a detailed description of the tangent viscosity concept we refer the reader to [25].

Modified Orr-Sommerfeld and Squire equations
Writing equation (23) in terms of the normal velocity v ′ (x, t) and the normal vorticity η(x, t) = ∂u ′ /∂z − ∂w ′ /∂x and assuming modal solutions of the form: we get the following eigenvalue problem for the frequency ω: where with k 2 = α 2 + β 2 , D = d/dy,∆ = D 2 − k 2 . Together with the boundary conditions v = Dv = η = 0 at y = ±1.

Bounds for the Squire and Orr-Sommerfeld modes and one-dimensional stability
In this section we present bounds for the eigenvalues of the modified Orr-Sommerfeld and the modified Squire operators, equations (29) and (32). We follow closely the work by Joseph [18] and Davies and Reid [8]. For the rest of the section instead of working with the frequency ω as the eigenvalue we consider the wave speed, c = ω/α unless otherwise stated.
As in [18] we define and take v ∈ H, where H is a complex-valued Hilbert space completed under the norm I 2 2 by the addition of limit points of sequences of functions in C 4 ([−1, 1]) satisfying (33).
The following isoperimetric inequalities hold in H: where λ 3 is the smallest eigenvalue of a vibrating rod with displacement v satisfying (33).
See A.1 for the proof.

3.
U ′′ min ≤ 0 : See A.2 for a proof. We should note that these results are not restrictive to regularized viscoplastic fluids, they hold for any shear-thinning (we are using the fact that µ is a decreasing function oḟ γ) viscosity function.

One-dimensional stability
From Theorem 3.2 we have that the Squire modes are always damped, thus for the rest of this section we consider only the modified Orr-Sommerfeld equation and α = 0. If instead of working with the perturbation of the normal velocity v we consider a stream function formulation, i.e. u = Ψ y and v = −Ψ x , and taking a modal solution of the form Ψ = φ(y) exp[i(αx − ωt)], then we can interchange v for φ in (28). Note that φ ∈ H and Theorem 3.3 holds for this formulation. Therefore using (89) we have a bound for the imaginary part of the frequency ω: Here we are interested in one-dimensional perturbations, i.e. α = 0. This is a special case in terms of boundary conditions. The fact that v ≡ 0 does not mean that φ = 0 at y = ±1. This boundary condition can be obtained if the perturbation to the flow rate 1 −1 u(y)dy is taken to be zero and a uniform pressure gradient in the x-direction is allowed, see [36]. Now, even for the case α = 0 we have φ ∈ H and (39) becomes: where we used the isoperimetric inequalities defined in Lemma 3.1. Therefore, onedimensional perturbations are always linearly stable. Note again that this result is valid for any decreasing function µ(γ).

Solution of modified Orr-Sommerfeld equation and code validation
There is no equivalent of Squire's theorem for generalized Newtonian fluids but some authors have performed several numerical tests for large range of axial and transverse wave numbers (α and β respectively) for different non-linear viscosity functions µ(γ), see [24][25][26][27]. Their numerical results show that the lowest critical Reynolds number is obtained for spanwise homogeneous perturbations (β = 0, α = 0). We should also note that in the situation of a homogeneous streamwise perturbation (α = 0, β = 0), the imaginary part of ω satisfies: which comes from letting α = 0, multiplying by the complex conjugate of v, say v * , integrating by parts and taking imaginary part in (28). Note that if α = 0 then L and S decouple, thus the flow is unconditionally stable. Using these facts we would only consider 2D perturbations in the streamwise direction for the rest of the section. We know from Theorem 3.2 that all the Squire modes are always damped, thus we are concerned only with the modified Orr-Sommerfeld equation: Together with the boundary conditions v = Dv = 0 at y = ±1.
We solve the eigenvalue problem (42)-(43) by discretizing using Chebyshev polynomials in the usual fashion, as described for example in [39]. In order to validate our code, we benchmark by solving the Newtonian problem m = 1 and B = 0, i.e. µ ≡ µ t ≡ 1 in (42) and compare with results obtained by Mack [21]. In Figure 6a we present the first 33 modes calculated using our code with N = 150 and the ones reported by Mack with the eigenvalues of the two computations overlaying one another. Figure 6b shows the error norm between our calculations and Mack's for N = 50 and N = 150. As expected the least stable eigenvalue has converged but the split of the S-branch is present when N is not large enough. In order to stably compute the spectrum of (42)-(43) we need substantially more nodes than for the Newtonian problem. In Figure 7 we present the calculations using N = 150, 200, 250. The first thing to note is that even with N = 150 the S-branch of the spectrum has not converged. Another important fact about the spectrum is the presence of an R-branch, which was first documented by Frigaard and Nouar in [11]. They found that this branch is physically spurious, meaning that these eigenvalues do not exist for the Orr-Sommerfeld operator corresponding to the viscoplastic model and their stability depends on ǫ. As ǫ → 0 we can find unstable modes from this R-branch even for moderate Reynolds numbers. We have found that by choosing a smooth regularized non-linear viscosity (either the one due to Becovier et. al. or the one due to Papanastasiu) increases the range of ǫ for which the R-branch remains stable and the least stable eigenvalue belongs to the A-branch as expected. We have fixed ǫ = 0.002 for all our results and we have checked that the critical eigenvalue is indeed in the A-branch. For the rest of the section we fix the number of nodes to be N = 250.

Numerical results
Stability results for the regularized Casson model are presented in the following section. First we introduce the "plasticity number " defined as: Note that P l only depends on the rheological properties of the fluid and the geometry of the problem, thus as we increase Re in our analysis P l will remained fixed. All the results that follow show the critical Reynolds number for the regularized model normalized with respect to Re N ewt = 5772.2. In Figure 8a we fix m = 2.5 and with the increase P l we enhance the stability. These results are not surprising, we are mimicking a viscoplastic fluid which is known, at least numerically, to be linearly stable [27] with a pseudoplastic fluid with a very large zero-shear rate viscosity. The use of this type of non-linear viscosity has also been shown to enhance stability [25] which is also represented in Figure  9a where we fix P l = 10000 and vary the power-law index m. Figures 8b and 9b show the corresponding critical wave number for increasing plasticity number and increasing power law index respectively. Note that these values remain in close proximity to α = 1.02 which is the critical wave number for the Newtonian case.

Linear stability analysis for Elastoviscoplastic model
In this section we study the linear stability for the elastoviscoplatic model presented in equations (13)-(16).

Baseflow profile
where T = 1/2 T 2 xx + 2T 2 xy . Substituting (46) into (44) we get: This equation is solved in exactly the same way as (21) in Section 3.1. Once we havė γ = |U ′ | from (49) we solve the following nonlinear equation for T xx : where we have used the definition of T xy in (46). We now solve equation (50) using fzero in MATLAB R for each point y in our discrete domain. Finally we can construct Φ using (48). In Figure 10 we present the solution of equations (44)-(48) for increasing Bingham numbers and fixed Weissenberg number W e = 1. As we can see the velocity profile presents the pseudo plastic behaviour which in this case, in addition of a intensely stratified viscosity, we have a viscoelastic plug, due to the presence of non-zero normal stresses inside this region. We should note the existence of a thin region (O(w)) where both solid (viscoelastic) and fluid phases coexist as it's clearly seen in the plot of Φ s .

Orr-Sommerfeld equation for Elastoviscoplastic model
As before, we consider an infinitesimal perturbation (εu ′ , εp ′ , ετ , εΦ ′ ) superimposed onto the base flow and the field equations are linearized around (U , P, T , Φ), thus we have: As there is no version of Squire Theorem for this type of models for simplicity we restrict ourselves to two-dimensional perturbations below. We introduce a stream function Ψ, such that: u ′ = Ψ y , v ′ = −Ψ x . We consider modal linear disturbances of the form: −ωt) , ij = xx, xy, yx, yy, Φ ′ = ϕ(y)e i(αx−ωt) , .
(55) Denoting D = d dy , the eigenvalue problem for ω is: where and with boundary conditions: As it is customarily done in linear stability analysis for viscoelastic flows we introduce the elasticity number defined as: which is the ratio of kinematic viscosity and relaxational diffusivity: ǫℓ depends only on the properties of the fluid and the flow geometry. For the rest of the paper we replace W e by ǫℓRe.

One-dimensional stability
In this section we consider one-dimensional disturbances of (56)-(60). In order for φ to satisfy boundary conditions (62) we proceed in the same way as for the regularized viscoplastic fluid, by letting the perturbation to the flow rate 1 −1 u(y)dy to be zero and we impose a uniform pressure gradient in the x-direction. Letting α = 0 in (56)-(60) we are left with the following eigenvalue problem: (1 − iǫℓReΦµω) τ xx = 2ǫℓReΦµDU τ xy + 2ǫℓReΦµT xy D 2 φ + 2ΦµDU T xy ϕ, (64) We should point out that from equations (64)-(67) a continuous spectrum exists and this consists of purely imaginary eigenvalues in the strip [−1/ǫℓRe max(µ), −∞) (note that for our choice of regularization parameter ǫ, we always have (K d + K r ) > 1/ǫℓRe max(µ)). One has to be careful when considering this part of the spectrum, note that as Φ → 0 then ω i → −∞ , clearly these eigenvalues are always damped, thus we consider only the case τ yy = 0. Also note that (64) and(67) decouple from (63) and (65), then we multiply (63) by the conjugate of φ, say φ * , (65) by τ * xy and integrate by parts to get: Expanding in real and imaginary parts and applying the mean value theorem for integrals as necessary we have the following equation for the imaginary part of ω, Therefore, (56)-(60) is linearly stable to one-dimensional perturbations.

Code validation
As a first test for our code we consider the case µ r = 0.5, Φ = µ = µ t ≡ 1 and ϕ ≡ 0 in (56)-(60), by doing this we recover the Oldroyd-B model and compare our results with the work by Sureshkumar and Beris in [40]. We set Re = 3960, W e = 3.96 (or ǫℓ = 10 −3 ), α = 1.15 and have used N = 128 and N = 256 to solve for the spectrum. We show the results in Figure 11, the least stable eigenvalue is ω = 0.3409 + 1.9888 × 10 −7 i, the value reported in [40] is ω = 0.3409 + 1.9696 × 10 −7 i. We have convergence of the discrete spectrum and the difference of the least stable eigenvalue between computations is 2.7 × 10 −12 .
We now turn our attention to the full model (56)-(60). For the rest of the section we consider the wave speed c = ω/α as our eigenvalue. The first thing to note is the location of the continuous spectrum, from (57)-(60). For simplicity let us define the Deborah function as De(y) = ǫℓΦµ. We can see that the real part of the continuous spectrum is in the line c r = U and the imaginary is in the line c i ∈ [a, ∞), where: Clearly, the continuous spectrum is always stable but one has to be careful when doing computations. Note that for y in the pseudo-plug region (where T < B), De(y) is very large (recall that µ ∼ B/ǫ), once we begin to approach the yield limit De(y) decays to zero very fast as Φ → 0 and µ → 1. This means that the continuous spectrum will approach zero when y is on the pseudo-plug region and it will be pulled down very fast towards minus infinity as T approaches B. This could cause severe numerical problems, Φ is actually zero when the fluid is fully yielded. We discretize the problem using the same technique as in Section 3.4.3. In order to prove the convergence of our code we compare the spectrum of a regularized version of (56)-(60) with the one of the full model using two different eigenvalue solvers in MATLAB R . Let where ϑ ≪ 1 is a small parameter. We replace Φ with Φ r in De(y) and solve the generalized eigenvalue problem using the function eig in MATLAB R , we choose ϑ such that we get finite results. We then solve the full problem (56)-(60) using eigs in MATLAB R , this function provides the reverse communication required by the Fortran library ARPACK which is based upon an algorithmic variant of the Arnoldi process called the Implicitly Restarted Arnoldi Method. This function can calculate k eigenvalues based on a user defined σ, thus we choose k such that we get finite results and σ is taken to be the least stable eigenvalue of the regularized Casson model calculated in Section 3.4.3. We present the results of these calculations with N = 250 in Figure 12a, the spectra of both problems overlap each other. Note the presence of the almost undisturbed spectrum of the generalized Newtonian problem (A-branch, P -branch, S-branch and R-branch), the continuous spectrum corresponding to equation (60) centered at c i = (K d + K r )/α 2 (we have assumed K d = K r for this example). Because the choice of parameters (see caption in figure), the pseudo-plug region is very small we have only a few eigenvalues corresponding to c r = U ∼ 1 and c i = −1/α 2 De(y)Re. Note that there exists a set of discrete eigenvalues which emanates from the region close to −1/α 2 De(0)Re which distorts the A, P, R branches of the modified Orr-Sommerfeld operator. The difference on the least stable eigenvalue between the problem using Φ r and the one using Φ is 5.1101 × 10 −6 . This result is not surprising due to the fact that this eigenvalue belongs to the wall modes (A-branch), there the fluid is fully yielded and it does not feel the effect of elasticity, when we use Φ r close to the wall we can assume that ǫℓ ∼ O(10 −8 ) for this choice of parameters. In view of these results we use the function eigs for the full discretized problem (56)-(60). In Figure 12b we present convergence results for N = 250, 350, 450, the change on the least stable eigenvalue is on the sixth significant figure. For the rest of the paper we fix N = 450 unless otherwise stated.

Numerical results
Next we present the results of our cacluations. Figure 13a shows the critical Reynolds number for the elastoviscoplastic model as a function of plasticity number, P l. As with the regularized Casson model the existence of a pseudo-plug region (stratified viscosity) is suffciente to greatly enhance the stability of the flow. The critical Reynolds number appears to be a monotone increasing function of plasticity number, just as with the regularized viscoplastic model. Clearly, in Figure 13b we can see that the behaviour of the critical wave number α c is very similar to the wave number of the Casson model. We should also note that the inclusion of a highly viscous viscoelastic fluid as a plug destabilises the flow when comparing it with the regularized model (shown as the dotted curve in Figure 13a). In relative terms, when P l = 1000 the critical Reynolds number for the elastoviscoplastic model is 2.66% smaller than the critical Reynolds number for regularized model. For P l = 10 5 this percentage increases to around 6%, this is due to the fact that the pseudo-plug and solid-fluid regions increase and are closer to the wall. As a way to understand better the effects of the elastoviscoplasticity on the instability mechanism we will turn our attention to the analysis of the distribution of the production and dissipation of disturbance kinetic energy across the channel. The linearized Reynolds-Orr equation is derived by multiplying equation (56)  model) by the conjugate of the stream function perturbation (or velocity perturbation in (42)) and integrating over a periodic domain. This can be written as: where j =regularized (reg) or elastoviscoplastic (ev) and · = 1 −1 ·dy . The left hand side of (72) is the temporal variation of the averaged kinetic energy, the first term in the right hand side is the exchange of energy between the base flow and the perturbation and the last term in the equation is the rate of energy dissipation. Explicitly we have for the regularized model and for the elastoviscoplastic model Following Govindarajan et al. [13], we compare the normalized space-averaged energy production, and normalized space-averaged energy dissipation At criticality these two quantities balance. Positive values of I 2,j indicate where in the flow domain energy is supplied from the base flow to the perturbed flow. It should be mentioned that even though that the energy production and energy dissipation balance at critical conditions, the mechanism of instability is governed by I 2,i and not by I 3,j . We refer the reader to [13] for an extensive discussion about Reynolds stress distribution analysis. In Figure 14a we can see how Γ + reg and Γ − reg balance at critical Reynolds number Re reg = 13364.41 for P l = 1000. The inclusion of an elastic plug an its effect on the stability of the flow is clearly seen in Figure 14b. We keep the same Reynolds number as for the regularized case. Note that both Γ + ev and Γ − ev are in the same orders of magnitude as before, but the balance is slightly broken. The energy production has increased but note that the energy dissipation has gone negative close to the centre line of the channel (insert figure) whereas for the regularized model is always positive due to the presence of a highly viscous pseudo-plug. This is due to the presence of elastic forces. Even though the elastic plug is confined to a very small region away from the wall, it is enough to break the balance between energy production and dissipation and stability is lost. Experimental evidence of the existence of an elastic plug can be found in [28]. Now we explore the effects that increasing elasticity will have on the stability of the flow. In Figure 15a we present the normalized critical Reynolds number for increasing ǫℓ. Note that the change is minimal, the same happens for the critical wave number shown in Figure 15b. This is not surprising either, elastic forces are confined to regions of very low shear-rates where the value of the nonlinear viscosity is very large, thus an increase on ǫℓ does not make any difference. This is also seen in Figure 17a, where ǫℓ = 10 −3 and Γ + ev and Γ − ev are almost unchanged with respect to the case ǫℓ = 10 −4 . In Figure 16a-b we show the effects of increasing the parameter w which in turns makes the solid-fluid region wider. Note that for small values of w things remain more less unchanged but as soon as w gets closer to 1 we see a steep decrease on the critical Reynolds number. This again is not too surprising, by increasing w the region where solid-fluid coexist stretches towards the wall. When w = 1 we have that both Φ and T xx are different from zero at the wall. This is clearly seen in Figure 17b where the normalized production of energy Γ + ev has increased in the critical layer.

Discussion
We have studied linear stability analysis of plane Poiseuille flow of an elastoviscoplastic fluid. Our results show that there is an increase in the critical Reynolds number for two-dimensional perturbations as the yield stress increases. The main difference between our study and previous works is that we consider a rather new yielding scenario recently suggested in [34] according to which the transition from solid to fluid regime is not direct but mediated by a solid-fluid phase coexistence regime within which the material behaves as a viscoelastic fluid. Thus, we were interested in understanding how the presence of elasticity modifies the hydrodynamic stability of the flow. The relevance of a particular yielding scenario to the stability problem might not be obvious at a first glance as the yielding typically occurs at low Reynolds numbers (the experiments described in [34] were performed at Re<1) whereas the loss of hydrodynamic stability emerges at significantly higher Re. However, it has been demonstrated experimentally in [14] that the transition to turbulence in the pipe flow of a viscoplastic fluid (Carbopol 940) is simultaneous with the breakdown of the unyielded plug initially located around the centerline of the pipe. These results suggested us that the loss of the hydrodynamic stability and the yielding are in fact connected, which motivated us to revisit the stability problem in the context of the yielding scenario suggested in [34]. These results are in a close agreement with the ones of a regularized viscoplastic model, this is not surprising due to the fact that the viscoelastic core is confined to a region away from the wall. We have also shown that the existence of a region where solid and fluid structure coexist destabilizes the flow.
We have chosen the regularized Casson model for the viscosity function but we suspect       that the stability behaviour is qualitatively similar for other models such as the regularized Hershel-Bulkley and Bingham models, and therefore that our results are quite generic.
There are practical limitations to be acknowledged. First, the domain is a geometric idealisation of a planar geometry of large aspect ratio. Depending on the actual geometry it becomes necessary to consider other effects, e.g. entry/start-up effects of the flow, curvature or imperfections of the walls, etc. On the other hand in the experimental setting, the point at which instability is actually observed is very sensitive to control of apparatus imperfections and the level of flow perturbations. For example, in Hagen-Poiseuille flow of Newtonian fluids one typically observes transition to turbulence starting for Re 2000. However, an experimental flow loop in Manchester UK produces stable laminar flows for Re ≈ 24, 000, [17,29], and stable flows have even been reported up to Re ≈ 100, 000, [30]. This suggests that enhanced stability may be achieved experimentally, where predicted by the linear theory.  c |η| 2 dy = U |η| 2 dy + i αRe µ(|Dη| 2 − k 2 |η| 2 ) + β 2 k 2 (µ t − µ)|Dη| 2 dy. (82)

A. Bounds for Orr-Somerfeld and Squire modes
Using the second mean value theorem for integrals and taking real and imaginary parts, we get: c r I 2 0 = U (y 1 )I 2 0 , with y i ∈ [−1, 1] to be the mean values. The first bound immediately comes from the fact that U min < U (y i ) < U max . Note that for shear thinning fluids, µ(y i ) > min µ > min µ t > 0. We also have µ t (y i ) > min µ t . Thus, for p = min µ t > 0, so all the Squire modes are damped since I 2 0 is positive. Dividing through by I 2 0 and applying our isoperimetric inequality I 2 1 /I 2 0 ≥ π 2 /4 given in Lemma 3.1 we get our result as required.

A.2. Proof of Theorem 3.3
Proof. Using the fact that the Squire modes are always damped, we multiply (42) by the conjugate v * , integrate by parts, use boundary conditions, take real and imaginary parts and use the second mean value theorem for integrals to get: c r (I 2 1 + α 2 I 2 0 ) = U |v ′ | 2 + (α 2 U + U ′′ /2)|v| 2 dy, where Q = (i/2) U ′ v * ′ v dy,μ = µ(y 1 ),μ t = µ t (y 2 ) with y 1 , y 2 ∈ [−1, 1]. We can expand the last integral of the right hand side of (86) and integrate by parts. This gives We will now proceed with the proof. As before, note thatμ > min µ > min µ t , thus we divide by (I 2 1 + α 2 I 0 ), using Cauchy-Schwartz inequality on the first term of the right hand side of (88) and our definitions for p and q, we have Clearly, equations (87) and (89) are the same as equations (3) and (4a,b) in [18] respectively, with the inclusion of the positive constant p in (89). Therefore the proof follows in the same way as in [18] and our results hold.