Franck Richard

Aditi Chakrabarti

Basile Audoly

Yves Pomeau

Serge Mora

Buckling of a spinning elastic cylinder: linear, weakly nonlinear and post-buckling analyses

An elastic cylinder spinning about a rigid axis buckles beyond a critical angular velocity, by an instability driven by the centrifugal force. This instability and the competition between the different buckling modes are investigated using analytical calculations in the linear and weakly nonlinear regimes, complemented by numerical simulations in the fully post-buckled regime. The weakly nonlinear analysis is carried out for a generic incompressible hyperelastic material. The key role played by the quadratic term in the expansion of the strain energy density is pointed out: this term has a strong effect on both the nature of the bifurcation, which can switch from supercritical to subcritical, and on the buckling amplitude. Given an arbitrary hyperelastic material, an equivalent shear modulus is proposed, allowing the main features of the instability to be captured by an equivalent neo-Hookean model.

Introduction

As is well known, a liquid sphere deforms into an oblate shape when spun, as the centrifugal forces are stronger at the equator than at the poles [START_REF] Poincaré | Sur l'équilibre d'une masse fluide animée d'un mouvement de rotation[END_REF][START_REF] Chandrasekhar | Ellipsoidal figures of equilibrium : an historical accounts[END_REF]. A spinning elastic body can deform by a similar mechanism. The magnitude of the elastic deformation depends on several factors including the angular velocity, mass density, elastic stiffness, the size and the initial shape [START_REF] Green | Theoretical Elasticity[END_REF]. Among the many possible body shapes, the case of a spinning right circular cylinder is special, as it yields a cylindrically symmetric distribution of centrifugal forces: the direct effect of the centrifugal force is a radial expansion, but the latter is typically limited as it has to work against the bulk modulus of the material. In this paper, we characterize the non cylindrically-symmetric deformations of an initially circular spinning cylinder. To this end, we use bifurcation methods as the symmetry breaking is induced by an elastic instability. Such deformations involve shear strain but no local change of volume, and are favorable for almost-incompressible materials such as typical soft solids.

The linear stability of rotating elastic cylinders has been investigated theoretically by Haughton and Ogden [START_REF] Haughton | a Bifurcation of finitely deformed rotating elastic cylinders[END_REF][START_REF] Haughton | Bifurcation of rotating thick-walled elastic tubes[END_REF]. These authors considered homogeneous cylinders of finite length, as well as hollow tubes, rotating about their axis with angular velocity ω and subjected to a prescribed axial loading. The unbuckled configuration involves a bi-axial stress that has been calculated in terms of the prescribed force and of ω. It has been found to be linearly unstable with respect to prismatic (invariant along the axis), axisymmetric (invariant under rotations) or asymmetric steady deformations when ω exceeds a threshold value depending on the axial force and on the elastic properties of the homogeneous cylinder [START_REF] Haughton | a Bifurcation of finitely deformed rotating elastic cylinders[END_REF][START_REF] Patterson | The stability of a solid rotating neo-Hookean cylinder[END_REF] or of the tube [START_REF] Haughton | Bifurcation of rotating thick-walled elastic tubes[END_REF]. In recent work, the influence of an internal pressure on the linear stability against axisymmetric deformations of rotating tubes has been considered [START_REF] Wang | Localized Bulging of Rotating Elastic Cylinders and Tubes[END_REF]. In the appendix of their book devoted to the mathematical analysis of elastic waves propagating at the surface of spinning elastic cylinders, Rabier and Oden [START_REF] Rabier | Bifurcation in rotating bodies[END_REF] provide a short account of the linear stability against steady prismatic deformations of spinning cylinders based on the plane-strain approximation, ignoring any deformation along the axis of the cylinder, in contrast with the systems studied by Haughton and Ogden [START_REF] Haughton | a Bifurcation of finitely deformed rotating elastic cylinders[END_REF][START_REF] Haughton | Bifurcation of rotating thick-walled elastic tubes[END_REF] where the prescribed load generate a base uniaxial deformation.

To the best of our knowledge, prior work on these instabilities is limited to linear bifurcation analyses. Hence it is not known whether they are supercritical (the buckled configuration is present beyond the linear bifurcation threshold only, and the bifurcation is continuous) or subcritical (the buckled configuration exists below the linear bifurcation threshold, and the bifurcation is discontinuous). In this paper, we carry out a nonlinear buckling analysis of a spinning cylinder. As we have in mind applications to hydrogels or elastomers such as silicone rubbers, the cylinder is modelled as an incompressible hyperelastic material. We consider the case of a slender elastic cylinder whose total length is prescribed, as happens for example if the cylinder is clamped to the rotating device at its ends. In addition, the axis of the cylinder is assumed to be undeformable; in experiments, a stiff core of this kind can be achieved by moulding the soft cylinder around a thin and stiff wire. This system can be viewed as the limit case of a rigid shaft coated by an elastic material whose thickness is far larger than the radius of the rigid shaft itself. In view of these assumptions, we carry out the stability analysis assuming plane strain; this assumption is relaxed and discussed in the appendix, where other modes of bifurcations are analyzed. We work in the frame rotating with the cylinder, and seek equilibrium (i.e. static) solutions in the presence of centrifugal forces. Transients and other dynamic effects are ignored.

The outline of this paper is as follow. We start by recalling the nonlinear equations governing the equilibrium of a cylinder undergoing finite deformations in Section 2. Next, we carry out its linear bifurcation analysis in Section 3. The first unstable mode occurs at a critical angular velocity ωc that is derived in terms of the shear modulus and of the mass density. Next, in Section 4, a nonlinear analysis based on a Koiter post-bifurcation expansion addresses weakly post-buckled solutions. We show that the instability can be either supercritical or subcritical, depending on the first nonlinear coefficient in the series expansion of the strain energy density function. Along the bifurcated branch, the amplitude is found to grow as |ω -ωc| 1/2 and the moment of inertia of the spinning cylinder as |ω -ωc|, with ω > ωc for a constitutive law corresponding to a supercritical bifurcation, and ω < ωc in case of a subcritical bifurcation. Finally, finite-amplitude solutions are studied numerically using the nonlinear finite elements method in Section 6. Bifurcated equilibria for cylinders made of a neo-Hookean elastic material [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF][START_REF] Macosko | Rheology : principles , measurements and applications[END_REF] and involving a strain up to 100% are presented. The approximation consisting of considering an apparent shear modulus (instead of the initial shear modulus) enables us to explain the subcritical-to-supercritical transition based on an equivalent neo-Hookean model.

Equilibrium equations based on a finite strain theory

In this section the nonlinear equations governing the equilibrium of a rotating elastic cylinder are derived, considering an arbitrary hyper-elastic incompressible isotropic material. Starting from the condition that the total potential energy of the system is stationary, the classical equilibrium equations relevant to finite deformations are derived, including the equilibrium conditions both in the bulk and along the lateral boundary.

Let r 0 denote the radius of the undeformed cylinder, ρ its mass density and µ its initial shear modulus, i.e. the shear modulus for infinitesimal strain. The cylinder is spun with an angular velocity ω about its rigid axis, as sketched in Fig. 1. The equilibrium of the system is governed by the only dimensionless parameter in the problem, namely α = ρ r 2 0 ω 2 µ . In the co-rotating frame, both the elastic force and the centrifugal force are conservative. The equilibrium can therefore be derived from the condition that the total potential energy is stationary. The position R of a material point in final configuration is given as a map R(r) in terms of the position r in undeformed configuration. For an isotropic and incompressible elastic material, the strain energy density is a function of the two first invariants, I 1 and I 2 , of Green's deformation tensor C = F T • F, where F = ∂R/∂r is the deformation gradient: I 1 = tr C and

I 2 = 1 2 (tr C) 2 -tr C 2 .
The strain energy density can be written as µ W (I 1 , I 2 ) where W is the dimensionless strain energy density, obtained by scaling out the initial shear modulus µ. The function W (I 1 , I 2 ) captures the response of the material to finite strain. For the strain energy µ W (I 1 , I 2 ) to be consistent with the initial shear modulus µ, the following normalization condition must be enforced,

∂W ∂I 1 (0, 0) + ∂W ∂I 2 (0, 0) = 1 2 . (2.1)
For instance, an incompressible neo-Hookean solid [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF][START_REF] Macosko | Rheology : principles , measurements and applications[END_REF] and an incompressible Mooney-Rivlin solid [START_REF] Mooney | A theory of large elastic deformation[END_REF][START_REF] Rivlin | Large elastic deformations of isotropic materials. IV. Further developments of the general theory[END_REF] correspond, respectively, to W = 1 2 (I 1 -3), and to

W = 1 2 (a 1 (I 1 -3) + (1 -a 1 )(I 2 -3)), with a 1 a material constant in the range [0; 1].
Incompressibility is represented as a constraint J(r) = 1, where J = det F is the Jacobian of the transformation. To characterize equilibrium configurations, we seek stationary points of the augmented energy

E = 0<r<r0; 0<z<h dr µ W (I 1 , I 2 ) - 1 2 ρ ω 2 R • R -(R • ez) 2 + µ q (J -1) . (2.2)
The terms in the integrand are the strain energy, the potential of the centrifugal force, and the Lagrange term taking care of the inextensibility constraint J = 1 by means of a multiplier q(r), respectively. We consider a slender cylinder (h r 0 ), whose terminal sections are blocked by the rotating device. The stiff central wire favors buckling modes that are plain strain, except in a small region of size ∼ r 0 near the ends, which we ignore. Accordingly, attention is limited to modes that are invariant along the axial coordinate z, i.e. ∂/∂z = 0. This assumption is justified in the appendix where the linear thresholds of the other modes (associated to non plane strain deformation) are found to be beyond the threshold of the first unstable plane strain mode. We use polar coordinates (r, θ), with r the distance to the axis and θ the angle, both taken in the reference configuration (Fig. 1). The coordinates (r, θ) are used as Lagrangian coordinates; they label material points in the system. In the deformed configuration where the centrifugal forces are in effect, the polar coordinates are R(r, θ) and Θ(r, θ), as sketched in Fig. 1. The determinant of the deformation gradient reads J = R r R,rΘ ,θ -R ,θ Θ,r , where a comma in subscript denotes a partial derivative. The invariants of Green's deformation tensor read

I 1 = R 2 ,r + 1 r 2 R 2 ,θ + R 2 Θ 2 ,r + R 2 r 2 Θ 2 ,θ + 1, (2.3)
I 2 = R 2 r 2 R 2 ,r Θ 2 ,θ + Θ 2 ,θ -2R,rR ,θ Θ,rΘ ,θ + R 2 ,θ Θ 2 ,r + Θ 2 ,r + Θ 2 r R 2 + R 2 ,θ r 2 + R 2 ,r . (2.4)
The equilibrium equations are derived from the condition that the first variation of Eq.2.2 with respect to the unknowns R(r, θ), Θ(r, θ) and q(r, θ) is zero. Let t = (R, Θ, q) denote the collection of unknowns, and δt = (δR, δΘ, δq) a virtual displacement that is kinematically admissible (abbreviated as 'k.a.'), i.e. such that δR(0, θ) = 0 as imposed by the stiff wire. The field t(r, θ) is a solution of the problem if it satisfies the kinematic boundary condition R(0, θ) = 0 and ∀δt k.a., DE(α, t) [δt] = 0.

(2.5)

Here, DE(α, t) [δt] denotes the first variation of the energy evaluated in the configuration t with an increment δt, also known as the first Gâteaux derivative of the functional E [START_REF] Gâteaux | Fonctions d'une infinité de variables indépendantes[END_REF]. Defining G = r (W (I 1 , I 2 ) + q(J -1)), where the factor r comes from the cylindrical volume element dr = 2 π r dr dz, and integrating by parts, we can derive the equations in the interior from Eqs.2.5 as

J -1 = 0, (2.6
)

∂G ∂R - ∂ ∂r ∂G ∂R,r - ∂ ∂θ ∂G ∂R ,θ = α r 2 0 rR, (2.7)
∂ ∂r ∂G ∂Θ,r + ∂ ∂θ ∂G ∂Θ ,θ = 0, (2.8)
where the first equation is the incompressibility constraint and the two other equations are the equilibrium in the radial and circumferential directions, respectively. These equations are complemented by the condition of zero traction at the lateral boundary r = r 0 , ∂G ∂R,r r=r0 = 0 and ∂G ∂Θ,r r=r0 = 0.

(2.9)

As mentioned above, the central wire imposes the condition R(r = 0, θ) = 0.

(2.10)

Linear bifurcation analysis

The linear bifurcation analysis under the assumption of plane strain has been done by Rabier and Oden [START_REF] Rabier | Bifurcation in rotating bodies[END_REF]. In this section, we carry out this linear analysis in a cylindrical coordinate system. Doing so, we shall recover the previous results of Rabier and Oden, and derive an explicit expression for the modes that will be useful for the nonlinear buckling analysis later in Section 4.

(a) Unbuckled solution

We start by analyzing the unbuckled configuration (base state), and label all quantities relevant to it using a subscript '0'. Because of the incompressibility assumption, the symmetric configuration remains undeformed, i.e. R 0 (r, θ) = r and Θ 0 (r, θ) = θ. The Lagrange multiplier q proportional to a hydrostatic stress is non-zero, however; its value can be found from the radial equilibrium (2.7) and (2.9) as q 0 = α 2 1 -(r/r 0 ) 2 -1. Altogether, the unbuckled solution of Eqs.2.6-2.10 is written as t 0 = (r, θ, q 0 ).

(b) Linearization of the equations about the unbuckled solution

A small perturbation is added to the unbuckled solution, and the equations of Section 2 are linearized with respect to the amplitude of the perturbation,

t = t 0 + εt 1 = (r + εu 1 (r, θ), θ + εΘ 1 (r, θ), q 0 + εq 1 (r, θ)) . (3.1)
In view of the rotational symmetry of the base solution, the perturbation can be assumed to depend on the angular variable θ as exp(i n θ) without loss of generality, where the circumferential wave number n is an integer parameter: the perturbation is sought in the form

         u 1 (r, θ) = Re fu(r)e inθ Θ 1 (r, θ) = Re -if Θ (r)e inθ q 1 (r, θ) = Re fq(r)e inθ (3.2)
where Re denotes the real part. The conventional complex factor (-i) has been included for convenience, anticipating on the fact that the phase of t 1 is shifted by π/2 compared to the phase of the two other unknowns. The goal of the linear bifurcation analysis is to calculate critical values of the loading α = αc(n) such that the linearized equations of equilibrium have non-trivial solutions: this is generically a sign that a branch is bifurcating from the unbuckled branch. Note that the buckling amplitude cannot be obtained from this linear bifurcation analysis; it will be derived later in Section 4. At linear order in ε, Eqs.2.6-2.8 yield

dfu dr + 1 r fu + nf Θ = 0, (3.3)
d 2 fu dr 2 + 1 r dfu dr - n 2 r 2 + 1 r 2 fu - 2n r + αnr r 2 0 f Θ + dfq dr = 0, (3.4)
n αr 2 r 2 + 2 fu -r 3 d 2 f Θ dr 2 -3r 2 df Θ dr + n 2 rf Θ + nrfq = 0. (3.5)
This set of differential equations is complemented with the boundary condition Eqs.2.9-2.10 which read, upon linearization,

dfu dr - 1 r 0 fu -nf Θ + fq = 0 at r = r 0 , (3.6)
df Θ dr - n r 2 0 fu = 0 at r = r 0 , (3.7
)

fu = 0 at r = 0. (3.8)
Upon elimination of f Θ and fq in Eqs.3.3-3.5 one gets an equation for fu:

n 2 -1 2 fu -2n 2 + 1 r dfu dr + 5 -2n 2 r 2 d 2 fu dr 2 + 6r 3 d 3 fu dr 3 + r 4 dfu dr 4 = 0. (3.9)
The general solution of Eq.3.9 that is regular at r = 0 can be expressed in terms of two constants A and B as fu = Ar n+1 + Br n-1 . The boundary conditions Eqs.3.6-3.7 then yield:

r 2 0 n 2 -n -2 - nρω 2 r 2 0 2µ A + n(n -1) - nρω 2 r 2 0 2µ B = 0 (3.10)
r 2 0 (n + 1)A + (n -1)B = 0. (3.11)
Non-trivial solutions correspond to A = 0, and/or B = 0. This is possible whenever the determinant of this linear system is zero. This condition yields the critical value of the control parameter α in terms of the circumferential wave number n as

αc(n) = 2 n - 1 n . (3.12)
We have recovered the result of Rabier and Oden [START_REF] Rabier | Bifurcation in rotating bodies[END_REF]. From Eq.3.11, if n = 1 then A = 0 and u = B is constant, corresponding to an in plane rigid translation, and thus one has to take B = 0 from Eq.3.8. Therefore, non-trivial solutions are possible for n ≥ 2 only. In view of the graphical representation of the dispersion relation (3.12) in Fig. 2(left), the first unstable mode encountered as the angular velocity is increased from zero is the ovalization mode corresponding to n = 2. The corresponding bifurcation takes place at αc(n = 2) = 3. The critical angular velocities predicted by Eq.3.12 are plotted in Fig. 2(right), using the typical mass density ρ ∼ 1 g/cm 3 for a soft material, and n = 2. A cylinder of radius 20 mm made in a soft gel [START_REF] Ferry | Viscoelastic Properties of Polymers[END_REF] having a shear modulus µ = 50 Pa bifurcates when the angular velocity is of order that of a typical car wheel. A cylinder of radius 20 mm made of rubber is predicted to be unstable for the characteristic angular velocity of a ultracentrifuge. These two examples show that the phenomenon we are studying can be encountered in real situations, in rubber joints in rotating parts for example.

In the Appendix, the assumptions of plain and cylindrical invariance are relaxed and an extension of this linear bifurcation analysis is presented. The modes are then characterized by an axial wave vector k and a circumferential wave number n. The first critical load αc is always found to be attained with cylindrically invariant modes, k = 0, see appendix. This validates our assumption of plane strain. For instance, the axisymmetric varicose mode n = 0 and k = 0 takes place with α equal to 6 (Fig. 6(left)). The mixed varicose/ovalization mode with n = 2 and k = 0 is in fact a long-wavelength one (k → 0) that is consistent with the prismatic mode just studied, as it has α = 3 (Fig. 6(right)).

For later reference, we provide the explicit expressions of the complex mode shape fu(r), f Θ (r) and fq(r) (as defined in Eqs.3.2). These expressions are valid for any circumferential wave number n:

fu(r) = 1 2 (n + 1) r r 0 n-1 + (1 -n) r r 0 n+1 , (3.13)
f Θ (r) = 1 r 0 - n 2 + 1 2 r r 0 n-2 + n 2 - 1 n + 1 2 r r 0 n , (3.14
)

fq(r) = 1 r 0 -n 2 + n + 1 - 1 n r r 0 n + n 2 -n -1 + 1 n r r 0 n+2 . (3.15)
In the expressions, the load parameter α has been eliminated by means of the critical condition α = αc(n). The buckling amplitude are normalized using the convention fu(r 0 ) = 1, which warrants that ξ is the radial buckling amplitude on the lateral surface, scaled by the parameter ε appearing in Eq.3.1

u 1 (r 0 , θ) = ξ cos(nθ). (3.16)

Weakly nonlinear analysis

In this section, we carry out a Koiter expansion and derive an expansion of the bifurcated solution in the vicinity of the bifurcation point. This involves pushing to the next order the formal expansion started earlier during the linear bifurcation analysis. Specifically, the displacement field and the Lagrange multiplier are expanded to order 2 in terms of an arc-length parameter ε that we will be able to express in terms of the increment of the load parameter α -αc(n). This method ultimately yields the buckling amplitude ε ξ in terms of the distance to threshold α -αc(n). The Koiter post-bifurcation expansion is a perturbation method that has been already used in various elastic systems [START_REF] Koiter | On the stability of an elastic equilibrium[END_REF][START_REF] Hutchinson | Imperfection sensitivity of externally pressurized spherical shells[END_REF][START_REF] Hutchinson | Postbuckling theory[END_REF][START_REF] Budiansky | Theory of buckling and post-buckling behavior of elastic structures[END_REF][START_REF] Peek | Worst shapes of imperfections for space trusses with many simultaneously buckling members[END_REF][START_REF] Peek | Postbuckling behavior and imperfection sensitivity of elastic structures by the Lyapunov-Schmidt-Koiter approach[END_REF][START_REF] Van Der | Koiter's elastic stability of solids and structures[END_REF]. In our recent paper dealing with the post-buckling analysis in the elastic Rayleigh-Taylor instability [START_REF] Chakrabarti | Selection of hexagonal buckling patterns in transversely isotropic structures: example of the elastic Rayleigh-Taylor instability[END_REF], Koiter's method has been rederived in a self-consistent way; the present work makes use of the same notations and of the same solution strategy as in this previous paper.

(a) Principle of the Koiter expansion

Both the load parameter α and the solution t = (R, Θ, q) are expanded in terms of the small arclength parameter ε as:

α = αc + α 2 ε 2 (4.1) t(α) = t 0 (α) + εt 1 + ε 2 t 2 + ε 3 t 3 + • • • (4.2)
where αc = αc(n) is the critical dimensionless angular velocity from the linear bifurcation analysis, see Eq.3.12, and t 1 is the corresponding linear mode. The load increment (α -αc) has been assumed to scale like ε 2 and not ε, as usual for a symmetric bifurcation: this makes the expansion invariant by the transformation ε ← (-ε), t 1 ← (-t 1 ) and t 3 ← (-t 3 ) corresponding to a rotation of the buckling pattern by half a wavelength in a plane perpendicular to the axis. Onwards, the dependence of αc on n will become implicit in our notation, for the sake of readability. Note that the base solution t 0 depends on the load α through the pressure parameter q 0 . The first-order correction t 1 is the linear mode from Section 3, normalized according to our convention that the radial amplitude on the lateral surface is ξ, u 1 (r 0 , θ) = ξ cos(nθ).

The Koiter method proceeds by inserting the expansion in Eqs.4.1-4.2 into the nonlinear equilibrium written earlier in Eq.2.5 as

∀δt, DE αc + α 2 ε 2 , t 0 (α) + εt 1 + ε 2 t 2 + ε 3 t 3 + • • • [δt] = 0, (4.3)
where δt(r, θ) is the set of virtual functions δR, δΘ, δq that represent infinitesimal increments of the displacements (including the Lagrange multiplier) satisfying the kinematic boundary conditions. Equation 4.3 is then expanded order by order in ε [START_REF] Van Der | Koiter's elastic stability of solids and structures[END_REF][START_REF] Triantafyllidis | Stability of solids: from structures to materials[END_REF].

At order ε, we recover the linear bifurcation problem solved earlier, now written in weak form,

∀δt, D 2 E(αc, t 0 (αc)) • [t 1 , δt] = 0 (4.4)
whose solution yields the critical value of the load αc and the linear mode t 1 ; At order ε 2 , we find an equation for the second-order correction t 2 , ∀δt,

D 2 E(αc, t 0 (αc)) • [t 2 , δt] + 1 2 D 3 E(αc, t 0 (αc)) • [t 1 , t 1 , δt] = 0. (4.5)
which can be solved modulo an unknown (and irrelevant) component proportional to t 1 ; At order ε 3 , one obtains the equation

∀δt, D 2 E(αc, t 0 (αc)) • [t 3 , δt] + D 3 E(αc, t 0 (αc)) • [t 2 , t 1 , δt] • • • + α 2 dD 2 E(α, t 0 (α)) dα α=αc • [t 1 , δt] + 1 6 D 4 E(αc, t 0 (αc)) • [t 1 , t 1 , t 1 , δt] = 0.
Upon insertion of the particular virtual displacement δt = t 1 , the first term cancels out by Eq.4.4 and we are left with

D 3 E(αc, t 0 (αc)) • [t 2 , t 1 , t 1 ] • • • + α 2 dD 2 E(α, t 0 (α)) dα α=αc • [t 1 , t 1 ] + 1 6 D 4 E(αc, t 0 (αc)) • [t 1 , t 1 , t 1 , t 1 ] = 0. (4.6)
This equation yields α 2 , since all the other quantities are known at this stage. In view of Eqs.4.1-4.2, the value of α 2 finally allows to connect the buckling amplitude ε to the increment of the load α -αc = α 2 ε 2 . In the equations above, D 2 E [δt; t 1 ] denotes the second Gâteaux derivative of E, which is a bi-linear symmetric form on the increment t 1 and on the virtual increment δt.

Similarly, D 3 E [δt; t 1 ; t 2 ] is the third Gâteaux derivative (a tri-linear symmetric form).

(b) Second-order correction to the displacement

The second order displacements t 2 = (u 2 , Θ 2 , q 2 ) results from the nonlinear interaction of the linear mode t 1 with itself, see Eq.4.5. As a result, it involves a superposition of Fourier modes with respect to the variable θ having wave numbers n -n = 0; n; and n + n = 2n. The component with wave number n corresponds to the linear mode t 1 with an amplitude that cannot be determined at this order: observe that it is possible to add an arbitrary contribution γ t 1 to the solutions t 2 of Eq.4.5, thanks to Eq.4.4. We will therefore omit this Fourier mode, and retain only the meaningful 0 and 2n components: we seek the second-order correction t 2 to the displacement as

u 2 = ḡu(r) + Re gu(r)e 2inθ , (4.7
)

Θ 2 = ḡΘ (r) + Re -ig Θ (r)e 2inθ , (4.8)
q 2 = ḡq(r) + Re gq(r)e 2inθ . (

In what follows, we focus on the branch emanating from the first unstable load, i.e. we take n = 2 and αc = 3. To calculate the Gâteaux derivative of the total potential energy, an expansion of the dimensionless strain energy W is needed in terms of I 1 -3 and I 2 -3. For the case of plane strain and incompressible materials considered here, I 1 = I 2 , and one can expand the strain energy as in [24]

W = 1 2 (I 1 -3) + β(I 1 -3) 2 + • • • , (4.10)
where β is a nonlinear constitutive parameter. For both a Neo-Hookean solid, and a Mooney-Rivlin solid in plane strain, β = 0. Note that β > 0 as well as β < 0 are compatible with material stability.

Inserting the expression of the linear mode in Eqs.3.13-3.15 and the second-order correction in Eqs.4.7-4.9 into the nonlinear equilibrium in Eqs.2.6-2.10 leads to a set of equations at order ε 2 . Each equation can be split into two equations : one equation for the terms depending on the θ variable as cos(4πθ) (or sin(4πθ)); the other one for the remaining terms. The incompressibility condition at order 2 writes: The radial equilibrium of the boundary in Eq.2.9 reads:

dgu dr + gu r + 4g Θ = 5r 4 2r 6
-(4β + 1) r 0 dgu dr (r 0 ) + (1 -4β) gu (r 0 ) + 4 (1 -4β) r 0 g Θ (r 0 ) -r 0 gq (r 0 ) = 1 2 -10β ξ 2 r 0 , (4.17)
and

-(4β + 1) r 0 dḡu dr (r 0 ) + (1 -4β) ḡu (r 0 ) -r 0 ḡq (r 0 ) = 6β + 1 2 ξ 2 r 0 . (4.18)
The equilibrium of the boundary in the orthoradial direction (see Eq.2.9) reads:

-

r 2 0 dg Θ dr (r 0 ) + 4gu (r 0 ) = ξ 2 r 0 , ( 4

(c) Amplitude equation

At this point, we have calculated the first and second-order displacements t 1 and t 2 . With the help of a symbolic calculation language, we can then determine the quantities appearing Eq.4.6 as The bifurcated branch is found above the critical load in the supercritical case, and below the critical load in the subcritical case. The amplitude equation captures the influence of an arbitrary material model through the single parameter β. For both the neo-Hookean and Mooney-Rivlin material, β = 0 and the bifurcation is supercritical. For the Gent hyper-elastic model, which captures the concept of limiting chain extensibility, the dimensionless strain energy density function is [START_REF] Gent | A new constitutive relation for rubber[END_REF]:

D 3 E(αc, t 0 (αc)) • [t 2 , t 1 , t 1 ] = - 329πβξ 4 40 - 325πξ 4 64 , D 4 E(αc, t 0 (αc)) • [t 1 , t 1 , t 1 , t 1 ] = 975πβξ 4 4 + 6387πξ 4
W Gent = - Jm 2 log 1 - I 1 -3 Jm ⇒ β = 1 2Jm , (4.31)
where the parameter Jm is related to the maximal stretch of the material under a uni-axial deformation. For a maximal stretch equal to 700%, Jm ∼ 50. From Eq.4.29, the ratio of the amplitude ζ of the ovalization mode n = 2 for a Gent elastic material versus a neo-Hookean solid (β = 0) is ζ Gent /ζ neo Hookean ∼ 0.9: even near the bifurcation threshold where the solid is almost undeformed and the chains are not extended, the limited chain extensibility accounts for a reduction of the buckling amplitude by 10%, which is not negligible. The ovalization of the cylinder induces an increase in the moment of inertia of the cylinder, an effect that could be probed by the rotating device. The moment of inertia being quadratic in the amplitude of the displacement, its calculation requires the use of the series expansion of R(r, θ) at order 2 (Eq.4.2). From Eqs.3.13, 4.25, 4.29 and 4.30, the moment of inertia per unit axial length is at order ε 2 (see Fig. 4):

I ∆ = ρ r0 0 2π 0 R 2 rdrdθ = π 2 ρr 4 0 1 + 80 (α -αc) 1296β + 63 , (4.32)
with α > αc if β > -7/144 and α < αc if β < -7/144. The predictions of the amplitude equation in Eq.4.28 are verified against numerical simulations: Fig. 3(left) shows as expected a perfect agreement asymptotically close to bifurcation for four particular material models.

A simple physical account of the influence of the material model on the buckling amplitude

In Section 4 we have shown that the post-buckling behavior is influenced by the constitutive law through the nonlinear elasticity coefficient β. As a result, the neo-Hookean model (β = 0) cannot predict the buckling amplitude of a different material model (β = 0), and does not capture the possibility of a subcritical bifurcation. In this section, we explain the dependence of the postbuckling behavior on β by a simple mean-field type of argument: we modify the Koiter solution for a neo-Hookean material (β = 0) by introducing an effective shear modulus depending on β, and obtain a simple prediction for the buckling amplitude ζ in terms of β, which turns out to be similar to the prediction based on the exact Koiter expansion. In particular, we recover the existence of a transition to a subcritical bifurcation for small values of β.

Our argument makes use of two features of the Koiter solution for a neo-Hookean material: for β = 0 the buckling amplitude in Eq.4.29 is

ζ = r 0 40 63 √ α -αc, (5.1)
and the average value of the invariant I 1 is, from its definition:

I 1 = 3 + 59 6 ζ 2 r 2 0 + O ζ 3 r 3 0 . (5.2)
This expression is accurate to order ε 2 , and has been found by setting β = 0 in the displacement found by the Koiter expansion in Eqs.3.13-3.15 and 4.22-4.27.

For an arbitrary elastic constitutive law we define an apparent nonlinear shear modulus µapp(I 1 ) for plane strain by

µapp(I 1 ) = 2µ W (I 1 ) I 1 -3 . (5.3)
In view of the series expansion of W in Eq.4.10, one obtains

µapp(I 1 ) = µ 1 + β(I 1 -3) + O (I 1 -3) 2 .
(5.4)

In particular, µapp = µ for a Mooney-Rivlin or a neo-Hookean solid (β = 0): one can think of the apparent modulus µapp as the modulus of an equivalent neo-Hookean model for a prescribed value of I 1 .

We postulate that the post-buckling behavior of a cylinder having an arbitrary constitutive law can be predicted by replacing the modulus µ in the prediction of the neo-Hookean model by the apparent modulus µapp( I 1 ) defined based on the average value of the invariant. Replacing α = ρ r 2 0 ω 2 /µ in equation 5.1 with the effective load parameter α eff ,

α eff = ρ r 2 0 ω 2 µapp( I 1 )
.

(5.5)

and combining with Eqs.5.4 and 5.2, we obtain

α eff = α 1 -β 59 6 ζ 2 r 2 0 , (5.6)
after dropping higher-order terms. Replacing α with α eff in the amplitude equation for β = 0 in Eq.5.1 as well, we find Using α ≈ αc = 3 in the left-hand side, we obtain the buckling amplitude as

ζ 2 = r 2 0 (α -αc) / 63 40 + 59 2 β . (5.8)
This approximate prediction is similar to the exact one based on the Koiter expansion in Eqs.4.29-4.30. The only difference is the coefficient of β, which is 162/5 = 32.4 in the exact theory and is now approximated as 59/2 = 29.5. In particular, the approximate formula in Eq.5.8 still predicts a change in the type of bifurcation, now at β = -(63/40)/(59/2) = -0.053.

In view of the above derivation, the dependence of the buckling amplitude on β can be interpreted as follows. If β > 0, the apparent modulus µapp is larger than the initial modulus µ and as a result the effective load parameter α eff is smaller than the true load parameter α, see Eq.5.6, and the buckling amplitude is smaller than in the neo-Hookean case, see Eq.5.8. Conversely, if β < 0 is moderately negative, the buckling amplitude is larger than for a neo-Hookean material; if β is more negative, α eff can be larger than αc even if α is smaller than αc and the bifurcation becomes subcritical as the bifurcated branch Qualitatively, the main features of the solution in the weakly nonlinear regime are well captured by the proposed approximation, which combines the nonlinear expansion for a neo-Hookean solid with an apparent modulus depending on the material non-linearity β and on the amplitude of the deformation.

Post-buckling analysis

To probe the equilibrium solutions in the fully nonlinear regime, we have set up numerical simulations of the finite-elasticity problem described in Section 2. The simulations make use of the finite element method, and of its open-source implementation in the FEniCS library [START_REF] Logg | Automated Solution of Differential Equations by the Finite Element Method[END_REF], and are based on the variational formulation written symbolically in Eq.2.5. We first use the numerical solution close to the bifurcation threshold to verify the predictions of the weakly nonlinear Koiter expansion (Section (a)). Next, proceeding farther from the bifurcation, we use it to probe the range of validity of the weakly nonlinear analysis (Sections (d) and (e)).

(a) Finite element implementation

The simulations make use of a set of units such that r 0 = 1, ρ = 1 and µ = 1, so that α = ω 2 . The domain is a 2-dimensional disk Ω of radius r 0 = 1. A Cartesian coordinate system (x, y) is used, whose origin is at the center of Ω. The orthonormal frame associated with the coordinate axes is denoted by (ex, ey). An isotropic two-dimensional elastic body initially occupying the domain Ω is set up, to represent an isotropic incompressible three-dimensional solid in plane strain. It is subjected to the action of the centrifugal volume force α (xex + yex)/r, where r = x 2 + y 2 is the radial coordinate.

The displacement vector u = uxex + uyey and the Lagrange multiplier p are discretized using Lagrange finite elements with a quadratic interpolation, on a triangular mesh (see chapter 20 of [START_REF] Logg | Automated Solution of Differential Equations by the Finite Element Method[END_REF]). The number of triangular mesh elements is written n 2 mesh , where n mesh is the typical linear mesh density in either direction. The nonlinear problem in the (u, p) variables is solved using a Newton algorithm based on a direct parallel solver (MUMPS). Full equilibrium branches are obtained by progressively incrementing the load parameter α, recording the displacement field, the Lagrange multiplier and the total potential energy of the system, and reaching convergence at each step.

(b) Verification using a neo-Hookean material Different material models have been used, corresponding to different expressions of the scaled strain energy density. We first consider the neo-Hookean case, for which W (I 1 ) = 1 2 (I 1 -3). Starting from the reference configuration (u = 0), α is gradually increased by increments δα = 1/4000. A random disturbance of amplitude ξ 0 = 10 -3 is systematically added to the solution from the previous step, and used as the initial condition for the next step.

For α lower than a threshold which depends on n mesh , no deformation is observed (u = 0). Beyond the threshold value of α the calculation converges to a non-zero solution u. The corresponding deformed shape is ellipsoidal, as shown in Fig. 4-(a). Upon increasing values of α, the ovalization is more and more pronounced. The difference of the total potential energy of the system and the energy of the undeformed cylinder, ∆E, is found to be negative, i.e. the ovalization is energetically favorable. A plot of ∆E as a function of the load α is shown in Fig. 3

(right).

For a given mesh size, one can estimate the buckling threshold α * (n mesh ) accurately by fitting ∆E as ∆E ≈ k(α -α * (n mesh )) 2 . The convergence of α * (n mesh ) with the mesh size, i.e. for 1/n mesh → 0 is shown in the inset. The limit α * (∞) is consistent with the analytical prediction αc = 3, with a residual error well below 10 -4 . In the following, the linear mesh density n mesh is always taken larger than 12, which warrants a residual error of the order (α * (n mesh = 12) -αc)/αc < 0.2%.

(c) Verification using a general material

Simulations have also been performed with a more nonlinear material model

W (I 1 ) = 1 2 (I 1 - 3) + β 2 (I 1 -3) 2 .
For all the tested values of β larger than -7/144, again, no deformation is observed if α < α * (n mesh ), the deformation grows quadratically above α * (n mesh ), and the buckling threshold convergences to the same buckling threshold αc = 3. For all the tested values of β smaller than -7/144, no deformation takes place when α passes α = 3. However, when α is decreased by small increments, the solution bifurcates from u = 0 as α becomes less that an critical value α * (n mesh ). We have checked that α * (n mesh ) converges to αc = 3 for finer and finer meshes in this subcritical case as well.

The buckling amplitude is measured from the simulations as ζ = 1 2 (Rmax/R min -1), where R min and Rmax are the smallest and largest distance from the deformed lateral boundary to the origin. This amplitude ζ is plotted as a function of α in Fig. 3(left) for different values of β, and compared to the prediction based on the Koiter expansion from Eqs.4.29-4.30. In this figure and the following figures, the horizontal axis is α -αc, with αc = 3 for the analytic curves for mode n = 2, and αc = α * (n mesh ) for the simulation data: this choice cancels the small shift in the instability threshold caused by the finite mesh size. The agreement is good in the limit of the small amplitudes, which indeed corresponds to the domain of validity of the Koiter expansion.

(d) Gent elastic cylinders

To assess the domain of validity of the Koiter expansion, which is in principle restricted to a small neighborhood of the bifurcation point, we have plotted in Fig. 4 the moment of inertia and the buckling amplitude of a neo-Hookean cylinder over a wide range of α. The buckling amplitude defined here by (Rmax -R min )/Rmax = 1 -R min /Rmax, is also one minus the aspect ratio for the ovalization mode n = 2 considered here. Comparison to the prediction of the Koiter postbifurcation expansion at order (α -αc) 2 , yields good agreement close to bifurcation.

In Section 5, we have shown that non neo-Hookean constitutive laws can be accounted for by using an effective load parameter α eff in the predictions of the neo-Hookean model. Here, we show how this approximation can be extended to the case of finite deformation. We start by capturing the average value of the invariant I 1 as obtained by the numerical simulation of a neo-Hookean material (black dots in Fig. 5) using a fitting function f nH ,

I 1 = f nH (α). (6.1)
This fitted function f nH serves as a replacement for Eq.5.2 which, being part of a Koiter expansion, is inapplicable beyond the neighborhood of the bifurcation point. Next, we define an effective load parameter based on I 1 as earlier in Section 5: inserting Gent's constitutive law in Eq.4.31 into the definition of the apparent modulus µapp in Eq.5.3, and the latter into the definition of the effective modulus in Eq.5.5, we find Finally, we replace the load parameter α appearing in Eq.6.1 by α eff , as we did earlier in Section 5.

α eff = α I 1 -3 -Jm log 1 - I1 -3 Jm . (6.2) (d) (c) (b) (a) . slo pe =1 /2 α -α * (n mesh ) R max /R min -1 1 10 -1 10 -2 10 -3 10 -4
The result is

I 1 = f nH (α eff ) . (6.3)
This implicit equation yields the two solid curves in Fig. 5, which agree well with the predictions of the numerical simulations. Fig. 5 confirms that the nonlinear solution depends mainly on the load α and on the constitutive parameter β through the combination α eff : for a fixed value of I 1 -3 or, equivalently, of α eff , the solutions corresponding to different values of the constitutive parameter Jm are almost indistinguishable (see snapshots in Fig. 5).

(e) Case of a subcritical bifurcation

We analyze the numerical solutions having a finite amplitude, obtained when the constitutive law produces a subcritical bifurcation (β < -7/144, see Eq.4.30). To this end, we consider the constitutive law

W = 1 2 (I 1 -3) -0.1(I 1 -3) 2 + 0.2(I 1 -3) 3 . (6.4)
The constitutive parameter β = -0.1 is indeed less than -7/144 and we expect a subcritical bifurcation. The higher-order ensures that the energy is bounded from below (material stability). Note that W is an increasing function of I 1 for the range of the deformations that are considered here.

Starting from α well above αc = 3 and decreasing α progressively, we find that the buckling indicator I 1 -3 departs from zero at α = α * (n mesh ), see Fig. 5(right).

The excess of energy ∆E is positive and increases as α decreases (inset of Fig. 5(right)), meaning that the energy of the buckled solution is larger than that in the unbuckled configuration. A fold point is then encountered: for a value αm of α, the slope of I 1 -3 as a function of α becomes infinite, and no solution is found below αm except the trivial one. Conversely, starting from the fold point at α = αm, and increasing α now, one finds a new branch of solutions with material with the dimensionless strain energy density as in Eq.6.4: approximate prediction based on the effective load parameter α eff , see Eqs.6.3 and 6.5, (line) and finite element simulations (dots). The stability of the branches is inferred from the Leray-Schauder degree theorem: (+1) indicates a stable branch, (-1) an unstable one. Inset: Dimensionless total potential energy as a function of ααc, using the unbuckled solution as a reference. larger buckling amplitude, as seen in Fig. 5(right). The excess of energy decreases and eventually becomes negative as one moves along this part of the branch, increasing α from the initial value αm. For α < αm the only solution found is the trivial one, u = 0; we call this the stable unbuckled branch; For αm < α < α * (n mesh ) there are three solutions, located on the stable unbuckled branch, on the lower bifurcated branch (enclosed between the initial bifurcation point and the fold point), and on the upper bifurcated branch; For α > α * (n mesh ) there are two solutions, corresponding to the so-called unstable unbuckled branch, and to the upper bifurcated branch. An approximation of this bifurcation diagram can be obtained without numerical calculation by the trick based on the effective load parameter. To do so, it suffices to insert into Eq.6.3 the expression of α eff relevant to the constitutive model in Eq.6.4, namely

α eff = α I 1 -3 1 -0.1 ( I 1 -3) + 0.2 ( I 1 -3) 2 . (6.5)
The result is shown by the continuous line in Fig. 5(right). While the agreement with the numerical simulations is not perfect, the approximate theory captures well the general aspect of the bifurcation diagram, including the presence of a fold point.

The stability of the branches can be addressed with the Leray-Schauder rule [START_REF] Normand | Convective instability: a physicist's approach[END_REF]. An index (+1) is assigned to a stable branch and (-1) to an unstable branch. From the Leray-Schauder degree theory the sum of the indexes has to be constant for any value of the control parameter α. The left unbuckled branch is stable, and its index is (+1). A small disturbance of a solution in the right unbuckled branch causes a negative ∆E: The index of this branch is therefore (-1). The buckled branches count double, since for a given value of I 1 the amplitude of the mode can be equally positive or negative. One concludes that the index of the upper branch has to be (+1) and the index of the lower branch is (-1) [START_REF] Cross | Pattern Formation and Dynamics in Nonequilibrium Systems[END_REF], as summarized Fig. 5(right)). We end up with the usual situation for a subcritical bifurcation: for α < αm the unique stable solution is the base state; two stable solutions coexist for αm > α > αc, namely the base state and the upper branch; going from αm to αc the lowest-energy solution transitions from the unbuckled one to the one on the upper branch; for α > αc, the only stable solution is the one on the upper branch.

Conclusion

We have studied the prismatic deformations induced by centrifugal forces in a hyperelastic cylinder spinning about its axis. These deformations are driven by the dimensionless control parameter α which is proportional to the square of the angular velocity, the radius of the cylinder, its mass density, and the inverse of the shear modulus.

A sinusoidal radial disturbance with a circumferential wave number n and an infinitesimal amplitude has been shown to become critical by a linear bifurcation analysis for a set of values α = αc(n) that depend on the wave number but not on the constitutive law of the material. The first critical value αc(n) corresponds to an ovalization of the cross-section, n = 2, and will therefore be observed when α is increased from zero as the cylinder is set into rotation.

The amplitude of the deformation of the cross-section has been characterized by a weakly nonlinear expansion. The first nonlinear elasticity modulus, β, has been shown to influence both the buckling amplitude and the nature of the bifurcation which, for the first mode n = 2, is subcritical when β < -7/144 and supercritical when β > -7/144.

Fully post-buckled configurations have been calculated numerically using the finite elements method, both in the supercritical and subcritical cases. The post-buckled shapes of a cylinder having an arbitrary constitutive law characterized by a material parameter β are well approximated by a neo-Hookean model (β = 0) making use of an apparent shear modulus depending on β and of the buckling amplitude. This interpretation in terms of an effective neo-Hookean model provides a simple and unified way to capture the nonlinear features of this instability, including the effects of the constitutive model on the buckling amplitude and on the nature of the bifurcation, as well as the general aspect of the bifurcation diagram, both in the subcritical and supercritical cases.

This instability causes significant modifications in the mechanical properties of rotating cylinders. It has therefore to be considered in the design of devices involving compliant elastic pieces that can rotate at high angular velocity, as joints or insulators placed in between rotating rigid pieces. Since the elastic properties of most of these soft materials are temperature dependent, these changes may occur-or disappear-while the device is operating for a long time, due to local heating. Another potential application is related to the development of patterning induced by external forces (here the centrifugal force): once a targeted deformation of a soft elastic cylinder is achieved by adjusting its angular velocity, a chemical reaction may be initiated in order to imprint the deformed configuration into the material permanently.

Ethics. This work did not involve experiments.

A. Linear stability analysis including variations along the axis of the cylinder

In this appendix, we go beyond the plane strain deformation. A linear stability analysis taking into account deformations in the radial, circumferential and axial directions is carried out. The first unstable mode is found to be the prismatic (plane-strain) one with the circumferential wave number n = 2. This result validates the relevance of the plane strain assumption for a post-buckling analysis.

(a) Linearization of the base equations

In the following we use cylindrical coordinates, with r the distance to the axis, θ the angle and z the height in the unperturbed state (Fig. 1). After deformation, the polar coordinates become R(r, θ, z), Θ(r, θ, z), and Z(r, θ, z). In this coordinate system, the deformation gradient F = ∇R(r) writes:

F =    R,r RΘ,r Z,r 1 r R ,θ R r Θ ,θ 1 r Z θ R,z RΘ,z Z,z    . ( A 1)
Expressions of J, I 1 and I 2 are obtained from Eq.A 1. Defining G = r (W (I 1 , I 2 ) + q(J -1)), Eqs.2.5 write in strong form: In a standard way, we assume for the linear stability analysis a harmonic θ and z dependence of any variation of the perturbation of R, Θ, Z and q: u 1 (r, θ, z) =

J -1 = 0, ( A 2
Re fu(r)e inθ+ikz , Θ 1 (r, θ, z) = Re -if Θ (r)e inθ+ikz , z 1 (r, θ, z) = Re -ifz(r)e inθ+ikz , q 1 (r, θ, z) = Re fq(r)e inθ+ikz , where n is the circumferential wave number and k is the axial wave number. Since R(r, θ, z) = R(r, θ + 2π, z), n has to be an integer. Furthermore, from the symmetry θ ↔ -θ, one can take n positive. One can take k > 0 as well. Eqs.A 2-A 6 lead at order ε to a differential equation for fu of order 6:

-r 6 d 6 fu dr 6 -9r In the next sections, the main steps to derive the linear solvability condition for asymmetric modes with circumferential numbers n = 0 (section (b)) and n = 2 (section (c)) are developed.

(b) Axisymmetric modes

We now consider deformations that are invariant by rotation but depend on z, giving to the cylinder a varicose shape, which formally amounts to imposing n = 0 and f Θ = 0. Eq.A 7 simplifies to L -k 2 2 fu = 0 with L = d 2 dr 2 + 1 r d dr -1 r 2 . Omitting terms diverging at r = 0, the general solution of this equation is found by using the Wronskian method and can be presented as u(r) = AI 1 (kr) + Bk 2 K 1 (kr) r0 r

Figure 1 .

 1 Figure 1. Sketch of the (a) reference and (b) deformed configurations. The axis of the cylinder is an undeformable wire,and the displacement is blocked at its two ends. We work under the assumption of plane strain. A material point with initial polar coordinates (r, θ) has final cylindrical coordinates R(r, θ), Θ(r, θ).

10 Figure 2 .

 102 Figure 2. Left: Buckling load α = ρω 2 r 2 0

160 , dD 2 E= -πξ 2 .

 22 (α, t 0 (α)) dα α=αc Equation 4.6 then yields the sought relation between the scaled load increment α 2 and the buckling amplitude ξ: α 2 = (63/40 + 162/5β) ξ 2 /r 2 0 . Multiplying both sides by ε 2 and identifying the load increment α -αc from Eq.4.1 and the true buckling amplitude ζ = ε ξ, we find the amplitude equation for the mode n = 2 as αnormalization convention in Eq.3.16, the buckling amplitude ζ is the amplitude of the radial displacement on the outer surface, i.e. R(r 0 ) = r 0 + Re(ζ e 2iθ ) + • • • Based on the amplitude equation in Eq.4.28, the bifurcation is supercritical (continuous) if the material parameter β is such that 63/40 + 162β/5 > 0, and subcritical (discontinuous) otherwise: ζ = r 0 (α -αc) 0486 (supercritical) (4.29) ζ = r 0 (αc -α) / 63 40 + 162 5 β , α < αc if β < -7 144 -0.0486 (subcritical) (4.30)

5 - 2 . 10 - 5 - 3 . 10 - 5 - 4 . 10 - 5 - 5 . 10 - 5 - 6 . 10 - 5 .Figure 3 .

 5210531054105510561053 Figure 3. left: Dimensionless amplitude of mode n = 2 as a function of the control parameter ααc for a spinningcylinder with different β in the elastic constitutive laws (Eq. 4.10), in the limit of the small deformations. The lines result from the Koiter post-bifurcation expansion (Eqs.4.29 and 4.30). The filled circles result from simulations using the finite elements method. right: Dimensionless energy ∆E as a function of the control parameter α computed from finite elements simulations with a linear mesh density equal to n mesh = 12. ∆E is zero below the critical value α * (n mesh ) of α, and starts to decrease quadratically beyond these value. Inset: Convergence of the critical load α * with the linear mesh size 1/n mesh towards the theoretical value α * = 3 for the ovalization mode (n = 2).

  = α -αc.

Figure 4 .

 4 Figure 4. Moment of inertia (top left) and aspect ratio (bottom left) from the finite element simulations (dots), for a neo-Hookean cylinder rotating about its axis, and comparison with the result of the Koiter expansion for β = 0, from Eq.4.29 and 4.32 (dotted line) (a)-(b) Snapshots corresponding to the labels placed above the curves, with a color map indicating the amplitude of the displacement, relative to the maximum displacement which respectively 0.05r 0 (snapshot (a) for α = 3.006), 0.13r 0 (snapshot (b) for α = 3.025), 0.20r 0 (snapshot (c) for α = 3.05) and 0.41r 0 (snapshot (d) for α = 3.15).

Figure 5 .

 5 Figure 5. Buckling indicator I 1 -3 as a function of the control parameter ααc. Note that this choice of the buckling indicator makes the tangent at the bifurcation point no longer vertical. Left: Case of spinning cylinders made of a neo-Hookean material (black disks), and a Gent material with Jm = 100 (dark grey) and Jm = 50 (light grey). The Gent materials are defined by Eq.4.31. The two black lines correspond to the predictions of the neo-Hookean model modified so as to use the effective load parameter, see Eq.6.2-6.3. Right: Case of a spinning cylinder made of an hyperelastic

5 d 5 fu dr 5 + 4 d 4 fu dr 4 + 6r 3 3 (kr) 2 + 1 + n 2 d 3 fu dr 3 - 2 - 4 0 d 4 fu dr 4 + 6r 0 3 d 3 fu dr 3 + r 2 0 5 - 6 0 d 6 fu dr 6 + 9r 5 0 d 5 fu dr 5 + r 4 0 12 -d 3 fu dr 3 +r 2 0n 4 +Figure 6 . 6 0 d 6 6 0 d 4 3 +r 2 0n 4 +

 54432132425654122466664324 Figure 6. Critical value of α beyond which the instability occurs, as a function of k. left: For axisymmetric modes. right: For asymmetric mode with n = 2

  Eliminating g Θ and gq from Eqs.4.11, 4.13 and 4.15 one gets a differential equation for gu:

	r 3 3	d 4 gu dr 4 + 2r 2 d 3 gu dr 3 -9r	d 2 gu dr 2 -11	dgu dr	+	75gu r	= -	36ξ 2 0 r 2	.	(4.21)
	Using the boundary conditions Eqs.4.17, 4.19 one obtains gu as
				gu(r) = -	9rξ 2 16r 2 0	+	5r 3 ξ 2 8r 4 0	-	r 5 ξ 2 0 16r 6	,	(4.22)
	from which we find g Θ and gq as
				g Θ (r) =			9ξ 2 8r 2 0	-	r 2 ξ 2 r 4 0	+	r 4 ξ 2 0 4r 6	,	(4.23)
				gq(r) = 0.		(4.24)
	Eqs.4.12, 4.14, 4.16, 4.18 and 4.20 can be solved by a similar method,
	ḡu(r) =	27r 2 4r 2 0	-	15r 4 2r 4 0	+	7r 6 4r 6 0		ξ 2 4r	,	(4.25)
	ḡΘ (r) = 0,								(4.26)
	ḡq(r) =	-(72β + 6) + (144β + 9)	r 2 r 2 0	+ (9 -72β)	r 4 r 4 0	-6	r 6 r 6 0	ξ 2 0 4r 2	.	(4.27)
											.19)
	and							dḡ Θ dr		(r 0 ) = 0.	(4.20)

Data Accessibility. This article has no additional data. Authors' Contributions. Theory and calculations have been developped by AC, BA, YP and SM. FR and SM implemented and performed the simulations.

Competing Interests. The authors have no competing interests.

Funding. No funding to report 1 (kr )r dr -I 1 (kr) r0 r r I 1 (kr )K 1 (kr )dr , (

where Ia and Ka are the modified Bessel functions of first and second kind of order a respectively, and A and B are integration constants. The boundary conditions Eqs.A 8 and A 10 writes 2A + B = 0 and -2kI 0 (kr 0 ) + 2 r0 I 1 (kr 0 ) + α r0 I 1 (kr 0 ) A + kI 0 (kr 0 )B = 0. The nontrivial solution (A = 0 and/or B = 0) exists if the determinant of this linear system is zero. This condition yields α = 4kr0I0(kr0) I1(kr0)

-2. The axisymmetric deformation first appears in the long scale modes (k = 0), beyond a threshold for α which is found to be equal to 6 (see Fig. 6(left)). This value being larger than the threshold for the prismatic mode with n = 2, one concludes that the mode invariant by rotation would occur far beyond the threshold for the prismatic mode n = 2.

(c) Asymmetric modes

In addition to the previous section and to section 3, we consider asymmetric modes with k = 0 and n = 0. In order to find the behavior of the solutions of the differential equation A 7 in the limit r → 0, let first assume that fu(r) ∼ r a as r → 0. Substituting this particular expression of fu in Eq.A 7 and considering the limit r → 0, one obtains the condition on exponent a:

The prismatic mode n = 2 has been found to corresponds to the first prismatic mode that appears upon increasing values of αc (Section 3 ). Aiming at finding the modes corresponding to the smallest critical value of α, one considers now asymmetric modes with n = 2. The roots of Eq.A 12 being 1,3,5,-1,-3, with 1 a double root, we take s 1 (r) ∼ r, s 2 (r) ∼ r 3 and s 3 (r) ∼ r 5 as r → 0, and we look for series expansion of s 1 (r), s 2 (r) and s 3 (r) in the form s i (r) = ∞ m=1 am(kr) m . The conditions for s i (r) to be a solution of Eq.A 7 is, for m ≥ 6:

a 1 , a 3 and a 5 are not fixed up to now. Coefficients with an even index have to be 0. We choose a 1 = 1, a 3 = a 5 = 0 for s 1 (r) ; a 1 = a 5 = 0 and a 3 = 1 for s 2 (r) ; a 1 = a 3 = 0 and a 5 = 1 for s 3 (r). Writing fu(r) as fu(r) = As 1 (r) + Bs 2 (r) + Cs 3 (r) and substituting this expression of fu(r) in the boundary conditions Eqs.A 8-A 10, one obtains the threshold value of αc as a function of k for n = 2 (See Fig. 6(right)). The mixed modes with the fixed circumferential number n = 2 first develops with an infinite wave length (k = 0), the instability threshold being the one of the prismatic mode n = 2.

Mixed modes with n > 2 are also first unstable in the long wave length limit kr 0 → 0. The corresponding instability threshold is equal to the one of the prismatic mode with the same circumferential wave number, and αc(n, k) increases as k increases.

To conclude, we have considered in the appendix 3d deformations. The first unstable mode is the prismatic mode with the circumferential wave number n = 2. Mixed modes with the same circumferential number (n = 2) are found to be linearly unstable as the control parameter is slightly increased beyond the threshold value of the prismatic mode. The axial wave number then continuously increases.