S T O P A Y E V

I F R I G A A

Taylor vortices in yield stress fluids with static wall layers on the outer cylinder

Here, we consider the stability of this flow with respect to a finite amplitude perturbation. We focus particularly on the case where the basic flow has an unyielded fluid layer on outer cylinder. A weakly nonlinear stability analysis is developed for a wide and a narrow gap. A third-order Ginzburg Landau equation is derived and the influence of the different nonlinearities on bifurcation features is investigated in detail. The results indicate that : (i) the nonlinear inertial terms act in favor of pitchfork supercritical bifurcation and the nonlinear yield stress terms promote a subcritical bifurcation; (ii) for a range of Bingham numbers, B, the extent of which depends on the radius ratio and outer Reynolds number, the nonlinear yield stress terms are dominant and the primary bifurcation is subcritical. The amplitude analysis indicates that, in the supercritical bifurcation regime, near the threshold, when the nonlinear inertial terms are dominant, the amplitude decreases slightly with increasing B. Once the nonlinear yield stress terms start to become significant, the equilibrium amplitude increases substantially with increasing B. Similar trends are observed for Taylor vortex strength. Finally, the erosion of the static layer is analyzed. It is shown that the nonlinear yield stress terms play a significant role.

Introduction

Viscoplastic material behaves as a "solid" from kinematic point of view when the applied stress is below a threshold value τy and flows like a viscous fluid for stresses higher than τy . The solid-like behavior is associated with elasticity, whereby the continuum deforms when subjected to a given stress and there is a complete strain recovery when the forcing is removed. In general, the critical strain before yielding is small [START_REF] Coussot | Saffman-Taylor instability in yield-stress fluids[END_REF], the the elastic properties may then reasonably be neglected when the stress is below τy .

In the present work, we follow this assumption. Many materials exhibit a yield-stress [START_REF] Bird | The rheology and flow of viscoplastic materials[END_REF], like drilling mud in the oil industry, cement, paints, cosmetic and pharmaceutical preparations as well as a large variety of food products. Although, many models have been proposed to describe the behavior of such materials, e.g. the Herschel-Bulkley, Casson or Robertson-Stiff models (see [START_REF] Agwu | A critical review of drilling mud rheological models[END_REF]), the Bingham model is the most well known and simple. Furthermore, it contains all the ingredients of viscoplastic materials namely a yield stress and a nonlinear variation of the effective viscosity with the shear-rate. In this model, the material is considered to be rigid below a yield criterion described by the von Mises criterion and is nonlinearly viscous above the yield criterion. The determination of these two regions is not a trivial task, especially in two-and three-dimensional flows. A review on yield-stress fluids can be found in [START_REF] Balmforth | Yielding to stress: recent developments in viscoplastic fluid mechanics[END_REF][START_REF] Coussot | Binghams heritage[END_REF][START_REF] Bonn | Yield stress materials in soft condensed matter[END_REF]. Although commonly used as industrial fluids, there are surprisingly few published works that focus on the stability of yield-stress fluid flows.

In the present work, we consider the stability of circular Couette flow for a Bingham fluid. The first linear stability analysis was done by [START_REF] Graebel | Stability of a Stokesian fluid in Couette flow[END_REF] using a narrow gap approximation. He found that the yield-stress has a stabilizing effect. This problem was later reconsidered by [START_REF] Peng | Linear stability of Bingham fluids in spiral Couette flow[END_REF] assuming axisymmetric disturbances. The most interesting feature of the results is the nonmonotonicity of the critical inner cylinder Reynolds number Re 1c for wide gap co-rotating cylinders as the Bingham number is increased. It is the only study that we know of where a yield stress fluid is less stable than the corresponding Newtonian fluid flow. It is explained by [START_REF] Landry | Stability and instability of Taylor-Couette flows of a Bingham fluid[END_REF] that in co-rotating cylinders (at small B), the decrease of the critical Reynolds number is due to an increase of strain rate of the basic flow which amplifies the production term in the linear energy equation. The production term provides the only means by which energy is exchanged between the base flow and the perturbation. This linear analysis is completed with the transient growth characteristics of both axisymmetric and non axisymmetric perturbations [START_REF] Agbessi | Linear stability of Taylor-Couette flow of shear-thinning fluids: modal and non-modal approaches[END_REF][START_REF] Chen | Transient growth in Taylor-Couette flow of a Bingham fluid[END_REF]. It is shown in particular that the yield stress reduces strongly the transient growth.

A numerical simulation of axisymmetric Taylor-Couette flow of Bingham fluids was performed by [START_REF] Jeng | Numerical simulation of Taylor Couette flow of Bingham fluids[END_REF], in the case of a wide gap with a radius ratio η = R 1 /R 2 = 0.5. To overcome the discontinuity in Bingham model, in the transition from solid-like to liquid-like, Papanastasiou regularization [START_REF] Papanastasiou | Flows of materials with yield[END_REF]) is used, which treats the whole material domain as a fluid of variable viscosity and locally assigns a large but finite value of viscosity to the unyielded region. A review of popular regularization models has been carried out by [START_REF] Frigaard | On the usage of viscosity regularisation methods for viscoplastic fluid flow computation[END_REF].

When the outer cylinder is at rest, and for a fixed inner cylinder Reynolds number Re 1 (not for a fixed relative distance to the onset of vortices, ǫ = (Re 1 -Re 1c )/Re 1c )) [START_REF] Jeng | Numerical simulation of Taylor Couette flow of Bingham fluids[END_REF] found that the intensity of the vortex flow is weaker than that obtained for a Newtonian fluid. For the co-rotation situation with fixed outer and inner Reynolds numbers, the authors found that first the intensity of the vortex flow is initially strengthened with increasing the yield stress, and then weakened as the yield stress is raised further.

Experimental results are sparse. Using Carbopol solution as a model fluid, [START_REF] Naimi | Etude dynamique et thermique de l'écoulement de Couette-Taylor-Poiseuille; cas d'un fluide présentant un seuil d'écoulement[END_REF] observed axisymmetric Taylor vortices. It is indicated that the yield stress has a stabilizing effect. Actually, the focus of this study was on heat transfer rather than on hydrodynamic stability.

A natural sequel to the linear theory developed by [START_REF] Peng | Linear stability of Bingham fluids in spiral Couette flow[END_REF] and [START_REF] Landry | Stability and instability of Taylor-Couette flows of a Bingham fluid[END_REF] is to consider the stability of a circular Couette flow of a Bingham fluid with respect to a finite amplitude perturbation. In this case nonlinear effects can no longer be neglected and the linear framework used previously becomes inapplicable.

In yield stress fluids equations of motion, in addition to the quadratic nonlinearity of the inertial terms, we have also a nonlinearity in the rheological law. In the research of nonlinear phenomena, much attention has been devoted to the weakly nonlinear phase around criticality, as analytical modelling of the weak nonlinearity is possible.

Multiple-scale expansion method is generally applied for the study of weak nonlinearity. This mathematical tool is employed in the present work. Alternatively, one can use the amplitude expansion method surveyed by [START_REF] Herbert | On perturbation methods in nonlinear stability theory[END_REF] or the center manifold reduction (Fujimura 1991). Fujimura (1989Fujimura ( , 1991) ) demonstrated that these theoretical tools are equivalent in the derivation of Ginzburg-Landau equation which models the temporal evolution of the disturbance amplitude. This equation allows to determine whether the nonlinearities saturate the linear instability and at what value of the disturbance amplitude. For a Newtonian Taylor-Couette flow with fixed outer cylinder, the primary bifurcation is supercritical, i.e. the nonlinear inertial terms saturate the amplitude of the vortices. Using the amplitude expansion method, Davey (1962) investigated the structure of the supercritical Taylor vortex flow. It is shown that the distortion of the mean azimuthal velocity profile by the Reynolds stress plays a significant role in the saturation process. For a purely viscous shear-thinning fluid described, for example by the Carreau model, Topayev et al. (2019) showed that the nonlinearity of the rheological model tends to accelerate the mean flow in the annular space due to the reduction of the viscous dissipation induced by the viscosity perturbation (see also Chekila et al. (2011) for another configuration). This reduction in the viscous dissipation is modest and does not alter the supercritical nature of the primary bifurcation, i.e. the nonlinear inertial terms remain dominant. This is perhaps not surprising as at sufficiently high shear-rate, the viscosity perturbation due to a shear-rate disturbance is very weak. In the case of a Bingham fluid with the presence of a static layer attached to the outer cylinder, the reduction of the viscous dissipation can be more substantial due to the strong nonlinear increase of the viscosity near the yield surface. The role of nonlinear yield stress terms may become sufficiently important comparatively to the nonlinear inertial terms to change the nature of the bifurcation.

The objective is to investigate the first effects of nonlinear yield stress terms on the nature of the primary bifurcation, the intensity of vortex flow and the erosion of the static layer. An approach based on weakly nonlinear analysis with multiple scales method is adopted. Very few people have carried out a weakly nonlinear analysis for these yield stress fluid models. To the best of our knowledge, there is only the work of Metivier et al. (2010) dealing with plane Poiseuille-Rayleigh-Bénard flow of a Bingham fluid.

The results of this study could have applications in Couette rheometry, especially with a wide gap [START_REF] Ovarlez | Wide-gap couette flows of dense emulsions: Local concentration measurements, and comparison between macroscopic and local constitutive law measurements through magnetic resonance imaging[END_REF]. They can also be considered as a first step in the study of rotating filtration systems involving yield stress fluids [START_REF] Martinand | Linear and weakly nonlinear analyses of cylindrical couette flow with axial and radial flows[END_REF].

This article is organized as follows. In §2 we formulate the physical problem, state the governing equations and define the dimensionless parameters. The velocity and viscosity profiles of the base state are discussed and the disturbance equations are derived.

Subsequently, the linearization of the disturbance equations and the eigenvalue problem derivation for the linear stability analysis are presented in §3. In §4, the weakly nonlinear scheme based on multiple scales method is described in detail as well as the method to derive the Ginzburg-Landau equation and to obtain the first Landau constant. We focus mainly on the situation where a static layer is attached to the outer wall. In §5, we present and discuss the numerical results dealing with the influence of Bingham number on the nature of primary bifurcation and the flow structure. Finally, in §6, the relevant results of the present study are summarized.

Problem formulation

We consider the flow between two infinitely long concentric cylinders, figure 1, with inner and outer radii, R1 and R2 that rotate independently with angular velocity speed Ω1 (inner) and Ω2 (outer). The scaled momentum and mass conservation equations are where U is the velocity, P the pressure, τ the deviatoric stress tensor, and Re 1 is the inner cylinder Reynolds number:

∂ t U + Re 1 (∇U ) • U = -∇P + ∇ • τ (2.1) ∇ • U = 0, (2.2)
Re 1 = ρ R1 Ω1 d μp .
(2.3)

Here ρ and μp are the density and the plastic viscosity. The velocity vector is of the form Using the von Mises yield criterion, the dimensionless constitutive equations for Bingham fluids are:

U = U e r + V e θ + W e z ,
τ ij = 1 + B γ γij ⇐⇒ τ > B (2.4) γ = 0 ⇐⇒ τ B, (2.5)
where γ = γij γij /2 and τ = τ ij τ ij /2 are the second invariant of the strain rate γ and deviatoric stress τ tensors respectively. The components of γ are γij = U i,j + U j,i .

The Bingham number B is defined as

B = τy d μp R1 Ω1 , (2.6)
which represents the ratio of the yield stress τy to a nominal viscous stress μp R1 Ω1 / d. In the regions where the yield stress is not exceeded, the rate of strain tensor is identically zero (i.e. no local deformation occurs) and the stress tensor is undetermined. The fluid within these regions is constrained to move as a rigid body and will be hereafter referred to as "plug region". In contact with a quiescent wall, the plug zone remains static.

Two further dimensionless parameters will be used: the outer Reynolds number Re 2 and the radius ratio η:

Re 2 = ρ R2 Ω2 d μp , η = R1 R2 .
(2.7)

Remark

Peng & Zhu (2004) used the Hedström number [START_REF] Hedström | Flow of plastic materials in pipes[END_REF][START_REF] Hanks | On the flow of Bingham plastic slurries in pipes and between parallel plates[END_REF] He = ρ d2 τy /μ 2 p = B Re 1 , rather than the Bingham number. It can be interpreted as the ratio of the yield stress τy to a viscous stress μp vd / d, where vd is the viscous diffusion velocity scale. Following [START_REF] Landry | Stability and instability of Taylor-Couette flows of a Bingham fluid[END_REF], the choice of B rather than He in the present work is driven by two considerations: firstly, as it will be shown in the next section, the basic Couette flows depends solely on B, η and Re 2 /Re 1 . Secondly, in using He, large values of He may correspond to modest values of B, leading to relatively small changes in the Couette flow. Finally, the use of either He or B can be considered as simply a matter of choice.

Basic flow

The base Couette flow velocity, U b = (0, V b (r), 0), is derived from

0 = 1 r 2 d dr r 2 τ rθ , r ∈ η 1 -η , 1 1 -η
(2.8)

τ rθ = 1 + B | γrθ | γrθ ⇐⇒ |τ rθ | > B (2.9) | γrθ | = 0 ⇐⇒ |τ rθ | B , (2.10) γrθ = dV b dr - V b r , (2.11)
with boundary conditions:

V b = 1 at r = R 1 = η 1 -η , (2.12
)

V b = Re 2 Re 1 at r = R 2 = 1 1 -η .
(2.13)

The basic flow equations are fully determined by the set of parameters B, η and Re 2 /Re 1 .

From equation (2.8) (2.15)

τ rθ = τ i η 2 r 2 (1 -η) 2 where τ i = τ rθ (r = R 1 ) . ( 2
The base solutions are of three types: (1) the inner and the outer cylinders rotate with the same angular velocity, the fluid is fully unyielded in the annular gap; (ii) there may be a layer of unyielded fluid attached to the outer wall; (iii) the fluid may be fully yielded through the annular gap. The regions where the three solutions may be found can be visualized in the plane (Re 1 , Re 2 ). Fully unyielded flows are found along the line

Re 1 = ηRe 2 .
(2.16)

Partially yielded flows are found in the domains bounded by the above line and the lines

Re 2 Re 1 = 1 η ± B f (η) , (2.17)
where f (η) is defined by [START_REF] Landry | Stability and instability of Taylor-Couette flows of a Bingham fluid[END_REF])

f (η) = 1 + η 2η 2 - ln(1/η) 1 -η .
(2.18)

As an example, figure 2 shows the domains where the three solutions hold in the plane (Re 2 , Re 1 ) in the case of a narrow gap η = 0.883 and a wide gap η = 0.4. The Bingham number is fixed B = 5. In domains of the plane (Re 2 , Re 1 ) where the fluid is fully yielded or partially yielded, the velocity profile V b (r) is given by (2.19) where R 0 = min(R y , R 2 ). An illustration of basic velocity profiles V b (r) at different values of B in the case of co-rotating (Re 1 = 1000, Re 2 = 100) and counter-rotating cylinders (Re 1 = 1000, Re 2 = -100) for a narrow and a wide gap is given by figure 3.

V b (r) =        Re 2 Re 1 r R 2 + τ i R 2 1 2 r 1 R 2 0 - 1 r 2 + B r ln R 0 r sgn(τ i ) R 1 r R 0 , Re 2 Re 1 r R 2 R 0 r R 2 ,
The position R y of the yield surface is indicated by a vertical dotted line. In the static layer (R y r R 2 ), V b varies linearly with r. The nonlinear variation of the effective viscosity µ = 1 + B/ γ in the yielded zone is shown in figure 4. According to the Bingham model, µ increases from the inner wall and tends to infinity near the yield surface. The degree of nonlinearity of the rheological behavior becomes stronger with increasing B.

Remarks

-It can be shown straightforwardly that there is a similarity mapping that maps the basic solution II onto a basic solution III with an outer radius equal to the yield surface radius. Mapping of the basic solution II consists simply on rescaling the length, by substituting r = r (R 0 -R 1 ), with η 1 -η r 1 1 -η , where η = R 1 /R 0 [START_REF] Landry | Stability and instability of Taylor-Couette flows of a Bingham fluid[END_REF]. This leads to the following scalings:

B = B(R 0 -R 1 ), τi = τ i (R 0 -R 1 ) and Re 2 Re 1 = Re 2 Re 1 1 - η R 2 .
-The terms "narrow gap" and "wide gap" are of course related to the radius ratio [START_REF] Chandrasekhar | The stability of viscous flow between rotating cylinders[END_REF]; [START_REF] Donnelly | Experiments on the stability of viscous flow between rotating cylinders i. Torque measurements[END_REF]; [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF]; [START_REF] Razzak | Numerical study on wide gap Taylor Couette flow with flow transition[END_REF]).

In the narrow gap problems, the annular gap is much smaller than the mean radius and usually indicates that radius ratio, η > 0.5. For the wide gap problems, η 0.5.

From physical point of view, the critical Reynolds number Re c at which Taylor vortices appear increases with increase in radius ratio in narrow gap problems, while the opposite behavior is observed for the wide gap problems.

Perturbation equations

A perturbation of a small amplitude (u, p) is imposed to the basic flow (U b , P b ). The perturbed flow is then given by

(U b + u, P b + p) = (u, V b + v, w, P b + p) .
(2.20) The perturbation equations in the yielded parts of the flow are straightforwardly derived:

(a) (b) (c) (d)
∇ • u = 0 , (2.21) ∂u ∂t = L(u, p) + N (u) . (2.22)
The operator L(u, p) denotes the linear part of the perturbation equations, which consists of 4 parts:

L(u, p) = Re 1 LI(u) -∇p + LV(u) + BLY(u) (2.23)
describing inertial, pressure, viscous and yield stress effects. The inertial, pressure and viscous operators are identical with those for a Newtonian fluid:

LI(u) = -[(U b • ∇) u + (u • ∇) U b ] ,
(2.24)

LV(u) = ∇ • γ . (2.25)
The yield stress term is addressed below along with the nonlinear parts. The terms in N are significant for nonlinear perturbations. This includes contribution from the inertial terms and from perturbations of shear-stress tensor, where only the yield stress contributions are nonlinear. This can be written as:

N (u) = Re 1 N I(u) + BN Y(u) , (2.26)
where the inertial operator N I(u) reads

N I (u) = -(u • ∇) u .
(2.27)

Yield stress terms

The components of the shear-stress tensor can be written as:

τ ij = γij + BM ij .
(2.28)

In a part of the fluid where the yield stress is not exceeded, τ ij are indeterminate.

Conversely, in a part of the fluid, where the yield stress is exceeded, the M ij is given by

M ij = γij γ .
(2.29)

The contributions to the yield stress in the perturbation equations all come from the expressions,

M ij (U b + u)-M ij (U b )
, which for small finite u, can be found by expanding about the base flow in a Taylor series:

M ij (U b + u) -M ij (U b ) = m 1,ij (u) + m 2,ij (u, u) + m 3,ij (u, u, u) , ... (2.30)
The terms in m k,ij contain products of k components of u. For k = 1, 2, 3, the general expression of m k,ij are provided in the Supplementary section S1. For the specific base flow that we have, where only γrθ (U b ) = γθr (U b ) = 0, some simplifications are made:

m 1,ij = 0, ij = rθ, θr (2.31) m 1,ij = γij (u) γ(U b ) ij = rθ, θr (2.32) m 2,ij = sgn( γrθ (U ) b ) 2 γ2 (U b ) γ2 rθ (u) -γ2 (u) ij = rθ, θr (2.33) m 2,ij = - sgn ( γrθ (U b )) γ2 (U b ) γij (u) γrθ (u) ij = rθ, θr (2.34) m 3,ij = γrθ (u) [ γ(U b )] 3 γ2 (u) -γ2 rθ (u) ij = rθ, θr (2.35) m 3,ij = 1 2 [ γ(U b )] 3 3 2 γ2 rθ (u) -γ2 (u) γij (u) ij = rθ, θr (2.36)
In terms of m k,ij , k = 1, 2, 3, we may now define the linear (k = 1) and nonlinear (k = 2, 3) Bingham perturbation terms as follows:

LY = ∇ • m 1 and N Y k = ∇ • m k (2.37)

Linear stability analysis

The basic flow is supposed to be perturbed by an infinitesimal disturbance. The linearized perturbation equations can be written formally as:

∇ • u = 0 , (3.1) ∂u ∂t = -∇p + Re 1 LI(u) + LV(u) + BLY(u) .
(3.2)

Boundary conditions

If the fluid is fully yielded in all the annular space (case III in Fig. 2), the no-slip conditions at both walls are imposed:

u = 0 at r = R 1 (3.3) u = 0 at r = R 2 . (3.4)
If the fluid is partially yielded (case II of base solutions), in the region where the yield stress is exceeded, the components of the deviatoric stress tensor are assumed linearly perturbed and u). Therefore, it can be assumed that the yield surface will be also linearly perturbed from its initial position R y :

M ij (U b + u) -M ij (U b ) = m 1,ij (
r = R y + h with h << R 2 -R y . (3.5)
In other words, it is assumed that the plug zone is able to withstand an infinitesimal perturbation without breaking up. The continuity of the velocity components at the yield surface R y + h gives

(U b + u) [R y + h] -, z, t = (U b + u) [R y + h] + , z, t , (3.6)
Where the superscripts ± indicate that the limit is taken from each side of the yield surface. Since u = 0 uniformly in the plug zone, the linearization of the boundary conditions at the yielded side reads

u = 0 at r = R y . (3.7)
Additional conditions arise from the von Mises yield criterion and the continuity of stress (i.e. traction) through the fluid domain, which demand that γ (U b + u) = 0 at the perturbed yield surface. [START_REF] Frigaard | On the stability of Poiseuille flow of a Bingham fluid[END_REF], [START_REF] Frigaard | On three-dimensional linear stability of Poiseuille flow of Bingham fluids[END_REF], [START_REF] Landry | Stability and instability of Taylor-Couette flows of a Bingham fluid[END_REF], [START_REF] Nouar | Modal and non-modal linear stability of the plane Bingham-Poiseuille flow[END_REF]). This results in :

γrr (u) = γθθ (u) = γzz (u) = 0 at r = R y , (3.8) γrz (u) = γθz (u) = 0 at r = R y , (3.9) γrθ (u) = 2 h B sgn(τ i ) R y at r = R y . ( 3
√ u 2 + w 2 = O((r -R y ) 2 ) as r → R - y .
Indeed, the compatibility conditions (3.8) & (3.9) reduce to ∂u/∂r = ∂w/∂r = 0 at r = R y .

Finally, note that there is no kinematic condition since the yield surface is not a material surface or interface.

Summary of subsequent analysis

In a classical way, the disturbance components of velocity as well as the pressure and the yield surface disturbances are put into normal mode form. An eigenvalue problem is then derived whose numerical solution allows to obtain the critical values of Reynolds number Re 1c , axial wavenumber k c and azimuthal wavenumber m c at the onset of the instability as a function of Bingham number and external Reynolds number. In particular, we have determined the ranges of B and Re 2 for which the instability is axisymmetric. A detailed analysis is given in the Supplementary section S2

Lastly, note for the case I basic flow, the finite plug that fills the annulus cannot be perturbed by an infinitesimal perturbation. Therefore, in this case there is no linear stability problem.

Weakly nonlinear analysis

From here on, we consider only axisymmetric disturbances. In this case, the continuity equation simplifies and is satisfied via introduction of a function φ such that

u = - ∂φ ∂z and w = 1 r ∂ ∂r (rφ) . (4.1)
The pressure is eliminated by cross-differentiating the radial and axial momentum equations. Finally, the perturbation equations are written in terms of f = (φ, v) T as follows

C ∂f ∂t = Lf + N (f ) , (4.2)
where again L and N denote linear and nonlinear operators, respectively. They are splitted into inertial, viscous and yield stress parts as:

L = Re 1 LI + LV + B LY , (4.3) N = Re 1 N I + B N Y . (4.4)

Multiple scales method

As the Reynolds number is increased above the onset Re 1c , the growth rate of the perturbation is positive for any wavenumber k within a band √ ǫ around the critical wavenumber, where ǫ = (Re 1 -Re 1c )/Re 1c is the distance from the onset. Taylor expansion of the dispersion curve near its maximum shows that s ∝ ǫ and (k

-k c ) ∝ √ ǫ.
For ǫ > 0, the solution of the nonlinear problem can be described by a sum of unstable modes of the form exp (s t/τ 0 ) exp (i

k c z) exp (i √ ǫz/ξ 0 ),
∂ ∂t → ∂ ∂t + δ 2 ∂ ∂T , ∂ ∂z → ∂ ∂z + δ ∂ ∂Z . (4.6)
The fast spatial variables vary on the order of a typical wavelength. The slow variables describe the temporal and the spatial modulations of these fast variables. Furthermore, as the marginal mode is stationary, then

∂ ∂t → δ 2 ∂ ∂T . (4.7)
In the neighborhood of the critical conditions, corresponding to the onset of convection, the solution is expanded in powers of δ:

f = δf 1 + δ 2 f 2 + δ 3 f 3 + O δ 4 (4.8) Re 1 = Re 1c + δ 2 Re (2) 1 + O(δ 4 ) . (4.9)
The Reynolds number is increased by increasing the rotation rate of the inner cylinder.

Therefore, the Bingham number will be perturbed as

B + δ 2 B (2) = B -δ 2 B Re 1c Re (2) 1 + O(δ 4 ) , (4.10)
where B is the Bingham number that is used in the determination of the critical conditions. Actually, in equation (4.10), we have written B = He/Re 1 , where He is the Hedström number. In the case (II), where the fluid is partially yielded, the yield surface at r = R y is disturbed to: Since there are temporal and spatial derivatives ∂/∂t and ∂/∂z in operators C, L and N as introduced in equations (4.2)-(4.4), these operators also need to be expanded as series of δ:

r = R y + δ h 1 + δ 2 h2 + δ 3 h 3 + O(δ 4 ) . ( 4
C = C 0 + δC 1 + O(δ 2 ) , (4.13) LI = LI 0 + δLI 1 + δ 2 LI 2 + O δ 3 (4.14) N I = δ 2 N I 2 + δ 3 N I 3 + O(δ 4 ) (4.15) N Y = δ 2 N Y 2 + δ 3 N Y 3 + O(δ 4 ) . (4.16)
Formally, Taylor expansion of LV and LY is similar to that of the inertial linear operator LI. The explicit expressions of C, LI, LV , LY , N I and their subscales are given in Supplementary section S3.

Derivation of Ginzburg-Landau equation

We substitute equations (4.6), (4.7) and (4.13)-(4.16) into equation (4.2) and collect the terms at the same order of δ. In the following, the first three orders in δ are developed and the Ginzburg-Landau equation is derived by applying the solvability condition to the equation written at order δ 3 .

For axisymmetric perturbations, the linear stability analysis suggests an exchange of stability as we traverse Re 1c . Thus in the classical way, we look for periodic solutions of form: 4.17) where the subscript refers to the order in δ and the superscript to the corresponding Fourier mode, with:

f = δ Af (1) 1 E 1 + c.c. + δ 2 |A| 2 f (0) 2 E 0 + A 2 f (2) 2 E 2 + c.c. +δ 3 |A| 2 Af (1) 3 E 1 + A 3 f (3) 3 E 3 + c.c. + ..., (
f (i) j = [F ij , V ij ] and E n = e inkcz . (4.18)
In equations (4.17), the amplitude A = A(Z, T ) of the perturbation depends on slow variables. Furthermore, we have assumed that the amplitudes of the higher frequency spatial modes, which require interactions between lower frequency modes will scale accordingly. Similarly, the perturbation h of the yield surface is assumed of the form

h = δ AH 11 E 1 + c.c. + δ 2 |A| 2 H02 E 0 + A 2 H 22 E 2 + c.c. +δ 3 |A| 2 AH 13 E 1 + A 3 H 33 E 3 + c.c. + ... (4.19)
As far as the boundary conditions are concerned, in the case III of base solution (fluid fully yielded), the no-slip condition is used at the inner and outer walls. In the case III of base solution (fluid partially yielded), the conditions of continuity of velocity and stress at the disturbed yield surface are considered. They are described in detail in the Supplementary section S4.

Finally, gathering powers of E at each order in δ, we derive the following hierarchical structure.

Solution at order δ

At the first order, we recover the linear problem

L 0 f (1) 1 = 0 . (4.20)
In the case III of base solutions (fully yielded) the no-slip conditions at the walls give: of velocity and yield conditions leads to:

F 11 = DF 11 = V 11 = 0, at r = R 1 , R 2 . ( 4
F 11 = DF 11 = D 2 F 11 = 0 at r = R y , (4.22)
V 11 = 0 at r = R y , (4.23)

DV 11 = -H 11 D 2 V b at r = R y . (4.24)
Since any multiple of an eigenfunction is also a solution of the linear problem (4.20) and in order to define a reference solution, F 11 and V 11 can be normalized such that

max(V 11 ) = 1 , (4.25)
which fixes the amplitude of the perturbation. Computations indicate that the eigenfunction F 11 is a pure imaginary valued while V 11 is real. As an example, figure 5 shows the structure of the critical eigenfunctions obtained for a wide gap, η = 0.4. The main features are as follows. For small values of 0 B 0.88, we have the type III basic solution: fully yielded fluid filling the annular gap. In this range of B, the maximum of V 11 and F 11 is shifted towards the inner wall with a significant decrease of the maximum of F 11 . As B increases, the critical eigenfunctions are non-zero only in a progressively smaller yielded layer of width (R y -R 1 ).

As for the basic flow, it can be shown that every case II stability problem maps to a case II-III stability problem. This procedure is quite common in linear stability problem of viscoplastic fluid [START_REF] Landry | Stability and instability of Taylor-Couette flows of a Bingham fluid[END_REF][START_REF] Nouar | Modal and non-modal linear stability of the plane Bingham-Poiseuille flow[END_REF]. Indeed, if we rescale the radial distance as for the basic flow and we use the following definitions:

k = k (R 0 -R 1 ), s = s (R 0 -R 1 ) 2 , Re 1 = Re 1 (R 0 -R 1 ) and Re 2 = Re 2 R 0 (1 -η) (R 0 -R 1 ), we recover
an identical eigenvalue problem. The mapping was made possible because the boundary conditions are homogeneous [START_REF] Landry | Stability and instability of Taylor-Couette flows of a Bingham fluid[END_REF].

In figure 6(a) we have represented the yield surface perturbations by the linear mode

H 11 E 1 ,

Linear adjoint mode

The adjoint mode is required to obtain the Ginzburg-Landau equation. Its definition is given as

f ad , Lf 1 = L ad f ad , f 1 , (4.26)
where

f ad = [F ad , V ad ]
T is the adjoint eigenfunction, L is the linear stability operator and L ad = Re 1 LI ad + LV ad + B LY ad is the corresponding adjoint operator. In this definition, the inner product is given as where f * is the complex conjugate of f . We find that C (ℓ) , LV (ℓ) and LY (ℓ) , (ℓ = 0, 1, 2), are real and self-adjoint whereas LI (ℓ) is not self-adjoint. We have:

f , g = R0 R1 f * • g r dr , ( 4
LI ad =     0 -ik c D * V b 2ik c V b /r 0     (4.28)
The linear adjoint problem, |A| 2 E 0 arises from the nonlinear interaction of the fundamental mode with its complex conjugate and the second one proportional to A 2 E 2 arises from the nonlinear interaction of the fundamental mode with itself.

L ad (f ad ) = 0 , ( 4
-Mode 0, factor

|A| 2 E 0
This harmonic is a correction at the second order of the base flow. It is obtained by solving the following system of equations

L 0 f (0) 2 = -Re 1c N I f (1) 1 |f (-1) 1 |A| 2 E 0 -B N Y f (1) 1 |f (-1) 1 |A| 2 E 0 , (4.30)
where f

(-1) 1 = f (1) * 1
is the complex conjugate of f

(1)

1 . As previously, the boundary conditions are of two types.

In the case (II) of base solutions, the no-slip condition at the walls gives

F 02 = DF 02 = V 02 = 0 at r = R 1 , R 2 . (4.31)
In the case (III) of base solutions, the boundary conditions and yield conditions at the yield surface lead to:

F 02 = DF 02 = 0 at r = R y , (4.32) D2F 02 = -H 11 D 3 F * 11 + H * 11 D 3 F 11 at r = R y , (4.33) V 02 = H 11 H * 11 D 2 V b at r = R y . (4.34)
As for the linear eigenfunction, the nonlinear corrections are computed numerically. The results show, as expected, that F 02 = 0, i.e. there is no radial or axial mean flow.

The correction at the second order of the azimuthal velocity profile of the base state is illustrated by the profiles of V 02 represented in figure 7 for different values of B.

Near the inner cylinder, V 02 < 0, i.e. the azimuthal velocity is reduced, and near the outer cylinder or near the yield surface, V 02 > 0, i.e. the azimuthal velocity is increased.

Actually, the profiles of V 02 can be related to the outward and inward radial flows.

The radial outward flow carries fluid particles with high azimuthal momentum from the inner cylinder, increasing the azimuthal velocity near the outer cylinder. The radial inward flow carries fluid particles with low azimuthal momentum from the outer cylinder, decreasing the azimuthal velocity near the inner cylinder. For a wide gap and in the case of a Newtonian fluid, the deficit of the azimuthal velocity is higher than the surplus.

With increasing the Bingham number, the opposite result is observed. The perturbation of the yield surface by the mode H02 E 0 is obtained from the yield condition. It is given by:

H02 D 2 V b = -DV 02 -H 11 D 2 V * 11 + H * 11 D 2 V 11 -H 11 H * 11 D 3 V b at r = R y ,(4.35)
with, H02 = H 02 + H0 . The first term arises from the interaction of the fundamental mode with its complex conjugate and the second one arises from the variation of the Bingham number as the Reynolds number is increased. The variation of H02 with Bingham number is shown in figure 9. Positive values of H02 mean that the width of the plug zone is reduced. The reduction is weaker with increasing B. Nevertheless, one can highlight the high values of H02 compared with H 11 . This result shows that in the presence of a static layer, the nonlinear terms quickly become significant. -Mode 2 factor of A 2 E 2

The first harmonic mode is solution of the following system of equations:

L 0 f (2) 2 = -Re 1c N I f (1) 1 |f (1) 1 A 2 E 2 -B N Y f (1) 1 |f (1) 1 A 2 E 2 , (4.36)
with the appropriate boundary conditions:

In the case (III) of base solutions, the no-slip condition at the walls gives

F 22 = DF 22 = V 22 = 0 at r = R 1 , R 2 .
(4.37)

In the case (II) of base solutions, the boundary conditions and yield conditions at the yield surface lead to: of the reduction of energy exchange between the fundamental and its harmonic via the nonlinear yield stress terms.

F 22 = DF 22 = 0 at r = R y , (4.38) D 2 F 22 = -H 11 D 3 F 11 at r = R y , (4.39) V 22 = 1 2 H 2 11 D 2 V b at r = R y .
The perturbation of the yield surface by the mode H 22 A 2 E 2 is obtained from yield conditions and is given by:

H 22 D 2 V b = -DV 22 -H 11 D 2 V 11 - H 2 11 2 D 3 V b at r = R y . (4.41)
It is represented in figure 11 as a function of B. One notices that: (i) H 22 decreases with increasing the width of the static layer and (ii) H 22 is one magnitude order larger than

H 11 .

Remark

As for the linear problem, there is a similarity mapping that maps V 02 , F 22 and V 22 obtained in the case where the fluid in the annular space is partially yielded, onto solutions Ṽ02 , F22 = F 22 /(R 0 -R 1 ) and Ṽ22 obtained for an outer cylinder of radius represents the feedback at order (δ 3 ) on the fundamental mode through the nonlinear interactions of the fundamental with the second harmonic and with the modification of the basic state. The second one is the second harmonic.

The first component, f

(3) 1 = [F 13 , V 13 ]
T satisfies the following non-homogeneous equation:

|A| 2 A L 0 f (1) 3 = C 0 ∂A ∂T f (1) 1 -Re (2) A LI 0 - B Re c LY 0 f (1) 1 -A L 2 f (1) 1 -Re 1c |A| 2 A N I f (1) 1 |f (0) 2 + N I f (-1) 1 |f (2) 2 -B |A| 2 A N Y f (1) 1 |f (0) 2 + N Y f (-1) 1 |f (2) 2 -B |A| 2 A N Y f (1) 1 , f (1) 1 |f (-1) 1 . (4.42)
Concerning the boundary conditions:

For the case III of base solutions, we have:

F 13 = DF 13 = 0 at r = R 1 , R 2 , (4.43) V 13 = 0 at r = R 1 , R 2 . (4.44)
For the case II of base solutions, the boundary conditions combined with the yield conditions at the yield surface are:

F 13 = 0 , (4.45) DF 13 = H 11 H * 11 D 3 F 11 + H 2 11 2 D 3 F * 11 , (4.46) D 2 F 13 = -H 11 H * 11 D 4 F 11 -H * 11 D 3 F 22 -H02 D 3 F 11 -H 22 D 3 F * 11 - H 2 11 2 D 4 F * 11 - H 11 R 0 H 11 D 3 F * 11 + H * 11 D 3 F 11 , (4.47) V 13 = H 11 H * 11 D 2 V 11 + H 2 11 2 D 2 V * 11 + H 2 11 H * 11 D 3 V b , (4.48) + H * 11 H 22 + H 11 H02 D 2 V b . (4.49)
The condition on DV 13 defines H 13 , but is lengthy and omitted for brevity.

The amplitude equation is found in the usual way, as a solvability condition on the 3rd order problem. In the case where the annular space is fully yielded, the boundary conditions for f

(3) 1 are homogeneous. The solvability condition leads to:

0 = f ad , C 0 f (1) 1 ∂A ∂T (4.50) -Re (2) f ad , LI 0 - B Re 1c LY 0 f (1) 1 A -f ad , L 2 f (1) 1 A -Re 1c f ad , N I f (0) 2 |f (1) 1 |A| 2 A -Re 1c f ad , N I f (2) 2 |f (-1) 1 |A| 2 A -B f ad , N Y f (0) 2 |f (1) 1 |A| 2 A -B f ad , N Y f (2) 2 |f (-1) 1 |A| 2 A -B f ad , N Y (f (1) 1 , f (1)
1 )|f

(-1) 1 |A| 2 A.
Using the departure from the linear threshold

ǫ = Re 1 -Re 1c Re 1c = δ 2 Re (2) Re 1c , (4.51)
and after returning to the fast variables δA(Z, T

) = A ′ (z, t), ∂ ∂Z = 1 δ ∂ ∂z and ∂ ∂T = 1 δ 2 ∂ ∂t ,
the following Ginzburg-Landau equation is derived:

∂A ∂T = ǫ τ 0 A + ξ 2 0 τ 0 ∂ 2 A ∂z 2 + g 1 |A| 2 A. (4.52)
In equation ( 4.52), we have dropped the prime in A ′ and we expect no confusion to the reader. In this equation, τ 0 is the characteristic time for the instability to grow,

τ -1 0 = Re 1c f ad , LI 0 - B Re 1c LY 0 f (1) 1 f ad , C 0 f (1) 1 (4.53)
ξ 0 the coherence length

ξ 2 0 = 1 Re 1c f ad , L 2 f (1) 1 f ad , LI 0 - B Re 1c LY 0 f (1) 1 (4.54)
and g 1 the first Landau coefficient

g 1 = Re 1c f ad , N I f (0) 2 |f (1) 1 f ad , C 0 f (1) 1 + Re 1c f ad , N I f (2) 2 |f (-1) 1 f ad , C 0 f (1) 1 (4.55) + B f ad , N Y f (0) 2 |f (1) 1 f ad , C 0 f (1) 1 + B f ad , N Y f (2) 2 |f (-1) 1 f ad , C 0 f (1) 1 + B f ad , N Y (f (1) 1 , f (1)
1 )|f

(-1) 1 f ad , C 0 f (1) 1
.

The integrals are evaluated numerically by means of Clenshaw and Curtis method [START_REF] Trefethen | Spectral methods in MATLAB[END_REF]. At critical conditions, and assuming that the amplitude does not vary with the axial position, f

(3) 1 = f (3)
1H satisfies equation of the form

L 0 f (3) 1H = g 1H C 0 f (1) 1 -Re 1C N I -B N Y . (4.56)
In the case, where there is a static layer on the outer wall, the boundary conditions at the yield surface (defined in (4.46)-(4.49)) are inhomogeneous. In order to derive the

Results and discussion

Characteristic time of instability growth and coherence length

The variation of the characteristic time τ 0 , for the instability to grow, as a function of Bingham number is shown in figure 12(a). As it can be observed, τ 0 decreases with increasing B in particular when there is a static layer attached to the outer cylinder, i.e.

for B > 0.88 at η = 0.4 and for B > 8.5 at η = 0.883. This result shows that the onset of Taylor vortex flow is faster with increasing B and even more in presence of a plug zone on the outer cylinder. To account for the reduction of the annular space in presence of a static layer on the outer cylinder, τ 0 and B can be rescaled as B = B(R y -R 1 ) and τ0 = τ 0 /(R y -R 1 ) 2 . In the representation of τ0 as a function of B, the slope of the curve is substantially reduced.

Similarly the coherence length which can be related to the curvature of the marginal stability curve, decreases with increasing B as shown in figure 12(b). This means that the marginal stability curve flattens with increasing B. The strong decrease of ξ 0 with increasing B, for B > 0.88 at η = 0.4 accounts for the reduction of the width of the annular space where the fluid is yielded. This can be taken into account by using the rescaled parameter B and ξ0 = ξ 0 /(R y -R 1 ) 2 .

Transition from supercritical to subcritical bifurcation

In figures 13(a) and (b) we have represented the variation of the first Landau constant g 1 as a function of Bingham number for different values of the radius ratio η, with

Re 2 = 0. The sign of g 1 determines whether we are dealing with a supercritical or a subcritical bifurcation. The following characteristics are noticed:

(i) For a given radius ratio and at B = 0 (Newtonian fluid), g 1 is negative and the bifurcation is supercritical. With increasing B, g 1 decreases, reaches a minimum then (ii) From η η ℓ , with η ℓ = 0.89 in the present case, (thick curve in figure 13(a)) g 1 is negative and the bifurcation supercritical for all the range of values of B considered.

(iii) For η < η ℓ , g 1 is positive and the bifurcation is subcritical for a range of Bingham number whose extent is all the larger as the annular gap is wide. This is illustrated by figure 13(c) where we have represented in the plane (η, B) the boundaries of the subcritical bifurcation domain. Outside this domain, the bifurcation is supercritical.

(iv) The numerical results show that at the upper boundary of the subcritical bifurcation domain (numbered (3) in figure 13(c)), the inner to the yield surface radius

R 1 /R y ≈ η ℓ .
In other words, when R 1 /R y > η ℓ , the bifurcation is supercritical.

(v) At the lower boundary of the subcritical bifurcation domain, (numbered (1) in figure 13(c)), R 1 /R y = 0.6. Hence, for η < 0.6, subcritical bifurcation occurs only in presence of a static layer on the outer wall.

(vi) The minimum value of B at which the primary bifurcation is subcritical is B = 1.8 obtained for η = 0.6.

To shed more light on the mechanism of supercritical/subcritical instability, the first Landau constant is decomposed as

g 1 = g 1I + g 1Y + g 1N H , (5.1)
where g 1I is the contribution of nonlinear inertial terms (terms where Re 1c is in factor in equation (4.55)), g 1Y is the contribution of nonlinear yield stress terms (terms where B is in factor in equation (4.55)) and g 1N H is the contribution of nonhomogeneous boundary conditions. In figure 14, we have represented the contribution of nonlinear inertial terms and nonlinear yield stress terms to the first Landau constant for different values of η. One can notice that: (i) g 1I is negative and g 1Y is positive, i.e. the nonlinear inertial terms promote a supercritical bifurcation leading to saturation of the stationary instability at finite amplitude, whereas the nonlinear yield stress terms promote a subcritical bifurcation. The position of the curve -g 1I relatively to g 1Y accounts for the nature of the primary bifurcation described in figure 13 (iv) As for the non-homogeneous boundary conditions, the obtained values of g 1N H

show that they promote a subcritical bifurcation. However, their contribution remains very weak, almost 10 -3 smaller than (-g 1I ) as it is shown in the Supplementary figure F7. Therefore, g 1N H does not play a significant role in the transition from supercritical to subcritical bifurcation.

Obviously, the increase of -g 1I with increasing B has to be due to changes with B, Regarding the variation of g 1Y with B, where g 1Y can be written formally as 1) , two situations must be distinguished depending on whether the annular space is fully or partially yielded. The first one, where the annular space is fully yielded, is as we have a purely viscous shear-thinning fluid. The numerical results show that the increase of g 1Y with increasing B can be fitted by a cubic polynom -For counter-rotating and co-rotating cylinders, variations of the first Landau constant g 1 , the contribution of nonlinear inertial terms g 1I and that of nonlinear yield stress terms g 1Y as a function of Bingham number are qualitatively similar to those obtained for fixed outer cylinder. They are described in Appendix A.0.1 and A.0.2.

in
g 1Y = B f ad , N Y / f ad , C 0 f (

Stationary amplitudes

A steady solution of equation (4.52) is

A = A 0 e iqz with A 0 = ǫ -ξ 2 0 q 2 τ 0 (-g 1 ) 1/2 and q = k -k c .
(5.2)

Flow structure

One has to note that the numerical values of the Landau constant and hence the values of the amplitude depend on the normalization condition used for the eigenfunctions in the linear theory. However, physical quantities such as velocity components, kinetic energy or torque are independent of the normalisation. In figure 16, we have represented the profiles of the radial velocity component The change in the flow structure when departing from the critical conditions modifies the second invariant γ of the strain rate tensor γ and thus the effective viscosity.

u = -ik c AF 11 E 1 + 2ik c A 2 F * 22 E 2 + c.c. + O(A 3 ), ( 5
In figure 18, from left to right, we show, the contours of the stream function ψ = rφ at ǫ = 0.5 × 10 -2 (figure 18a), which will serve as a guide for the description of the results, the contours of γ at ǫ = 0 (figure 18b) and at ǫ = 0.5 × 10 -2 (figure 18c). The case of a wide gap η = 0.5 is considered with stationary outer cylinder and B = 1 (fluid partially

Perturbation of the yield surface

For a sufficiently large Bingham number a static layer forms and is adheres to the outer cylinder. Near the onset of the TVF regime, the disturbed yield surface is given by

r = R y + AH 11 E 1 + c.c. + |A| 2 H02 + A 2 H 22 E 2 + c.c. + O(A 3 ) .
(5.5)

In figure 19, we have represented contours of the stream function (figure 19a), the deformation of the yield surface (figure 19b) and contours of the azimuthal component of the velocity disturbance (figure 19c). The radius ratio is η = 0.4, the Bingham number B = 1 and the outer cylinder is stationary. In figure 19 (b) a zoom is made on the modification of the yield surface. At critical conditions, the yield surface is represented in dashed line and its modification near the onset of vortices in continuous line. On average, over one wavelength and at second order in amplitude, the thickness of the static layer is reduced by |A| 2 H02 . This reduction is due on the one hand to the increase in Re 1 , i.e.

increase of the speed of the inner cylinder, and on the other hand to the interaction of the fundamental mode with its complex conjugate. The numerical results show that less than 7% of the reduction in the static layer thickness is due to the increase in Re 1 and the remainder is due to the interaction of the fundamental mode with its complex conjugate.

More precisely, it is observed that the erosion of the static layer arises mainly from the nonlinear yield stress terms. Indeed, if these nonlinear terms are cancelled artificially, the erosion will be substantially reduced as it is illustrated by the curves (1) and ( 2

Conclusion

The present work focuses on the first principles understanding the influence of the yield stress on the stability of a circular Couette flow with respect to a finite amplitude perturbation. A weakly nonlinear analysis based on the multiple scales method is used as a first approach to take into account nonlinear effects. A Bingham model is used the primary bifurcation results from the competition between nonlinear inertia terms which favor a supercritical bifurcation and nonlinear yield stress terms which favor a subcritical bifurcation. The nonlinear inertia terms are dominant either in the case where the Bingham number is close to zero or in the case where the annular space is sufficiently narrow, i.e. the radius ratio is greater than a limit value η ℓ , or equivalently when the yielded fluid domain is sufficiently narrow, R 1 /R y > η ℓ . However, when η < η ℓ , the nonlinear yield stress terms become dominant and the bifurcation becomes subcritical for a range of B, the extent of which depends on Re 2 and η. To our knowledge, such variations in the nature of the primary bifurcation with B and η were not observed previously in the literature. Concerning the flow structure in the supercritical regime, it is shown that at low values of B, i.e. when the nonlinear inertial terms are dominant, the strength of the vortices decreases slightly with increasing B. However, once the nonlinear yield stress terms which account for the viscosity perturbation start to become significant, the strength of the vortices increases strongly with increasing B. This effect is more pronounced in presence of a static layer.

Our analysis provides information about the first stages in the evolution of the yield surface when departing from the critical conditions. It is shown that the static layer is reduced due to the increase in the velocity of the inner cylinder and the interaction of the fundamental mode with its complex conjugate involving nonlinear yield stress terms. The contribution of the second effect is much greater than that of the first one.

Furthermore, the yield surface is sinusoidally disturbed by the fundamental mode and the second harmonic.

The obvious next step of our analysis is to use a direct numerical simulation to study more deeply the subcritical regime by determining the stable branch and the critical Reynolds number at the onset of vortices. Another point of interest concerns the flow structure and the evolution of the static layer in a strongly nonlinear supercritical regime. Finally, it is worth pointing out the lack of experimental results in the literature.

Supplementary material. Supplementary materials are available at https://doi.org /10.1017/jfm.2023... The thin dash-dotted lines are the boundaries of the subcritical domain when the outer cylinder is fixed.

Figure 1 :

 1 Figure 1: Geometry sketch and parameters of the Couette flow of a yield stress fluid with a plug zone represented by a dashed area. The yield surface radius Ry is defined by equation (2.15).

  where U, V, W are the velocity components and e r , e θ , e z are unit vectors in the radial (r), azimuthal (θ) and axial (z) directions. Lengths are scaled with the annular gap d = R2 -R1 . The velocities are scaled with Ω1 R1 (velocity of the inner cylinder). The time is scaled with the viscous diffusion time ρ d2 /μ p . The pressure and stresses are scaled with μp R1 Ω1 / d. By convention, we take Ω1 > 0.

  .14) Therefore τ rθ does not change sign in the annulus and |τ rθ | decreases with r. Consequently, if there is an unyielded plug region in the annulus, it must be bounded inside by a yield surface, say at r = R y , and must extend to the outer wall. The position of R y is defined by |τ rθ | = B:

Figure 2 :

 2 Figure 2: Flow regimes in the plane (Re 2 , Re 1 ) at B = 5: (a) narrow gap η = 0.883 and (b) wide gap η = 0.4 .

Figure 3 :

 3 Figure 3: Basic velocity profiles, V b (r). Case of co-rotating cylinders with Re 2 = 100 and Re 1 = 1000: (a) η = 0.883 and (b) η = 0.4. Case of counter-rotating cylinders with Re 2 = -100 and Re 1 = 1000: (c) η = 0.883 and (d) η = 0.4. (1) B = 0, (2) B = 1, (3) B = 5, (4) B = 10.

Figure 4 :

 4 Figure 4: Basic viscosity profiles, µ b (r). Case of co-rotating cylinders with Re 2 = 100; Re 1 = 1000 for (a) η = 0.883 and (b) η = 0.4. Case of counter-rotating cylinders with Re 2 = -100 and Re 1 = 1000: (c) η = 0.883 and (d) η = 0.4. (1) B = 0, (2) B = 1, (3) B = 5, (4) B = 10.

  .10) These conditions are not strictly boundary conditions. Instead (3.8) & (3.9) are compatibility conditions. Indeed, each γij (u) in (3.8) & (3.9) also appears in the Bingham terms (2.32) divided by γ (U b ). With these conditions satisfied, all the terms m 1,ij (equation 2.32) remain bounded and the linear problem is thus well defined. Equation (3.10) is not required for compatibility, since m 1,rθ = 0 (equation 2.31). Instead it defines the perturbation of the yield surface h. Using the normal-mode ansatz of the perturbation, it can be shown straightforwardly that γrθ (u) = ∂v/∂r at r = R y . In other words, the yield surface perturbation depends only on the azimuthal component of the velocity disturbance. The strength of Taylor vortices

  where τ 0 is the characteristic time for the instability to grow and ξ 0 the coherence length, i.e. a characteristic length scale that governs the spatial modulation of the solution. In a similar vein,[START_REF] Manneville | Aspect-ratio dependence of transient Taylor vortices close to threshold[END_REF] interpret ξ 0 as a measure how easily the unstable mode can accommodate modulations. For small ǫ, we can separate the dynamics into fast eigenmodes and slow modulation of the form exp (s ǫ/τ 0 ) and exp (i √ ǫz/ξ 0 ). We denote δ = √ ǫ. The multiplescales approach is used to obtain the amplitude equation, which describes the slow temporal and spatial variation of the variables. The slow scales Z = δ z and T = δ 2 t (4.5) are treated as independent of the fast scales z and t. The derivatives with respect to the new variables are

  .11) It is worth noting that the 2nd order perturbation term above h2 includes the change in the Bingham number due to the change in Re from Re c to Re c + δ 2 Re (2) , as the angular rotation of the inner cylinder is increased. This second order change in Reynolds number effectively leads to a change in the dimensionless yield stress, which in turn shifts the position of the yield surface: h2 = h 2 -B Re 1c Re (2) 1 dR y dB c . (4.12)

Figure 5 :

 5 Figure 5: Eigenfunctions, (a) Imag(F 11 ) and (b) V 11 , associated with the critical mode at η = 0.4 and for different values of B: 0, 0.2, 0.4, 0.6, 0.8, 1, 2, 3, 4, 5 and 10.

  along the reduced axial position z/λ z , where λ z = 2π/k c is the axial wavelength. The case II of the base flow is considered with different values of Bingham number. The minimum of Real H 11 E 1 corresponds to the jet impingement near the yield surface and the maximum of to the radial flow from the outer to the inner cylinder. It can also be noticed in figure 6(b) that H 11 decreases significantly with increasing Bingham number.

Figure 6 :

 6 Figure 6: (a) Perturbation of the yield surface associated with the critical mode at η = 0.4 and for three different values of Bingham number: (1) B = 1, (2) B = 5 and (3) B = 10. (b) Variation of H 11 as a function of B for three different values of radius ratio: (1) η = 0.4, (2) η = 0.6 and (3) η = 0.8. Vertical dashed lines indicate the value of B from which a static layer appears on the outer cylinder.

  Furthermore, the fluid is yielded in a smaller domain of width (R y -R 1 ). Canceling artificially the nonlinear yield stress terms in (4.30) allows to highlight the contribution of the nonlinear inertial terms (figure 8(a)) and vice-versa, to highlight the contribution of the nonlinear yield stress terms (figure 8(b)) on the modification of the basic flow. The interaction of the fundamental mode with its complex conjugate through nonlinear yield stress terms accelerate the fluid in the whole yielded fluid domain. A positive correction of the basic azimuthal flow leads to a destabilizing effect and therefore may be considered as a precursor to the emergence of a subcritical bifurcation.

Figure 7 :Figure 8 :

 78 Figure 7: Wide gap geometry with η = 0.4. Modification of the base flow at the second order in δ for B = 0, 0.2, 0.4, 0.6, 0.8, 1, 2, 3, 4, 5 and 10.

Figure 9 :

 9 Figure 9: Modification at the second order in δ of the width of the static layer as a function of the Bingham number: (1) η = 0.4, (2) η = 0.6 and (3) η = 0.8. Vertical dashed line indicate the value of B from which a static layer appears on the outer cylinder.

Figure 10 :

 10 Figure 10: Correction of the linear mode at the second order in δ at critical conditions at η = 0.4 and for different values of B: 0, 0.2, 0.6, 0.8, 1, 2, 3, 4, 5 and 10. (a) Imag(F 22 (r)) and (b) V 22 (r).

Figure 11 :

 11 Figure 11: (a) Yield surface perturbation of mode 2 for η = 0.4 and three different values of B: (1) B = 1, (2) B = 5 and (3) B = 10. (b) Amplitude H 22 of the second harmonic of the yield surface perturbation as a function of Bingham number: (1) η = 0.4, (2) η = 0.6 and (3) η = 0.8.

Figure 12 :

 12 Figure 12: Variations of the characteristic time of instability (left frame) and the coherence length (right frame) as a function of Bingham number for a narrow and a wide gap: (1) η = 0.883, (2) η = 0.4

  (a): For η η ℓ = 0.89, (figure 14(a)-(b)), (-g 1I ) > g 1Y and the bifurcation is supercritical. For η < η ℓ , (figure 14(c)-(d)), (-g 1I ) > g 1Y for a range of B close to zero and for large values of B. For intermediate values of B, i.e. between the two crossings of the curves -g 1I (B) and g 1Y (B) (figure 14(c)), the bifurcation is subcritical. (ii) (-g 1I ) and g 1Y increase with B but differently depending on whether the annular space is fully or partially yielded. (iii) The abrupt change in the convexity of g 1 (B) observed in figures 13(a) and (b) arises from the change in the curvature of g 1Y (B) when a static layer appears on the outer cylinder.

Figure 13 :

 13 Figure 13: Variation of the first Landau constant as a function of the Bingham number for different values of the radius ratios: (a) from η = 0.92 until η = 0.81 by step of 0.01 and then η = 0.79, 0.77, 0.75, 0.73, 0.7, 0.66 and 0.6. (b) η = 0.6, 0.5 and 0.4. (c) Boundaries of the subcritical domain in the plane (η, B): (1) and (2) are the lower boundaries and (3) is the upper boundary.

  either Re 1c or in the factor term of Re 1c in equation (4.55). This last term contains quadratic products of the eigenfunctions and their derivatives which are non-zero in a progressively smaller yielded fluid domain. The numerical results show that more than 99% of the variation of g 1I arises from the variation of Re 1c as B increases. Furthermore, it is observed that -g 1I can be fitted by an affine function of Re 1c , with different coefficients, in partially and fully yielded domains. The factor term of Re 1c in (4.55) varies mainly with η.

  order. The nonlinear yield stress term N Y derived from mk ij expressions (equations (2.31)-(2.36)) contains products of eigenfunctions and their derivatives divided by γk (U b ) (k = 1, 2, 3). It is precisely these terms, i.e. 1/ γk (U b ) that are responsible for the strong increase in g 1Y with increasing B. Indeed, at the outer wall and at critical conditions, γ(U b )| R2 decreases linearly with increasing B and tends to zero as B approaches B sℓ the value of B at which a static layer appears at the outer wall, i.e. γ(U b )| R2 ∝ (B sℓ -B). Therefore, 1/ γ(U b )| R2 increases sharply allowing the nonlinear yield stress term to become dominant as B → B sℓ . In the second situation, where the annular space is partially yielded, g 1Y increases linearly with B. The quantity f ad , N Y / f ad , C 0 f (1) remains almost constant as B increases. Here, the characteristics of the flow just after yielding, where γ(U b ) → 0, plays the dominant role. Recall that the compatibility conditions (Supplementary equations (4.8)-(4.13)) ensure that the limit of N Y is finite when γ(U b ) tends to zero. Remarks -A detailed study of the contribution of the different terms that intervene in the g 1I and g 1Y expressions (equation 4.55) is provided in section S5. The data show that the feedback of the mean flow correction plays an important role.

Figure 15 Figure 14 :Figure 15 :Figure 16 :

 15141516 Figure15shows the bifurcation diagram obtained in the case where k = k c , i.e. q = 0, at different Bingham numbers. For a wide gap η = 0.4, near the critical conditions, the stable stationary equilibrium amplitude in the supercritical regime is shown in continuous line. At low values of B (B 0.5), the nonlinearities dominated by the nonlinear inertial

  .3) at z = 0, ǫ = 0.01, in a wide and a narrow gap and for different Bingham numbers. These profiles represent the outward radial flow between the inner and outer cylinders, near the onset of the Taylor vortex flow regime. For values of B small enough comparatively to that for which a plug zone appears on the outer cylinder, the increase of |g 1I | with B inner to the outer cylinder and low momentum fluxes from near the outer to the inner cylinder, respectively. These result in the formation of positive and negative azimuthal streaks. They are determined by the azimuthal velocity expressionv = AV 11 + A 2 V 22 + |A| 2 V 02 + c.c. + O(A 3 ) (5.4)Again, at low values of B, the nonlinear inertial terms are dominant and have a stabilising effect. Once, the nonlinear yield stess terms become significant the intensity of the azimuthal streaks increases strongly with B and even more in presence of a static layer as it is shown in figure17

  yielded). At the onset of Taylor vortices (ǫ = 0), γ = γb (r) decreases monotonically from the inner wall to the yield surface represented by a dashed line. At ǫ = 0.5 × 10 -2 , in the outward radial flow region, i.e. around z/λ z = 1, γ < γb near the inner wall, γ = γb around the center of vortices and γ > γb at further radial positions. An opposite behavior is observed in the inward radial flow region, i.e. around z/λ z = 0.5 and 1.5. Actually, an analysis of the different γij (U b + u) shows that γ is dominated by | γrθ (U b + u)|. In the outward radial flow region, γrθ (u) is positive and maximum at the inner wall ( γ < γb ), it decreases and vanishes around the center of Taylor vortices ( γ = γb ) and then becomes negative ( γ > γb ). Of course, in the regions where γ > γb , we have µ < µ b and vice versa.

Figure 18 :

 18 Figure 18: Wide gap, η = 0.4 with B = 1 and stationary outer cylinder. (a) Contours of the stream function ψ = rφ at ǫ = 0.5 × 10 -2 : (+) clockwise rotation, (-) anticlockwise rotation. The dashed line indicates the position of the yield surface at critical conditions. (b) Contours of the second invariant γ of the stain rate tensor γ at the onset of Taylor vortices, ǫ = 0. (c) Contours of γ at ǫ = 0.5 × 10 * 2 .

  figure 19 (b).Around the average radius of the yield surface R y + |A| 2 H02 , a disturbance of small amplitude occurs due to the fundamental mode AH 11 E 1 +c.c. and to the second harmonic (interaction of the fundamental with itself) A 2 H 22 E 2 + c.c.. This disturbance of the yield surface shows a 'hump' in the region where there is an inward radial flow with a negative azimuthal streak (blue contours in Figure19 (b)) and a trough in the region where there is an outward radial flow with a positive azimuthal streak (yellow contours). However, this disturbance around the average radius is quite weak comparatively to the reduction of the width of the plug layer.Qualitatively similar results are obtained for a higher Bingham number, as shown in figure20, for B = 1.4. However, for the same relative distance to the onset of Taylor vortices, the erosion of the static layer is less pronounced for a larger static layer width.

Figure 19 :

 19 Figure 19: Wide gap, η = 0.4 with stationary outer cylinder, B = 1 and ǫ = 0.5 × 10 -2 . (a) Contours of the stream function. The dashed line indicates the position of the yield surface at critical conditions. (b) Zoom on the modification of the yield surface. (1) Nonlinear yield stress terms are artificially cancelled. (2) Nonlinear yield stress terms are taken into account. (c) Contours of azimuthal velocity perturbations : (blue) negative azimuthal streaks, (yellow) positive azimuthal sreaks.

Figure 20 :

 20 Figure 20: Wide gap, η = 0.4 with stationary outer cylinder, B = 1.4 and ǫ = 0.5 × 10 -2 . (a) Contours of the stream function. (b) Zoom on the modification of the yield surface. (1) Nonlinear yield stress terms are artificially cancelled. (2) Nonlinear yield stress terms are taken into account. (c) Contours of azimuthal velocity perturbations : (blue) negative azimuthal streaks, (yellow) positive azimuthal sreaks.

Figure 21 :Figure 22 :

 2122 Figure 21: Counter-rotation of the outer cylinder. Variation of the first Landau constant as a function of the Bingham number for different values of the radius ratios at Re 2 = -50: (a) from η = 0.9 until η = 0.75 by step of 0.01 and then η = 0.73. (b) η = 0.6, 0.6, 0.5 and 0.4. (c) Boundaries of the subcritical domain in the plane (η, B: (1) and (2) are the lower boundaries, (3) is the upper boundary. The thin dash-dotted lines are the boundaries of the subcritical domain in the case where the outer cylinder is fixed. (d) Influence of the outer Reynolds number on the first Landau constant for η = 0.7.

  Solution at order δ 2 : Quadratic modes At order δ 2 , the solution has two components. The first component proportional to

	4.2.3.
	.29)
	is subject to appropriate boundary conditions, matching those of the linear problem,
	i.e. (4.21), (4.22) and (4.23). As the adjoint is linear and can be arbitrarily scaled, we
	normalize the adjoint eigenfuctions so that the maximum of the adjoint azimuthal velocity
	max (V ad ) = 1.

Acknowledgments. This research is funded by the Science Committee of the Ministry of Education and Science of the Republic of Kazakhstan (Grant No. AP14972680).

Declaration of interests.

The authors report no conflict of interest.

Author ORCIDs.

C. Nouar https://orcid.org/0000-0002-0457-6694; I. Frigaard https://orcid.org/0000-0002-0970-240X.

Author contributions. All authors contributed equally to analysing data and reaching conclusions, and in writing the paper.

solvability condition, we decompose, as in [START_REF] Sen | On the stability of plane Poiseuille flow to finiteamplitude disturbances, considering the higher-order landau coefficients[END_REF] and [START_REF] Bouteraa | Weakly nonlinear analysis of Rayleigh-Bénard convection in shear-thinning fluids: nature of the bifurcation and pattern selection[END_REF], f

(3) 1 into homogeneous f

(3) 1H and inhomogeneous parts f (4.57) where, f

1N H is a correction term that accounts for the non-homogeneity of the boundary conditions. Furthermore at critical conditions f

(3) 1

satisfies an equation that can be written as:

Substituting (4.57) into (4.58), we obtain

By applying the solvability condition to (4.59), we get

The technique of solution adopted is to iterate a few times between (4.57) and (4.60).

At the start, f

1N H is assumed to be identically zero in (4.60). A first approximation of g 1 is then obtained: g

(1) 1 = g 1H . This is put into (4.58), which is solved at the critical conditions, with non-homogeneous boundary conditions, to obtain a first approximation of f

(3) 1 . Using (4.57) a first approximation of f

1N H is put into (4.60). This process is repeated until it converges to a desired level of accuracy. Note that equation (4.58) is solved with an additional condition

where r max the radial position at which V 11 (r max ) = 1. This normalization was suggested by [START_REF] Herbert | Nonlinear stability of parallel flows by high-ordered amplitude expansions[END_REF][START_REF] Herbert | On perturbation methods in nonlinear stability theory[END_REF]; [START_REF] Sen | On the stability of plane Poiseuille flow to finiteamplitude disturbances, considering the higher-order landau coefficients[END_REF]; [START_REF] Generalis | Range of validity of weakly nonlinear theory in the Rayleigh-Bénard problem[END_REF].

Without this normalization, f

(3) 1

is defined up to an arbitrary multiple of the solution

1 . A validation of the procedure used is provided in the Supplementary section S5.

A.0.2. Co-rotation of the outer cylinder

The variation of g 1 as a function of B at Re 2 = 100 is displayed in figure 22(a)

for different radius ratios. Probably, the most interesting points that can be noted are: The radius ratio limit η ℓ , above which the primary bifurcation remains supercritical, increases slightly with Re 2 . All the results obtained for the variation of η ℓ as a function of Re 2 , (in co-and counter-rotation) are summarised in Supplementary figure F8.