Ana Cristina

Moreira Freitas
email: amoreira@fep.up.pt

Jorge Milhazes Freitas
email: jmfreita@fc.up.pt

Mike Todd
email: m.todd@st-andrews.ac.uk

Sandro Vaienti
email: vaienti@cpt.univ-mrs.fr

RARE EVENTS FOR THE MANNEVILLE-POMEAU MAP

We prove a dichotomy for Manneville-Pomeau maps f : [0, 1] → [0, 1]: given any point ζ ∈ [0, 1], either the Rare Events Point Processes (REPP), counting the number of exceedances, which correspond to entrances in balls around ζ, converge in distribution to a Poisson process; or the point ζ is periodic and the REPP converge in distribution to a compound Poisson process. Our method is to use inducing techniques for all points except 0 and its preimages, extending a recent result [HWZ14], and then to deal with the remaining points separately. The preimages of 0 are dealt with applying recent results in [AFV14]. The point ζ = 0 is studied separately because the tangency with the identity map at this point creates too much dependence, which causes severe clustering of exceedances. The Extremal Index, which measures the intensity of clustering, is equal to 0 at ζ = 0, which ultimately leads to a degenerate limit distribution for the partial maxima of stochastic processes arising from the dynamics and for the usual normalising sequences. We prove that using adapted normalising sequences we can still obtain non-degenerate limit distributions at ζ = 0.

Introduction

One of the standard ways to investigate the statistical properties of a dynamical system f : X → X with respect to some measure P is to look at its recurrence to certain points ζ in the system. This can be connected to Extreme Value theory: by taking a suitable observable ϕ : X → R taking its unique maximum u F at ζ, one can look at the behaviour of the iterates x, f (x), f 2 (x), . . . via the observations

X i = X i (x) = ϕ • f n (x).
If P is an f -invariant probability measure then X 0 , X 1 , . . . is a stationary stochastic process. We furthermore assume that P is ergodic in order to isolate specific statistical behaviour. So Birkhoff's Ergodic Theorem implies that these random variables will satisfy the law of large numbers. We can now consider the random variable given by the maximum of this process:

M n = M n (x) := max{X 0 (x), . . . , X n-1 (x)}.

Date: March 5, 2015. 2000 Mathematics Subject Classification. 37A50, 60G55, 60G70, 37B20, 60G10, 37C25. ACMF was partially supported by FCT grant SFRH/BPD/66174/2009. JMF was partially supported by FCT grant SFRH/BPD/66040/2009. Both these grants are financially supported by the program POPH/FSE. ACMF, JMF and SV are partially supported by FCT (Portugal) project PTDC/MAT/120346/2010, which is funded by national and European structural funds through the programs FEDER and COMPETE. ACMF, JMF and MT are also supported by CMUP (UID/MAT/00144/2013), which is funded by FCT (Portugal) with national (MEC) and European structural funds through the programs FEDER, under the partnership agreement PT2020. SV has been supported by the ANR grant PERTURBATIONS and by the project MOD TER COM of the french Region PACA.

Again by the ergodic theorem, if ζ is in the support of P, we expect M n → u F , so to obtain a non-trivial limit law, we need to rescale {M n } n . Indeed, we say that we have an Extreme Value Law (EVL) for M n if there is a non-degenerate distribution function H : R → [0, 1] with H(0) = 0 and for every τ > 0 there exists a sequence of levels u n = u n (τ ) such that

nP(X 0 > u n ) → τ as n → ∞, (1)
and for which the following holds:

P(M n u n ) → H(τ ) = 1 -H(τ ) as n → ∞,
where the convergence is meant at the continuity points of H(τ ).

In recent years, there has been a great deal of work on EVLs in the context of dynamical systems (see for example [Col01, FF08a, VHF09, FFT10, GHN11, HNT12, LFW12, Kel12, FHN14, AFV14]), the standard form of the observable ϕ being a function of the distance to ζ, for example ϕ(x) = -d(x, ζ) for d a metric on X . In many cases it has been shown that for P-a.e. ζ ∈ X , this setup gives an EVL with H = e -τ . More recently it has been shown that if ζ is a periodic point of period p then H = e -θτ where θ ∈ (0, 1) depends on the Jacobean of the measure for f p , and is referred to as the Extremal Index (EI). The EI is known to measure the intensity of clustering of exceedances of the levels u n . In fact, in many cases, the EI is equal to the inverse of the average cluster size, so that the EI is equal to 1 when there is no clustering. In the case of a class of uniformly hyperbolic dynamical systems, a stronger property, a dichotomy, has been shown: either ζ is periodic and we have an EVL with some extremal index θ ∈ (0, 1), or there is an EVL H = e -τ . This was shown for f some uniformly expanding interval maps with a finite number of branches in [START_REF] Ferguson | Escape rates for gibbs measures, Ergodic Theory[END_REF] (see also [START_REF]The extremal index, hitting time statistics and periodicity[END_REF]Section 6]) and with a countable number of branches in [START_REF] Aytaç | Laws of rare events for deterministic and random dynamical systems[END_REF]; here, depending on the precise form of the map, the measure can be absolutely continuous with respect to Lebesgue (acip), or an equilibrium state for some Hölder potential. Inducing methods have been used to extend some of these results to non-uniformly hyperbolic dynamical systems (see [START_REF]The compound Poisson limit ruling periodic extreme behaviour of non-uniformly hyperbolic dynamics[END_REF] which built on [START_REF] Bruin | Return time statistics via inducing[END_REF]), but the results have not thus far extended to such a complete dichotomy.

We can further enrich our process by considering the point process formed by entries into the regions {X > u n }, which in good cases gives rise to a Poisson process. An analogous dichotomy can often be shown there also: in the case of a periodic point ζ, we obtain a compound Poisson process. We leave the details of this construction to later.

In this note, we extend the dichotomy to a simple non-uniformly hyperbolic dynamical system, the Manneville-Pomeau (MP) map equipped with an absolutely continuous invariant probability measure. The form for such maps given in [LSV99, BSTV03] is, for α ∈ (0, 1), the corresponding MP map is

f = f α (x) = x(1 + 2 α x α ) for x ∈ [0, 1/2) 2x -1 for x ∈ [1/2, 1]
Members of this family of maps are often referred to as Liverani-Saussol-Vaienti maps since their actual equation was first introduced in [START_REF] Liverani | A probabilistic approach to intermittency[END_REF]. As it can be seen for example in [LSV99, You99, Hu04], these maps have polynomial decay of correlations:

ϕ • (ψ • f t )dµ α -ϕdµ α ψdµ α ≤ C ϕ H β ψ ∞ 1 t 1 α -1 , (2)
where H β denotes the space of Hölder continuous functions ϕ with exponent β equipped with the norm ϕ

H β = ϕ ∞ + |ϕ| H β ,
where

|ϕ| H β = sup x =y |ϕ(x) -ϕ(y)| |x -y| β .
Let h α = dµα dx . From [START_REF] Hu | Decay of correlations for piecewise smooth maps with indifferent fixed points[END_REF] we know that lim x→0 h(x)

x -α = C 0 > 0. Hence, for small s > 0 we have that

µ α ([0, s)) ∼ c s 1-α , (3)
where the notation ∼ c is used in the sense that there is c > 0 such that lim s→∞ µα([0,s))

s 1-α = c.
In this case there are canonical induced maps which capture all but a countable number of points in the phase space, so with some extra consideration for those points not captured, we can prove the full dichotomy, where for ζ a periodic point of period p, the extremal index is

θ = 1 -1/|Df p (ζ)|.
For the special case in which ζ is the indifferent fixed point, we prove that there exists an EI equal to zero, which corresponds to a degenerate limit law, when the usual normalising sequences are used. Moreover, we show that by changing the definition of (u n ) n given by (1) in a suitable way, we recover a non-degenerate EVL. This latter result relies on information on the transfer operator in [START_REF] Hirata | Statistics of return times: a general framework and new applications[END_REF] as well as a refinement of the techniques for proving EVLs at periodic points developed in [START_REF]Speed of convergence for laws of rare events and escape rates[END_REF].

1.1. Point process of hitting times. We will use our observations on our dynamical system to generate point processes. Here we adopt the approach and notation of [START_REF] Zweimüller | Mixing limit theorems for ergodic transformations[END_REF]. Let M p ([0, ∞)) be the space of counting measures on ([0, ∞), B [0,∞) ). We equip this space with the vague topology, i.e., ν n → ν in M p ([0, ∞)) whenever ν n (ψ) → ν(ψ) for any continuous function ψ : [0, ∞) → R with compact support. A point process N on [0, ∞) is a random element of M p ([0, ∞)). We will be interested in point processes N n : X → M p ([0, ∞)). If we have a fixed measure µ on X, we say that

(N n ) n converges in distribution to N if µ • N -1 n converges weakly to µ • N -1 . We write N n µ =⇒ N .
So given X 0 , X 1 , X 2 , . . . and some u ∈ R, we begin the construction of our point process R → M p ([0, ∞)) as follows. Given A ⊂ R we define

N u (A) := i∈A∩N 0 ½ X i >u .
So N u [0, n) counts the number of exceedances of the parameter u among the first n observations of the process X 0 , X 1 , . . . , X n-1 or, in other words, the number of entrances in U (u) up to time n.

We next re-scale time using the factor v := 1/P(X > u) given by Kac's Theorem. However, before we give the definition, we need some formalism. Let S denote the semi-ring of subsets of R + 0 whose elements are intervals of the type [a, b), for a, b ∈ R + 0 . Let R denote the ring generated by S. Recall that for every J ∈ R there are k ∈ N and k intervals I 1 , . . . , I k ∈ S such that J = ∪ k i=1 I j . In order to fix notation, let a j , b j ∈ R + 0 be such that

I j = [a j , b j ) ∈ S.
For I = [a, b) ∈ S and α ∈ R, we denote αI := [αa, αb) and

I + α := [a + α, b + α). Similarly, for J ∈ R define αJ := αI 1 ∪ • • • ∪ αI k and J + α := (I 1 + α) ∪ • • • ∪ (I k + α).
We suppose that we are given τ > 0 such that P(X > u n ) = τ n . We let U (u n ) = {X > u n } and let v n be the corresponding scaling factor defined above.

Definition. We define the rare event point process (REPP) by counting the number of exceedances (or hits to U (u n )) during the (re-scaled) time period v n J ∈ R, where J ∈ R. To be more precise, for every J ∈ R, set

N n (J) := N un (v n J) = j∈vnJ∩N 0 ½ X j >un .
π κ = θ(1 -θ) κ-1 , for every κ ∈ N 0 .
Theorem 2. For ζ = 0, consider the maximum function M n = M n (x) defined above.

(1) Let (u n ) n = (u n (τ )) n be chosen as in (1), then

P(M n u n ) → 1 as n → ∞ for any τ > 0. (2) If α ∈ (0,
√ 5 -2), then for each τ > 0, there exists a sequence of thresholds (u n ) n = (u n (τ )) n so that P(M n u n ) → e -τ as n → ∞.

Remark 1. We note that in the case ζ = 0, while it is possible to rescale the thresholds to recover an EVL as in Theorem 2 (2), the corresponding REPP remains degenerate. This result will form part of a forthcoming work [FFR].

1.3. Comments on history and strategy. Before discussing our approach we introduce some notation. For a dynamical system f : X → X and a subset A ⊂ X , for x ∈ X define

r A (x) := inf{n ∈ N : f n (x) ∈ A},
the first hitting time to A. Note that there is a connection with the behaviour of the variable r U (un) and our REPP since we can break that process down into a sequence of first hits to U (u n ). This gives a connection with our REPP and the asymptotics of r U (un) , the Hitting Time Statistics (HTS). One basic difference is that here we are concerned with all hits to U (u n ), not just the first.

Our main result for the case of P-typical points and for periodic points in (0, 1) follows quickly from previous works, including works already mentioned above, and indeed in some of these papers mention MP explicitly. We also remark that some of the earliest works on HTS for dynamical systems considered the case of MP maps with α 1, see for example [CG93, CGS92, CI95], with a focus on the behaviour at 0. In these cases, the sets A n considered were formed from dynamically defined cylinder sets and the analysis was done at 1/2, the preimage of 0, so that finite measure sets could be used. In this paper we consider the case α ∈ (0, 1), so f has an acip, and we also consider more general points and sets A n .

We will first consider all points in (0, 1), using inducing methods. This will require us to generalise the already very flexible result of [START_REF] Haydn | Return-time statistics, hitting-time statistics and inducing, Ergodic theory, open dynamics, and coherent structures[END_REF] to point processes. Finally we use the approach which goes back to Leadbetter [START_REF] Leadbetter | On extreme values in stationary sequences[END_REF] of proving some short range and long range recurrence conditions to prove that we have a degenerate law at 0 (the extremal index is 0).

Induced point processes

Here we aim to generalise [START_REF] Haydn | Return-time statistics, hitting-time statistics and inducing, Ergodic theory, open dynamics, and coherent structures[END_REF] to point processes. In that paper, they use [Zwe07, Corollary 5] as one of their key tools. In our, fairly analogous, setting we use [Zwe07, Corollary 6] instead. Note that previous results here include [BSTV03, Theorem 2.1], where they proved that for balls around typical points, the HTS of first return maps are the same as that for the original map -they also remarked, without details, that this can be extended to successive return times. Also in [START_REF]The compound Poisson limit ruling periodic extreme behaviour of non-uniformly hyperbolic dynamics[END_REF], we extended this idea to periodic points. The strengths of the approach in [START_REF] Haydn | Return-time statistics, hitting-time statistics and inducing, Ergodic theory, open dynamics, and coherent structures[END_REF] to HTS are that it covers all points, and that the proof is rather short.

We will give our result comparing the point process of the induced system to that coming from the original system in a general setting and then later apply this to our MP example. In this section, we take a dynamical system f : X → X with an ergodic f -invariant probability measure µ, choose a subset Y ⊂ X and consider F Y : Y → Y to be the first return map f r Y to Y (note that F may be undefined at a zero Lebesgue measure set of points which do not return to Y , but most of these points are not important, so we will abuse notation here). Let

µ Y (•) = µ(•∩Y ) µ(Y ) be the conditional measure on Y . By Kac's Theorem µ Y is F Y -invariant. Setting v Y n = 1/µ Y (X > u n ), for the induced process X Y i , N Y n (J) := N Y un (v Y n J) = j∈v Y n J∩N 0 ½ X Y j >un .
In keeping with [START_REF] Haydn | Return-time statistics, hitting-time statistics and inducing, Ergodic theory, open dynamics, and coherent structures[END_REF], we denote our inducing domain by Y . Denote the speeded up return time r A by r A,Y and the induced measure on Y by µ Y .

Theorem 3. For η > 0, setting J η := ∪ s∈J B η (s), we assume that N (J η ) is continuous in η, for all small η.

N Y n µ Y =⇒ N as n → ∞ implies N n µ =⇒ N as n → ∞.
Proof. By [Zwe07, Corollary 6], for a general sequence of point processes (N n ) n and an er-

godic reference measure m, if P ≪ m then N n P =⇒ N in M p ([0, ∞)) implies N n Q =⇒ N in M p ([0, ∞)) for any Q ≪ m.
So replacing both m and Q with µ and replacing P with µ Y we see that for our sequence of processes, if

N n µ Y =⇒ N in M p ([0, ∞)), then N n µ =⇒ N in M p ([0, ∞)).
Thus it suffices to show that for every J ∈ R and all k ∈ N,

µ Y (N Y n (J) k) n→∞ -→ µ Y (N (J) k) implies µ Y (N n (J) k) n→∞ -→ µ Y (N (J) k).
For δ > 0, let

E M = E δ M := 1 -δ µ(Y ) j r j Y 1 + δ µ(Y ) j for all j M and F N := {r U (un),Y N }.
As in [START_REF] Haydn | Return-time statistics, hitting-time statistics and inducing, Ergodic theory, open dynamics, and coherent structures[END_REF],

µ Y (F c N ) N µ Y (U (u n )) → 0 as n → ∞. Also the ergodic theorem says that µ Y ((E δ M ) c ) → 0 as M → ∞. For x ∈ E δ M , r k U (un),Y (x) 1 -δ µ(Y ) r k U (un) (x) = r k U (un),Y (x)-1 j=0 r Y •F j Y (x) = r r k U (un),Y (x) Y (x) r k U (un),Y (x) 1 + δ µ(Y )
We can deduce that for

x ∈ F N ∩ E δ M , µ(Y )r k U (un) (x) ∈ B δr k U (un),Y (x) (r k U (un),Y (x)). So if r k U (un),Y (x) ∈ v Y n J then µ(Y )r k U (un) (x) ∈ v Y n J δ and so r k U (un) (x) ∈ v n J δ . Therefore, µ Y ({N n (J δ ) k} ∩ (E N ∩ F M )) µ Y N Y n (J) k ∩ (E N ∩ F M ) .
Setting δ ′ := δ 1+δ , we also obtain that

1 µ(Y ) r k U (un),Y (x) ∈ B δ ′ r k U (un) (x) (r k U (un) (x))
for x ∈ F N ∩ E M . Analogously to above, this leads us to

µ Y N Y n (J δ ′ ) k ∩ (E N ∩ F M ) µ Y ({N n (J) k} ∩ (E N ∩ F M )) .
So since ε, δ > 0 were arbitrary, we are finished.

Application of inducing to Manneville-Pomeau

In this section we prove our main theorem for all points ζ ∈ (0, 1).

Let P be the renewal partition, that is the partition defined inductively by

Z ∈ P if Z = [1/2, 1)
or f (Z) ∈ P. Now let Y ∈ P and let F Y be the first return map to Y and µ Y be the conditional measure on Y . It is well-known that (Y, F Y , µ Y ) is a Bernoulli map and hence, in particular, a Rychlik system (see [START_REF] Rychlik | Bounded variation and invariant measures[END_REF] or [AFV14, Section 3.2.1] for the essential information about such systems) and so the REPP is understood as in [FFT13, Corollary 3]. Hence by [START_REF] Aytaç | Laws of rare events for deterministic and random dynamical systems[END_REF] we have the following theorem. (2) ζ is periodic with period p and N Y n converges in distribution to a compound Poisson process N with intensity θ = 1 -D(F -p Y )(ζ) and multiplicity d.f. π given by 1 π κ = θ(1 -θ) κ-1 , for every κ ∈ N 0 .

For points in Y \ ∪ n 1 f -n (0), this theorem is Proposition 3.2 of [START_REF] Aytaç | Laws of rare events for deterministic and random dynamical systems[END_REF]. For the boundary points ∪ n 1 f -n (0), in the language of [START_REF] Aytaç | Laws of rare events for deterministic and random dynamical systems[END_REF], any such point is called aperiodic non-simple. Hence by Proposition 3.4(1) of that paper, we have a standard extremal index of 1 at all such points. Varying Y means that we have considered all points in (0, 1). So combining Theorems 3 and 4 completes the proof of Theorem 1 for ζ = 0.

Analysis of the indifferent fixed point

The tangency of the graph of the MP map with the identity map, creates an intensive clustering of exceedances of levels (u n ) n∈N , when they are chosen as in (1), that leads to the existence of an EI equal to 0, which leads to a degenerate limit distribution for M n . However, if we choose the levels (u n ) n∈N not in the classical way, but rather a sequence of lower thresholds, so that the exceedances that escape the clustering effect have more weight, then we can recover the existence of a non-degenerate distribution for the maxima.

The proof of an EI equal to 0 for the usual normalising sequences follows easily from the existing connections between Return Times Statistics (RTS), Hitting Times Statistics (HTS) and EVL, which we will recall in the next subsection. The proof of the existence of a non degenerate limit, under a different normalising sequence of thresholds, is more complicated and requires some new results from [START_REF]Speed of convergence for laws of rare events and escape rates[END_REF], which we will recall below.

The usual normalising sequences case. For any

ζ ∈ [0, 1], let B ε (ζ) = (ζ -ε, ζ + ε) ∩ [0, 1].
If there exists a non degenerate d.f. G such that for all t 0,

lim ε→0 µ α r Bε(ζ) ≤ t µ α (B ε (ζ)) = G(t),
then we say we have Hitting Time Statistics (HTS) G for balls, at ζ. Similarly, we can restrict our observations to U (u): if there exists a non degenerate (d.f.) G such that for all t 0, lim

ε→0 µ α r Bε(ζ) ≤ t µ α (B ε (ζ)) B ε (ζ) = G(t),
then we say we have Return Time Statistics (RTS) G for balls, at ζ.

The normalising term in the definition of HTS/RTS is inspired by Kac's Theorem which states that the expected amount of time you have to wait before you return to B ε (ζ) is exactly

1 µα(Bε(ζ)) .
The existence of exponential HTS (G(t) = 1-e -t ) is equivalent to the existence of exponential RTS ( G(t) = 1 -e -t ). In fact, according to the Main Theorem in [START_REF] Haydn | Hitting and return times in ergodic dynamical systems[END_REF], a system has HTS 1 We note that there is an error in [FFT13, Theorem 1], propagated throughout the main results there: the κ should be replaced by κ -1.

G if and only if it has RTS G and

G(t) = t 0 (1 -G(s)) ds. (4)
This formula was later generalised to obtain a relation between the distributional properties of Hitting Times and Return Times Points Processes in [START_REF] Haydn | Averaged number of visits[END_REF].

Moreover, according to the main theorem in [START_REF] Cristina | Hitting time statistics and extreme value theory[END_REF] the existence of HTS G for balls at ζ is equivalent to the existence of an EVL H for the process generated dynamically by an observable with a maximum at ζ and, in fact, G = H.

Let U = [0, b) and A = [a, b), where a is such that f (a) = b, i.e., b = a + 2 α a 1+α . Using (3) we easily get µ α (U ) ∼ c a 1-α + (1 -α)2 α a + o(a) and µ α ([0, a)) ∼ c a 1-α .
Next we compute the RTS distribution, which we denote by G(s). For s ≤ 0, we easily have that G(s) = 0, since r U ≥ 1, by definition of hitting time. Let s > 0 then

G(s) = lim b→0 µ U r U ≤ s µ(U ) = lim b→0 1 µ(U ) µ r U ≤ s µ(U ) ∩ U ≥ lim b→0 µ(U \ A) µ(U ) = lim b→0 µ([0, a)) µ([0, b)) = 1
Using the [HLV05] formula, we get G(t) = t 0 1 -F (s)ds = 0, which , by [START_REF] Cristina | Hitting time statistics and extreme value theory[END_REF], corresponds to an EI equal to 0. Recall that H(τ ) = e -θτ = 1, which means that, in this case, H(τ ) = 0.

Adjusted choice of thresholds.

In order to prove the existence of EVLs in a dynamical systems context, there are a couple of conditions on the dependence structure of the stochastic process that if verified allow us to obtain such distributional limits. These conditions are motivated by the conditions D(u n ) and D ′ (u n ) of Leadbetter but were adapted to the dynamical setting and further developed both in the absence of clustering, such as in [START_REF] Collet | Statistics of closest return for some non-uniformly hyperbolic systems[END_REF][START_REF]On the link between dependence and independence in extreme value theory for dynamical systems[END_REF][START_REF] Holland | Extreme value theory for non-uniformly expanding dynamical systems[END_REF], and in the presence of clustering in [START_REF]The extremal index, hitting time statistics and periodicity[END_REF]. Very recently, in [START_REF]Speed of convergence for laws of rare events and escape rates[END_REF], the authors provided some more general conditions, called Д(u n ) and Д ′ q (u n ), which subsumed the previous ones and allowed them to address both the presence (q ≥ 1) and the absence (q = 0) of clustering. To distinguish these conditions the authors used a Cyrillic D to denote them. We recall these conditions here.

Given a sequence (u n ) n∈N of real numbers satisfying (1) and q ∈ N 0 , set

A (q) n := {X 0 > u n , X 1 ≤ u n , . . . , X q ≤ u n }.
For s, ℓ ∈ N and an event B, let

W s,ℓ (B) = s+ℓ-1 i=s f -i (B c ). (5)
Condition (Д q (u n )). We say that Д(u n ) holds for the sequence X 0 , X 1 , . . . if, for every ℓ, t, n ∈ N

P A (q) n ∩ W t,ℓ A (q) n -P A (q) n P W 0,ℓ A (q) n ≤ γ(q, n, t), (6)
where γ(q, n, t) is decreasing in t and there exists a sequence (t n ) n∈N such that t n = o(n) and nγ(q, n, t n ) → 0 when n → ∞.

For some fixed q ∈ N 0 , consider the sequence (t n ) n∈N given by condition Д q (u n ) and let (k n ) n∈N be another sequence of integers such that

k n → ∞ and k n t n = o(n). (7)
Condition (Д ′ q (u n )). We say that Д ′ q (u n ) holds for the sequence X 0 , X 1 , . . . if there exists a sequence (k n ) n∈N satisfying (7) and such that

lim n→∞ n ⌊n/kn⌋ j=1 P A (q) n ∩ f -j A (q) n = 0. (8)
We note that, when q = 0, condition

Д ′ q (u n ) corresponds to condition D ′ (u n ) from [Lea74]. Now let ϑ = lim n→∞ ϑ n = lim n→∞ P(A (q) n ) P(U n ) . (9)
From [FFT15, Corollary 2.4], it follows that if the stochastic process X 0 , X 1 , . . . satisfies both conditions Д q (u n ) and Д ′ q (u n ) and the limit in (9) exists then

lim n→∞ P(M n ≤ u n ) = e -ϑτ .
Now, we consider the fixed point ζ = 0. For every n ∈ N, let b n be such that

U n = {X 0 > u n } = [0, b n ) and µ α (U n ) ∼ τ /n. Also set a n ∈ U n so that f α (a n ) = b n , i.e., b n = a n + 2 α a 1+α n . Using (3) we easily get µ α (U n ) ∼ c a 1-α n + (1 -α)2 α a n + o(a n ) (10) µ α ([0, a n )) ∼ c a 1-α n (11) µ α ([a n , b n )) ∼ c (1 -α)2 α a n + o(a n ) (12)
Now, since we are assuming that

µ α (U n ) ∼ τ /n, then a n ∼ c 1/n 1/(1-α) . Observe that µ α (U n ∩ f -1 α (U n )) = µ α ([0, a n )) ∼ c a 1-α n ∼ c 1/n.
Hence, if we consider q = 0, the periodicity of ζ implies that Д ′ q (u n ) does not hold since

n ⌊n/kn⌋ j=1 P U n ∩ f -j (U n ) ≥ nµ α (U n ∩ f -1 α (U n )) > 0,
for all n ∈ N. Hence, here, given that ζ is a periodic point of period 1 the natural candidate for q is q = 1. From here on we always assume that q = 1.

In this case, A

n = [a n , b n ) =: Q n . (q)
However, if we plug (12) and (10) into (9), we obtain that ϑ = 0, which means that the natural candidate for a limit distribution for µ

α (M n ≤ u n ) is degenerate.
The problem is that the indifferent fixed point creates too much dependence. In [START_REF]The extremal index, hitting time statistics and periodicity[END_REF], under a condition called SP , we have seen that when ζ is periodic, the probability of having no entrances in U n , among the first n observations, is asymptotically equal to the probability of having no entrances in Q n , among the first n observations, i.e.,

lim n→∞ P(M n ≤ u n ) = lim n→∞ P(W 0,n (U n )) = lim n→∞ P(W 0,n (Q n )).
In [START_REF]Speed of convergence for laws of rare events and escape rates[END_REF], it was shown that it is possible to replace U n by Q n even without the SP condition (see [START_REF]Speed of convergence for laws of rare events and escape rates[END_REF]Proposition 2.7]). Making use of this upgraded result, we can now change the normalising sequence of levels (u n ) n∈N so that we can still obtain a non-degenerate limit for P(M n ≤ u n ). To understand the need to change the normalising sequence in order to obtain a non-degenerate limit, recall that condition (1) guaranteed that M n was normalised by a sequence of levels that kept the average of exceedances among the first n observations at an (almost) constant value τ > 0. When ϑ > 0, condition (1) also guarantees that the average number of entrances in Q n among the first n observations is kept at an (almost) constant value θτ > 0. Here, since ϑ = 0, we need to change u n so that the average number of entrances in

Q n is controlled, i.e., lim n→∞ nP(A (q) n ) = τ > 0. ( 13
)
From equations (2.15) and (2.16) from [START_REF]Speed of convergence for laws of rare events and escape rates[END_REF] one gets:

P(W 0,n (A (q) n )) -1 - n k n P(A (q) n ) kn ≤2k n t n P(U n ) + 2n ⌊n/kn⌋-1 j=1 P A (q) n ∩ f -j A (q) n + γ(q, n, t n ) (14)
Note that since by (13) we have lim n→∞ 1 -n kn P(A

n )

kn = e -τ , then if both conditions Д q (u n ) and Д ′ q (u n ) hold, then all the terms on the left of (14) vanish, as n → ∞, and consequently:

lim n→∞ P(M n ≤ u n ) = lim n→∞ P(W 0,n (A (q) n )) = e -τ . (15)
Hence, in order to show that we can still obtain a non-degenerate limiting law for the distribution of M n when ζ = 0, we start by taking a sequence (u n ) n∈N so that (13) holds. Note that this implies that by ( 12) and (10) we have that a n ∼ c 1/n and µ α (U n ) ∼ c 1/n 1-α . In particular, this means that lim n→∞ nµ α (U n ) = ∞, which contrasts with the usual case where condition (1) holds.

To prove the existence of the limit in (15) we need to verify conditions Д q (u n ) and Д ′ q (u n ), where q = 1. We start by the latter, which is more complicated. 4.2.1. Proof of Д ′ q (u n ). We will next focus on the proof of Д ′ q (u n ) in the case of part (2) of Theorem 2. That is, (a n ) n will be chosen so that a n ∼ c 1/n, as described above. Later we will note that we can change (a n ) n to recover a degenerate law as in part (1) of that theorem.

We have to estimate the quantity

∆ ′ n := n [n/kn] j=1 µ α (Q n ∩ f -j Q n ) where Q n = [a n , b n ), for a n ∼ c 1 n and b n = f (a n ).
We follow the proof of [HSV99, Lemma 3.5]. By denoting by P the transfer operator and by τ n ∈ N the first return time of the set Q n into itself, we have:

∆ ′ n n [n/k n ]µ a (Q n ) sup j=τn,...,[n/kn] sup Qn P j (1 Qn h) h
where h is the density of µ α . In order to compute P τn (1 Qn h) we need to know how many branches of f τn have their domain intersecting Q n . If ξ 0 is the original partition into the sets [0, 1/2), [1/2, 1], we denote with ξ k the join ξ k := ξ 0 ∨ f -1 ξ 0 , . . . , ∨f k-1 ξ 0 .

We begin to observe that Q n contains at most one boundary point of the partition ξ τn-1 , otherwise one point of Q n should be sent into the same set, being f τn-1 onto on each domain of injectivity. Then when we move to ξ τn , the interval Q n will be crossed by at most 4 cylinders of monotonicity of the partition ξ τn . By denoting them with C τn,1 , . . . C τn,4 , we have

P τn (1 Qn h) = i=1,4 h • f -τn i 1 f τn i Qn Df τn • f -τn i
where f τn i denotes the branch of f τn restricted to C τn,i . Notice that the density is computed in Q n whose left boundary point is 1/n, so h is bounded from above by a constant times n α . We have now to estimate the derivative Df τn

• f -τn i on the sets Q n ∩ C τn,i . Let us define r m as the m-left preimage of 1, r m := f -m 1 (1) and define m(n) as r m(n) a n r m(n)-1 . Then the interval [a n , b n ) will intersect the two cylinders (r m(n) , r m(n)-1 ) and (r m(n)-1 , r m(n)-2
) and the first return of Q n will be larger than the first returns of those two cylinders; on the other hand the first return of (r m(n) , r m(n)-1 ) is m(n). The derivative Df τn will be computed at a point ι n which will be in one of those two cylinders; suppose without any restriction that ι n ∈ (r m(n) , r m(n)-1 ). Since we need to bound from below the derivatives, we begin to replace Df τn (ι n ) with Df m(n) (ι n ); then we observe that the map f m(n) : [r m(n) , r m(n)-1 ] → [0, 1] is onto and we use the distortion bound given, for instance, in [LSY, Lemma 5] which states that there exists a constant C such that for any m 1 and any x, y ∈ [r m , r m-1 ] we have

Df m (x) Df m (y)
C. We finally note that m(n) ∼ c n α . This implies immediately that

1 Df m(n) (ι n ) C|r m(n)-1 -r m(n) | ∼ c C 1 m(n) 1 α +1 ∼ c C 1 n 1+α .
Consequently (C will continue to denote a constant which could vary from one bound to another) P τn (1 Qn h) ∼ c 1 n We now continue as in [HSV] by getting for the other powers of the transfer operator:

P j (
∆ ′ n n [n/k n ]µ a (Q n ) C inf h 1 n We now know that µ a (Q n ) ∼ c 1 n ; hence ∆ ′ n Cn [n/k n ]µ a (Q n ) 2 1 µ a (Q n ) 1 n ∼ c [n 2 µ a (Q n ) 2 ] 1 k n .
So letting n → ∞, we see that Д ′ q (u n ) holds.

4.2.2. Proof of Д q (u n ). This follows since, as in (2), we have decay of correlations of Hölder functions against bounded measurable functions and condition Д q (u n ) was designed to follow from sufficiently fast decay of correlations, as shown in [Fre13, Proposition 5.2]. In order to compute the required rate of decay of correlations, which will impose a restriction on the domain of the parameter α, we recall here the above-mentioned result so that we can follow the computations closely.

Proposition 1 ([Fre13, Proposition 5.2]). Assume that X is a compact subset of R d and f : X → X is a system with an acip P, such that dP Leb ∈ L 1+ǫ . Assume, moreover, that the system has decay of correlations for all ϕ ∈ H β against any ψ ∈ L ∞ so that there exists some C > 0 independent of ϕ, ψ and t, and a rate function

̺ : N → R such that ϕ • (ψ • f t )dP -ϕdP ψdP ≤ C ϕ H β ψ ∞ ̺(t), (16)
and n 1+β(1+max{0,(ǫ+1)/ǫ-d}+δ) ̺(t n ) → 0, as n → ∞ for some δ > 0 and t n = o(n). Then condition Д q (u n ) holds.

Remark 2. We note that during the proof, in order to obtain the condition on the rate of decay of correlations, it is assumed that P(A

(q) n ) ∼ c 1/n.
Observe that since we are working in dimension 1, which means d = 1, then max{0, (ǫ + 1)/ǫ -d} = 1/ǫ. Also, from [START_REF] Hu | Decay of correlations for piecewise smooth maps with indifferent fixed points[END_REF], we may assume that the decay of correlations is written for Lipschitz functions, which allows us to take β = 1. Hence, for condition Д q (u n ) hold, we need that the rate of decay of correlations ̺ is sufficiently fast so that there exists some δ > 0 such that lim

n→∞ n 2+1/ǫ+δ ̺(t n ) = 0, (17)
where t n = o(n). From (3), in order that the density h α ∈ L 1+ǫ , we need that ǫ < 1/α -1.

Since by (2), we have that ̺(t) = t -(1/α-1) , then by (17) it is obvious that we must have α < 1/2, which implies that 1/ǫ < α + 2α 2 . Taking t n = n 1-α , we obtain: n 2+1/ǫ+δ n 1-α -(1/α-1) = n 2+1/ǫ+δ n -1/α+2-α < n 4+2α 2 +δ-1/α .

Hence, if α < √ 5 -2 we can always find δ > 0 so that (17) holds and consequently condition Д q (u n ) is verified.

1. 2 .

 2 Main results. Theorem 1. Given ζ ∈ (0, 1], consider the REPP N n defined above. Then either (1) ζ is not periodic and N n converges in distribution to a Poisson process N with intensity 1. (2) ζ is periodic with period p and N n converges in distribution to a compound Poisson process N with intensity θ = 1 -|det D(f -p )(ζ)| and multiplicity distribution function π given by

Theorem 4 .

 4 Given ζ ∈ Y , consider the REPP N Y n defined above. Then either (1) ζ is not periodic and N Y n converges in distribution to a Poisson process N with intensity 1.


Riccardo Artoni
email: riccardo.artoni@ifsttar.fr

Patrick Richard
email: patrick.richard@ifsttar.fr

Average balance equations, scale dependence, and energy cascade for granular materials

Keywords: numbers: 45.70.-n, 47.57.Gc, 83.10.Ff

A new averaging method linking discrete to continuum variables of granular materials is developed and used to derive average balance equations. Its novelty lies in the choice of the decomposition between mean values and fluctuations of properties which takes into account the effect of gradients. Thanks to a local homogeneity hypothesis, whose validity is discussed, simplified balance equations are obtained. This original approach solves the problem of dependence of some variables on the size of the averaging domain obtained in previous approaches which can lead to huge relative errors (several hundred percentages). It also clearly separates affine and nonaffine fields in the balance equations. The resulting energy cascade picture is discussed, with a particular focus on unidirectional steady and fully developed flows for which it appears that the contact terms are dissipated locally unlike the kinetic terms which contribute to a nonlocal balance. Application of the method is demonstrated in the determination of the macroscopic properties such as volume fraction, velocity, stress, and energy of a simple shear flow, where the discrete results are generated by means of discrete particle simulation.

I. INTRODUCTION

The behavior of granular materials under deformation and flow is often studied by means of numerical methods that solve the equations of motion for every single body and account for the interactions between particles through appropriate models (see, e.g., Refs. [1][2][3]). These methods, commonly grouped under the name of discrete element methods (DEM) [START_REF] Smit | Understanding Molecular Simulation[END_REF][START_REF] Schäfer | [END_REF][6], can be extremely useful for simulating discrete media with the purpose of establishing micro-macro relations, and in this perspective they need consistent averaging procedures [7][8][9][10][11][12]. An averaging method computes continuum-like variables (e.g., a velocity field) starting from their discrete, particle-scale counterparts (e.g., particle velocities). When using such an approach, two issues have to be considered with utmost attention: the representativity of the average and the physical meaning of the obtained estimates. We can say that an average performed at a given point is representative if a sufficiently large domain in space and time can be defined, where particles share the same properties. Such a volume, usually called "representative volume element," may be defined with respect to only some variables.

An important step in the development of averaging techniques was made by Babic [7]. In that work the author developed a general framework for weighted space-time averages, applicable to a wide range of conditions, and giving a large set of self-consistent continuum balance equations. In Babic's method, the definition of a representative volume is intrinsically that of a zone where affine velocity fields are locally uniform, i.e., no gradients occur. The assumption that such a volume can be defined is called by Babic the "continuum assumption." Now, the lack of scale separation typical of granular flows (the scale of spatial variation of variables has the same order of magnitude of the particle size) implies that in a three-dimensional (3D) flow we cannot in practice define such a volume. As a consequence, the physical meaning of averages obtained using Babic's approach in a 3D flow is generally questionable. At least some average estimates will then probably depend on the size of the averaging volume [8,13].

We can see a simple example of this by considering a uniform shear flow; if all particles follow the mean flow with no frustration, we will tend to say that fluctuations are negligible. But this depends on the definition of fluctuation. If we calculate the fluctuational kinetic energy ε T [which is related to granular temperature [14]; see also Eq. [START_REF] Woods | The Thermodynamics of fluid Systems[END_REF]] by considering fluctuations as Babic does, we find ε T ∼ ( γD) 2 , where D is the averaging scale. Thus the fluctuations depend on the averaging scale which casts doubt on such a method. Also the kinetic contribution to the stress tensor will suffer the same problem. Thus, as we will show below, such an approach may lead to errors that can reach several hundred percentages.

In this work we focus on this problem, with a particular attention on the mechanical energy equation. We discuss both the issue of dependence on the size of the averaging domain and the energy cascade picture associated with the averaging method.

It is clear that the dependence of some average variables on the averaging domain size is intrinsic to Babic's definition of fluctuations. In the following we show how a different fluctuation decomposition taking into account gradients yields a new derivation of continuum balance equations. These balance equations have the merit of clearly separating affine and nonaffine fields and identifying the terms responsible for averaging scale dependence. With this decomposition it is possible to attribute a more sound physical meaning to the terms appearing in the final balance equations. We show how these balance equations may be greatly simplified by a "local homogeneity assumption," which is strongly related to the fluctuation decomposition itself and which yields less strict requirements than Babic's "continuum assumption." Applying this method to a simple class of flows, we find that two separate paths for fluctuating energy dissipation are implied: while kinetic stress power enters a nonlocal balance, the contact stress power seems to be dissipated locally. This result, as well as the robustness of the method, is verified with discrete element simulations of a simple shear flow.

Given a particle property ψ p (x p , t ′ ), Babic defined the averaged property ψ at point x and time t as a weighted spacetime average:

ρ ψ = ∞ -∞ p w(x p -x, t ′ -t)m p ψ p dt ′ , (1)
where w(x px, t ′ -t) = w p is a normalized weighting function, ρ the average density, and m p the mass of a grain. The fluctuation was simply defined by Babic with respect to the center of the averaging volume:

ψp = ψ p -ψ(x). (2)
This may seem harmless, but in fact-as shown in Ref. [8]-the dependence of some terms of the continuum averaged equations on the size of the averaging domain is strictly related to this choice for the fluctuation decomposition.

We propose to define the fluctuation not with respect to the average value at point x (as was done in Ref. [7]), but with respect to the particle center, ψp = ψ p -ψ(x p ).

(

In the following we will redevelop continuum-averaged equations in the spirit of Babic's approach, but using this new fluctuation decomposition.

In order to simplify the result, we will make a constitutive assumption. We suppose that a scale exists where the gradients of ψ are smooth and that if this scale is the averaging scale intrinsic to the weighting function w p , we can approximate ∇ ψ as a constant near the averaging point. We call this assumption the local homogeneity assumption. Together with this assumption we also assume a homogeneous distribution of particle centers near the averaging point and the decorrelation of positions and velocities. These three assumptions are strictly related.

While the absence of scale separation usually prevents identification of a scale where gradients are zero, such a first-order approximation is reasonable and-as we will see-very fruitful.

The local homogeneity assumption, involving gradients instead of variables, is less strict than Babic's continuum hypothesis and therefore more likely to be valid under a proper choice of the weighting function. Based on this assumption, we can therefore approximate ψ(x p ) by

ψ(x p ) = ψ(x) + (x p -x) • ∇ ψ(x). ( 4
)
Following Babic's approach, it is possible to develop a generic balance equation for the particle property ψ(x p ). Details are given in Appendix A; here we show the main findings, with a particular regard to the continuity, momentum, and translational energy equations. When performing the derivation, many new terms appear (with respect to Babic's treatment), but most of them vanish due to the local homogeneity hypothesis.

B. Continuum balances

The local homogeneity hypothesis leaves the continuity equation unchanged:

∂ ∂t ρ + ∇ • ρv = 0, (5)
where average density and velocity are defined by:

ρ = ∞ -∞ p w p m p dt ′ , (6) ρv
= ∞ -∞ p w p m p v p dt ′ . (7)
As for the momentum balance, we end up with the classical equation:

∂ ∂t ρv + ∇ • ρvv = ∇ • T + ρg, (8)
where the total stress tensor is composed by three contributions:

T = T c + T k + T γ .
The first two terms are the same as in Babic. The first contribution is due to contact forces (both long-lasting contacts and collisions) and is given by

T c = ∞ -∞ p q>p w F pq l pq f pq dt ′ , (9)
where w F pq is a weight function related to the fraction of the branch vector joining the two particles p and q lying within the averaging function. Due to the new decomposition, Babic's kinetic stress tensor is decomposed in a truly kinetic part,

T k = - ∞ -∞ p w p m p ṽp ṽp dt ′ , (10)
and a new term,

T γ = -ρ(D • ∇v)(D • ∇v), (11)
which contains all the dependence on averaging domain size through the vector D. It is useful to recall that the emergence of a part of the stress tensor depending on velocity gradients has nothing to do with constitutive relations but is the joint product of coarse graining and scale dependence. Due to their definition, both T γ and T k are symmetric. Evidently, in absence of velocity gradients, T γ will vanish; in that case, the kinetic stress tensor as defined by Babic will no longer depend on the coarse-graining length. The components of D, related to the characteristic size of the domain along each direction, are defined by the following equation:

ρD 2 i = ∞ -∞ p w p m p (x pi -x i ) 2 dt ′ . ( 12
)
where x i is the i-th coordinate of the averaging point.

If the local homogeneity hypothesis holds, D i can be shown to scale with D m,i , the size of the averaging domain in the i direction (see Appendix B). The scaling will depend also on the shape of the weighting function.

Moreover, if the assumption of local homogeneity holds, T γ contains all the dependence on the coarse-graining width of the stress tensor while the kinetic stress tensor T k defined by means of the present fluctuation decomposition appears to be independent on averaging domain size.

The translational kinetic energy balance is only slightly more complicated than the one obtained by Babic. Once more, most of the new terms are negligible under the local homogeneity assumption.

The total translational kinetic energy is found to be the sum of three contributions, ρE T + ρE γ T + ρε T , which are respectively the translational kinetic energy of the mean flow, the translational kinetic energy related to the velocity gradients, and the truly fluctuational translational kinetic energy. These three variables are defined as

ρE T = ρ 1 2 v • v, (13)
ρE γ T = 1 2 ρ(D • ∇v) • (D • ∇v), (14)
ρε T = 1 2 ∞ -∞ p w p m p ṽp • ṽp dt ′ . ( 15
)
The translational kinetic energy balance equation is derived as:

∂ ∂t (ρE T + ρE γ T + ρε T ) + ∇ • (ρE T + ρE γ T + ρε T )v = ∇ • (q k + q c ) +∇ • ((T k + T c + T γ ) • v) +ρg • v -γ k -γ γ , ( 16
)
where q k and q c are Babic's energy fluxes (see Appendix A).

The last two terms correspond to the rate of conversion of kinetic energy to other forms of energy (dissipation due to friction, collisions, plasticity, but also to storage and restitution of elastic energy). In the rigid limit they represent dissipation through friction and inelastic collisions and are therefore positive. The first one contains the contribution of fluctuations:

γ k = - ∞ -∞ p q>p f pq • ( ṽp -ṽq )w S pq dt ′ , (17)
and the second one the rate of conversion of mechanical energy due to affine deformations:

γ γ = - ∞ -∞ p q>p f pq • ((x p -x q ) • ∇v)w S pq dt ′ , (18)
where w S pq is the mean of w p and w q . Both terms are independent of the size of the averaging domain. We can see here again the benefit of our approach: On one hand the domain size-dependent terms are clearly identified in the equation (ρE γ T and T γ ); on the other hand, the contribution of affine deformations is separated from the contribution of the fluctuations.

In fact, by superficially interpreting Babic's approach, one could conclude that all the energy is dissipated by means of fluctuations. This is incorrect, since affine deformations do also dissipate energy; the reason for this misunderstanding is that affine deformations are considered as fluctuations in Babic's method.

We can remove the mean flow kinetic energy terms by manipulating the momentum equation in a classical way (for example, see Ref. [START_REF] Woods | The Thermodynamics of fluid Systems[END_REF]). The result is an equation expressing the balance for the fluctuating energy related to fluctuations and to affine deformations:

∂ ∂t (ρE γ T + ρε T ) + ∇ • (ρE γ T + ρε T )v = (T c + T γ + T k ) † : ∇v + ∇ • (q k + q c ) -γ k -γ γ , ( 19
)
where † stands for transposition. Let us do some more manipulation. By looking at the term accounting for the rate of conversion of mechanical energy due to affine deformations, manipulation of Eq. ( 18), making use of l pq = x qx p , allows us to rewrite γ γ in the following form:

γ γ = T c * † : ∇v, (20)
where

T c * = ∞ -∞ p q>p w S pq l pq f pq dt ′ , (21)
is a new contact stress tensor, differing from the T c by the weight used. Indeed, as mentioned above, w S pq is the simple mean of w p and w q , whereas w F pq is related to the fraction of the branch vector l pq lying within the averaging region.

III. ONE DIRECTIONAL, STEADY, FULLY DEVELOPED FLOW

In this section we will apply our method to a simplified flow situation in order to gather some information on the effect of the new terms on force balance and on energy cascade.

Let us consider a one directional, steady, fully developed flow [defined by v = (v x , 0, 0), ∂ x = 0, ∂ t = 0]. In such a case, the continuity equation is identically zero. The momentum equation becomes:

∇ • (T c + T γ + T k ) + ρg = 0, ( 22
)
and the fluctuating energy equation

(T c + T γ + T k ) † : ∇v + ∇ • (q k + q c ) -γ k -γ γ = 0. ( 23
)
According to the definition of T γ , it is easy to see that its only nonzero component is:

T γ xx = -ρ(D y ∂ y v x + D z ∂ z v x ) 2 , ( 24
)
and it is therefore straightforward to see that

(T γ ) † : ∇v = 0, (25)
that is, the contribution of affine deformations to the stress tensor does not perform any work. Moreover, it is obvious that T γ directly affects only the x component of the force balance. Given that the flow is fully developed, we can safely remove it out from the balances, obtaining:

∇ • (T c + T k ) + ρg = 0, (26)
(T c + T k ) † : ∇v + ∇ • (q k + q c ) -γ k -γ γ = 0. ( 27
)
According to Eq. ( 20):

(T c ) † : ∇v -γ γ = (T c -T c * ) : ∇v. ( 28
)
It is clear that T c and T c * are not in principle the same, particularly when the averaging domain is smaller than a particle diameter. However, if the averaging domain spans more than one particle diameter, we can assume that T c ≈ T c * . Equation (27) becomes therefore:

T k : ∇v + ∇ • (q k + q c ) -γ k = 0. ( 29
)
This means that the work provided by the contact part of the stress tensor is converted locally by affine deformations, while only the kinetic part of the stress tensor enters into a nonlocal energy balance. This result has important consequences, since it shows that separate paths may exist for the dissipation of the work of kinetic and contact stresses. However, it does not cast doubt on the validity of hydrodynamic theories [START_REF] Bocquet | [END_REF][17][18][19] since our results do not deal with constitutive relations but suggest that a revision of the balances used for such theories could lead to a simplified description in some cases. Concerning the size of the averaging domain, it may seem strange to talk about domains smaller than a particle diameter, because it may seem obvious that such domains are inappropriate: Few particles are indeed contained in a snapshot of such a domain. However, we are always dealing with space-time averages, so even if the averaging domain is small, in many cases (e.g., slowly evolving or steady-state processes) the domain indeed contains a lot of particles since many of them traveled through it during their history without affecting local homogeneity. In other cases, periodic boundary conditions may help too, allowing use of small domains along some directions. That is one of the reasons why in the next section we will employ the averaging method to some DEM simulation results pertaining to this class of simplified flows.

IV. DEM SIMULATIONS A. Method & parameters

We use our own 2D implementation [2,20] of the classical discrete element method where Newton's equations of motion for a system of N "soft" disks are integrated. Such a technique, which is able to reproduce successfully the experimental results in many configurations (e.g., gravity-driven flows [1,2,21,22], sheared systems [3,23], granular materials close to jamming [24], silos [25], and rotating drums [26,27]), requires giving an explicit expression for the inter-particle forces. The discrete element method is classical and well known and can found in the aforementioned references. Therefore, we just present here the forces used in this work. For the normal force between two overlapping disks we use a standard linear spring-dashpot interaction model [6],

F n = k n δ n -γ n v n ,
where δ n is the normal overlap, k n is the spring constant, γ n the damping coefficient, and v n the normal relative velocity. The damping models the dissipation characteristic of granular materials. Likewise, the tangential force is modeled as a linear elastic and linear dissipative force in the tangential direction F t = k t δ t -γ t v t , where k t is the tangential spring constant, δ t the tangential overlap, γ t the tangential damping, and v t the tangential velocity at the contact point. The magnitude of δ t is truncated as necessary to satisfy the Coulomb law, |F t | ≤ µ |F n |, where µ is the grain-grain friction coefficient. The simulated system is two dimensional (Fig. 1). The granular material is a dense assembly of N dissipative disks of average diameter d and average mass m. A small polydispersity of ±20% is considered to prevent crystallization.

The granular material is submitted to a plane shear, without gravity, leading to a uniform stress distribution. The material is sheared between two parallel rough walls, separated by a distance H. One of the walls is fixed, while the other moves at the prescribed velocity V . The flow and transverse directions are respectively called x and y. Periodic boundary conditions are applied along the flow direction, and the length of the simulation box, L x , is set to 60 grain sizes. The wall roughness is made of disks sharing the characteristics of the flowing grains (same polydispersity and mechanical properties, no rotation). Their centers are equally spaced by a distance equal to the largest disk diameter 1.2d. Along the y axis, the position y = 0 corresponds to the center of the glued grains on the fixed wall. The normal stress applied on the moving wall, P 0 , is controlled. The vertical position of the moving wall is thus not fixed and, using the method described in Ref. [28], the height of the system H obeys Ḣ = (P 0 -P w )L x /g p , where g p is a viscous damping parameter and P w is the normal stress exerted by the grains on the moving wall.

The following values of the parameters are used: k n /P 0 = 10 5 , k t = 2k n /7, γ t = 0, and µ = 0.5 and g p = 100 √ k n m. The value of γ n is adjusted to obtain a normal restitution coefficient e n = 0.5 [START_REF] Schäfer | [END_REF]. The equations of motion for the translational and rotational degrees of freedom are integrated with a velocity-Verlet scheme with a time step 10 -4 m/P 0 .

The number of grains, N , is adjusted so the size of the system in the y direction, H, is roughly equal to 80d. The initial state of the system is a randomly diluted hexagonal lattice with disorder both on grain positions and velocities. The attainment of a steady state is verified by observing timeinvariant total kinetic energy and distance between the two walls.

B. Averaging procedure

For the simulation set up described above, we implemented an averaging procedure following the method developed in the earlier sections of the paper.

Several snapshots of the state of the system (positions, velocities, contact forces) at given times were captured. For each time a space average was performed; given that the system was periodic, averages were computed at different y values. A Heaviside step function was chosen as the weighting function. The only parameter of the weighting function was therefore the averaging window size. Some authors suggested [9,10] smoother weighting functions. Here we preferred a step function for its simplicity and to avoid smoothing the profiles too much. Given that the system was simple and periodic, we felt that a nonsmooth function was the best choice, so as not to lose information about layering, localization, and so on.

Given that the simulations displayed a stationary state, the time average was performed as a simple mean over the results obtained for each snapshot. Standard errors of the time averages were estimated by the blocking method [29] in order to take into account time correlation of the data.

Note that, in order to calculate the full set of variables, two steps were necessary: a first step to calculate average values of solid fraction, velocity, and contact stresses and a second step to calculate terms containing velocity fluctuations. Two steps are needed since velocity fluctuations are defined with respect to the mean field.

C. Results

In the following we will apply the method to a simulation with a velocity of the top wall V top = 100d P/m; as stated above, the flow is at constant pressure, so the vertical position of the top wall evolved freely to a stationary value, y top ≈ 83.9d.

It is convenient to characterize such a shear using the socalled inertial number I that compares the typical time scale of microscopic rearrangements with the typical time scale of macroscopic deformations: I = γd/ P/ρ S . Note that I corresponds to the square root of the Savage number [14] and is also called the Coulomb number [30]. Our system corresponds to a rapid flow with a global value of the inertial number corresponding to I > 1; such a flow regime was chosen because kinetic and contact stresses have the same order of magnitude, allowing us to test the results of the previous sections. A snapshot of the system at steady state is given in Fig. 1. All profiles shown in the following are calculated using an averaging window of 2d unless otherwise stated.

Figure 2(a) shows the velocity profile at steady state. As highlighted in the figure, the profile is clearly linear except from two zones (∼ 10d wide) near the bumpy walls, where the shear rate increases approaching the boundary.

The solid fraction profile [Fig. 2(b)] displays a central zone with a constant value, and a decreasing behavior approaching the walls. The state of the system, characterized by a lower compaction near the walls, can be appreciated also from Fig. 1. In fact, all the variables suggest that the system is in a simple shear condition (uniform shear, stresses, no energy flux) in the center and deviates from this state only in two narrow bands (∼ 10d wide) near the bumpy walls. Another feature of the system that can be appreciated in Fig. 2(b) is a slight asymmetry of the profiles near the walls. As it could be seen also in Fig. 1, in the shear zone near the upper wall the material is denser and less agitated than in the shear zone near the bottom wall. This slight difference is related to the asymmetry of boundary conditions: while the bottom wall is fixed, the upper one is free to move vertically in order to impose a constant normal stress. This yields slightly asymmetric profiles for all the variables. Figure 3(a) displays the fluctuating energy profile calculated with the new procedure developed in this work for an averaging window of 2d. To compare the two approaches (the present one and Babic's one) we report in Fig. 3(b) the quantity k -1, where k is the ratio of the fluctuating energy profile calculated using one of the methods with different averaging sizes to a reference profile obtained with the present method and an averaging size equal to D = 1d. It clearly shows that the Babic's method leads to an estimate of the fluctuating energy which increases with the domain size. On the contrary, the present method leads to an estimate which does not depend on the size of the averaging window. However, if the latter size is too large it influences the averages close to the walls due to the inhomogeneity of the shear in such an averaging window. To illustrate this point, we have reported in Fig. 3(b) results obtained with D = 5d. Note also that, using the present method, the fluctuating energy profiles for D = 1d and D = 2d are almost equal for any position y. Figure 4(a) displays the profiles of the components of the contact part of the stress tensor (here mainly due to collisions, since the flow is quite rapid). Again, a plateau is found at the center; when approaching the walls all the components decrease in absolute value.

At the same time [Fig. 4(b)], the components of the kinetic part of the stress tensor increase when approaching the walls. This is a well-known direct consequence of the decreasing of solid fraction. In our simulation, near the walls, kinetic stresses are larger than collisional ones.

The sum of the two parts gives the total stress tensor, shown in Fig. 4(c). As predicted by the stress balance, the (xy) and (yy) components of the total stress tensor are uniform along y. Concerning the third part of the stress tensor, T γ , it is clear that according to Eq. ( 24) its only nonzero component is the (xx) one. For the reference simulation it was verified that, apart for a small difference in a thin zone (∼ 2d) close to the walls (where slight deviations from local homogeneity appear due to the shape of the wall), D y ≈ D m,y / √ 12, confirming the scaling given above.

Figure 5 collects contact and kinetic flux profiles. All the components of these fluxes differ from zero only when approaching the wall; as regards the y component of the fluxes, fluxes are negative when the fluctuating energy decreases with y and positive otherwise. This agrees with the sign convention used for q in the derivation, and supports a Fourier-like expression for the flux dependence on fluctuating energy gradient. On the other hand, the horizontal components of the fluxes are not negligible near the walls. This contrasts with a linear transport theory since there is no gradient of fluctuating energy in the x direction. Figure 6 displays the estimates of the stress power related to contact and kinetic contributions and of the dissipative terms related to affine and non affine motions. Firstly, it can be seen that far from the walls, the algebraic relations hold

T c † : ∇v = γ γ , T k † : ∇v = γ k , (30)
therefore justifying the separate dissipation paths for kinetic and contact stress powers. Then, near the walls, the energy fluxes seem to influence only the kinetic part, therefore supporting the result obtained in the previous section: only kinetic terms enter in a nonlocal energy balance, while contact stress power is dissipated locally.

V. DISCUSSION

In the previous sections we developed an averaging procedure for granular materials which was based on a weighted average plus an original fluctuation decomposition. Through numerical simulation results we proved the strength of the approach which gives consistent predictions for kinetic terms, and-more importantly-gives a new picture of the energy cascade in granular flows, at least for simple flows. In this section we discuss the validity of such an approach. It is clear that the assumption of local homogeneity holds only when the solid fraction and the shear rate profiles do not display large gradients with respect to the averaging domain size. For the simulations presented above, two zones were identified near the walls where such gradients exist. The independence of the fluctuating energy profile on the size of the averaging window and the respect of the stress balance are two proofs that the assumption is still acceptable.

Figure 7 allows us to discuss another important point: the physical meaning of the variables obtained by the method, with a particular regard to solid fraction. In the figure, the solid fraction profile was computed by two other methods: estimating the fraction of (1) a horizontal line and (2) of a horizontal slice (1d wide) intersecting the particles. These two methods have a clear geometrical meaning, the second being a moving average of the first.

For the rapid flow simulation which was taken as the reference [Fig. 7(a)], where no layering occurs, the three methods agree very well. Figure7(b) displays the results of the three methods for a denser flow (V top = 10, y top ≈ 50,⇒ I ≈ 0.2). Please note that the quite ordered bumpy wall was explicitly prepared to display layering. In this case the estimates on lines strongly fluctuate near the wall, and the average on slices is a smoother profile which still fluctuates, while the estimate produced by the present averaging method seems to be intermediate between the two. Far from the wall results from the methods coincide. Therefore, in the presence of layering, the solid fraction calculated by use of the present method may not have the direct geometrical meaning that its name implies. This problem may be reduced by choosing a smoother weighting function (in this paper we chose on purpose a very abrupt weighting function); we also suspect that in three dimensions this rather pathological issue may be less important.

Another issue caused by layering is that the strong spatial fluctuation of variables may imply that some of the assumptions of the local homogeneity hypothesis do not hold: The distribution of particles around the averaging point is not in principle homogeneous. Fluctuations are evidently associated with strong local gradients of solid fraction and shear rate. Therefore, when strong layering occurs, care has to be taken in applying the present method due to the approximations therein.

The inhomogeneous distribution of particles inside the averaging domain yields also a deviation of the virtual center of mass X from the center of the domain. This is not, in principle, a problem, because the averaging procedure is related to a volume and not to a point; average variables may be harmlessly referred to the virtual center of mass instead of the simple geometric center of the averaging domain. In other words the relation X = x which seemed to be a consequence of the local homogeneity assumption can be taken as a definition of x, the point where averages are supposed to be related.

In addition, results from simulations at V top = 1 and V top = 10 (which show quite strong layering), not shown here for the sake of brevity, show that stress balance is respected and confirm locality of the dissipation of the contact stress power. However, if layering occurs, local homogeneity may be difficult to assume, and therefore the simplified balances obtained in this work have always to be critically analyzed when dealing with regions very close to the walls.

VI. CONCLUSIONS

In this work, we propose a new derivation of continuum balance equations for granular materials which is inspired by the work of Babic [7]. We introduce a new decomposition of particle scale velocities into mean and fluctuations, taking into account velocity gradients. A very reasonable local homogeneity hypothesis allows us to simplify the balances, ensuring at the same time representativity and physical meaning of the continuum estimates obtained.

Two important results are obtained. First, the method solves the problem of dependence of some variables (e.g., fluctuating energy, kinetic stresses) on the averaging domain size which is intrinsic to previous methods [7] by clearly identifying scale-dependent terms. Second, the development suggests a new physical picture for the energy cascade in granular flows: kinetic terms enter in a nonlocal energy balance, while contact terms (containing both collisions and long-lasting contacts) seem to be dissipated locally. Both results are verified against discrete element simulations of granular flows, with respect to which the method is shown at work. This point is particularly important with regards to the hope of obtaining a statistical mechanics for granular flows whose predictions can be faithfully tested in physical experiments and numerical simulations. Note also that the question of the size of the averaging window close to a boundary [11] is problematic since it is bounded by the distance between the center of a grain and that boundary.

A last note has to be added on the domain of application of the present method. Discrete element methods are sometimes divided into smooth and nonsmooth ones [START_REF] Jean | Micromechanics of Granular Materials[END_REF]. The present averaging method (as Babic's one) applies rigorously only to smooth methods, since the idea of smooth (differentiable) particle fields was implicit in the expression of the equation of motion. The adaptation of such a method to non smooth dynamics is a very important point which will be the subject of future works.

FIG. 1 :

 1 FIG.1: Snapshot of the simulated system. A 2D granular material is sheared between two bumpy walls. The lower wall is fixed, the upper one moves horizontally at a constant velocity Vtop = 100 and is submitted to a constant vertical stress P = 1.

FIG. 2 :FIG. 3 :

 23 FIG. 2: (Color online) (a) x velocity and (b) solid fraction profiles along the y direction (Vtop = 100). The system appears to be in a simple shear condition (linear velocity profile, constant solid fraction) in the center of the cell. The size of the averaging window is 2d. (a) (b)

FIG. 4 :

 4 FIG. 4: (Color online) Components of the stress tensor [contact part (a), kinetic part (b), and total tensor (c)] along the y direction. The upper wall velocity is Vtop = 100 and the size of the averaging window is 2d.

FIG. 5 :

 5 FIG. 5: (Color online) Kinetic (a) and contact (b) fluxes of fluctuating energy along the y direction (Vtop = 100). The size of the averaging window is 2d. (a) (b)

FIG. 7 :

 7 FIG. 7: (Color online) Layering close to a bumpy wall: solid fraction computed using the present method with an averaging window of 1d, on lines and on slices (1d wide) for (a) Vtop = 100, and (b) Vtop = 10.

Acknowledgments

We are indebted to M. Y. Louge for fruitful discussions and to S. McNamara for a critical reading of the manuscript. The numerical simulations were carried out at the CCIPL (Centre de Calcul Intensif des Pays de la Loire) under the project MTEEGD.

Appendix A: Derivation of balance equations

The derivation follows closely that of Babic, except for the choice of the fluctuation decomposition. The evolution equation for a generic particle property ψ p is:

We multiply by w p , sum on p, and integrate on time:

The left hand side term can be developed into [7]:

As regards the right hand side of Eq. (A2), Babic also showed that the interaction term could be developed as:

where w F pq and w S pq are two weighting functions respectively defined as:

and

Now let us see the effect of the fluctuation decomposition on Eq. (A3). Decomposing the velocity field with the new decomposition leads to:

The divergence term on Eq. (A3) can be written as:

Decomposing also ψ p yields for the second term in Eq. (A8):

where X is the average center of mass of the averaging domain, the integral in the third term of Eq. (A8) becomes

Following the homogeneity assumption, we neglect correlations of fluctuations and of particle positions, and we consider that the average center of mass coincides with the averaging point. This means that terms under the divergence operator in Eq. (A8) simplify to

Developing the k component of the third term in Eq. (A12) [using the definition ∆x pi = (x pi -x i )] yields:

Assumption of local homogeneity also implies decorrelation of the components of particle position with respect to the averaging point. This means that:

) where D i is related to the scale of the averaging volume in the i direction. If D is the vector with components D i , the third term of Eq. (A12) is

This term renders explicit the scale dependence which was otherwise implicit in Babic's method. Equation (A3) is finally reduced to:

As regards to the effect of the fluctuation decomposition on the right hand side of Eq. (A2), the result depends on the particular property ψ p considered. Thus in the following we will analyze each property separately. If mass conservation is considered (ψ p = 1, P pq = 0, g p = 0), the right-hand-side term of Eq. ( A2) is zero and it is straightforward to see that the classical continuity equation is obtained [Eq. [START_REF] Schäfer | [END_REF].

If force balances are considered (ψ p = v p , P pq = f pq , g p = g), Eq. (A3) becomes

The last two terms can be considered as related to two components of an effective stress tensor as defined by Eqs. (11) and (10). On the other hand, the right-hand side of Eq. (A2) is not affected by the nature of the fluctuation decomposition. Substitution of the P pq = f pq and g p = g yields to the standard momentum equation, Eq. ( 8), with the third component of the effective stress tensor given by Eq. ( 9).

Let us derive now the translational kinetic energy equation. In this case the reference particle property is

The average translational kinetic energy is

and the fluctuation decomposition allows us to split this quantity into six terms, three of which disappear due to the local homogeneity assumption, yielding:

where, with respect to Babic's development, ρε T is likely to be scale independent, and a new term appears, which is the translational kinetic energy related to affine deformation within the average volume ρE γ T . The divergence term in Eq. (A3) corresponds, for the kinetic energy, to:

When developing this term, the fluctuation decomposition gives 27 terms, most of which can be deleted as in the previous development through (1) the postulate X = x, (2) the assumption of decorrelation among components of x px, and (3) the assumption of decorrelation between components of x px and velocity fluctuations. It can be shown that the irreducible terms are

where

Let us treat the terms coming from the right hand side term of Eq. (A2). First, we must identify the particle scale interaction terms. Taking the product of v p and Newton's equation,

Therefore, the particle interaction term is

and the particle source term is

Concerning the source term, local homogeneity implies:

(A27) Developing the interaction terms yields:

Using the fluctuation decomposition, the relation f pq = -f qp , and the assumption of local homogeneity this becomes:

which correspond respectively to:

where the first two terms are the conversion rates of translational kinetic energy into other forms of energy respectively due to velocity fluctuations (-γ k ) and to affine deformations (-γ γ ); q c can be seen as a translational kinetic energy flux related to contact forces. The translational kinetic energy balance equation is therefore given by Eq. ( 16), which can be further simplified to Eq. ( 19) when subtracting the mean flow energy balance.

Appendix B: Derivation of a scaling for D

In the following, a simplified scaling for the vector D which enters in the definition of T γ is derived using the local homogeneity hypothesis. Let us assume equal-sized spheres, and a step function on space and time as a weighting function. Let us consider variables computed at a point (x, t) in space and time. The average density is therefore:

where w p = 1/(V T ) for particles residing in the averaging volume V during the period (t -T /2, t + T /2) and zero otherwise. If we call N p the time average of the number of particles lying in the averaging volume, it is clear that ρ = m p N p /V . Developing the definition of the vector D, we find:

where x pi is the i component of the particle position. The local homogeneity hypothesis says that X = x and that x pi is uniformly distributed around x i . This yields for large N p T to:

Therefore if the local homogeneity holds, D i is likely to scale on the size of the averaging domain in the i direction. The proportionality will then depend on the choice of the weighting function, on polydispersity, on shape, and so on.