M Bouteraa

C Nouar

C Metivier

E Plaut

A Kalck

Weakly nonlinear analysis of Rayleigh-Bénard convection in shear-thinning fluids: nature of the bifurcation and pattern selection

A linear and weakly nonlinear analysis of convection in a layer of shear-thinning fluids between two horizontal plates heated from below is performed. The objective is to examine the effects of the nonlinear variation of the viscosity with the shear rate on the nature of the bifurcation, the planform selection problem between rolls, squares and hexagons, and the consequences on the heat transfer coefficient. Navier's slip boundary conditions are used at the top and bottom walls. The shear-thinning behavior of the fluid is described by the Carreau model. By considering an infinitesimal perturbation, the critical conditions, corresponding to the onset of convection, are determined. At this stage, non-Newtonian effects do not play. The critical Rayleigh number decreases and the critical wave number increases when the slip increases. For a finite amplitude perturbation, nonlinear effects enter in the dynamic. Analysis of the saturation coefficients at cubic order in the amplitude equations shows that the nature of the bifurcation depends on the rheological

Introduction

When a thin horizontal fluid layer is heated from below and cooled from above, a density stratification appears because of the thermal expansion of the fluid. This stratification is potentially unstable: when the temperature difference between the bottom and the top exceeds a threshold value controlled by the viscosity and heat diffusivity, by a small amount, convection sets in various forms of ordered regular patterns. Since the pioneering studies of Bénard (1900) and Rayleigh (1916), a large number of theoretical and experimental investigations were devoted to the study of this buoyancy-driven instability.

Reviews can be found in [START_REF] Getling | Rayleigh-Bénard convection. Structures and dynamics[END_REF] and [START_REF] Bodenschatz | Recent developments in Rayleigh-Bénard convection[END_REF]. Some of these studies were concerned with the nonlinear competition between different structures that develop above the linear convective threshold. It is found that, under Boussinesq approximations with a linear variation of the density with the temperature, rolls are stable right above onset [START_REF] Schluter | On the stability of steady finite amplitude convection[END_REF]. If the Boussinesq approximation is invalid, hexagons are preferred to rolls [START_REF] Busse | The stability of finite amplitude cellular convection and its relation to an extremum principle[END_REF][START_REF] Busse | Nonlinear properties of thermal convection[END_REF] due to triad wavevector resonance. Compara-tively to the Newtonian fluids, very few studies were devoted to non-Newtonian fluids.

In the following, a literature review on convection in a horizontal layer of non-Newtonian fluid heated from below and cooled from above is presented. Most non-Newtonian fluids have two common properties: viscoelasticity and shear-thinning. Polymer and colloïd solutions as well as particulate dispersions display this behavior above a certain concentration threshold. The influence of an elastic response, particularly, the possibility of oscillatory convection due to the elastic restoring forces has been discussed in the literature. According to [START_REF] Vest | Overstability of a viscoelastic fluid layer heated from below[END_REF], [START_REF] Sokolov | Convective stability of a general viscoelastic fluid heated from below[END_REF], [START_REF] Shenoy | Advances in heat transfer[END_REF] and [START_REF] Larson | Instabilities in viscoelastic flows[END_REF], viscoelastic effects may in principle produce an oscillatory instability at a lower Rayleigh number than the Newtonian stationary mode. However, the observation of the oscillating cells requires a very high temperature difference incompatible with realistic experimental conditions [START_REF] Larson | Instabilities in viscoelastic flows[END_REF]. Oscillatory convection was observed for binary viscoelastic fluids, when the binary fluid aspects are significant compared to the thermal diffusion, such as in DNA suspension [START_REF] Kolodner | Oscillatory convection in viscoelastic dna suspensions[END_REF]. The problem of pattern selection in viscoelastic fluids has also been considered in the literature, e.g., by [START_REF] Li | Finite-amplitude Rayleigh-Bénard convection and pattern selction for viscoelastic fluids[END_REF]. Using an Oldroyd-B model, they found that, near onset, rolls or hexagons can be stable, depending on secondary parameters.

Hereafter, we neglect the elastic response. We focus only on the shear-thinning effects, i.e., the influence of nonlinear decrease of the viscosity with the shear-rate.

Review on Rayleigh-Bénard convection in shear-thinning fluids

To our knowledge, the first experimental investigation of convection in a shear-thinning fluid layer confined between two horizontal plates was carried out by [START_REF] Pierre | Experimental investigation of natural convection heat transfer in confined space for non-newtonian fluid[END_REF].

The fluids used were aqueous solutions of Methocel (1w%) and Carbopol 934 (0.5, 0.75 and 1w%). The rheological behavior of these fluids was described by a power-law model, with a shear-thinning index ranging between 0.4 and 1. The results were presented in terms of a correlation relating the Nusselt number N u to Rayleigh and Prandtl numbers, for 10 5 Ra 10 6 . Later on [START_REF] Tsuei | Free convection heat transfer in a horizontal layer of non newtonian fluid[END_REF] extended this correlation to a wider range, 10 3 Ra 10 6 . For power-law fluids, Rayleigh and Prandtl numbers are defined with a viscosity calculated at a characteristic shear-rate, which is the inverse of the thermal diffusion time. [START_REF] Tien | Thermal instability of a horizontal layer of non newtonian fluid heated from below[END_REF] attempted to establish a stability criterion for shear-thinning fluids described by a power-law model. As indicated by the authors, the linear marginal stability curve cannot be determined, because of the unphysical infinite viscosity, at zero shear-rate, introduced by the rheological model. An approximate method was used for the determination of the critical Rayleigh number. It was based on the energy principle of [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF]: Instability occurs at the minimum temperature gradient at which a balance can be steadily maintained between the kinetic energy dissipated by viscosity and the internal energy released by the buoyancy force. The solution in the limit of Newtonian fluids was used for power-law model at zero shear-rate. The critical Rayleigh number was determined for two convective patterns: rolls and hexagons.

The authors found that the critical Rayleigh number decreases when the shear-thinning index decreases. This evolution was also found in their experimental measurements using aqueous solutions of carboxymethylcellulose, described by a power-law model with a shear-thinning index 0.75 n p 1. [START_REF] Liang | Experiments on buoyancy-driven convection in non-newtonian fluids[END_REF] conducted an experimental study of the buoyancy-driven convection in horizontal layers of dilute aqueous solutions of polyacrylamide (Separan AP 30 at 0.5 and 1w%). These fluids are shear-thinning with approximately constant viscosity at low shear-rate. The variation of the viscosity was about one order of magnitude over a strain-rate variation of three orders of magnitude. [START_REF] Liang | Experiments on buoyancy-driven convection in non-newtonian fluids[END_REF] found that the critical Rayleigh number is practically the same as for a Newtonian fluid and that the shear-thinning behavior tends to increase the heat transfer.

The first numerical computation of the onset of convection in a horizontal layer of a shear-thinning fluid was done by [START_REF] Ozoe | Hydrodynamic stability and natural convection in Ostwaldde-Waele and Ellis fluids: The development of a numerical solution[END_REF]. Two rheological models were considered: power-law and Ellis model. This later model has the advantage to converge to the Newtonian behavior in the limit of zero rate of strain. However, for certain range of rheological parameters, the viscosity in the Ellis model is not differentiable at zero shear-stress. The computations were carried out for roll-cells with both rigid and dragless vertical boundaries. This later case corresponds to a roll inside a periodic row of counter-rotating rolls. The critical Rayleigh number was obtained by extrapolating the Nusselt (N u) curve to N u = 1. Qualitatively, the influence of the shear-thinning on the critical Rayleigh number and on the Nusselt number are similar to those obtained by [START_REF] Tien | Thermal instability of a horizontal layer of non newtonian fluid heated from below[END_REF]. Nevertheless, the critical values found by Ozoe and Churchill are higher than those given by [START_REF] Tien | Thermal instability of a horizontal layer of non newtonian fluid heated from below[END_REF]. The same trends were observed for rigid and dragless vertical boundaries. In a companion paper, [START_REF] Ozoe | Hydrodynamic stability and natural convection in newtonian and non-newtonian fluids heated from below[END_REF] presented the computed results in terms of a correlation relating the Nusselt number to the shear-thinning index:

N u N u newt = 0.87n 2 p -2.28n p + 2.41, (1.1)
for Ra c Ra 2 Ra c and 0.5 n p < 1.

The case of very viscous fluids (infinite Prandtl number) with a power-law model and 0.11 n p 1 was considered by [START_REF] Parmentier | A study of thermal convection in non-newtonian fluids[END_REF]. Numerical solutions were obtained in two-dimensional periodic convective modes. [START_REF] Parmentier | A study of thermal convection in non-newtonian fluids[END_REF] shows that when the Rayleigh number is based on a strain-rate squared averaged viscosity, a good correlation for the heat transfer is obtained over a wide range of Rayleigh numbers. Three decades later, the two dimensional Rayleigh-Bénard convection for a power-law fluid in a rectangular cavity with adiabatic vertical walls, was investigated numerically by [START_REF] Lamsaadi | Natural convection of power law fluids in a shallow horizontal rectangular cavity uniformly heated from below[END_REF] and [START_REF] Alloui | Kalck convection of power-law fluids in a shallow cavity heated from below by a constant heat flux[END_REF]. Their findings are qualitatively in agreement with the literature. The decrease of power-law index, n induces a precocious onset of convection and enhances the rate of heat transfer. [START_REF] Alloui | Kalck convection of power-law fluids in a shallow cavity heated from below by a constant heat flux[END_REF] explain that for power law fluids, the system is unconditionally stable to infinitesimal disturbances. Note that [START_REF] Lamsaadi | Natural convection of power law fluids in a shallow horizontal rectangular cavity uniformly heated from below[END_REF] and [START_REF] Alloui | Kalck convection of power-law fluids in a shallow cavity heated from below by a constant heat flux[END_REF] considered a shallow rectangular cavity heated and cooled with uniform heat fluxes. In this case, the system convects with one cell. Relying on a parallel flow concept for infinite aspect ratio, [START_REF] Alloui | Kalck convection of power-law fluids in a shallow cavity heated from below by a constant heat flux[END_REF] show that the onset of convection occurs at subcritical Rayleigh number.

Weakly nonlinear stability analysis of thermal convection for shear-thinning fluid between two plates maintained at different temperatures was performed by [START_REF] Balmforth | Weakly nonlinear viscoplastic convection[END_REF]. Assuming a two-dimensional situation, with stress-free boundary conditions, the authors found that when the degree of shear-thinning α = |dµ/dΓ| Γ=0 is greater than 24/(601 π 4 ) the bifurcation becomes subcritical. In the previous expression, the viscosity µ and the second invariant of the strain rate deformation Γ (defined by equation 2.7) are rendered dimensionless using the zero shear-rate viscosity and thermal-diffusion time as characteristic scales. Recently, a systematic weakly nonlinear analysis for Carreau fluid in two and three-dimensional situation with stress-free boundary conditions was carried out by [START_REF] Albaalbaki | Pattern selection in the thermal convection of nonnewtonian fluids[END_REF]. When the convection takes place in the form of rolls, the threshold value of α for a subcritical bifurcation found was 14×10 -4 , in disagreement with [START_REF] Balmforth | Weakly nonlinear viscoplastic convection[END_REF]. [START_REF] Albaalbaki | Pattern selection in the thermal convection of nonnewtonian fluids[END_REF] found also, that depending on the degree of shear-thinning, the fluid can convect in the form of rolls, squares or hexagons. This result is surprising. Usually, near the onset of convection, hexagons are observed in convection systems lacking midplane reflection symmetry such as in fluids with strongly temperature-dependent viscosity [START_REF] Palm | On the tendency towards hexagonal cells in steady convection[END_REF][START_REF] Golubitsky | Symmetries and pattern selection in rayleigh-bénard convection[END_REF] or in Bénard-Marangoni convection [START_REF] Thess | Planform selection in Bénard-Marangoni convection: l hexagons versus g hexagons[END_REF]. Indeed, in such systems, quadratic terms present in amplitude equations enable triadic resonant wavevector in-teractions that can explain the occurence of hexagons near onset. In §5.5.4, it will be shown that under Boussines approximations, with identical boundary conditions at the two horizontal plates, the non-Newtonian terms do not break the midplane reflection.

1.2. Objectives, methodology and outline of the paper

Here, we consider shear-thinning fluids with a finite zero shear-rate viscosity µ 0 . The [START_REF] Carreau | Rheological equations from molecular network theories[END_REF] model (equation 2.8) is adopted to describe the nonlinear variation of the viscosity µ with the second invariant of the strain-rate tensor Γ. This model is chosen because it has a sound theoretical basis, and is C ∞ with respect to Γ, unlike the power-law model or the general Carreau-Yasuda model, which are singular at Γ = 0.

Interestingly, the Carreau model approaches the power-law model, as the viscosity µ 0 or the characteristic time λ of the fluid become large.

In light of the previous works, there are two points that need to be re-examined. The first one concerns the critical value of the degree of shear thinning, α, above which the bifurcation becomes subcritical. This value was determined only for stress-free boundary conditions (SFBC) which are not quite physical, and moreover the existing results are contradictory. The second point concerns the competition between different patterns of convection near the criticality. Only [START_REF] Albaalbaki | Pattern selection in the thermal convection of nonnewtonian fluids[END_REF] dealt with this problem, using SFBC.

The purpose of the present work is to revisit the Rayleigh-Bénard problem for shearthinning fluids using more general boundary conditions, with both slip and stress, i.e., Navier-type boundary conditions with a slip parameter. Note that, for Newtonian fluids with a Navier slip boundary conditions, only the linear stability of the Rayleigh-Bénard problem has been studied by [START_REF] Webber | The destabilizing effect of boundary slip on Bénard convection[END_REF], [START_REF] Kuo | Effects of slip boundary conditions on Rayleigh-Bénard convection[END_REF]. Here, a general weakly nonlinear analysis is performed. The calculation of the saturation coefficient at the cubic order allows to determine the nature of the bifurcation depending on the slip and rheolog- we examine the relevance of the principle of maximum heat transfer for non-Newtonian fluids. Calculation at higher order allows one to characterize the convection for a significant departure from the critical conditions, in particular, through a correlation for the Nusselt number using a generalized Rayleigh number as suggested by [START_REF] Parmentier | A study of thermal convection in non-newtonian fluids[END_REF].

The article is organized as follows. In section 2, the governing equations of mass, momentum and energy are presented in dimensionless form. Section 3 deals with the linear stability theory. The influence of the slip parameter on the critical Rayleigh number and critical wavenumber is examined. In section 4, the main steps of the weakly nonlinear analysis are outlined. The results are presented and discussed in section 5. The nature of the bifurcation is determined and the competition between different patterns near onset is analyzed. It is found that only rolls are stable. This result is confirmed in section 6, by computing higher-order Landau constants. The flow structure, the modification of the viscosity field and the heat transfer for steady rolls are described in section 7. Finally, section 8 is devoted to a concluding discussion.

Physical and mathematical model

General equations and parameters

Hereafter, quantities with hats are dimensional quantities. We consider a layer of shearthinning fluid of depth d confined between two horizontal plates, infinite in extent, which are perfect heat conductors. The bottom and top plates are kept at constant temperatures, respectively T0 + δ T /2 and T0 -δ T /2, with δ T > 0. The fluid has density ρ, thermal diffusivity κ, thermal expansion coefficient β and viscosity μ0 at zero shear-rate.

Because of the thermal expansion, the temperature difference between the two plates induces a vertical density stratification. Heavy cold fluid is above a light warm fluid.

For small δ T , the fluid remains at rest and the heat is transferred by conduction. The hydrostatic solution for the pressure P and the temperature profile are:

d P dẑ = -ρ0 ĝ 1 -β T -T0 and Tcond -T0 = δ T 2 1 - 2ẑ d , (2.1)
with ρ0 the fluid density at the reference temperature and ĝ the acceleration due to gravity. The z-axis is directed upwards, with its origin located at the bottom plate.

The stability of the hydrostatic solution is considered by introducing temperature and pressure perturbation as well as a fluid motion. Boussinesq approximation is adopted, i.e., the temperature dependence of the fluid properties can be neglected except for the temperature-induced density difference in the buoyancy force. The heat production due to viscosity is neglected. Using the units d2 /κ, d, κ/ d and ∆ T for time, length, velocity and temperature, the dimensionless perturbation equations are:

∇ • u = 0 , (2.2) 1 P r ∂u ∂t + (u • ∇) u = -∇p + Ra θ e z + ∇ • τ , (2.3) ∂θ ∂t + u • ∇θ = u • e z + ∇ 2 θ .
(2.4)

Here, e z denotes the unit vector in the vertical direction, u(x, t) the fluid velocity, p(x, t) and θ(x, t) represent the pressure and temperature deviations from their values in the conductive state. We denote (x, y, z) as the components of the position vector x, and (u, v, w) the components of the velocity vector u. The Rayleigh number Ra and the Prandtl number P r are Ra = ρ0 ĝ β δ T d3 κ μ0 ; P r = μ0 ρ0 κ .

(2.5)

Generally, for non-Newtonian fluids, P r ≫ 1, i.e., the viscous diffusion time is shorter than the thermal diffusion time.

Rheological model and parameters

The fluid is assumed to be purely viscous and shear-thinning. The viscous stress-tensor

τ = µ (Γ) γ with γ = ∇u + (∇u) T (2.6)
the rate-of-strain tensor, of second invariant

Γ = 1 2 γij γij . (2.7)
The Carreau model is given by

μ -μ∞ μ0 -μ∞ = 1 + λ2 Γ nc-1 2 , (2.8)
with μ0 and μ∞ the viscosities at low and high shear rate, (n c < 1) the shear-thinning index, λ the characteristic time of the fluid. The location of the transition from the Newtonian plateau to the shear-thinning regime is determined by λ, since 1/ λ defines the characteristic shear rate for the onset of shear-thinning. Increasing λ reduces the Newtonian plateau to lower shear rates. The infinite shear viscosity, μ∞ , is generally associated with a breakdown of the fluid, and is frequently significantly smaller (10 -3 to 10 -4 times smaller) than μ0 , see [START_REF] Bird | Dynamics of polymeric liquids[END_REF] and [START_REF] Tanner | Engineering rheology[END_REF]. The ratio μ∞ /μ 0 will be thus neglected in the following. The dimensionless effective viscosity is then

µ = μ μ0 = 1 + λ 2 Γ nc-1 2 with λ = λ d2 /κ .
(2.9)

The Newtonian behavior, μ = μ0 , is obtained by setting n c = 1 or λ = 0.

For a small amplitude disturbance, the viscosity can be expanded about the hydrostatic solution,

µ = 1 + n c -1 2 λ 2 Γ + 1 2 n c -1 2 n c -3 2 λ 4 Γ 2 + ... (2.10)
At lowest nonlinear order, a relevant rheological parameter is the 'degree of shearthinning'

α = dµ dΓ Γ=0 = 1 -n c 2 λ 2 .
(2.11)

Boundary conditions with slip

The plates are not permeable, i.e., u • n = 0, (2.12) n being the unit vector normal to the wall, pointing towards the fluid. Concerning the component of the fluid velocity tangent to the plates, it can be significantly affected by liquid-surface wall interactions. Polymer melts and solutions usually slip at a plane wall [START_REF] Denn | Extrusion instabilities and wall slip[END_REF]. This slip may result from an adhesive failure of the polymer chains at the solid surface or from disentanglement between chains adsorbed to the wall and those in the polymer bulk [START_REF] Brochard | Shear-dependent slippage at a polymer/solid interface[END_REF][START_REF] Baljon | Adhesion and friction of thin films[END_REF]. Another class of complex fluids prone to wall slip are colloidal suspensions and emulsions. In this case, slip arises from a depletion of particles adjacent to the shearing surfaces [START_REF] Barnes | A review of the slip (wall depletion) of polymer solutions, emulsions and particle suspensions in viscometers: its cause, character, and cure[END_REF]. The wall slip is often modeled macroscopically using Navier's slip law. This law is adopted in the present study, to take into account of a possible wall slip. For purely viscous non-Newtonian fluids, the tangent velocity u t is proportional to the tangent wall shear stress τ t via an empirical coefficient L s , called slip parameter [START_REF] Ferras | Analytical solutions for newtonian and inelastic non-newtonian flows with wall slip[END_REF]: This symmetry plays an essential role in the pattern selection.

u t = L s τ t , (2.13) with u t = u -(u • n) n, τ t = τ • n -τ n n and τ n = n • τ • n. No-

Reduction: elimination of the pressure

The pressure field is eliminated by applying the curl to (2.3). Then, we take curl curl of (2.3). Using the continuity equation, and projecting onto e z , we get the following evolution equations for the vertical vorticity ζ = ∂v/∂x -∂u/∂y and the vertical velocity w:

1 P r ∂ζ ∂t + e z • ∇ × [(u • ∇) u] = ∆ζ + e z • ∇ × [∇ • (µ -1) γ] ,
(2.15)

1 P r ∂∆w ∂t -e z • [∇ × ∇ × ((u.∇) u)] = ∆ 2 w + Ra ∆ H θ - (2.16) [∇ × ∇ × [∇ • (µ -1) γ]] • e z , ∂θ ∂t + u • ∇θ = w + ∆θ, (2.17)
where the 'horizontal Laplacian'

∆ H = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 .
From the continuity equation and the vertical vorticity definition, one deduces the horizontal velocity components:

∆ H u = - ∂ 2 w ∂x∂z - ∂ζ ∂y ; ∆ H v = - ∂ 2 w ∂y∂z + ∂ζ ∂x .
(2.18)

The boundary conditions are

w = 0, θ = 0, u = L s µ ∂u ∂z , v = L s µ ∂v ∂z , at z = 0, (2.19a) w = 0, θ = 0, u = -L s µ ∂u ∂z , v = -L s µ ∂v ∂z , at z = 1. (2.19b)
For horizontal Fourier modes to be used below, it is interesting to combine the two boundary conditions at each plane by taking their derivatives with respect to x and y, to obtain the equivalent conditions

∂w ∂z = L s µ ∂ 2 w ∂z 2 - ∂µ ∂x ∂u ∂z - ∂µ ∂y ∂v ∂z , ζ = L s µ ∂ζ ∂z + ∂µ ∂x ∂v ∂z - ∂µ ∂y ∂u ∂z , at z = 0, ∂w ∂z = -L s µ ∂ 2 w ∂z 2 - ∂µ ∂x ∂u ∂z - ∂µ ∂y ∂v ∂z , ζ = -L s µ ∂ζ ∂z + ∂µ ∂x ∂v ∂z - ∂µ ∂y ∂u ∂z , at z = 1.

Linear stability analysis

Direct eigenvalue problem: critical conditions

In the linear theory, u and θ are assumed infinitesimal. By neglecting the nonlinear terms in (2.15)-(2.17), one obtains the linear problem:

1 P r ∂ζ ∂t = ∆ζ, (3.1) 1 P r ∂∆w ∂t = ∆ 2 w + Ra∆ H θ, (3.2) ∂θ ∂t = w + ∆θ. (3.3)
No non-Newtonian effects enter the problem at this order. The vertical vorticity decouples and obeys a diffusion equation and thus can be ignored in the linear theory. For equations

(3.2)-(3.3), we seek a normal mode solution     w(x, y, z, t) θ(x, y, z, t)     =     F 11 (z) G 11 (z)     f (x, y) exp (st) , (3.4)
with f (x, y) = exp (ik x x + ik y y), k = (k x , k y , 0) the horizontal wavenumber and s = s r + is i a complex number. This leads to the differential equations

s P r -1 D 2 -k 2 F 11 = -k 2 RaG 11 + D 2 -k 2 2 F 11 , (3.5) s G 11 = F 11 + (D 2 -k 2 )G 11 , (3.6)
with D the derivative with respect to z and k the norm of the vector k. The boundary conditions are

F 11 = 0 , DF 11 -L s D 2 F 11 = 0 , G 11 = 0 at z = 0, (3.7a)
F 11 = 0 , DF 11 + L s D 2 F 11 = 0 , G 11 = 0 at z = 1. (3.7b)
It is easy to show that the principle of exchange of stability still holds, i.e. It can be written

s M • X 11 = L • X 11 . (3.8)
Since any multiple of the eigenvector X 11 is also a solution of (3.8), and for symmetry reasons, X 11 can be normalized such that

G 11 (z = 1/2) = 1. (3.9)
A spectral Chebyshev method is used. The eigenfunctions F 11 and G 11 are expanded in terms of the Chebyshev polynomials, T j , and G 11 (z) are real-valued functions. The critical mode is such that w and θ are even with respect to the midplane reflection symmetry. Hence, according to (2.18), u and v are odd. For SFBC, the critical mode is obtained analytically:

F 11 (z) = N j=0 a j T j (2 z -1) , G 11 (z) = N j=0 b j T j (2 z -1) . ( 3
w = 3π 2 2 sin (πz) f (x, y) , θ = sin(πz)f (x, y), (3.11) u = 3π cos(πz) ∂f ∂x , v = 3π cos(πz) ∂f ∂y .
(3.12)

Adjoint eigenvalue problem: Adjoint mode

For vectors fields f and g, one defines a scalar product by To the direct problem (3.8) corresponds the adjoint problem

f , g = 1 0 f • g dz . ( 3
s M + • X ad = L + • X ad with X ad = (F ad , G ad ) , (3.14)
where the adjoint operators M + and L + are defined by

X ad , M • X = M + • X ad , X , X ad , L • X = L + • X ad , X , (3.15)
where X fulfills the 'linear' boundary conditions (3.7). By integrating by part we get the linear adjoint problem and the corresponding boundary conditions with

s P r -1 D 2 -k 2 F ad = D 2 -k 2 2 F ad + G ad , (3.16) s G ad = -k 2 Ra F ad + D 2 -k 2 G ad , (3.17) 0 0.2 0.4 0.6 0.8 1 -5 -4 -3 -2 -1 0 x 10 -3 z F ad (1) (2) (3) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 z G ad (1) (3)
F ad = 0 , DF ad -L s D 2 F ad = 0 , G ad = 0 at z = 0, (3.18a) F ad = 0 , DF ad + L s D 2 F ad = 0 , G ad = 0 at z = 1. (3.18b)
The solution of these equations is obtained using the same method as for the direct problem. Similarly, the normalization adopted for the adjoint mode is

G ad (z = 1/2) = 1. (3.19)
At Ra = Ra c , the so-called adjoint critical mode does not depend on the Prandtl number.

It is displayed in figure 3 for three values of the slip parameter. For SFBC, the critical adjoint mode is given by

w + = - 4 9π 4 sin (πz) f (x, y) , θ + = sin(πz)f (x, y) , (3.20) u + = - 8 9π 5 cos(πz) ∂f ∂x , v + = - 8 9π 5 cos(πz) ∂f ∂y . (3.21)
Note that F ad is three order of magnitude smaller than G ad . This indicates that the system is more receptive to thermal perturbations than to velocity perturbations. (1) P r = 1; (2) P r = 1.5; (3) P r = 2; (4) P r = 3; (5) P r = 5; (6) P r = 10; (7) P r = 100.

Characteristic time of the instability

In slightly supercritical conditions, the growth rate s can be approximated using Taylor expansion,

s = ǫ τ 0 + O(ǫ 2 ) with ǫ = Ra -Ra c Ra c . (3.22)
The determination of the characteristic time τ 0 of the instability follows the methodology described in [START_REF] Cross | Derivation of the amplitude equation at the Rayleigh-Bénard instability[END_REF]. The details are given in Appendix A. Figure 4 shows the variation of τ 0 as a function of L s for different values of P r. It increases from

τ 0 = 1 + 1.9544 P r 38.4429 P r (3.23)
for NSBC [START_REF] Cross | Derivation of the amplitude equation at the Rayleigh-Bénard instability[END_REF][START_REF] Segel | Distant side-walls cause slow amplitude modulation of cellular convection[END_REF][START_REF] Daniels | Nonlinear convection in a rigid channel uniformly heated from below[END_REF] to

τ 0 = 2 3 π 2 1 + P r P r (3.24)
for SFBC [START_REF] Newell | Finite bandwidth, finite amplitude convection[END_REF]. The effect of L s is all the more significant when the Prandtl number is low.

Weakly nonlinear stability analysis: Formulation and procedure

For given boundary conditions, the critical Rayleigh number for the onset of convection, determined from the linear stability analysis, depends only on the norm k c of the wavevector. Because of the isotropy of the extended horizontal plane, the direction of the wavevector is arbitrary. In addition, any linear combination of modes

c p exp (i k p .r) (F 11 (z), G 11 (z))
, where r = (x, y) and k p = (k px , k py ), |k p | = k c and c p 's are constant coefficients, is a solution of the linear problem, i.e. there is also a pattern degeneracy. Hereafter, we study the existence and stability of simple regular patterns, rolls made of a single wavevector ±k c r, squares made of two pairs of wavevectors at right angles, or hexagons made of three pairs of wavevectors at 2π/3 angles apart. A weakly nonlinear analysis using the amplitude expansion method is adopted as a first approach.

At leading order, one writes

w(x, y, z, t) = f (x, y, t) F 11 (z) + c.c. , θ(x, y, z, t) = f (x, y, t) G 11 (z) + c.c. , (4.1)
with f (x, y, t) = Configuration with N = 1 corresponds to rolls, N = 2 to squares and N = 3 to hexagons.

Omitting the temporal dependence, the planform function

f (x, y) = 2A cos(k c x) for rolls, f (x, y) = 2 [A 1 cos(k c x) + A 2 cos(k c y)] for squares, f (x, y) = 2 A 1 cos(k c x) + A 2 cos(- 1 2 k c x + √ 3 2 k c y) + A 3 cos(- 1 2 k c x - √ 3 2 k c y)
for hexagons.

The weakly nonlinear analysis is applied to each of these three patterns. To avoid overloading the article, the details of the method are presented only for rolls.

Principles of the amplitude expansion method: Case of rolls

The amplitude expansion method was introduced by [START_REF] Stuart | On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. part 1. the basic behaviour in plane Poiseuille flow[END_REF] and [START_REF] Watson | On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. part 1. the development of a solution for plane Poiseuille flow and for plane Couette flow[END_REF] and later modified by [START_REF] Reynolds | Finite amplitude instability of parallel shear flows[END_REF]. It was surveyed by [START_REF] Herbert | Nonlinear stability of parallel flows by high order amplitude expansions[END_REF][START_REF] Herbert | On perturbation methods in nonlinear stability theory[END_REF].

The amplitude expansion method was shown to be equivalent to the center manifold reduction, which is another technique for deriving the Landau equation [START_REF] Fujimura | Methods of centre manifold and mutiple scales in the theory of weakly nonlinear stability of fluid motions[END_REF][START_REF] Fujimura | Centre manifold reduction and the stuart-landau equation for fluid motions[END_REF]. For a roll pattern, the problem is two-dimensional: ∂/∂y = 0, v = 0 and ζ = 0.

The interaction of the fundamental with itself and with its complex conjugate generates higher harmonics and a modification of the basic state. It is natural to write the nonlinear perturbation as the Fourier series

[u(x, z, t), w(x, z, t), θ(x, z, t)] = +∞ n=-∞ [u n (z, t), w n (z, t), θ n (z, t)] E n , (4.2) with E n = e inkcx and ink c u n = -Dw n . (4.3)
The growth of the disturbance or transitory evolutions are taken into account by the temporal evolution of the Fourier coefficients u n , w n and θ n . Because w and θ are real, we have w -n = w * n and θ -n = θ * n , where the star denotes complex conjugation. Substituting (4.2) and (4.3) into (2.16) and (2.17) and separating out the coefficients of like exponentials, we obtain an infinite set of nonlinear partial differential equations for the Fourier components w n and θ n :

1 P r ∂ ∂t S n w n = S n 2 w n -n 2 k 2 c Ra θ n + [N I w ] E n + [N V ] E n , (4.4) ∂ ∂t θ n = w n + S n θ n + [N I θ ] E n , (4.5)
with

S n = D 2 -n 2 k 2 c , (4.6) [N I w ] E n , [N I θ ] E n and [N V ]
E n the coefficients of E n in the nonlinear inertial and viscous terms respectively. The nonlinearity and coupling of the infinite set of partial differential equations (4.4), (4.5) make its solution difficult. However, if the amplitude A(t) of the fundamental mode (w 1 , θ 1 ) is small, the Fourier components w n and θ n can be sought using a perturbation method expanding around the solution of the linear problem:

(w 1 (z, t), θ 1 (z, t)) = A(t) (F 1 (z, t), G 1 (z, t)) . (4.7)
The amplitude of the perturbation is defined by setting

A(t) = θ 1 (z = 1/2, t). (4.8)
It is clear that, if the fundamental mode is O(A) at leading order, then the leading term of (w 2 , θ 2 ) is O(A 2 ), due to the nonlinear forcing terms. The same reasoning applied for higher harmonics indicates that

(w n (z, t), θ n (z, t)) = A n (t) (F n (z, t), G n (z, t)) if n > 0 , (4.9)
and (w 0 (z, t), θ 0 (z, t)) = A 2 (t) (F 0 (z, t), G 0 (z, t)) .

(4.10) Substituting (4.9) and (4.10) into (4.4) and (4.5) and equating like powers of A(t), the following set of equations is obtained for F n and G n :

1 P r n g + ∂ ∂t S n F n = S n 2 F n -n 2 k 2 c Ra G n + [N I w ] E n A n + [N V ] E n A n , (4.11) n g + ∂ ∂t G n = F n + S n G n + [N I θ ] E n A n , (4.12)
where the subscript E n A n means the coefficient of E n A n and g = 1/A dA/dt. The time evolution of the amplitude A(t) is given by the Stuart-Landau equation

g = 1 A dA dt = +∞ m=0 g m A 2m , (4.13)
where in particular g 0 = s, the linear eigenvalue. Since

F n (G n ) is O(1) or O(A 2 ) as
A → 0, the nonlinearities generate terms in ascending powers of A 2 . Hence, F n and G n are expanded as follows:

F n (z, t) = +∞ m=0 F n,2m+n (z)A 2m , G n (z, t) = +∞ m=0 G n,2m+n (z)A 2m . (4.14)
Substitution of (4.14) into (4.11) and (4.12) yields the differential equations for F n,2m+n

and G n,2m+n ,

L 1nm F n,2m+n + L 2nm G n,2m+n = [N I w ] E n A 2m+n + [N V ] E n A 2m+n -(1/P r) m j=1 (2(m -j) + n) g j S n F n,2(m-j)+n , (4.15) -F n,2m+n + L 3nm G n,2m+n = [N I θ ] E n A 2m+n - m j=1 (2(m -j) + n) g j G n,2(m-j)+n , (4.16)
with

L 1nm = 1 P r (2m + n) sS n -S 2 n , L 2nm = n 2 k 2 c Ra , L 3nm = (2m + n) s -S n .

Solution procedure

The A 7 in amplitude for rolls and order A 5 for squares and hexagons. Note that, according to (2.10), the influence of the non-linearity of the rheological behavior appears only at orders A 3 , A 5 and A 7 .

Results and discussion

This section is divided into five subsections. The first three subsections are devoted to the modification of the base state, the generation of the first harmonic as well as coupling modes for squares and hexagons, induced by the interaction of the fundamental mode with itself and its complex conjugate. These elements are necessary for the calculation of the first Landau constant that determines the nature of the bifurcation. This is done in the fourth subsection. The fifth subsection deals with the competition between different patterns of convection.

Modification of the conductive temperature profile

The interaction of the fundamental (4.1) with itself through the nonlinear quadratic terms produces a correction of the basic state: As for the linear problem, equation (5.2) with the boundary conditions (5.3) is solved numerically using a spectral Chebyshev collocation method. Figure 5 The boundary conditions on F 22 and G 22 are identical to the ones on F 11 and G 11 , equation (3.7). The results are shown in figure 6. For SFBC, we have

F 22 = G 22 = 0
(5.7) in agreement with the numerical results obtained at large L s .

Quadratic interaction between Fourier modes with different wavevectors

The quadratic interaction of the fundamental mode with itself generates the first harmonic mode described in the above section, but also modes resulting from the interaction between modes with wavevectors k p and k q (p = q). In the present study, the wavevectors lie on either a square or hexagonal lattice in the wavevector plane.

Square lattice

With k 1 = k c e x and k 2 = k c e y , modes generated at order A 1 A 2 are

A 1 A 2     F A1A2 (z) G A1A2 (z)     exp [i (k 1 + k 2 ) • r] + c.c. , (5.8)
and

A 1 A 2     FA1A2 (z) ḠA1A2 (z)     exp [i (k 1 -k 2 ) • r] + c.c.
(5.9)

0 0.2 0.4 0.6 0.8 1 -5 0 5 F A 1 A 2 (1) (2)
(3) Since F 11 and G 11 are real, (F A1A2 , G A1A2 ) and FA1A2 , ḠA1A2 are identical and real.

z 0 0.2 0.4 0.6 0.8 1 -1.2 -0.8 -0.4 0 0.4 0.8 1.2 G A 1 A 2 z (3) (1)
They satisfy

D 2 -2k 2 c 2 -2(s/P r) (D 2 -2k 2 c ) F A1A2 -2k 2 c Ra G A1A2 = (2/P r) F 11 D 3 F 11 + (DF 11 ) D 2 F 11 -(4/P r) k 2 c F 11 (DF 11 ] ,
(5.10)

D 2 -2k 2 c -2s G A1A2 + F A1A2 = 2F 11 (DG 11 ) .
(5.11)

Boundary conditions on (F A1A2 , G A1A2 ) are the same as the ones on (F 11 , G 11 ). The functions F A1A2 and G A1A2 are shown in figure 7 for different values of L s . The coupling between the modes exp(ik c x) and exp(ik c y) is significant and it is more important in the situation of wall adhesion than for SFBC. In this latter case, at the critical conditions,

F A1A2 (z) = - 27 100 π 2 π P r + 1 sin (2πz) ,
(5.12)

G A1A2 = -π 27 + 150P r 473P r sin (2πz) .
(5.13)

Note that, at order A 1 A 2 , the nonlinear inertial terms in (2.15) cancel, and the vertical vorticity ζ obeys a diffusion equation: it thus can be ignored.

Hexagonal lattice

With

k 1 = k c e x , k 2 = k c - 1 2 e x + √ 3 2 e y and k 3 = k c - 1 2 e x - √
3 2 e y , modes generated at order A p A q with 1 p 3, 1 q 3 and p = q are

A p A q     F ApAq (z) G ApAq (z)     exp [i (k p + k q ) • r] + c.c.,
(5.14) and

A p A q     FApAq (z) ḠApAq (z)     exp [i (k p -k q ) • r] + c.c.
(5.15)

The functions F ApAq and G ApAq satisfy

D 2 -k 2 c 2 -2(s/P r) D 2 -k 2 c F ApAq -k 2 c Ra G ApAq = (1/P r) F 11 D 3 F 11 + 2 (DF 11 ) D 2 F 11 -3k 2 c F 11 (DF 11 ) , (5.16) D 2 -k 2 c -2s G ApAq + F ApAq = 2F 11 (DG 11 ) + G 11 (DF 11 ) .
(5.17)

The functions FApAq and ḠApAq satisfy (5.20)

D 2 -3k 2 c 2 -(2s/P r) D 2 -3k 2 c FApAq -3 k 2 c Ra ḠApAq = (5.18) (3/P r) F 11 D 3 F 11 -k 2 c F 11 (DF 11 ) , D 2 -3k 2 c -
0 0.2 0.4 0.6 0.8 1 -10 -5 0 5 10 z F A p A q (1) (2) (3) 0 0.2 0.4 0.6 0.8 1 -3 -2 -1 0 1 2 3 z G A p A q (1) (3) (a) (b) 0 0.2 0.4 0.6 0.8 1 -4 -2 0 2 4 z F _ A P A q (1) (2)
(3) .

0 0.2 0.4 0.6 0.8 1 -0.5 -0.25 0 0.25 0.5 z G _ A p A q (1) (2) (3) (c) (d)
As for squares, the amplitude of modes arising from the quadratic coupling between modes with vector k p and k q is more important than that of the first harmonic.

Modification of the fundamental at cubic order: nature of the bifurcation

The nonlinear interactions between the fundamental, its first harmonic, the modification of the conductive temperature profile and modes generated through different couplings, lead to a cubic correction O(A 3 p ) to the fundamental mode. The first Landau coefficient accounts for the feedback of these nonlinear interactions on the fundamental mode. It is determined for the three convective patterns and the nature of the bifurcation is deduced.

The nonlinearity of the rheological law intervenes through the term (µ -1) γ in (2.15) and (2.16). At cubic order, because of (2.10), it reduces to:

(µ -1) γ = -α Γ γ.

(5.21)

Bifurcation to rolls

The modification of the fundamental at order A 3 is governed by (4.15) and (4.16) with

m = n = 1, i.e., S 2 1 -3 s P r S 1 F 13 -k 2 c Ra G 13 = g 1 P r S 1 F 11 -[N I w ] E 1 A 3 -[N V ] E 1 A 3 , (5.22) (S 1 -3s) G 13 + F 13 = g 1 G 11 -[N I θ ] E 1 A 3 .
(5.23)

The boundary conditions are

F 13 = 0 , DF 13 -L s D 2 F 13 = αL s ∂ ∂x Γ ∂u ∂z E 1 A 3
, G 13 = 0 at z = 0, (5.24a)

F 13 = 0 , DF 13 + L s D 2 F 13 = -αL s ∂ ∂x Γ ∂u ∂z E 1 A 3 , G 13 = 0 at z = 1. (5.24b)
Generally, these boundary conditions are inhomogeneous. They are homogeneous only for SF and NSBC, since F 13 = DF 13 = 0 at z = 0, 1 for NSBC, (5.25)

F 13 = D 2 F 13 = 0 at z = 0, 1 for SFBC.

(5.26)

The system (5.22), (5.23) can be written

L • X 13 = g 1 M • X 11 -N I -N V + 3 s M • X 13 with X 13 = (F 13 , G 13 ) . (5.27)
The nature of the bifurcation is determined at the critical conditions, i.e., s = 0.

In order to calculate the first Landau coefficient g 1 from the Fredholm solvability condition, i.e., orthogonality of the inhomogeneous part of (5.27) to the null-space of the adjoint operator of L, we decompose X 13 into a homogeneous X 13H and inhomogeneous X 13N H parts:

X 13 = X 13H + X 13N H .
(5.28)

X 13H satisfies L • X 13H = g 1H M • X 11 -N I -N V ,
(5.29) with homogeneous boundary conditions, i.e. F 13H = DF 13H -L s D 2 F 13H = 0 at z = 0 and similarly at z = 1. By applying the solvability condition to (5.29), we obtain

g 1H = N I + N V , Xad , Xad = X ad M • X 11 , X ad .
(5.30)

X 13H is then determined by solving (5.29). X 13N H is a correction term that accounts for the non-homogeneity of the boundary conditions. Substituting (5.28) into (5.27), we obtain at the critical conditions

L • X 13H = g 1 M • X 11 -N I -N V -L • X 13N H .
(5.31)

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

g 1 = g 1H + L • X 13N H , Xad .
(5.32)

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

At the start, X 13N H is assumed to be identically zero in (5.32). A first approximation to g 1 is then obtained: g

(1)

1 = g 1H . This is put into (5.27). which is solved, at the critical conditions, with non-homogeneous boundary conditions, to obtain a first approximation of X 13 . Using (5.28), a first approximation of X 13N H is deduced. Then X 13N H is put into (5.32). This process is repeated until it converges to a desired level of accuracy. Note that (5.27) and (5.29) are solved with an additional condition (5.33) as suggested by [START_REF] Herbert | Nonlinear stability of parallel flows by high order 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, X 13 is defined up to an arbitrary multiple of the solution X 11 of the linear problem (3.8). Finally, g 1 can be written as the sum of three contributions. The first one arises from the nonlinear inertial terms N I, the second one from the nonlinear viscous terms N V and the third one from the inhomogeneity of the boundary conditions,

X 13 = X 13H = 0 at z = 1/2,
g 1 = g I 1 + g V 1 + g N H 1 (5.34) with g I 1 = N I, Xad , g V 1 = N V , Xad , g N H 1 = L • X 13N H , Xad .
(5.35) Using (5.21), the contribution of the nonlinear viscous term g V 1 can be written as

g V 1 = -αg N N 1 , (5.36)
with α defined by (2.11) and g N N 1 that does not depend on the rheological parameters.

Similarly, it can be shown that X 13N H = -α X13NH , where X13NH does not depend on the rheological parameters. Hence

g N H 1 = -αg BC 1 with g BC 1 = L • X13NH , Xad (5.37)
Note that, for NS and SFBC,

X 13N H = 0 thus g BC 1 = 0 .
(5.38)

The first Landau constant g 1 as well as the different contributions g I 1 , g V 1 and g N H In figure 9(a), g I 1 in plotted as a function of L s . As expected, g I 1 is negative, i.e., the bifurcation is supercritical for a Newtonian fluid. The absolute value of g I 1 decreases with increasing L s , i.e., slip promotes the development of the convection. For SFBC, an analytical expression of g I 1 is obtained,

g I 1 = - 9 π 4 8 P r 1 + P r .
(5.39)

The contribution of the nonlinear inertial terms gathers the feedback from thermal advection terms and velocity terms:

g I 1 = N I w (F 22 |F 11 ) , Fad + N I θ (G 22 |F 11 ) , Gad + N I θ (F 22 |G 11 ) , Gad + N I θ (G 02 |F 11 ) , Gad .
(5.40)

For P r 1, the analysis of the different terms shows that g I 1 is dominated by the nonlinear thermal convection terms, more precisely the term involving G 02 , the modification of the conductive temperature. The contributions of the other nonlinear inertial term is practically negligible.

In figure 9(b), we plot g N N 1 as a function of L s for different values of P r. For SFBC, an analytical expression is obtained,

g N N 1 = - 1803 π 8 64 P r 1 + P r .
(5.41)

As can be observed, g N N 1 is negative, and

g V 1 = -αg N N 1 > 0.
Therefore, shear-thinning effects promote a subcritical bifurcation, which is understandable since the viscosity, which damps convection, is reduced.

Concerning g N H

1

, its absolute value does not exceed g 1 /100. The maximum value is reached at L s ≈ 0.1.

In figure 10 we plot g 1 as a function of L s at P r = 10 and different values of α. For low shear thinning effects, g 1 < 0 and the bifurcation is supercritical, while for sufficiently high shear-thinning effects, g 1 > 0 and the bifurcation is subcritical. Using (5.34) and

(5.36), the critical degree of shear-thinning α c above which the bifurcation changes from supercritical to subcritical is given by

α c = g I 1 g N N 1 + g BC 1 .
(5.42)

Figure 11 shows the variation of α c as a function of P r for NS and SFBC. For this latter (5.43)

This result agrees with [START_REF] Balmforth | Weakly nonlinear viscoplastic convection[END_REF] but disagrees with the value α c = 14×10 -4 found by [START_REF] Albaalbaki | Pattern selection in the thermal convection of nonnewtonian fluids[END_REF]. For NSBC, and P r 1, α c = 2.15×10 -4 .

The plot of α c vs L s for P r = 10 is visible in figure 14, which also displays some results

for two other patterns.

Bifurcation to squares

As compared to rolls, for square patterns, the modification of the fundamental comprises an additional term, which arises from the interaction of modes with wavevectors

k 1 = k c e x and k 2 = k c e y .
Hence, the modification of the fundamental

w = F 13 (z) A 3 1 e ikcx + A 3 2 e ikcy + c.c. + F13 (z) A 2 2 A 1 e ikcx + A 2 1 A 2 e ikcy + c.c. ,
and similarly for θ by replacing F by G. The amplitude equations for A 1 and A 2 are obtained using symmetries introduced by the square lattice. Employing crystallographic terminology [START_REF] Mckenzie | The symmetry of convective transitions in space and time[END_REF][START_REF] Silber | Pattern selection in ferrofluids[END_REF][START_REF] Golubitsky | Symmetries and pattern selection in rayleigh-bénard convection[END_REF], these are the symmetries of a square D 4 in addition to the two-torus T 2 of translation in the two horizontal directions. Requiring equivariance with respect to the group D 4 × T 2 leads to the amplitude equations

dA 1 dt = s A 1 + g 1 A 2 1 + λ 1 A 2 2 A 1 , (5.44) dA 2 dt = s A 2 + g 1 A 2 2 + λ 1 A 2 1 A 2 .
(5.45)

Equivariance under the midplane reflection implies that the self-saturation coefficient g 1 and the coupling coefficient λ 1 are real.

Equations satisfied by F13 (z) and G13 (z) are similar to that satisfied by F 13 (z) and G 13 (z), i.e. (5.22) and (5.23) respectively, by replacing g 1 by λ 1 and evaluating the nonlinear forcing terms that are factor of A 2 2 A 1 exp (ik c x). At criticality, s = 0, the coupling coefficient λ 1 is obtained with a solvability condition. As for g 1 , it can be written

λ 1 = λ I 1 -αλ N N 1 -αλ BC 1 .
(5.46)

For SFBC,

λ I 1 = - 9π 4 120 + 72P r + 673P r 2 3784 P r (1 + P r) , λ N N 1 = - 921π 8 P r 32(1 + P r) , λ BC 1 = 0. (5.47)
The critical value of α above which the bifurcation to squares becomes subcritical,

α c = g I 1 + λ I 1 g N N 1 + λ N N 1 + g BC 1 + λ BC 1 .
(5.48)

Figure 12 shows the variation of α c as a function of P r for NS and SFBC. In this latter case, α c = 16 20 + 12 P r + 191 P r 2 63855 π 4 P r 2 .

(5.49)

Bifurcation to hexagons

For hexagonal patterns, the modification of the fundamental modes reads

w = F 13 (z) A 3 1 e ik1•r + A 3 2 e ik2•r + A 3 3 e ik3•r + c.c. (5.50) + F13 (z) A 2 2 + A 2 3 A 1 e ik1•r + A 2 3 + A 2 1 A 2 e ik2•r + A 2 1 + A 2 2 A 3 e ik3•r + c.c. ,
and similarly for θ by replacing F by G. The amplitude equations for A 1 , A 2 and A 3 are obtained using symmetries introduced by the hexagonal lattice. The symmetry group that leaves the hexagonal lattice invariant is D 6 × T 2 ; D 6 represents the rotational (2π/3) and reflection symmetries. The amplitude equations are [START_REF] Golubitsky | Symmetries and pattern selection in rayleigh-bénard convection[END_REF]:

dA 1 dt = sA 1 + g 1 A 2 1 + δ 1 A 2 2 + A 2 3 A 1 .
(5.51)

Equations for dA 2 /dt, dA 3 /dt are obtained by cyclic permutation of A 1 , A 2 and A 3 . The quadratic term proportional to A 2 A 3 in (5.51) is equivariant with respect to the group D 6 × T 2 , but vanishes here because of the up and down reflection symmetry. Equations satisfied by F13 (z) and G13 (z) are similar to that satisfied by F 13 (z) and G 13 (z), i.e. (5.22) and (5.23) respectively, by replacing g 1 by δ 1 and evaluating the nonlinear forcing terms that are factor of A 2 2 A 1 exp (ik c x). The determination of δ 1 the same procedure as for g 1 and λ 1 . As previously, the coupling coefficient δ 1 can be written

δ 1 = δ I 1 -αδ N N 1 -αδ BC 1 .
(5.52)

For SFBC

δ I 1 = - 9π 4 1728 + 1113P r + 1273P r 2 65000 P r (1 + P r) , δ N N 1 = - 2175 π 8 P r 64 (1 + P r) , δ BC 1 = 0. (5.53)
The critical value of the degree of shear-thinning α above which the bifurcation to hexagons becomes subcritical The variation of α c as function of P r is plotted in figure 13 for NS and SFBC. In this latter case,

α c = g I 1 + 2δ I 1 g N N 1 + 2δ N N 1 + g N H 1 + 2δ BC 1 . ( 5
α c = 24 3456 + 2226P r + 33599P r 2 16664375π 4 P r 2 .
(5.55)

Finally, the figure 14 shows α c vs L s for rolls, squares and hexagons, at P r = 10.

From now on, for the sake of simplicity, in order to not overload this article, we restrict ourselves to the two limiting cases L s = 0, NSBC, or L s = +∞, SFBC.

Remark

M. Bouteraa, C. Nouar, C. Metivier, E. Plaut and A. Kalck In this study, the midplane reflection symmetry is not broken (see Appendix D). The amplitude equation at cubic order does not contain a quadratic term. Therefore, the bifurcation is of pitchfork type. It is either supercritical or subcritical depending on the sign of the cubic term. In the presence of non-Boussinesq effects, for instance, in fluids with thermodependent viscosity, the midplane reflection symmetry is broken and a quadratic term appears in the amplitude equation. In general, the presence of the quadratic term can lead to a bifurcation of transcritical type with an hexagonal planform.

Pattern selection

In this section, we investigate the competition between rolls and squares, and between rolls and hexagons, when we depart from the critical conditions. The calculation proceeds in two stages. First, the possible steady states solutions of the amplitude equations are determined. Then, their linear stability is considered. However, this study requires the prior evaluation of the saturation and coupling coefficients in the amplitude equations, outside the critical conditions.

Saturation and coupling coefficients outside the critical conditions

Departing sufficiently from the critical conditions, we can no longer assume s = 0 in the equations that describe the modification of the fundamental at cubic order. Hence, these equations become unconditionally solvable. To calculate g 1 , λ 1 or δ 1 , an iterative process is used as suggested by [START_REF] Sen | On the stability of plane Poiseuille flow to finiteamplitude disturbances, considering the higher-order Landau coefficients[END_REF]. Details of the procedure are given in Appendix B.

Figure 15 shows the variation of the ratios g I 1 (ǫ)/g I 1 (ǫ = 0) and

g N N 1 (ǫ)/g N N 1 (ǫ = 0) vs
the reduced Rayleigh number ǫ = Ra/Ra c -1. The contribution to g 1 of the non Newtonian viscous terms increases more rapidly than that of the inertial terms. At ǫ = 0, the

saturation coefficient g 1 = g I 1 -αg N N 1 changes sign for α = α c . It is clear from figure 15 0 0.1 0.2 0.3 1 1.2 1.4 1.6 1.8 2 ε g 1 I /g 1 I (ε = 0) , g 1 NN /g 1 NN (ε = 0) (1) (2) (3) (4)
(1) 0 0.05 0.1 0.15 0.2 0.25 0.3 10 that, for ǫ = 0, g 1 changes sign for α < α c . This result will be used when determining the domain of existence of the stationary solutions of the amplitude equations at cubic order.

0 10 1 ε δ 1 / g 1 (1) (4) (5) (6) (7) (8)
Ratios between the coupling and saturation coefficients play a fundamental role in the pattern selection. They are shown in figure 16 vs the reduced Rayleigh number, for different values of α, P r = 10 and NSBC; similar curves are obtained for SFBC. For α 0.9 α c , one has: (i) λ 1 , δ 1 and g 1 are of the same, negative, sign. (ii)

|λ 1 | > |g 1 | and |δ 1 | > |g 1
| and (iii) λ 1 /g 1 and δ 1 /g 1 increase with increasing shear-thinning effects.

Pattern selection on a square lattice

The amplitude equations satisfying the symmetry requirements of a square lattice are (5.44) and (5.45). The coefficients s, g 1 and λ 1 in these equations depend on α and on the reduced Rayleigh number ǫ. The stationary solutions are obtained by setting

f i (A 1 , A 2 ) = 0
, where f i is the right hand side of the amplitude equations. Their stability is determined by the sign of the eigenvalues χ i of the Jacobian matrix J ij = ∂f i ∂A j evaluated at the steady states. In the following, the stability of the stationary solutions is examined in details.

(i) Conduction state, A 1 = A 2 = 0. The stability eigenvalues associated to this state are

χ 1 = χ 2 = s.
The conductive state is stable if ǫ < 0 and unstable if ǫ > 0.

(ii) Steady convection with rolls, A 1 = -s/g 1 , A 2 = 0 or A 1 = 0, A 2 = -s/g 1 . The associated eigenvalues are: χ 1 = -2s and χ 2 = s g 1 -λ 1 g 1 . In the supercritical regime, according to figures 10 and 16(a), g 1 and λ 1 are negative and |λ 1 | > |g 1 |. Thus stable rolls appear at a supercritical bifurcation from ǫ = 0. In the neighborhood of the points (ǫ, α) where g 1 changes sign, the stability of rolls cannot be determined and a higher order analysis is required. In the subcritical regime, roll solutions exist but are unstable.

(iii) Steady convection with square patterns, A 1 = A 2 = -s/(g 1 + λ 1 ). The associated eigenvalues are χ 1 = -2s and χ 2 = 2 s (λ 1 -g 1 ) λ 1 + g 1 . In the supercritical regime, g 1 and λ 1 are negative and |λ 1 | > |g 1 |, therefore χ 2 > 0 and squares are unstable. As explained by [START_REF] Fauve | Pattern forming instabilities[END_REF], when |λ 1 | > |g 1 |, the interaction between the two sets of rolls is too strong and one of the two sets of rolls nonlinearly damps out the other. As indicated previously, this interaction is more stronger with increasing shear-thinning effects (figure 16a). At cubic order, the stability of the squares cannot be determined in the region where g 1 + λ 1 changes sign. In the subcritical regime, square solutions exist but are unstable.

Pattern selection on a hexagonal lattice

The system of amplitude equations satisfying the symmetry requirements of a hexagonal lattice is given by (5.51). The stationary solutions are determined and their linear stability is investigated.

(i) Steady convection with rolls, A 1 = -s/g 1 , A 2 = A 3 = 0. The eigenvalues associated to this state are χ 1 = -2 s and χ 2 = χ 3 = s g 1 -δ 1 g 1 . In the supercritical regime, according to figures 10 and 16(b), δ 1 and g 1 are negative and |δ 1 | > |g 1 |, thus χ 3 < 0. As on the square lattice, stable rolls appear supercritically from ǫ = 0. This analysis is not valid in the vicinity of the points (ǫ, α) where g 1 changes sign.

(ii) Steady convection with hexagons, A 1 = A 2 = A 3 = -s/ (g 1 + 2δ 1 ). The associated eigenvalues are χ 1 = -2s, χ 2 = χ 3 = 2s δ 1 -g 1 g 1 + 2δ 1 . In the supercritical regime, g 1 and δ 1 are negative and |δ 1 | > |g 1 |, thus χ 3 > 0 and hexagons are unstable. This analysis is not valid in the vicinity of the points (α, ǫ) where g 1 + 2δ 1 changes sign. In the subcritical regime, hexagon solutions exist but are unstable.

Comparison with Albaalbaki & Khayat (2011)

Finally, our results show that only rolls are stable near critical conditions. Furthermore, shear-thinning effects reinforce convection in the form of rolls. This contradicts [START_REF] Albaalbaki | Pattern selection in the thermal convection of nonnewtonian fluids[END_REF] findings, where stable squares and hexagons were found near critical conditions. For an hexagonal lattice, [START_REF] Albaalbaki | Pattern selection in the thermal convection of nonnewtonian fluids[END_REF] wrote the amplitude equations up to O(A 4 ), where we denote O(A) the common order of magnitude of the amplitudes A 1 , A 2 , A 3 of all modes. They indicated that the terms of O(A 2 ) vanish, but that terms of O(A 4 ) -the terms proportional to a 5 , a 6 , a 7 in their equation (3.14) or to a 5 + a 6 + 2a 7 in their equation (5.6) -play an important role, because "non-Newtonian effects are symmetry breaking" (their §5.2 p. 527).

In fact, non-Newtonian effects do not break the midplane reflection symmetry. Therefore all the even terms in the amplitude equations vanish. In the Appendix D, detailed calculations are provided to show that the coefficients a 5 , a 6 , a 7 in equation (3.14) of (Albaalbaki & Khayat 2011) vanish.

Maximum heat transport principle

For Newtonian fluids, [START_REF] Malkus | Finite amplitude cellular convection[END_REF] introduced the maximum heat transport principle: "the only stable solution is the one of maximum heat transport". [START_REF] Schluter | On the stability of steady finite amplitude convection[END_REF] confirmed this principle for slightly supercritical conditions.

Here, we show that the maximum heat transport principle is valid for shear-thinning fluids, i.e., the maximum heat transport is obtained for the only stable solution: rolls.

The heat transfer through the horizontal fluid layer is described by the Nusselt number, N u, the ratio of the total heat flux to the purely conductive heat flux in the absence of fluid flow. It can be calculated either at the lower or upper plate. At the lower plate, we have (5.56) where the overbar denotes the horizontal average over one wavelength, N = 1 corresponds to rolls, N = 2 to squares and N = 3 to hexagons. The unperturbed solution, N u = 1, corresponds to the hydrostatic solution. The second term of N u refers to the convective transfer. As it can be seen in figure 5, (DG 02 ) z=0 < 0. Using the stationary solutions found in the previous subsection, we obtain:

N u = 1 - ∂ θ ∂z z=0 = 1 - N p=1 A 2 p (DG 02 ) z=0 ,
N u r = 1 + s g 1 (DG 02 ) z=0 for rolls, (5.57)

N u s = 1 + 2 s g 1 + λ 1
(DG 02 ) z=0 for squares, (5.58)

N u h = 1 + 3 s g 1 + 2δ 1
(DG 02 ) z=0 for hexagons.

(5.59)

As g 1 , λ 1 , δ 1 are all negative, and

|λ 1 | > |g 1 |, |δ 1 | > |g 1 |
, the maximum heat transfer occurs for rolls.

Solutions at higher order -Range of validity

So far, the results were obtained using amplitude equations truncated at O A 3 . For a significant deviation from the critical conditions, terms of higher order become large and should be taken into account. A weakly nonlinear expansion was carried out up to O A 5 ; for rolls, the expansion was extended to O A 7 . Amplitude equations at O A 5

for squares and hexagons, stationary solutions and the eigenvalues of the associated Jacobian matrices are given in the Appendix C. The analysis of the stability of the various patterns shows once again that, for supercritical conditions, only rolls are stable for NS and SFBC. For subcritical conditions, the solutions for the three patterns are unstable.

The range of validity of these results may be estimated roughly as 0 ǫ ǫ 0 , where ǫ 0 is the reduced Rayleigh number at which the viscosity perturbation, µ -1, reaches 30%.

Using this criterion, the figure 17 shows the curve (α, ǫ 0 ) for rolls calculated at O A 7 .

Below this curve, rolls are stable. With increasing shear-thinning effects, nonlinear terms increases rapidly, reducing ǫ 0 , i.e., the range of validity of the weakly nonlinear theory.

Heat transfer, flow structure and viscosity field in roll solutions

It is found that only rolls are stable near onset. In the present section, information on the heat transfer, the flow structure and viscosity field in rolls are provided. The influence of shear-thinning will be emphasized.

Heat transfer

Using (4.2) combined with (4.3), (4.10) and (4.14), one obtains:

N u = 1 - M m=1
A 2m (DG 0,2m ) z=0 .

(7.1)

Figure 18 shows, for a Newtonian fluid with NSBC, the evolution of the Nusselt number as a function of ǫ, depending on the order of truncation of the computations. Our results are in good agreement with the numerical solution of (2.16) and (2.17) obtained using the spectral code of [START_REF] Plaut | Low-Prandtl-number convection in a rotating cylindrical annulus[END_REF], at least up to ǫ = 0.35, when series are truncated at O A 7 . A good agreement is also observed with the numerical results of [START_REF] Plows | Some numerical results for two-dimensional steady laminar Bénard convection[END_REF].

For indication, the correlation

N u = 1 + ǫ 1 + ǫ 1 0.6992 -0.00472P r -1 + 0.00382P r -2 , (7.2)
proposed by [START_REF] Schluter | On the stability of steady finite amplitude convection[END_REF] using a weakly nonlinear approach is also represented.

The influence of the shear-thinning behavior is illustrated by figure 19 the rates of deformation in the roll solutions between NS and SFBC is, in fact, a complex process, which is continuous but not monotonous as suggested by the figure 14.

Conclusion

We studied the influence of shear-thinning effects on the convection in a horizontal layer of an inelastic shear-thinning fluid using the Carreau model as a typical rheological model. To take into account the possibility of wall slip, Navier's slip law has been used at the top and bottom walls. As the fluid has a finite viscosity at zero shear-rate, the critical Rayleigh number and wavenumber remain unchanged with respect to the Newtonian case. In agreement with [START_REF] Webber | The destabilizing effect of boundary slip on Bénard convection[END_REF] and [START_REF] Kuo | Effects of slip boundary conditions on Rayleigh-Bénard convection[END_REF], it is found that slip boundary condition has a destabilizing effect.

The nature of the bifurcation to rolls, squares and hexagons has been determined us-ing a three-dimensional weakly nonlinear analysis. The bifurcation is supercritical for moderately shear-thinning effects and becomes subcritical for strongly shear-thinning effects. In the case of a bifurcation to steady rolls with stress-free boundary conditions, the critical value α c of the degree of shear-thinning above which the bifurcation becomes subcritical found here confirms that given by [START_REF] Balmforth | Weakly nonlinear viscoplastic convection[END_REF], but contradicts the computations done by [START_REF] Albaalbaki | Pattern selection in the thermal convection of nonnewtonian fluids[END_REF]. In the case of slip with friction, new results were provided, where α c is given as a function of L s for the three convective patterns (figure 14). Except the case of rolls with stress-free boundary conditions, in all other situations, α c depends on P r. Nevertheless, for P r 10, this dependency is no longer significant (figures 11, 12, 13 ).

The analysis of the amplitude equations at cubic order for squares or hexagons has shown that the ratios between the coupling coefficients and the self saturation coefficient are greater than 1. Moreover, these ratios increase with increasing degree of shear-thinning and Rayleigh number. This reflects a strong interaction between sets of rolls that constitute squares or hexagons and this interaction is more stronger with increasing α. As a consequence, only rolls are stable, and this stability is reinforced by shear-thinning effects. These results, which contradict [START_REF] Albaalbaki | Pattern selection in the thermal convection of nonnewtonian fluids[END_REF], as explained in §5.5.4, are consistent with the maximum heat-transfer principle: the Nusselt number for rolls is larger than that for squares or hexagons. The influence of the shear-thinning on the pattern selection was confirmed by considering amplitude equation at the fifth order. A rough estimation of the range of validity of the weakly nonlinear approach used has been proposed. This range of validity decreases with increasing shear-thinning effects (figure 17).

Finally, roll solutions have been computed with an amplitude equation at the seventh order in the supercritical regime. The relevance of a "Parmentier-Rayleigh number", as defined by [START_REF] Parmentier | A study of thermal convection in non-newtonian fluids[END_REF] for power-law fluids has been confirmed for the Carreau fluids studied here, for mildly nonlinear convection above onset (figure 19).

In the present study, the temperature dependence of the viscosity has not been taken into account. The inclusion of temperature-dependence of viscosity will affect the onset of instability, the symmetry properties of the system, the value of α c and the preferred mode of convection. For Newtonian fluids, [START_REF] Busse | Square-pattern convection in fluids with strongly temperaturedependent viscosity[END_REF]; [START_REF] Jenkins | Rolls versus squares in thermal convection of fluids with temperaturedependent viscosity[END_REF] showed that the roll planform is preferred for low ratios r of the viscosities at the top and bottom boundaries. For larger values of r, the square planform is the preferred mode of convection. At still larger values of r, subcritical bifurcation is predicted. The experiments of [START_REF] White | The planforms and the onset of convection with a temperature dependentviscosity[END_REF] confirmed these tendencies. Based on our analysis of the pattern selection, -the stability of rolls are reinforced by the shear-thinning behavior -we think that the limit value of the ratio r, below which the roll planform is the preferred mode of convection, will be larger for shear-thinning fluids. The study of the effects of the temperature-dependent viscosity in this system is under way.

Coming back to a Boussinesq model with a viscosity that does not depend on the temperature, it would also be interesting to study the subcritical regime, and investigate strongly subcritical convection.

Finally, we stress that there are very few experimental data in the literature. Additional experiments are needed, particularly, for shear-thinning fluids for which the bifurcation should be subcritical, according to the present model. coefficients a 5 , a 6 and a 7 in their equation (3.14) vanish. For this purpose, we use SFBC like [START_REF] Albaalbaki | Pattern selection in the thermal convection of nonnewtonian fluids[END_REF], though this result is also fulfilled with the more general boundary conditions (2.13). The coefficients a 5 , a 6 , a With SFBC, equations (5.1) and (5.7) show that a fundamental mode coupled with itself does not generate velocity modes of order A 2 . Therefore, focusing on the calculation of a 5 for the sake of simplicity, the contributions proportional to (A 1 ) pq is even, i.e., Γ [1] (-z) = Γ [1] (z).

t ij = -α Γ [1] γ[2] ij + Γ [1,2] γ[1] ij , Γ [1] = 1 2 γ[1] pq γ[1] pq , Γ [1,2] = 1 2 γ[1]
-For the modes arising from the quadratic interaction between fundamental modes with different wavevectors, according to equations (5.19) and (5. pq is odd, i.e.

Γ [1,2] (z) = -Γ [1,2] (-z). Consequently, γ[2] ij and Γ [1] γ[2] ij have the same symmetry properties as Γ [1,2] Therefore, [N V ] given by equation (D 4) is odd, i.e. The proof that a 6 , a 7 vanish is quite similar, since when one applies the complex conjugate, the functions of z implied in all modes are unchanged.

  s i = 0, when Navier's slip boundary conditions are used. The set of differential equations (3.5)-(3.6) is an eigenvalue problem where s is the eigenvalue and X 11 = (F 11 , G 11 ) the eigenvector.

  .10) By canceling the residual (sM • X 11 -L • X 11 ) at the (N + 1) collocation points (for j = 0, 1, ..., N , one obtains a matrix eigenvalue problem solved using the QZ algorithm with Matlab. The marginal stability curve Ra(k) is obtained by the condition s(Ra, k) = 0. Using 20 Chebyshev polynomials, the first eigenvalue, i.e. that for which the real part is the largest, is calculated with an accuracy of 10 -4 . The minimum of the marginal stability curves gives the critical Rayleigh number Ra c and critical wave number k c . Figure1displays the variation of Ra c and k c as a function of the dimensionless slip parameter. These results are in very good quantitative agreement with those obtained by[START_REF] Webber | The destabilizing effect of boundary slip on Bénard convection[END_REF];[START_REF] Kuo | Effects of slip boundary conditions on Rayleigh-Bénard convection[END_REF].The critical Rayleigh number decreases with increasing slip parameter L s , from 1707.7 (NSBC) to 27 π 4 /4 = 657.5 (SFBC). The slip has therefore a destabilizing effect. The critical wavenumber decreases with increasing L s , from 3.116 (NSBC) to π/ √ 2 = 2.221 (SFBC). Additional properties of the critical mode are given by F 11 and G 11 at the critical conditions. They are displayed in figure2for different values of L s . Here, F 11 (z)

Figure 1 .Figure 2 .

 12 Figure 1. Critical Rayleigh number (a) and critical wavenumber (b) as function of the slip parameter.

Figure 3 .

 3 Figure 3. Adjoint functions at critical conditions, and different values of Ls: (1) Ls = 0 (NSBC); (2) Ls = 0.1 and (3) Ls = 10 4 very close to SFBC.

Figure 4 .

 4 Figure 4. Characteristic time of instability as a function of Ls for different Prandtl numbers.

A

  p (t) exp (i k p • r), |k p | = k c , and A p (t) the amplitude of the perturbation. According to the normalization of the eigenfunctions used in the linear theory, A p (t) represents the amplitude of the thermal perturbation measured at the midplane.

  set of differential equations (4.15), (4.16) is solved sequentially beginning from n = 1 and m = 0. The problem n = 1, m = 0 is the linear problem (3.5), (3.6), which gives the critical point around which the harmonic-amplitude expansion is carried out. The problem n = 0, m = 1 yields the first correction of the conductive temperature profile. The problem n = 2, m = 0 yields the first harmonic of the fundamental mode. The problem n = 1, m = 1 yields the feedback coefficient g 1 of the fundamental mode. More precisely, g 1 is determined using the condition for the solvability of the equation corresponding to the modification of the fundamental mode. The calculations were continued up to order

FFigure 5 .

 5 Figure 5. Modification of the conductive temperature profile at the critical conditions for P r = 10 and different values of Ls: (1) Ls = 0 NSBC; (2) Ls = 0.1 and (3) Ls = 10 4 very close to SFBC.

FFigure 6 .

 6 Figure 6. First harmonic of the fundamental at the critical conditions for P r = 10 and different values of Ls: (1) Ls = 0 NSBC; (2) Ls = 0.1 and (3) Ls = 10 4 very close to SFBC.

Figure 7 .

 7 Figure 7. Modes generated in a square lattice at order A1A2, at the critical conditions, for P r = 10, and different values of Ls: (1) Ls = 0 (NSBC); (2) Ls = 0.1 and (3) Ls = 10 4 very close to SFBC.

Figure 8 .

 8 Figure 8. Modes generated in an hexagonal lattice at the critical conditions for P r = 10: factor of ApAq exp [i (kp + kq) .r] in (a) and (b) and ApAq exp [i (kpkq) .r] in (c) and (d). (1) Ls = 0, (2) Ls = 0.1, (3) Ls = 10 4 .

  critical sets (Ra c , k c , L s , P r) and different values of the degree of shear-thinning α. The integrals in (5.35) are evaluated numerically by means of Clenshaw and Curtis method.

Figure 9 .Figure 10 .

 910 Figure 9. Contribution of the nonlinear inertial terms (a) and nonlinear viscous terms (b) to the first cubic Landau constant, vs the slip parameter Ls for different values of the Prandtl number: (1) P r = 1, (2) P r = 5, (3) P r = 10.

Figure 11 .

 11 Figure 11. Critical value of the degree of shear-thinning αc above which the bifurcation to rolls becomes subcritical, vs Prandtl number.

Figure 12 .

 12 Figure 12. Critical value of the degree of shear-thinning αc above which the bifurcation to squares becomes subcritical, vs Prandtl number.

Figure 13 .

 13 Figure 13. Critical value of the degree of shear-thinning αc above which the bifurcation to hexagons becomes subcritical, vs Prandtl number.

Figure 14 .

 14 Figure 14. For rolls (R), squares (S) and hexagons (H), critical value of the degree of shear-thinning αc, vs the slip parameter Ls for P r = 10.

Figure 15 .

 15 Figure 15. For rolls, ratios g I 1 /g I 1 (ǫ = 0) and g NN 1 /g NN 1 (ǫ = 0) vs the reduced Rayleigh number, ǫ. (1) g I 1 /g I 1 (ǫ = 0) for SFBC ; (2) g I 1 /g I 1 (ǫ = 0) for NSBC ; (3) g NN 1 /g NN 1 (ǫ = 0) for SFBC and

Figure 16 .

 16 Figure 16. For squares and hexagons, ratios (a) λ1/g1 and (b) δ1/g1 vs the reduced Rayleigh number, at P r = 10 and different values of α in the case of NSBC. (1) α = 0, Newtonian fluid,(2) α = 0.1 αc, (3) α = 0.2 αc, (4) α = 0.3 αc, (5) α = 0.4 αc, (6) α = 0.5 αc, (7) α = 0.6 αc and (8) α = 0.7 αc.

Figure 17 .

 17 Figure 17. Curves, (α, ǫ0), that delimit roughly the domain of validity of the weakly nonlinear approach at O A 7 , for rolls.

Figure 18 .

 18 Figure 20(a) shows the viscosity distribution for a Carreau fluid with α = 0.2 α c and ǫ = 0.2, in the case of NSBC. For comparison, the case of SFBC, discussed by Albaalbaki & Khayat (2011) is shown in figure 20(b). With NSBC, the interior of the roll is practically isoviscous with µ ≈ 1. The viscosity is minimal at the wall where the shear-rate is maximal. It is weakly reduced at the four corners of the roll because of the elongational rate γzz = -γxx , displayed in figure 20(d). The computation of the viscous dissipationindicates that this one is weak at the interior of a roll where the viscosity is high. In agree-

Figure 19 .Figure 20 .

 1920 Figure 19. For rolls at P r = 10, nc = 0.5 and NSBC. (a) Nusselt number vs the reduced Rayleigh number for different values of α. (1) α = 0, Newtonian fluid, (2) α = 0.1 αc, (3) α = 0.2 αc, (4) α = 0.3 αc, (5) α = 0.4 αc, (6) α = 0.5 αc, (7) α = 0.6 αc and (8) α = 0.7 αc. (b) Nusselt number vs the reduced Parmentier Rayleigh number ǭ. (continuous line) Newtonian fluid, (o) α = 0.1 αc, (⋄) α = 0.2 αc, (△) α = 0.3 αc, ( ) α = 0.4 αc, (+) α = 0.5 αc, (×) α = 0.6 αc.

  7 of the equation (3.14) of[START_REF] Albaalbaki | Pattern selection in the thermal convection of nonnewtonian fluids[END_REF] are given with our notations bya 5 , a 6 , a 7 = 1 0 [N V ] F ad dz 1 0 M X 11 • X ad dz , (D 1)with F ad ∝ sin (πz) according to our equation (3.20), M , X 11 and X ad defined in §3, N V the nonlinear viscous terms in the equation for the vertical velocity (2.16),N V = [∇ × ∇ × [∇ • (µ -1) γ]] • e z . (D 2) [N V ] designates the contributions in N V proportional to (A 1 ) 2 A 2 A 3 for a 5 , |A 1 | 2 A * 2 A * 3 for a 6 , |A 2 | 2 A * 2 A * 3 + |A 3 | 2 A * 2 A * 3 for a 7 .In (D 2), because of equation (2.10), the nonlinear viscous terms that can give contributions of order A 4 readN V = -α [∇ × ∇ × (∇ • Γ γ)] • e z = -α 2 [∇ × ∇ × [∇ • ( γpq γpq γij e i ⊗ e j )]] • e z .(D 3)

  z ∈ [0, 1] into z ∈ [-1/2, 1/2] by defining z = z -1/2. Let us demonstrate that [N V ] (-z) = -[N V ](z), i.e., that [N V ] is odd.-For the fundamental mode, according to equations (3.11) and (3.12), one hasu ∝ cos (πz) = -sin(πz) ; u(-z) = -u(z), v ∝ cos (πz) = -sin(πz) ; v(-z) = -v(z),w ∝ sin (πz) = cos(πz) ; w(-z) = w(z),

  20) combined with (2.18), one hasu ∝ cos (2πz) = -cos(2πz) ; u(-z) = u(z), v ∝ cos (2πz) = -cos(2πz) ; v(-z) = v(z), w ∝ sin (2πz) = -sin(2πz) ; w(-z) = -w(z),

  [N V ](-z) = -[N V ](z) (D 6)On the other hand, F ad (z), given by equation (3.20),F ad (z) ∝ sin(πz) = cos(πz), (D 7)is even under z → -z. Finally, the product [N V ]F ad is odd under z → -z, hence ]F ad dz = 0, i.e., a 5 = 0. (D 8)

  M. Bouteraa, C. Nouar, C. Metivier, E. Plaut and A. Kalck ical parameters. The study of the stability of fixed points of amplitude equations allows one to analyze the competition between different convection patterns near onset. Then,

  ApAq and FApAq , ḠApAq are the same as the ones on (F 11 , G 11 ). In figure 8, the functions F ApAq , G ApAq , FApAq and ḠApAq are plotted at the critical conditions, P r = 10 and different values of the slip parameter. For SFBC, at

	the critical conditions,
	F ApAq = -	9 π 3 104	1 +	3 P r	sin(2πz) , G ApAq = -	3π 52	9 +	1 P r	sin(2πz),	(5.19)
	FApAq = -	81π 3 5000	3 +	11 P r	sin(2πz) , ḠApAq = -	3π 2500	121 +	27 P r	sin(2πz).

2s ḠApAq + FApAq = 2F 11 (DG 11 ) -G 11 (DF 11 ) , Boundary conditions on F ApAq , G

  2 A 2 A 3 in (D 3) are, first, those due to γpq = A 1 γ[1] pq multiplied by itself and γij = A 2 A 3 γ[2] ij , second those due to γpq = A 1 γ[1] pq multiplied by γpq = A 2 A 3 γ[2] pq and γij = A 1 γ[1] ij . The superscripts [1] and[2] refer to the fundamental mode and to modes arising from the quadratic interaction between Fourier modes with different wavevectors, as calculated in §5.3.2. The nonlinear terms of (D 3) proportional to (A 1 ) 2 A 2 A 3 can thus be developed as

	with
	[N V ] =	∂ ∂z	∂ 2 ∂x 2 t xx +	∂ 2 ∂y 2 t yy -	∂ 2 ∂x 2 t zz -	∂ 2 ∂y 2 t zz + 2	∂ 2 ∂x∂y	t xy
	+	∂ 2 ∂z 2 -	∂ 2 ∂x 2	∂t xz ∂x	+	∂ 2 ∂z 2 -	∂ 2 ∂y 2	∂t yz ∂y	(D 4)

  Therefore, under z → -z, γ[2] xx , γ[2] yy , γ[2] zz , γ[2] xy are even and γ[2] xz , γ[2] yz are odd. The symme-

	try properties of	γ[2] ij are opposite to that of	γ[1] ij . Hence, Γ [1,2] =	1 2	γ[1]

pq γ[2]

  γ[1] ij .As a result, t ij given by equation(D 5) has the symmetry properties of γ[2] ij . Hence, t xx , t yy , t zz , t xy are even, xy are odd and t xz , t yz are odd.

	∂ ∂ z t xx ,	∂ ∂ z t yy ,	∂ ∂ z t zz ,	∂ ∂ z t

The same iterative procedure is used to evaluate the coupling coefficients λ 1 and δ 1 , by replacing X 13 in (B 1)-(B 3) by X13 and X13 respectively.

Acknowledgement

The authors gratefully acknowledge the financial support of the French National Research Agency, under the grant ANR-10-BLAN-925-01. The authors also thank M. Jenny for interesting comparisons with his Freefem code.

which do not depend on the Prandtl number. The integrals involved are evaluated numerically using Clenshaw-Curtis quadrature method [START_REF] Trefethen | Spectral Methods in Matlab[END_REF].

Appendix B. Calculation of the saturation and coupling coefficients outside the critical conditions (i) A first approximation of g 1 is given by:

(ii) A first approximation of X 13 is solution of :

which has to be solved with the boundary conditions (5.24). (iii) A correction of g 1 is given by

Iterations between (B 2) and (B 3) are then performed until convergence.

Appendix C. Amplitude equations, sationary solutions, eigenvalues of the Jacobian matrices Amplitude equations at quintic order as given by Hoyle ( 2006), [START_REF] Fujimura | Hexagons and triangles in the Rayleigh-Bénard problem: quintic-order equations on a hexagonal lattice[END_REF], are recalled, and their stationary solutions and stability in the supercritical regime, in the range of validity of our analysis, is studied.

C.1. Amplitude equation on a square lattice

The stationary solutions and the eigenvalues of the Jacobian matrix are:

-Steady convection with rolls:

The eigenvalues associated to this state are

We have checked numerically that , χ 1 , χ 2 are negative, i.e., rolls are stable.

-Steady convection with squares:

The eigenvalues assciated to this state are

We have checked numerically that χ 1 and χ 2 are positive, i.e., squares are unstable.

C.2. Amplitude equation on an hexagonal lattice

Equations for dA 2 /dt, dA 1 /dt are obtained by cyclic permutations of A 1 , A 2 , A 3 .

The stationary solutions and the eigenvalues of the Jacobian matrix are:

-Steady convection with rolls:

The eigenvalues associated to this state are

We have checked numerically that χ 1 , χ 2 are negative, i.e. , rolls are stable.

-Steady convection with hexagons:

The eigenvalues associated to this state are

We have checked that these eigenvalues are positive, i.e., hexagons are unstable.

Appendix D. Symmetry properties under the midplane reflection -Comparison with [START_REF] Albaalbaki | Pattern selection in the thermal convection of nonnewtonian fluids[END_REF] For an hexagonal lattice, [START_REF] Albaalbaki | Pattern selection in the thermal convection of nonnewtonian fluids[END_REF] wrote the amplitude equations up to O(A 4 ), where we denote O(A) the common order of magnitude of the amplitudes A 1 , A 2 , A 3 of the fundamental modes. They indicated (in their §5.2 p. 527) that the terms of O(A 2 ) vanish, but that terms of O(A 4 ) exist, because 'non-Newtonian effects are symmetry breaking'. Let us show that this does not hold, i.e., that the corresponding