Vivien Lecomte
email: vivien.lecomte@univ-paris-diderot.fr

Juan P Garrahan

Frédéric Van Wijland

Inactive dynamical phase of a symmetric exclusion process on a ring

We investigate the nature of the dynamically inactive phase of a simple symmetric exclusion process on a ring. We find that as the system's activity is tuned to a lower-than-average value the particles progressively lump into a single cluster, thereby forming a kink in the density profile. All dynamical regimes, and their finite size range of validity, are explicitly determined.

Introduction

The study of temporal large deviations in simple systems such as those of interacting diffusing particles has revealed an unexpected wealth of behaviors, ranging from puzzling correspondences between driven systems and their equilibrium counterparts [START_REF] Tailleur | Mapping Nonequilibrium onto Equilibrium: The Macroscopic Fluctuations of Simp[END_REF][START_REF] Tailleur | Mapping out-of-equilibrium into equilibrium in one-dimensional transport models[END_REF][START_REF] Imparato | Equilibriumlike fluctuations in some boundary-driven open diffusive systems[END_REF][START_REF] Lecomte | Current Fluctuations in Systems with Diffusive Dynamics, in and out of Equilibrium[END_REF], to intrinsically dynamic phase transitions [START_REF] Bodineau | Distribution of current in nonequilibrium diffusive systems and phase transitions[END_REF][START_REF] Bodineau | Cumulants and large deviations of the current through non-equilibrium steady states[END_REF][START_REF] Garrahan | Dynamical First-Order Phase Transition in Kinetically Constrained Models of Glasses[END_REF][START_REF] Garrahan | First-order dynamical phase transition in models of glasses: an approach based on ensembles of histories[END_REF][START_REF] Bertini | Current Fluctuations in Stochastic Lattice Gases[END_REF][START_REF] Bertini | Non Equilibrium Current Fluctuations in Stochastic Lattice Gases[END_REF]. The dynamical phase transition terminology dates back to the eighties [START_REF] Beck | Thermodynamics of chaotic systems: an introduction[END_REF], in an era when the dynamical systems community exploited the so-called thermodynamic formalism of Ruelle and others to study and characterize the variety of dynamical regimes displayed by simple iterated maps. There have been many efforts in the recent past to classify the various phase transitions that can be found between distinct dynamical regimes, but there is a great scarcity of studies attempting an even coarse description of the properties of these dynamical regimes. Exceptions are Bodineau and Derrida's study of a weakly asymmetric exclusion process [START_REF] Bodineau | Distribution of current in nonequilibrium diffusive systems and phase transitions[END_REF], Toninelli and Bodineau's work on kinetically constrained models [START_REF] Bodineau | Activity phase transition for constrained dynamics[END_REF] and the finite size study of Bodineau, Lecomte and Toninelli [START_REF] Bodineau | Finite size scaling of the dynamical free-energy in a kinetically constrained model[END_REF]. In the generic dynamical phase transition scenario, one of the dynamical phases is disordered and its physics is easy to grasp, while the ordered symmetry-breaking phase is usually hard to characterize, for in dynamics there is no such thing as an easy-to-implement free-energy based variational principle.

In the present work, we would like to devote our own efforts to the understanding of the nature of the various dynamical phases of the Simple Symmetric Exclusion Process (SSEP) on a ring, a system of mutually excluding particles hopping (if the target site is empty) with equal rates to either of their nearest neighbor sites on a one-dimensional lattice with periodic boundary conditions. The SSEP belongs to the broader class of diffusive systems whose dynamics can be described by fluctuating hydrodynamics, meaning that at the scale fixed by the system size, the field of occupation numbers becomes a smoothly varying function evolving through a Langevin equation with a vanishingly small noise (as the system size increases). We choose to classify the various time realizations the system can follow according to their activity level defined as the number of particle moves having taken place over a given time interval. In practice, however, we introduce a Lagrange multiplier s biasing trajectories towards a given average activity [START_REF] Lecomte | Thermodynamic Formalism for Systems with Markov Dynamics[END_REF]. The s > 0 (resp. s < 0, s ∼ 0) regime allows one to explore the less-than-average-activity (resp. larger-thanaverage, typical) trajectories. The latter activity is a physical observable that has received considerable attention in the past [START_REF] Garrahan | Dynamical First-Order Phase Transition in Kinetically Constrained Models of Glasses[END_REF][START_REF] Lecomte | Thermodynamic Formalism for Systems with Markov Dynamics[END_REF][START_REF] Merolle | Space-time thermodynamics of the glass transition[END_REF], not only in diffusive systems but also in a variety of systems with slow dynamics [START_REF] Bodineau | Large Deviations of Lattice Hamiltonian Dynamics Coupled to Stochastic Thermostats[END_REF][START_REF] Lester | Dynamic Order-Disorder in Atomistic Models of Structural Glass Formers[END_REF][START_REF] Pitard | Dynamic transition in an atomic glass former: a molecular dynamics evidence[END_REF], and which has strong connections with entropy production in nonequilibrium systems [START_REF] Maes | On and beyond Entropy Production: the Case of Markov Jump Processes[END_REF]. This is the simplest trajectory-dependent time-reversal symmetric observable one can think of. For realizations showing typical or larger than typical activity, the particles remain uniformly distributed, and the density profile remains flat. This is the disordered regime. By contrast, if we focus on realizations displaying a smaller than typical activity (s < 0), the uniform profile becomes unstable. This scenario was identified in [START_REF] Appert-Rolland | Universal cumulants of the current in diffusive systems on a ring[END_REF], and presents features similar to the transition at high enough values of the current in the weakly asymmetric exclusion process [START_REF] Bodineau | Cumulants and large deviations of the current through non-equilibrium steady states[END_REF]. The transition point, occurring right at the typical (or average) activity is characterized by a universal scaling function that does not depend on the specifics of the SSEP, in the large size limit. The purpose of our work is precisely to investigate the properties of the ordered phase in which translation invariance, as our intuition tells us, will be broken. The extreme lowactivity behavior is indeed going to be dominated by realizations in which particles have clustered into a single lump, thus permitting particle moves at the borders of the cluster, leading to a strongly subextensive activity.

Before we proceed, we would like to recall that studying the statistics of the activity in the SSEP is equivalent to investigating the ground-state properties of a ferromagnetic XXZ chain [START_REF] Appert-Rolland | Universal cumulants of the current in diffusive systems on a ring[END_REF], which allows one to connect our work to the integrable systems literature. We are indeed interested in the smallest eigenvalue -ψ(s) of the operator

H(s) = L 2 - e -s 2 L k=1 [σ x k σ x k+1 + σ y k σ y k+1 + e s σ z k σ z k+1 ] (1)
where the σ α i 's are the Pauli matrices. In standard notations (see Baxter [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF] for instance) the anisotropy parameter ∆ is related to the Lagrange multiplier s enforcing a given average activity by ∆ = e s . Another related problem [START_REF] Elliott | Exact Analysis of an Interacting Bose Gas. I. The General Solution and the Ground State[END_REF][START_REF] Gaudin | La fonction d'onde de Bethe[END_REF] is that of determining the ground state of a gas of hard-core bosons with interaction strength -s and density ρ 0 , which, in the present s > 0 case, would thus correspond to attractive interactions (this is the no-saturation case discarded by Lieb and Liniger at the beginning of their seminal work [START_REF] Elliott | Exact Analysis of an Interacting Bose Gas. I. The General Solution and the Ground State[END_REF]). In the regime of interest here, we have that ∆ > 1. The same problem of attractive bosons with hardcore interactions is of direct interest to study the directed polymer in the replica approach [START_REF] Kardar | Replica Bethe Ansatz studies of two-dimensional interfaces with quenched random impurities[END_REF], for which the problem was fully diagonalized [START_REF] Calabrese | Free-energy distribution of the directed polymer at high temperature[END_REF][START_REF] Dotsenko | Bethe ansatz derivation of the Tracy-Widom distribution for one-dimensional directed polymers[END_REF], see [START_REF] Corwin | The Kardar-Parisi-Zhang equation and universality class[END_REF] for a review. It was shown [START_REF] Albertini | XXZ model as an effective Hamiltonian for generalized Hubbard models with broken eta-symmetry[END_REF] that such a model and such a regime for ∆ were relevant for the description of the superconducting phase of some generalized Hubbard models, when working at fixed magnetization. The latter constraint, which renders the problem nontrivial (without the constraint, see section 8.8 in Baxter [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF]) turns out to be our working framework, since particles are conserved in our closed ring, and magnetization per site is 2ρ 0 -1. While the active dynamical phase lent itself moderately easily to a Bethe Ansatz approach [START_REF] Appert-Rolland | Universal cumulants of the current in diffusive systems on a ring[END_REF], which was shown to be fully equivalent to the fluctuating hydrodynamics approach, it can be understood from [START_REF] Albertini | XXZ model as an effective Hamiltonian for generalized Hubbard models with broken eta-symmetry[END_REF] that things only get worse for ∆ > 1. It therefore appears to be desirable to resort an alternative method. It seems to us that the phase separation mechanism conjectured in [START_REF] Albertini | XXZ model as an effective Hamiltonian for generalized Hubbard models with broken eta-symmetry[END_REF] has a simple interpretation in our approach, as will become clear later in the paper. Fluctuating hydrodynamics, in the integrable systems language, can be viewed as the correct effective field theory able to describe low lying excitations of the system. The questions we wish to answer can be phrased as follows. As the system's activity is lowered below its average, how does the flat density profile become unstable? Is there any limiting profile deep into this inactive phase? Which are the activity scales, and what are their system size dependence, governing these various dynamical inactive regimes? Given that these are generic questions that can be raised for other comparable systems, it is of interest to present a model with a dynamic phase transition for which they can be answered in an exact manner. We will begin with a brief overview of existing results related to large deviations in the SSEP, as predicted by fluctuating hydrodynamics-also termed Macroscopic Fluctuation Theory by Bertini et al. [START_REF] Bertini | Fluctuations in Stationary Nonequilibrium States of Irreversible Processes[END_REF][START_REF] Bertini | Macroscopic Fluctuation Theory for Stationary Non-Equilibrium States[END_REF][START_REF] Bertini | Large Deviations for the Boundary Driven Symmetric Simple Exclusion Process[END_REF][START_REF] Bertini | Current Fluctuations in Stochastic Lattice Gases[END_REF][START_REF] Bertini | Non Equilibrium Current Fluctuations in Stochastic Lattice Gases[END_REF]. This will allow us, in Section 3, to present explicit results on the density profile in the inactive regime. Three scaling regimes exist, which we will investigate separately in a mathematically controlled fashion. Open questions will be gathered in the final section.

The symmetric exclusion process and its fluctuating hydrodynamics description

Fluctuating hydrodynamics

We consider a Simple Symmetric Exclusion Process (SSEP) on an L-site onedimensional lattice with periodic boundary conditions. Each of the mutually excluding N particles can hop with rate 1 to either of its nearest neighboring site, if empty. We consider time realizations of the process over the time interval [0, t] and we choose to work in space and time units scaled by the system size L and the diffusion relaxation time L 2 , respectively. We thus introduce a smoothly varying field ρ(x, τ ) of occupation numbers, defined by n j (t ′ ) = ρ(j/L, t ′ /L 2 ), whose existence we assume. The field ρ evolves according to a conserving Langevin equation

∂ τ ρ + ∂ x j = 0 (2) with j = -D∂ x ρ -1 √ L
√ σξ, and where D(ρ) = 1 and σ(ρ) = 2ρ(1 -ρ) are functions reflecting the SSEP dynamics at a macroscopic scale. The Gaussian white noise,

ξ(x, τ ), has unit correlations, ξ(x, τ )ξ(x ′ , τ ′ ) = δ(x -x ′ )δ(τ -τ ′ ). The Green-Kubo relation ensures that 2D σ = βf ′′
, where f (ρ) is the free energy per unit length, which, in the SSEP, reduces to the purely entropic expression

-βf = -ρ ln ρ-(1-ρ) ln(1-ρ) [up to ρ 2 , f ′′ is the inverse isothermal compressibility, κ T = (ρ 2 f ′′ ) -1 ].
The conditions under which fluctuating hydrodynamics is valid, and the reasons why it applies to the SSEP, have been discussed by many authors. We refer the reader to Kipnis and Landim's book [START_REF] Kipnis | Scaling Limits of Interacting Particle Systems[END_REF] which specifically addresses the SSEP, and to the series of papers by Bertini et al. [START_REF] Bertini | Fluctuations in Stationary Nonequilibrium States of Irreversible Processes[END_REF][START_REF] Bertini | Macroscopic Fluctuation Theory for Stationary Non-Equilibrium States[END_REF][START_REF] Bertini | Large Deviations for the Boundary Driven Symmetric Simple Exclusion Process[END_REF][START_REF] Bertini | Current Fluctuations in Stochastic Lattice Gases[END_REF][START_REF] Bertini | Non Equilibrium Current Fluctuations in Stochastic Lattice Gases[END_REF], the latter employing the macroscopic fluctuation theory vocabulary. See also [START_REF] Tailleur | Mapping out-of-equilibrium into equilibrium in one-dimensional transport models[END_REF] for a field theory approach using coherent state construction of path integrals to represent the e -sK = e -tH(s) [the operator H(s) is defined in [START_REF] Tailleur | Mapping Nonequilibrium onto Equilibrium: The Macroscopic Fluctuations of Simp[END_REF]].

We will focus on the so-called activity K(t), which, in the SSEP, is the number of particle hops that have taken place over a given time window [0, t], over the whole ring. Up to irrelevant finite size corrections, we may view the activity as a functional of the local occupation numbers given by K

(t) = d 0 dt ′ j,± n j (t ′ )(1-n j±1 (t ′ )).
In terms of the field of occupation numbers, it is expressed by

K(t) = L 3 1 0 dx t/L 2 0
dτ σ(ρ(x, τ )). The equilibrium distribution in the SSEP is a simple Bernoulli distribution with parameter ρ 0 . Denoting by ρ 0 = N/L the average density, we thus find that, to leading order in the system size, K /(tL) = 2 N (L-N +1)

L(L-1)
≃ σ(ρ 0 ). For time realizations with activity greater than K we expect that the system remains as homogeneous as possible to favor particle hops able to accommodate a high value of K. We also expect, as we confirm below, that for low activities the system will group particles into a single cluster. We now briefly recall what the existing literature has established about these opposite regimes.

Universal fluctuations of the activity

The key to the calculation of large deviation properties in systems described by a Langevin equation ( 2) is that, since the noise becomes vanishingly small as the system size increases, a WKB-like saddle point expansion can be exploited. There are many ways to implement the saddle point expansion. One is to cast the generating function Z(s, t, L, ρ 0 ) = e -sK of the activity in the form of a path-integral over a pair of fields ρ and ρ (the latter is the actual occupation number field, while the former is conjugate to the noise ξ). The net result of that procedure that has been described many times in the past (see e.g. [START_REF] Tailleur | Mapping out-of-equilibrium into equilibrium in one-dimensional transport models[END_REF]) is that

Z = D ρDρ e -LS[ ρ,ρ;s] (3)
where the action S[ρ, ρ; s] reads

S = t/L 2 0 dτ 1 0 dx ρ∂ τ ρ + ∂ x ρ∂ x ρ - σ(ρ) 2 (∂ x ρ) 2 + sL 2 σ(ρ) (4)
Once cast in the above path-integral formulation, it is clear that, as L → ∞, the leading behavior to Z will be given by a saddle point evaluation of the path-integral.

One must then minimize the action S with respect to ρ and ρ. The latter gives us the optimal density profile able to accommodate a given value of the average activity, as tuned by s. In Appert et al. [START_REF] Appert-Rolland | Universal cumulants of the current in diffusive systems on a ring[END_REF] it was shown that the large deviation function of the activity, defined by

ψ(s) = lim t→∞ ln Z t (5)
could be written in the form

ψ(s) = -s K t + D L 2 F - σσ ′′ 8D sL 2 (6)
where F is a scaling function independent of the functions D and σ, and its argument reduces to F σ 2 sL 2 in the SSEP. The functions D and σ appearing in this formula are evaluated at ρ 0 , as in the sequel of the paper when no argument is made explicit. The hypotheses underlying this result are that the optimal profile in the path-integral formulation for the generating function Z is both stationary and uniform at the value ρ 0 . The function F (u) exhibits a branch cut along the real axis for u → (π 2 /2) -, which signals that for s > λc L 2 , λ c = π 2 /σ(ρ 0 ), the uniform profile ceases to be the one minimizing the action. The new results of the present work are devoted to a study of this sL 2 > λ c regime. We now define the parameter λ by λ = sL 2 .

It is piquant to note that F has different integral representations according to whether λ < 0 or 0 < λ < λ c . We need it on the positive side, namely

F (x) = 2x -4 √ 2x 3/2 1 -1 dy y 2 cot √ 2x 1 -y 2 .
The function F has the following limiting behaviors

F (x) ≃ 2 7/2 3π (-x) 3/2 for x → -∞ (7)
F (x) = C 1 -C 2 (π 2 /2 -x) 1/2 + . . . for x → π 2 2 - ( 8
)
where C 1 ≃ 51.61351... and C 2 ≃ 55.83091... > 0. The latter expansion explicitly displays the π 2 /2 -x singularity. The x → -∞ limiting behavior was first found by Lieb and Liniger in their study of the one-dimensional Bose gas with repulsive interactions [START_REF] Elliott | Exact Analysis of an Interacting Bose Gas. I. The General Solution and the Ground State[END_REF]. We now turn to our analysis of the inactive phase.

Inactive phase

The technical problem

Inspired by the recent work of Bodineau and Toninelli [START_REF] Bodineau | Activity phase transition for constrained dynamics[END_REF], we define a rescaled large deviation function φ(λ) ≡ lim L→∞ L lim t→∞ ln e -λK/L 2 t = lim L→∞ Lψ(λ/L 2 ). From (6) we get that φ(λ) = -λσ(ρ 0 ) which has finite size corrections of the order 1/L given by F (σλ/2). While F appears only through finite-size corrections, we note that its singularities signal the end of the regime over which the uniform saddle profile remains the optimal one. Mathematically, the quadratic expansion of the action around the uniform saddle ceases to be positive definite at λ > λ c . We therefore search for another saddle to Euler-Lagrange equations. Due to the conservation of the number of particles, we introduce a Lagrange multiplier µ enforcing the total density to be ρ 0 , by adding a µ (ρ -ρ 0 ) contribution to the action. Assuming a stationary saddle, the Euler-Lagrange equations of motion read

j = -∂ρ + σ∂ ρ, ∂ 2 ρ + σ ′ 2 (∂ ρ) 2 -λσ ′ -µ = 0 (9)
We see a posteriori that in the homogeneous regime we must have µ = -λσ ′ (ρ 0 ). The only way to find out the stability regime of the homogeneous profile is to study the quadratic action expanded in the vicinity of that saddle. Since biasing the trajectories with the activity does not break the left-right symmetry nor the time-reversal one (K is even by time-reversal), we search for a solution to [START_REF] Bertini | Current Fluctuations in Stochastic Lattice Gases[END_REF] in which the current j = 0. Then it is easy to verify that the equations in ( 9) are equivalent to Euler-Lagrange equations deduced from a Lagrangian L = (∂xρ) 2

2σ

+ λσ + (ρ -ρ 0 )µ, and were the momentum conjugate to ρ is π = ∂L ∂∂xρ = ∂xρ σ = ∂ x ρ. The corresponding Hamiltonian H = π∂ x ρ -L is a constant of motion, which means, when written in terms of the Lagrangian variables, that H given by

H = 1 2σ (∂ρ) 2 -λσ -(ρ -ρ 0 )µ ( 10
)
is independent of x. A simple rewriting of the above equation tells us that 1

2 (∂ x ρ) 2 + E P (ρ) = 0, E P (ρ) = -λσ(ρ) 2 -(ρ -ρ 0 )µ + H σ(ρ) (11)
Again, [START_REF] Beck | Thermodynamics of chaotic systems: an introduction[END_REF] has a simple mechanical interpretation: ρ is the position of a particle evolving in the force field given by the potential energy E P (ρ), with kinetic energy 1 2 (∂ x ρ) 2 , while maintaining a zero total mechanical energy. Time coordinate is x.

The two constants H and µ are for the moment undetermined. Periodic boundary conditions and particle conservation will lead to λ and ρ 0 dependent expressions for these parameters. In the homogeneous regime (λ < λ c ), we must have E P = 0, H = -λσ(ρ 0 ) and µ = -λσ ′ (ρ 0 ). In the mechanical analogy, searching for a periodic profile means that we are after a periodically oscillating solution between two extreme positions. For λ < λ c and the corresponding values of H and µ there is only a single value of ρ such that E P (ρ) = 0, which means that the oscillations are reduced to a standstill at this very density (which is of course ρ 0 ). Note that, mathematically speaking, a similar set of equations were obtained by Bodineau and Derrida [START_REF] Bodineau | Current Fluctuations in Nonequilibrium Diffusive Systems: An Additivity Principle[END_REF] in the context of the additivity principle for the large deviations of the current in a system in contact with two reservoirs (a situation in which particle number is not conserved and no chemical potential is required). We now set out to find the optimal profile ρ(x) along with H and µ in the regime where λ → λ + c .

Ρ Ρ

Close to the critical point

We use this formulation to find the optimal profile for λ = λ c + ε, with ε → 0 + . In the spirit of [START_REF] Bodineau | Distribution of current in nonequilibrium diffusive systems and phase transitions[END_REF], we assume that a small perturbation around the flat profile will develop, which leads us to write that ρ(x

) = ρ 0 + √ εφ 1 (x) + εφ 2 (x) + . . ., µ = µ c + √ εµ 1 + εµ 2 + . . ., with µ c = -π 2 σ ′ σ .
We write that H is independent of x order by order in ε, which, after a few manipulations, leads to

φ 1 (x) = 4 π [ρ 0 (1 -ρ 0 )] 3/2 cos(2πx) (12)
and

φ 2 (x) = - 1 2π 2 σ ′ (ρ 0 )σ(ρ 0 ) 2 cos(4πx) ( 13
)
together with µ 1 = 0 and µ 2 = -σ ′ (ρ 0 ) 3 /4. We then find that as λ = λ c + ε, with

ε → 0 + , φ(λ) = -π 2 -σ(ρ 0 )ε + 3 σ(ρ 0 ) 3 π 2 ε 2 + . . . (14)
We find that H = -π 2 -σ(1 -4σ)ε + . . .. In other words, the instability sets in via the largest wave-length Fourier mode, just as was the case in the system studied by Bodineau and Derrida [START_REF] Bodineau | Distribution of current in nonequilibrium diffusive systems and phase transitions[END_REF]. However, in our case, we can go deeper into the broken symmetry phase, as will now become clear.

In the inactive regime

In order to pursue our mechanical analogy, we plot E P (ρ) as a function of ρ and have the space variable x play the role of a time. A typical plot is shown in figure 1. There one can see that E P (ρ) vanishes at two densities ρ ± given by

ρ ± = 2λ + µ ± 8λ(H -µρ 0 ) + (2λ + µ) 2 4λ (15)
The particle thus oscillates between these two extreme positions, and it takes a full period

1 0 dx = 1 = 2 ρ+ ρ- dρ 1 -2E P (ρ) (16)
to do the round-trip between ρ -and ρ + . Besides, we must have that the total density is ρ 0 , which in turn imposes

ρ 0 = 1 0 dxρ(x) = 2 ρ+ ρ- dρ ρ -2E P (ρ) (17)
The two equations ( 16) together with [START_REF] Lester | Dynamic Order-Disorder in Atomistic Models of Structural Glass Formers[END_REF] thus completely fix the unknowns H and µ.

At this stage, it is more convenient to deal with the ρ 0 = 1/2 case, for which analytic expressions are somewhat simpler.

At half-filling

Particle-hole symmetry at ρ 0 = 1/2 imposes that E P (ρ) = E P (1 -ρ), which in turn forces µ = 0 (this can explicitly be verified from ( 16) and ( 17)). Then we have that

0 < ρ -= λ -λ(2H + λ) 2λ < ρ + = λ + λ(2H + λ) 2λ < 1 (18)
The density profile thus evolves between ρ -and ρ + given above, with the constant H being given by the implicit solution to

1 2 = ( √ 2H + λ - √ λ) K 2 √ λ (2H+λ) H+λ+ √ λ(2H+λ)
√ 2H (19)
where K is the complete elliptic integral of the first kind. Using that as u → 0

K(1 -u) = 2 ln 2 - 1 2 ln u + u 8 (-ln u -2 + 4 ln 2) + O(u 2 ) ( 20
)
we can solve for H as a function of λ from [START_REF] Maes | On and beyond Entropy Production: the Case of Markov Jump Processes[END_REF], which leads to

H = 8λe - √ λ 2 + 16e -2 √ λ 2 λ √ 2λ -4 + 16e -3 √ λ 2 λ(6λ -19 √ 2λ + 22) + O λ 3/2 e -4 √ λ 2 (21)
The large deviation function φ(λ) is then obtained from

φ(λ) = - 1 0 dx (∂ x ρ) 2 2σ + λσ + (ρ -ρ 0 )µ (22) = 1 0 dx H - (∂ x ρ) 2 σ (23)
Given that we have

1 0 dx (∂ x ρ) 2 σ = 2 ρ+ ρ- dρ √ -2E P σ (24)
we arrive, after some manipulations, at

φ(λ) = - √ 2 √ 2H + λ + √ λ E 2 λ(2H + λ) H + λ + λ(2H + λ) (25)
where H as a function of λ is extracted from the implicit formula [START_REF] Maes | On and beyond Entropy Production: the Case of Markov Jump Processes[END_REF]. The notations we arrive at the following large λ asymptotics for φ(λ):

φ(λ) = -2 √ 2λ -8 √ 2λe - √ λ 2 + 8 2λ -3 √ 2λ e -2 √ λ 2 + O λ 3 2 e -3 √ λ 2 (27)
This is the final result for the ρ 0 = 1/2 particular case. Translated into the language of quantum spin chains, φ(λ) in ( 25) is the ground-state energy per site of a ferromagnetic XXZ chain with anisotropy parameter ∆ = 1 + λ L 2 for λ > λ c , and at half-filling. To our knowledge, this is the first appearance of this expression in the literature. The optimal profile corresponding to that value of φ(λ) takes the form of a smooth kink with area equal to 1/2, as shown in figure 2 for increasing values of λ. Its expression reads (for x as function of ρ for 0 ≤ x ≤ 1 2 )

x(ρ) = 1 2 - √ 2 F 1 2 arccos (ρ-1)λ-H (1-ρ) √ λ(2H+λ) 2 √ λ(2H+λ) H+λ+ √ λ(2H+λ) √ 2H + λ + √ λ ( 28
)
where F is the incomplete elliptic integral of the first kind.

3.3.2. At arbitrary density ρ 0 For general ρ 0 , by exactly the same methods, though with somewhat less elegant simplifications (no symmetry enforces the Lagrange multiplier µ to vanish), the two equations ( 16) and ( 17) now take the following forms

1 2 ρ 0 = 2λ K µ+2λ(1-2 ρ) (µ+2λ(1-ρ))(1-ρ) + (µ -2λρ)Π µ+2λ(1-2 ρ) 2λ(1-ρ) µ+2λ(1-ρ) (µ+2λ(1-ρ))(1-ρ) 2λ (1 -ρ)(µ + 2λ(1 -ρ)) ( 29
)
where Π is the complete elliptic integral of the third kind, along with

1 2 = 1 (1 -ρ)(2λ(1 -ρ) + µ) K µ + 2λ(1 -2 ρ) (µ + 2λ(1 -ρ))(1 -ρ) (30)
Λ c where ρ is either equal to the minimum or the maximum density ρ ± [defined in [START_REF] Merolle | Space-time thermodynamics of the glass transition[END_REF]]. Once µ and H have been extracted from ( 29) and [START_REF] Bertini | Macroscopic Fluctuation Theory for Stationary Non-Equilibrium States[END_REF], one can determine the large deviation function φ(λ) as a function of λ and ρ 0 from the formula (see figures 3 and4)

φ(λ) = -2 (1 -ρ)(2λ(1 -ρ) + µ) E µ + 2λ(1 -2 ρ) (µ + 2λ(1 -ρ))(1 -ρ) (31)
where λ > λ c .The expansion close to λ c can be checked from those results. We see in particular that 1

L 3 t K 2 c = φ ′′ (λ) presents a jump φ ′′ (λ + c ) -φ ′′ (λ - c ) = 6σ(ρ0) 3 π 2
. We do not intend to labour beyond the leading order of φ(λ) as λ ≫ 1, which is independent of ρ 0 ,

φ(λ) ≃ -2 √ 2λ + H + (1 -ρ 0 )µ + . . . (32)
That the leading contribution is independent of ρ 0 could have been expected, as the main contribution to φ as given in [START_REF] Gaudin | La fonction d'onde de Bethe[END_REF] is the integral [START_REF] Kardar | Replica Bethe Ansatz studies of two-dimensional interfaces with quenched random impurities[END_REF], which is itself dominated by a small region of size 1/ √ λ where the derivative of ρ varies the most steeply between ρ -and ρ + , while density is close to ρ-+ρ+ 2 ≃ 1/2, and thus Lt K s is shown for a uniform profile at λ > λc is shown for comparison (dashed line). As expected, above λc, the activity for histories with a non-uniform profile is lower than for uniform profiles.

φ(λ) ∼ -1/ √ λ × 1 1/ √ λ 2 ∼ - √ λ (33)
Together with the implicit formulas ( 29) and ( 30), equation ( 31) is the ground state energy of an XXZ ferromagnetic Heisenberg chain with anisotropy parameter ∆ = 1 + λ L 2 at fixed magnetization L(2ρ 0 -1). The density profile is, for x as function of ρ for 0 ≤ x ≤ 1 2 (see figure 5),

The infinitely inactive regime

In the infinitely inactive regime, that is for λ ≫ L 2 , the density profile described in the previous subsections becomes infinitely steep and the assumptions underlying fluctuating hydrodynamics fail. Since the typical slope of the rise of the profile towards its plateau value increases at √ λ, fluctuating hydrodynamics must fail when the slope exceeds L, that is when profile variations become significant at the lattice scale, which is consistent with the λ ≫ L 2 condition. For large s it is easy to see that the ground state of the s-modified evolution operator (1) reads

ψ(s) = -2Lσ(ρ 0 )(1 -e -s ) + O(e -2s ) (35)
The above regime consistently connects to the calculation of the previous subsection where it was shown that φ(λ) ∼ -√ λ for large values of λ. The validity of that previous result holds until φ(λ) ∼ -2Lσ(ρ 0 ), when it connects with the infinitely inactive regime. This occurs as λ ∼ L 2 and confirms once more that the crossover is at a typical λ of the order of L 2 .

Prospects

The SSEP on a ring displays a phase transition from a homogeneous state to a kinklike profile as the trajectories' activity is lowered below its average value. We have described both the details of the kink profile that develops and the large deviation function of the activity itself in this nonuniform dynamical state. In the course of our analysis, we have identified three distinct finite-size scaling regimes. The system dealt with in the present work exhibits a second order continuous phase transition. It would certainly be quite desirable to find similar solvable examples for first-order transitions. The mechanism by which the symmetry breaking sets in may be more subtle. This is work in progress.

Figure 1 .

 1 Figure 1. The potential energy profile, E P (ρ), is shown for a set of typical values of the parameters λ, H, µ and ρ 0 .

Figure 2 .

 2 Figure 2. The density profile is shown for λ = λc + 1 10 (red, large dashes), 25 (blue), 200 (purple, small dashes) at average density ρ 0 = 1 2 . The kink becomes steeper as λ is increased. The slope of the rise to the plateau increases as √ λ for large λ.

K

  and E refer to the complete elliptic integrals of the first and second kind, respectively. This expression is valid for any λ > λ c . Using the following asymptotic formula for E as u → 0 E(1 -u) = 1 + uu -13 + 24 ln 2) + O(u 3 )[START_REF] Dotsenko | Bethe ansatz derivation of the Tracy-Widom distribution for one-dimensional directed polymers[END_REF]

Figure 3 .

 3 Figure 3. Large deviation function φ(λ) (thick line) as a function of λ, at average density ρ 0 = 23 . The function is singular at λ = λc = π 2 σ(ρ 0 ) (the second derivative is discontinuous), marking the entrance in the regime where the density profile is non-uniform (λ > λc). The value of the large deviation function for a uniform profile at λ > λc is shown for comparison (dashed line).

Figure 4 .

 4 Figure 4. Mean activity 1Lt K λ = -φ ′ (λ) (thick blue line) as a function of λ, at average density ρ 0 = 2 3 . The function is singular at λ = λc = π 2 σ(ρ 0 ) . The value of 1

  2λρ-)(ρ-ρ-) (2(2ρ--1)λ-µ)(1-ρ) 2λ(1-2ρ-)+µ (µ-2λρ-)ρρ -(2λρ --µ)(34)

Figure 5 .

 5 Figure 5. The density profile is shown for λ = λc + 1 10 (red, large dashes), 25 (blue), 200 (purple, small dashes) at average density ρ 0 = 2 3 (thin dashes). The kink becomes steeper as λ is increased.

Acknowledgements

We would like to thank Thierry Giamarchi for interesting discussions on XXZ related matters. This work was supported in part by British Council France Alliance Project 09.013.


Jérôme Didier

Béatrice Faverjon

Jean-Jacques Sinou
email: jean-jacques.sinou@ec-lyon.fr

Analyzing the dynamic response of a rotor system under uncertain parameters by Polynomial Chaos Expansion

Keywords: dynamics, rotor, uncertainties

In this paper, the quantification of uncertainty effects on response variability in rotor systems is investigated. To avoid the use of Monte Carlo simulation (MCS), one of the most straightforward but computationally expensive tools, an alternative procedure is proposed. Monte Carlo Simulation builds statistics from responses obtained from sampling uncertain inputs by using a large number of runs. However, the method proposed here is based on the stochastic finite element method (SFEM) using polynomial chaos expansion (PCE).

The efficiency and robustness of the method proposed is demonstrated through different numerical simulations in order to analyze the random response against uncertain parameters and random excitation to assess its accuracy and calculation time.

Introduction

In rotordynamics, two types of uncertainty on dynamic systems are of particular interest. The first of these can derive from variations in mechanical properties (such as mass, stiffness and geometrical imperfections) due to manufacturing errors [Lalanne and Ferraris (1990); [START_REF] Friswell | Finite Element Model Updating in Structural Dynamics[END_REF]; [START_REF] Erich | Handbook of Rotordynamics[END_REF]; [START_REF] Childs | Turbomachinery Rotordynamics: Phenomena, Modeling, and Analysis[END_REF]; Yamamoto and Ishida (2001)]. Besides this type of structural uncertainties, external and internal forcing functions can also be random.

Numerous methods have been used to quantify physical uncertainties in a variety of computational problems like the perturbation method, the Monte Carlo Simulations and the Polynomial Chaos Expansion [Ghanem and Spanos (1991)]. The perturbation method based on the expansion of random quantities into Taylor series [Nayfeh (1973)] and the Neumann method based on Neumann series expansion [START_REF] Benaroya | Finite element methods in probabilistic structural analysis: A selective review[END_REF]; Yamazaki et al. (1988)] provide acceptable results for small random fluctuations, then they are not adapted here for solving a dynamic problem in the case of an excitation frequency close to a resonance frequency. The direct Monte Carlo Simulations, which is the most straightforward and frequently used approach, is adapted to include uncertainties in a deterministic finite element model, by generating n independent samples of the random parameter. Then it leads to solve the deterministic problem n times in order to obtain n samples of the response vector and so the statistics of the response. Due to the slow convergence rate of this method, a very high number of samples is necessary then if solving the deterministic problem is already computationally intensive, the computational costs of the method can become impractical. Particularly, rotordynamics problems are quite complex to solve in a deterministic sense. The polynomial chaos expansion associated with a Galerkin projection so-called stochastic finite element method [Ghanem and Spanos (1991)] has shown to be a successful approach to solve uncertainty quantification problems. It represents the stochastic processes and variables in a set of orthogonal bases of random variables. The polynomial chaoses are from the homogeneous chaos theory of Wiener [Wiener (1938)] and the original polynomial chaos expansion [Ghanem and Spanos (1991)] used a mean-square convergent expansion as multidimensional Hermite polynomials of normalized Gaussian variables. Since the Hermite polynomials are orthogonal with respect to the Gaussian measure, the homogeneous polynomial chaos can achieve optimal exponential convergence for Gaussian inputs [START_REF] Ghanem | Stochastic finite elements with multiple random non-gaussian properties[END_REF]]. This last method then seems the most adapted to study the influence of the uncertainties on the parameters of the rotor structure on the response.

The present paper is organized as follows: firstly, we present the rotor system after which a brief explanation is given of the Stochastic Finite Element Method [Ghanem and Spanos (1991)] for the solution of mechanical problems with several random characteristics. Secondly, expansions of the operator of random material properties and of the random external forcing function on the chaos basis are explained and studied for application to the rotordynamics problem. The Polynomial Chaos Expansion procedure is illustrated by different numerical examples that include the most common sources of randomness in a rotordynamics problem (such as physical and geometric parameters). Thirdly, the results obtained by applying the Polynomial Chaos Expansion (PCE) procedure are compared with those evaluated by Monte Carlo Simulation (MCS) whose costs become prohibitive for large finite element models with large numbers of design parameters. Finally, the efficiency and robustness of the method proposed is demonstrated through several numerical simulations of the effects of uncertainties and orders of polynomial chaos.

Rotor Model

The system under study is illustrated in Figure 1. The rotor consists of a rotor shaft with two discs. The shaft is discretized into 10 Timoshenko beam finite elements with four degrees of freedom at each node [Lalanne and Ferraris (1990); [START_REF] Friswell | Finite Element Model Updating in Structural Dynamics[END_REF]] and a constant circular section. All the values of the parameters are given in Table 1. The beam element model is given by

[M e ]{ ẍe } + ([C e ] + ω[G e ]){ ẋe } + [K e ]{x e } = {q e } (1)
where the vector {q e } defines forces applied on the shaft and ω is the rotational speed of the shaft.

[M e ] and [G e ] are the mass and gyroscopic matrices of the shaft element.

[K e ] and [C e ] are the elementary stiffness and damping matrices. These matrices are described in A.

The model of the rigid discs is given by

[M d ]{ ẍd } + ω[G d ]{ ẋd } = {q d } (2)
where [M d ] and [G d ] are the mass and gyroscopic matrices of the disc. These matrices will be described in the following part of the paper. The vector {q d } defines the unbalance forces due to an eccentric mass. For the degree-of-freedom [v w θ ψ] T (see Figure 1), the unbalance forces are given by

q d = [m e d e ω 2 cos(ωt + φ) m e d e ω 2 sin(ωt + φ) 0 0] T (3)
where m e and d e are the mass unbalance and the eccentricity respectively. φ and ω define the initial angular position in relation to the z-axis and the rotational speed of the rotor. Finally, discrete stiffness components are located at either end of the shaft. After assembling the shaft elements and the rigid discs, the equation of motion for the complete rotor system is defined as follows:

[M] { x} + ([C] + ω[G]) {x} + ([K] + [K b ]){x} = {q} (4)
where [M] and [G] are the mass and gyroscopic matrices of the shaft and the two discs.

[K] and [C] are the stiffness and damping matrices of the shaft. [K b ] is the stiffness matrix of the bearings. {q} is the unbalance forces of the complete rotor system. Considering that the unbalance force can be written as {q} = {Q}e iωt , the response vector may be assumed to be {x} = {X}e iωt . By using Equation (4), the system governing the equation in the frequency domain is given by

[A(ω)]{X(ω)} = {Q(ω)} (5)
where

[A(ω)] = -ω 2 [M] + iω([C] + ω[G]) + [K] + [K b ] (6)
In the following part of the paper, the frequency dependence will be omitted to simplify the notation and [A(ω)] will be noted as [A].

Stochastic model

In rotordynamics, uncertainties on dynamic responses can occur due to manufacturing inaccuracies related to mechanical properties such as mass, stiffness, damping and geometrical imperfections. Moreover, forcing functions (external and internal) also lead to considerable uncertainties, so that they have to be taken as random quantities. Therefore, considering in this part [A], {X} and {Q} as random processes, where argument τ denotes the random character, Equation (5) can then be rewritten in a random way such that 6) becomes

[A(τ )]{X(τ )} = {Q(τ )} (7)
[A(τ )] = -ω 2 [M(τ )] + iω([C(τ )] + ω[G(τ )]) + [K(τ )] + [K b (τ )] (8)

The system equation on the Polynomial Chaos basis

The Polynomial Chaos basis is a set of orthogonal bases of random variables, represented in a meansquare convergent expansion by multidimensional Hermite polynomials of normalized Gaussian variables. A rapid overview on the construction of the basis is given below and more details can be found in [Ghanem and Spanos (1991)]. In the following of the paper, the random behavior of each physical or geometrical quantity M (scalar or matrix) considered would be sufficiently modeled by using the Karhunen-Loeve expansion implemented in the Galerkin formulation of the finite element method [Ghanem and Spanos (1991)], then we can expand M such as

M = M + L l=1 ξ l M l (9)
where {ξ 1 , ...ξ L } is a set of orthonormal random variables, M is the mean of M and M l is its lth expansion term. The response can be expanded on the Polynomial Chaos basis such that where Ψ j (ξ(τ )) refers to a rearrangement of the p-order finite dimensional orthogonal polynomials in relation to the Gaussian function that forms a complete basis in the space of second-order random variables ; {X} j is the unknown deterministic jth vector associated with Ψ j (ξ(τ )) and ξ = {ξ r } [Ghanem and Spanos (1991)]. Finally, the system to be solved, when expanded on the polynomial chaos basis, is

{X(τ )} = ∞ j=0 {X} j Ψ j (ξ(τ )) (10)
∞ j=0 [A(τ )]{X} j Ψ j (ξ(τ )) = {Q(τ )} (11)
with random quantities [A(τ )] and [Q(τ )] defined by

[A(τ )] = ∞ i=0 [A] i Ψ i (ξ(τ )) (12)
with

[A] i = -ω 2 [ M] i + iω([ Ĉ] i + ω[ Ĝ] i ) + [ K] i + [ Kb ] i i = 0, 1, ..., ∞ (13)
and

[Q(τ )] = ∞ k=0 { Q} k Ψ k (ξ(τ )) (14)
where quantity Q denotes the rearrangement of quantity Q in the polynomial chaos basis. Details of the rearrangement will be given in the following section of the paper. It should be noted that each random expansion will be described in Subsection 3.2. Finally, the system to be solved to the subspace spanned by

{Ψ k } ∞ k=0 , is ∞ i=0 [A] i Ψ i (ξ(τ )) ∞ j=0 {X} j Ψ j (ξ(τ )) = ∞ k=0 { Q} k Ψ k (ξ(τ )) (15)
which, after doing a Galerkin projection on the polynomial chaos basis, can also be rewritten as

∞ i=0 ∞ j=0 E{Ψ i Ψ j Ψ k }[A] i {X} j = { Q} k E{Ψ 2 k } k = 0, 1, ..., ∞ (16)
where E{} denotes the operation of mathematical expectation. It should be noted that coefficients E{Ψ i Ψ j Ψ k } and E{Ψ 2 k } only have to be calculated once. In practice, the expansion can be truncated after the P th term where P is the total number of polynomial chaoses used in the expansion excluding the 0th order term and can be determined by

P = 1 + p s=1 1 s! s-1 r=0 (L + r) (17)
and in which p is the order of homogeneous chaos used.

The random quantities in the stochastic rotor

There are different sources of randomness in the rotordynamics problem studied due to geometric and material parameters. This paper deals with uncertainties modeled by Gaussian random variables that represent the random character of the parameters, such as the Young modulus of the shaft, bearing stiffness, disc diameter and density and the amplitude of the unbalance force. All these random quantities are modeled by using expansion defined in Equation ( 9) and, for physically strictly positive parameters, the random variables of negative values have been removed.

Stiffness of the shaft

From the random character of the Young modulus modeled by one truncated Gaussian random variable ξ 1 , thus, from expansion described in Equation ( 9), we obtain the relation

E(τ ) = E(1 + δ E ξ 1 ) (18)
where E and δ E are respectively the mean value and the variation coefficient of the Young modulus.

The detailed deterministic expression of the elementary stiffness matrix of the shaft [K e ] as a function of the Young modulus is described in A. By introducing the random Young modulus defined above and after assembling the elementary stiffness matrices of the shaft, it is easy to find the random expansion of

[K] [K] = [K] 0 + ξ 1 [K] 1 (19) (20) j p Ψ j E{Ψ 2 j } 0 0 1 1 1 1 ξ 1 1 2 ξ 2 1 3 ξ 3 1 4 ξ 4 1 5 ξ 5 1 6 ξ 6 1 7 2 ξ 2 1 -1 2 8 ξ 1 ξ 2 1 9 ξ 1 ξ 3 1 10 ξ 1 ξ 4 1 11 ξ 1 ξ 5 1 12 ξ 1 ξ 6 1 13 ξ 2 2 -1 2 14 ξ 2 ξ 3 1 15 ξ 2 ξ 4 1 16 ξ 2 ξ 5 1 17 ξ 2 ξ 6 1 18 ξ 2 3 -1 2 19 ξ 3 ξ 4 1 20 ξ 3 ξ 5 1 21 ξ 3 ξ 6 1 22 ξ 2 4 -1 2 23 ξ 4 ξ 5 1 24 ξ 4 ξ 6 1 25 ξ 2 5 -1 2 26 ξ 5 ξ 6 1 27 ξ 2 6 -1 2
Table 2: Six-Dimensionnal Polynomial Chaoses and their variance i j

j E{Ψ i Ψ j Ψ k } 1 1 7 2 1 2 8 1 1 3 9 1 1 4 10 1 1 5 11 1 1
E{Ψ i Ψ j Ψ k }, E{Ψ i Ψ j Ψ k } = E{Ψ j Ψ i Ψ k } = E{Ψ i Ψ k Ψ j }, Six-Dimensional Polynomial Chaoses, Chaos order 2
Similarly, the hysteretic damping is defined by

[C e ] = η ω [K e ] ( 21
)
where η is the hysteretical damping factor. Assembling elementary damping matrices of the shaft yields:

[C] = [C] 0 + ξ 1 [C] 1 (22)
Bearing stiffness For the shaft corresponding to the degree-of-freedom [v w θ ψ] T , the deterministic elementary stiffness matrix of the bearing is defined as

[K e b ] =     k 1x 0 0 0 k 1y 0 0 0 0 sym 0     (23)
where k 1x and k 1y are the stiffnesses of the bearing in directions x and y. In this case, it has been chosen to only investigate the randomness of the stiffness on the first bearing in direction x modeled by the truncated Gaussian random variable ξ 2 , which yields:

k 1x (τ ) = k 1x (1 + δ k 1x ξ 2 ) (24)
in which k 1x and δ k 1x are the mean value and the variation coefficient of the stiffness. Finally, the expression of the assembled stiffness matrix [K b ] is written as

[K b ] = [K b ] 0 + ξ 2 [K b ] 1 (25)

Disc parameters

The parameters of the discs should be random. Here, we consider the randomness on the outer diameter D(τ ) and the density ρ(τ ) of one of the discs, which are the geometric and material parameters of the model. Describing them by using two truncated Gaussian random variables ξ 3 and ξ 4 yields

D(τ ) = D(1 + δ D ξ 3 ) (26) ρ(τ ) = ρ(1 + δ ρ ξ 4 ) (27)
where D and ρ are the mean values, δ D and δ ρ are the variation coefficients of the diameter and the density of the disc respectively. These quantities appear in the definition of the mass and gyroscopic elementary matrices [M d ] and [G d ] that are expressed for the disk relative to the degree of freedom [v w θ ψ] T such that

[M d ] =     m d (τ ) 0 0 0 0 m d (τ ) 0 0 0 0 I d (τ ) 0 0 0 0 I d (τ )     , [G d ] =     0 0 0 0 0 0 0 0 0 0 0 -I p (τ ) 0 0 I p (τ ) 0     (28) with m d (τ ) = 1 4 ρ(τ )πh(D(τ ) 2 -d 2 ) (29) I d (τ ) = 1 64 ρ(τ )πh(D(τ ) 4 -d 4 ) + 1 48 ρ(τ )πh 3 (D(τ ) 2 -d 2 ) (30) I p (τ ) = 1 32 ρ(τ )πh(D(τ ) 4 -d 4 ) (31)
in which h and d are the thickness and the inner diameter of the disc respectively. In addition, m d is the mass of the disk. I d and I p are the diametral moment of inertia about any axis perpendicular to the rotor axis and the polar moment of inertia about the rotor axis. Substiting Equations (26-27) in Equations ( 29) to (31) leads to the expression of the components of the mass and gyroscopic elementary matrices such that

m d (τ ) = 4 j=0 1 i=0 m dij ξ j 3 ξ i 4 I d (τ ) = 4 j=0 1 i=0 I dij ξ j 3 ξ i 4 I p (τ ) = 4 j=0 1 i=0 I p ij ξ j 3 ξ i 4 ( 32
)
where m dij , I dij and I p ij are given in B. The Polynomial Chaos Expansion for m d (τ ), I d (τ ) and I p (τ ), constructed for two random variables ξ 3 and ξ 4 is given by

m d (τ ) = N j=0 m dj Ψ j (ξ 3 , ξ 4 ) I d (τ ) = N j=0 I dj Ψ j (ξ 3 , ξ 4 ) I p (τ ) = N j=0 I p j Ψ j (ξ 3 , ξ 4 ) (33)
in which m dij , I dij and I p ij are determined after identification between Equations ( 32) and (33) using Tables 4 and5. The number of polynomial chaoses N is deduced from Equation ( 17) by two random variables: L = 2. It should be noted that this identification yields a minimum order p = 5. Then, expressing mass and gyroscopic elementary matrices of a disc [M d e ] and [G d e ] in the polynomial chaos basis yields

[M d ] = P j=0 [M d ] j Ψ j (ξ 3 , ξ 4 ) [G d ] = P j=0 [G d ] j Ψ j (ξ 3 , ξ 4 ) (34)
where

[M d ] j =     m dj 0 0 0 0 m dj 0 0 0 0 I dj 0 0 0 0 I dj     , [G d ] j =     0 0 0 0 0 0 0 0 0 0 0 -I p j 0 0 I p j 0     (35)
Finally, the assembled mass and gyroscopic matrices are given by

[M] = P j=0 [M] j Ψ j (ξ 3 , ξ 4 ) and [G] = P j=0 [G] j Ψ j (ξ 3 , ξ 4 ) (36) j p Ψ j E{Ψ 2 j } 0 0 1 1 1 1 ξ 3 1 2 ξ 4 1 3 2 ξ 2 3 -1 2 4 ξ 3 ξ 4 1 5 ξ 2 4 -1 2 6 3 ξ 3 3 -3ξ 3 6 7 ξ 2 3 ξ 4 -ξ 4 2 8 ξ 1 3 ξ 2 4 -ξ 3 2 9 ξ 3 4 -3ξ 4 6 10 4 ξ 4 3 -6ξ 2 3 + 3 24 11 ξ 3 3 ξ 4 -3ξ 3 ξ 4 6 12 ξ 2 3 ξ 2 4 -ξ 2 4 -ξ 2 3 + 1 4 13 ξ 3 ξ 3 4 -3ξ 3 ξ 4 6 14 ξ 4 4 -6ξ 2 4 + 3 24
Table 4: Two-Dimensionnal Polynomial Chaoses and their variance

1 Ψ 0 ξ 3 Ψ 1 ξ 4 Ψ 2 ξ 2 3 Ψ 3 + Ψ 0 ξ 3 ξ 4 Ψ 4 ξ 2 4 Ψ 5 + Ψ 0 ξ 3 3 Ψ 6 + 3Ψ 1 ξ 2 3 ξ 4 Ψ 7 + Ψ 2 ξ 3 ξ 2 4 Ψ 8 + Ψ 1 ξ 3 4 Ψ 9 + 3Ψ 2 ξ 4 3 Ψ 10 + 6Ψ 3 + 3Ψ 0 ξ 3 3 ξ 4 Ψ 11 + 3Ψ 0 ξ 2 3 ξ 2 4 Ψ 12 + Ψ 5 + Ψ 3 + Ψ 1 ξ 3 ξ 3 4 Ψ 13 + 3Ψ 4 ξ 4 4 Ψ 14 + 6Ψ 5 + 3Ψ 0 Table 5: Identification

Excitation characteristics

The unbalance force due to an eccentric mass on a disk can be written on the degree of freedom

[v w θ ψ] T {Q} = m b r b ω 2 e iφ [1 -i 0 0] T (37)
where m b and r b are the unbalance mass and the eccentricity respectively. Furthermore, φ defines the initial angular position. Parameters m b and φ are considered as random Gaussian type quantities and are defined as

m b (τ ) = m b (1 + δ m ξ 5 ) (38) φ(τ ) = σ φ ξ 6 ( 39
)
with m b and δ m being the mean value and the variation coefficient of the unbalance mass, and σ φ the standard deviation of the angular position of the force. For the reader comprehension, the angular position is illustrated in Figure 2. In addition, ξ 5 and ξ 6 are Gaussian random variables. (iσ φ ξ 6 ) j j! (41)

Thus the unbalance force due to an eccentric mass on a disk can be given by

{Q} = m b r b ω 2 (1 + δ m ξ 5 ) M j=0 (iσ φ ξ 6 ) j j! [1 -i 0 0] T (42)
Equation ( 42) can be rewritten as

{Q} = 1 k=0 M j=0 {Q} kj ξ k 5 ξ j 6 ( 43
)
where

{Q} kj = m b r b ω 2 (δ m ) k (iσ φ ) j j! [1 -i 0 0] T (44)
Finally, the random loading {Q} can be expanded on the polynomial chaos basis as follow

{Q} = R j=0 {Q} j Ψ j (ξ 5 , ξ 6 ) (45)
where the deterministic coefficients {Q} j are given by using the same identification process described for the mass and gyroscopic matrices, adapted here to coefficients ξ i 5 ξ j 6 and Ψ j (ξ 5 , ξ 6 ). R is the number of polynomial chaoses given by Equation ( 17) for L = 2 and depending on the equivalent M to the p = M + 1 order. Synthesis Finally, the system to be solved, given by Equation ( 16), is expanded on six-dimensional polynomial chaoses Ψ j (ξ) with ξ = {ξ 1 , ..., ξ 6 }, j = 0 to P where P is defined by Equation ( 17) with

L = 6. Thus [A] i is a function of [ M] i , [ Ĉ] i , [ Ĝ] i , [ K] i and [ Kb ] i (see Equation (13)) which refer to a rearrangement of [M] j , [C] j , [G] j , [K] j , [K b ] j
and {Q} j on a six-dimensional polynomial chaos basis.

Numerical results

In this section, the quantification of the uncertainty effects on the response variability of the rotor under study are presented using the Polynomial Chaos Expansion method. To show lower and higher dynamic responses of the rotor system under uncertain parameters, the stochastic response of the rotor system is proposed via the mean value and the variance of the random response, also represented graphically by an envelope. The envelope of the stochastic response is constructed by calculating the maximum and the minimum of all the responses computed by the PCE approach for samples generated by the MCS method. Then, the MCS method generates n values of Hermite polynomials and consequently n samples of the Frequency Response Functions.

Finally, the envelope is built by considering the maximum and the minimum of all the samples. In the following, the sections are organized as follows: firstly, the main dynamic characteristics are investigated in the deterministic case. The efficiency and robustness of the Polynomial Chaos Expansion method is then discussed for the dynamic response of the rotor system under uncertain parameters. Finally, the mean and the variance of the Frequency Response Function obtained via the PCE approach, and the envelope are compared with results obtained by using the Monte Carlo simulation.

δ E δ k 1x δ ρ δ D δ m σ φ Order Case 1 5% 5% - - - - 2 Case 2 - 10% - - - - 2 Case 3 - 10% - - - - 10 Case 4 - - - - 1%

Deterministic case

Before discussing the effects of uncertainties on the dynamic of the rotor system, a brief summary is given of the main dynamic characteristics of the deterministic rotor system. Considering the model parameters given in Table 1, Figure 3 shows the horizontal steady-state responses of the rotor for each of its transversal nodes.

It can be seen that the horizontal displacements indicate the presence of the first, second and third forward critical speeds around 28.3 Hz, 97.2 Hz and 240 Hz, respectively. To facilitate comprehension, the first, second and third backward critical speeds do not appear on the unbalance responses due to the fact that the bearing stiffnesses are identical in the vertical and horizontal directions. Table 7 summarizes the values of the three first forward and backward critical speeds of the rotor system.

Comparisons between the Polynomial Chaos Expansion approach and Monte Carlo simulation

In this part of the paper, the results of the Monte Carlo simulation and those of the Polynomial Chaos Expansion method are compared in order to validate the efficiency of the second approach.

Comparisons are given for the set of parameters defined by Case 1 in Table 6: to simulate the variation The Monte Carlo analysis is carried out to obtain a statistical sample of the random response. As explained previously, this method requires a large number of samples to provide a reference solution. In this study, to obtain convergence of the Monte Carlo method, 1000 samples were used. For this first case, the order of chaos equals 2. A convergence study, presented soon after, will justify the choice of this truncation.

Figures 4 and5 show the mean and the variance of the Frequency Response Function (at the node 2 in the direction x) obtained by the two methods. The two methods yield quasi-identical results for both quantities (a very low discrepancy can only be seen on peak around 100 Hz) which validates the PCE method.

Figure 6 shows the results for both the Monte Carlo simulations and the PCE method. All the Frequency Response Functions samples obtained by using the Monte Carlo simulations can be seen at node 2 in the horizontal direction, with their mean value. We can see that the mean of the Frequency Responses Function and so the envelopes built from the Monte Carlo simulations and the Polynomial Chaos Expansion are very close one to the other. It should be noted that the same samples have been used for both direct Monte Carlo method and PCE approach. Moreover, for this example, if we compare the CPU time, it appears that PCE approach is eight times faster than the direct Monte Carlo approach.

The variations of the Young modulus E of the shaft and the horizontal bearing stiffness k 1x can be seen to cause small changes in the critical speeds. Moreover, increases (and decreases) of the maximum amplitudes when the rotor passes through the forward critical speeds are also indicated.

Finally, it can be seen that backward critical speeds (at 61.8 Hz and 128.1 Hz) can occur due to the randomness of the bearing stiffness k 1x which introduces dissymmetry in the rotor system.

A convergence study of the PCE with the order of chaos is performed through two choices of orders : 2 and 10. Figures 7-8 and 9-10 present the results for cases 2 and 3 (10% of the variation for the horizontal bearing stiffness k 1x on the right side of the rotor system. See Table 6) for the Frequency Response Functions at node 2 in the horizontal direction. As expected, the order of chaos improves modeling. However, the effect of this discrepancy between both expansions does not harm the quality of the model, especially when taking into account the increase of computation costs (which depends on the size of the problem as a function of the order of chaos) obtained subsequently. Even if the effect of damping is not on the scope of the study, it can be mentionned that decreasing the damping factor needs to increase the PCE order, especially close to the resonances. For more details, the reader is referred to the research of [START_REF] Dessombz | Analyse dynamique de structures comportant des paramètres incertains[END_REF]].

Effects of uncertainties in mechanical properties and external forces

In this part of the paper, the effects of uncertainties from stiffness properties of the rotor, geometric parameters of the disc and external forces are investigated. In the next paragraphs, only the mean value and the envelope of the FRF are studied in order to highlight the results more clearly. It should be noted that in the following parts all the results are given at node 2 in direction x.

Excitation

The effects of uncertainties in the external forcing functions are now studied. Variations on the mass unbalance and the angular position are considered. Two cases are studied (cases 4 and 5): the first and second cases deal with 1% and 5% of uncertainties for both the mass unbalance and the angular position (m ± δ m and φ ± δ φ ).

Figures 11 and12 illustrate the mean values (using the Monte Carlo simulation and the Polynomial Chaos method) and the envelope. In this particular case, the random quantities are only located at the loading {Q}. Therefore, theoretically speaking, the mean of the response and the values of the critical speed obtained by the Polynomial Chaos method must be identical to the reference mean response. In Figure 11, the mean of the Polynomial Chaos approach is close to the reference mean response: the error between these two results is only due to the truncated expansion of the loading expression defined in Equation (41).

It can then be seen that the variations of the maximum amplitudes due to uncertainties for all the critical speeds are not very great. The amplitudes of the first, second and third critical speeds increase from only 2.814 × 10 -3 m, 2.803 × 10 -3 m and 3.76 × 10 -3 m to 2.954 × 10 -3 m, 3.534 × 10 -3 m and 3.987 × 10 -3 m respectively.

Uncertainties in disc properties

In this paragraph, variations for the properties of the disc located at the left side of the rotor system are considered. Two cases are investigated (cases 6 and 7): firstly, the density and the diameter of the disc are allowed to undergo 1% variations (ρ ± δ ρ and D ± δ D ). Secondly, 5% variations on the same set of parameters are introduced. As explained previously in Section 3.2, these random parameters affect the mass and gyroscopic matrices (see Equation ( 28)). The order of chaos has been chosen as equal to 2. Figures 13 and14 give the mean value of the Frequency Response Function and the envelope at node 2 in the horizontal direction. It appears that the mean values calculated by applying the Monte Carlo simulations and the Polynomial Chaos method are very similar. It should be noted that an order 2 gives accurate results in spite of the fact that the development given in Section 3.2 shows that an order 5 is needed to take all terms into account.

It can be seen that increasing uncertainties on the density and the diameter of the disc can dras- tically affect the values of the critical speeds and the associated maximum amplitudes. Even if the maximum variations of amplitudes are located at the critical speed, the evolutions of the rotor response far from the critical speed are significant. When comparing Figures 6 and 14, it may be concluded that the effects of uncertainties on disc properties are greater than those on shaft properties of the rotor under study.

Uncertainties in both mechanical properties and external forces

In order to demonstrate the efficiency and accuracy of the Polynomial Chaos procedure described above, this last part of the paper treats the cases in which uncertain quantities come from all the parameters studied previously (i.e. stiffness properties of the rotor, geometric parameters of the disc and external forces). Numerical simulations are given by considering the variations of mechanical properties of the shaft (i.e. the Young modulus E and the horizontal bearing stiffness k 1x ), the properties of the disc (i.e. the density ρ and the diameter D), and the excitation forces (i.e. the mass unbalance m and and the angular position φ), as indicated in Table 6 for cases 8 and 9. We recall that in this case the cost of calculation may be high since it is directly linked to the number of polynomials and consequently to the order of chaos and the number of random parameters. for cases 8 and 9 through the mean and the variance obtained from the Polynomial Chaos and the Monte Carlo approaches. The results from both methods are in very good agreement for both the mean and for the variance. Figures 17 and 20 illustrate all the Frequency Response Function samples (at node 2 in the horizontal direction) obtained by using the Monte Carlo simulations, and the lower and upper envelopes built with the Polynomial Chaos method. It appears that increasing uncertainties affects the maximum amplitudes of the dynamic response and the value of the critical speeds. Then, as explained previously in section 4.2, the dissymmetry due to the variations in the bearing stiffness k 1x leads to increases in the dynamic response of the rotor system around the backward critical speeds (at 61.8Hz and 128.1Hz). Finally, it can be seen that the Polynomial Chaos method may over-estimate vibrational amplitude (see for example the dynamic response of the rotor around the third critical speed, at 240Hz). However, whatever the levels and different kinds of uncertainty (such as material, geometrical and loading characteristics) presented here, the Polynomial Chaos method agrees very well with the Monte Carlo simulation, thereby demonstrating the robustness of the method.

Conclusion

This paper described a numerical procedure using the Chaos Polynomial approach to evaluate the stochastic response of a rotor system with uncertain mechanical parameters and uncertain external forces. It explained how this kind of problem can be solved with the Spectral Finite Element Method and how the random parameters can be modeled by random variables through a Karhunen-Loeve expansion. The results obtained by applying the Polynomial Chaos Expansion (PCE) procedure were compared with those evaluated by the Monte Carlo Simulation (MCS).

The stochastic response of the rotor system was proposed via the mean value and the variance of the random response and also represented graphically by an envelope. This envelope could be useful for designing rotor systems and predicting their lower and higher dynamic responses under uncertain parameters.

The efficiency and robustness of the Polynomial Chaos method were tested and validated through numerical simulations of the effects of uncertainties and orders of polynomial chaos.

Figure 1 :

 1 Figure 1: Rotor system and Shaft finite element

Figure

  Figure 2: Unbalance model

Figure 3 :

 3 Figure 3: Frequency Response Functions for the determinism case

Figure 4 :

 4 Figure 4: Mean of Frequency Response Functions (Case 1); Polynomial Chaos method (red dotteddashed line); Monte Carlo Simulation (black line)

Figure 5 :

 5 Figure 5: Variance of Frequency Response Functions (Case 1); Polynomial Chaos method (red dotteddashed line); Monte Carlo Simulation (black line)

Figure 6 :

 6 Figure 6: Frequency Response Functions (Case 1); Mean of the FRF with Polynomial Chaos method (red dotted-dashed line); Lower and upper envelopes (red dashed line); Mean of the FRF with the Monte Carlo Simulation (black solid line); Monte Carlo samples (grey solid line)

Figure 7 :

 7 Figure 7: Study with chaos order 2 : Mean of Frequency Response Functions (Case 2); Polynomial Chaos method (red dotted-dashed line); Monte Carlo Simulation (black line)

Figure 8 :

 8 Figure 8: Study with chaos order 2 : Variance of Frequency Response Functions (Case 2); Polynomial Chaos method (red dotted-dashed line); Monte Carlo Simulation (black line)

Figure 9 :

 9 Figure 9: Study with chaos order 10 : Mean of Frequency Response Functions (Case 3); Polynomial Chaos method (red dotted-dashed line); Monte Carlo Simulation (black line)

Figure 10 :Figure 11 :Figure 12 :Figure 13 :Figure 14 :Figure 17 :Figure 19 :Figure 20 :

 1011121314171920 Figure 10: Study with chaos order 10 : Variance of Frequency Response Functions (Case 3); Polynomial Chaos method (red dotted-dashed line); Monte Carlo Simulation (black line)

Table 1 :

 1 Model parameters

	Parameters	dimension
	Length of shaft	1 m
	Diameter of shaft	0.04 m
	Position of disc 1	0.6 m
	Position of disc 2	0.8 m
	Outer diameter of disc 1	0.2 m
	Outer diameter of disc 2	0.4 m
	Inner diameter of discs 1 and 2	0.04 m
	Thickness of discs 1 and 2	0.02 m
	Young modulus of elasticity E	2.1 10 11 Nm 2
	Shear modulus G	8.0 10 10 Nm 2
	Poisson ratio ν	0.3
	Density ρ	7800 kg m -3
	Mass unbalance	0.05 g
	Phase unbalance	0 •
	Eccentricity of the mass unbalance 0.02 m
	Damping factor η	0.03

Table 3 :

 3 Coefficient

			i	j	j E{Ψ i Ψ j Ψ k }	i	j	j E{Ψ i Ψ j Ψ k }
			7 12 12	2	15 16 23	1
			8 8 13	2	15 17 24	1
			8 9 14	1	16 16 25	2
			8 10 15	1	16 17 26	1
			8 11 16	1	17 17 27	2
	6 12	1	8 12 17	1	18 18 18	8
	2 2 13	2	9 9 18	2	18 19 19	2
	2 3 14	1	9 10 19	1	18 20 20	2
	2 4 15	1	9 11 20	1	18 21 21	2
	2 5 16	1	9 12 21	1	19 19 22	2
	2 6 17	1	10 10 22	2	19 20 23	1
	3 3 18	2	10 11 23	1	19 21 24	1
	3 4 19	1	10 12 24	1	20 20 25	2
	3 5 20	1	11 11 25	2	20 21 26	1
	3 6 21	1	11 12 26	1	21 21 27	2
	4 4 22	2	12 12 27	2	22 22 22	8
	4 5 23	1	13 13 13	8	22 23 23	2
	4 6 24	1	13 14 14	2	22 24 24	2
	5 5 25	2	13 15 15	2	23 23 25	2
	5 6 26	1	13 16 16	2	23 24 26	1
	6 6 27	2	13 17 17	2	24 24 27	2
	7 7 7	8	14 14 18	2	25 25 25	8
	7 8 8	2	14 15 19	1	25 26 26	2
	7 9 9	2	14 16 20	1	26 26 27	2
	7 10 10	2	14 17 21	1	27 27 27	8
	7 11 11	2	15 15 22	2

Table 6 :

 6 Sets of parameters

	0.05 rad 2

Table 7 :

 7 Critical speeds of the rotor system

	Critical speed Value (Hz)
	1st backward 27.9
	1st forward	28.3
	2nd backward 61.8
	2nd forward	97.2
	3rd backward 128.1
	3rd forward	240

Acknowledgements

Jean-Jacques Sinou gratefully acknowledges the financial support of the French National Research Agency through the Young Researcher program ANR-07-JCJC-0059-01-CSD 2.

A Model of the shaft

As illustrated in Figure 1, the nodal displacement of a beam element is defined by Ghanem, R., Spanos, P., (1991). Stochastic Finite Elements: A Spectral Approach. Springer-Verlag.

B Mass and gyroscopic matrices components

Lalanne, M., Ferraris, G., (1990)