Jean Jacques Moreau

NUMERICAL DYNAMICS OF GRANULAR MATERIALS

Contact Dynamics is a numerical method, suitable for computing the dynamical motion of large collections of rigid bodies, with Coulomb friction taken into account in the event of contact. The principles of the method are sketched, in particular the way possible collisions or other nonsmooth features of the evolution are handled. As applications to granular dynamics, the construction of dry deposits and banks is simulated, in order to investigate their microstructure: force chains, geometrical anisotropy, Cauchy stress and some unexpected features of force transmission.

Introduction 1.Addressed problems

In diverse domains, computation methods are needed for the statics or the dynamics of collections of rigid bodies subject to the constraints of non-interpenetrability, with friction taken into account in the event of contact. This includes the dynamics of machines, in particular robots, the dynamics of masonry works submitted to transient actions (earthquakes, gusts of wind or impacts), animated computer graphics and numerical simulation in granular mechanics. Possibly, some variables are added in order to also account for a certain deformability of the bodies [START_REF] Jean | Computational Modeling of Contact and Friction[END_REF] without essentially changing the computational strategies. Of course, for the handling of non-interpenetrability and friction, much may be learned from the rich literature devoted to the treatment of contact between deformable media, discretized through finite elements; see e.g. [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like methods[END_REF][5] [START_REF] Christensen | Formulation and comparison of algorithms for frictional contact problems[END_REF][7] [START_REF] Vola | Consistent time discretization for a dynamical frictional contact problem and complementarity techniques[END_REF].

The main part of the lecture at the Symposium consisted of projected animations, intended to show what sort of knowledge may be gained from the numerical simulation of granular processes. The present paper begins with an overview of the 'Contact Dynamics' time-stepping technique used in these simulations. More detailed accounts of the method, originated in [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], may be found in [START_REF] Moreau | Contact et frottement en dynamique des systèmes de corps rigides[END_REF] [START_REF] Moreau | An introduction to unilateral dynamics[END_REF]. Finally, some drawings created in the course of computations are commented.

Analytical setting

The analytical description of a multibody system begins with the choice of a parametrization -at least a local one -of the set of possible configurations through an element q := (q 1 , . . . , q n ) of R n . In terms of this parametrization, the geometric effect of the non-interpenetration of the system members, or of their confinement by external obstacles with prescribed motion, is assumed expressed by a finite set of inequalities g α (t, q) ≥ 0, α ∈ {1, . . . , κ},

where the functions g α (t, q) are measures of the gaps between bodies, counted as negative if the configuration q at time t involves overlap. The event of equality in [START_REF] Abadie | Dynamic simulation of rigid bodies : modelling of frictional contact[END_REF] for some value of α, corresponds to contact between a member B of the system and a body B which may be another member of the system or an external obstacle. Non-interpenetrability is a unilateral constraint; as always in Mechanics, its description cannot reduce to merely specifying the corresponding set of feasible states (such, in contrast, is the meaning of the word 'constraint' in Optimization and similar contexts). Some information must be added about the mechanical process through which the geometrical restriction is enforced. For instance the use of servomechanisms to secure (1) should result in other motions than those investigated in what follows.

Generically, let us conceive a contact law as a relationship involving the contact force R α experienced by B from B at the contact point M α and the local velocity U α of B relative to B , relationship which a priori depends on time and on the configuration attained by the system, say

law α (t, q, U α , R α ) = true. ( 2
)
On numerical and analytical purposes, the definitions of the elements U α and R α of E 3 are extended, at least locally, to configurations with nonzero g α (t, q). Saying that (2) models a contact phenomenon involves that this law yields R α = 0 whenever g α (t, q) > 0. Of course, more complex situations may be addressed, for instance with contact actions not reducing to single forces. The system motion on a time-interval I, with origin t 0 , is described by a function t → q(t) of I into R n . Classical Kinematics requires of this function to be locally absolutely continuous, so that its derivative q makes an element of L 1 loc (I, R n ) from which q may be retrieved through Lebesgue integration. The components of q (t) are commonly referred to as the velocity components of the system with regard to the parametrization in use.

When dealing with stereodynamical problems one may also characterize the system velocity by some otherwise defined set of components, making an element u of R n whose dependence on t may generically be called the velocity function. Commonly, one attaches to each rigid body an orthonormal frame emanating from its mass center and principal with regard to body inertia. The components relative to these axes of the spin vector of the body are entered as elements of the column u(t). Retrieving from the function u the evolution of q then rests on the integration of elementary kinematical relationships. The moments, about the same axes, of the forces experienced by the body are concomitantly entered at the corresponding places in the row of the covariant components of the forces acting on the system.

The classical framework of 'smooth' dynamics requires in turn of u to be locally absolutely continuous on the interval I. By Lagrange's method or any stereodynamical practice, one writes the dynamical equations as the following equality of elements of R n

A(t, q) du dt = F (t, q, u) + α r α , (3)
to hold for almost every t. Here A denotes in general the n×n inertia matrix , diagonal and constant in (t, q) if u has been constructed according to the technique of mass centers and principal axes mentioned above. Expression F comprises certain standard terms, sometimes referred to as 'centrifugal' and 'gyroscopic', and also the covariant components of possible applied forces, supposed given as functions of time, the position of the system and its velocity. The elements r α , α ∈ {1, 2, . . . , κ}, are made of the covariant components of the respective contact forces.

In order to take the contact laws (2) into account, one has to connect the elements r α and u of R n with the vectors R α and U α of E 3 . By the kinematical analysis of the way the parametrization (q) has been constructed, one obtains an expression of the form

U α = G α u + W α , (4)
where G α : R n → E 3 denotes a linear mapping, depending on t and q. In case the contact labelled α involves an external obstacle, the term W α ∈ E 3 , a known function of t and q, accounts for the given motion of this obstacle. No attention is paid at this stage to the imagined motion preserving contact or not. By applying the standard virtual power machinery, one correlatively obtains the following expression for the covariant components of the contact forces acting on the system at M α (a pair of opposite forces if B is a member of the system, a single force if it is an external obstacle)

r α = G * α R α , (5)
with G * α : E 3 → R n denoting the transpose of G α . As said before, the elements associated with each α, such as G α here or a common normal unit n α , directed toward B, invoked in the sequel, are conventionally extended in a smooth way to neighbouring configurations without requiring g α = 0.

Approximation

Relationships (1) to ( 5) convey all the retained information about the considered mechanical system. For the numerical treatment of evolution problems, it is easy to conceive time-stepping procedures of approximation to the differential equation (3), under the kinematical relationship connecting q an u. Through (4) and ( 5), the contact laws (2) are entered into the management of each time-step. The difficulty comes from the necessity of taking the non-interpenetration inequalities (1) also into account.

In the majority of the techniques described in literature and applied in commercial software packages, this is achieved through the classical trick of penalizing inequalities. Mechanically, this means that the strict non-interpenetrability of a pair of bodies is approximately replaced by elastic repulsion forces which become effective when the bodies come close to each other. Such a procedure amounts to approximate the joint conditions (1) and (2) by relationships connecting the same variables, but smoother. In principle, the resulting system of differential equations should be sufficiently regular for standard time-stepping procedures to apply. Similar interaction laws are used in the computer simulation of molecular motions, hence the name of Molecular Dynamics (abbr.: MD) commonly given to this approach.

The drawback is that the need of precision requires of the artificial repulsion laws to be very stiff. Numerical stability in integration then calls for the use of very short time-steps and frequently also for the introduction of some artificial damping or artificial increase of inertia. When treating dynamical applications, the effect of such alterations of the mechanical data may blur the picture. Significant simulations of loose, collisional, flows of granular materials have been obtained in that way, but when dense collections of bodies are concerned (masonry works or compact granulates) the method is mainly applied to quasi-static evolutions in which only a succession of equilibrium states is looked for.

For these reasons, other numerical strategies may be preferred in unilateral multibody dynamics.

In the line of the traditional stereodynamic practice are approaches qualified as event-driven (abbr: ED). Starting from a state in which certain contacts are in effect, one attempts to calculate the subsequent motion under the provisional assumption that these contacts persist and that their status -sticking or sliding -in regard to Coulomb law is invariant. Calculation is then the same as with bilateral constraints, but requires to watch the evolution of some indicators. If, for instance, the normal components of some contact forces take directions incompatible with unilaterality, one concludes that the subsequent motion has to be calculated otherwise. Finding the further status of contacts without having to try all combinations is a nontrivial matter. It is usually reduced to Linear Complementarity Problems [START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF][29][32] [START_REF] Abadie | Dynamic simulation of rigid bodies : modelling of frictional contact[END_REF].

The ED strategy becomes unpractical if the number of contacts in presence is too large. Some timestepping integration procedures with preselected time-intervals have been developed instead [START_REF] Pang | A unified approach to discrete frictional contact problems[END_REF] [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of Painlevé's problem[END_REF]. The following one belongs to this class.

The 'Contact Dynamics' approach 2.1 Handling non-interpenetration in terms of velocity

For every label α, it has been agreed to extend, at least locally, the definition of contact descriptors such as the normal unit n α , to configurations with nonzero g α . Then put

K α (t, q) := {V ∈ E 3 | V.n α ≥ 0} if g α (t, q) ≤ 0 E 3 otherwise.
If g α (t, q) ≥ 0, this is the set of the values of U + α , the contact right-velocity, which are kinematically compatible with non-interpenetration.

The elements for the proof of the following may be found in [START_REF] Moreau | Some basics of unilateral dynamics[END_REF]:

Lemma. Let a motion be defined by a locally integrable velocity function u : I → R n . If the corresponding expression of U α : I → E 3 , as it results from (4), satisfies U α ∈ K α (t, q) for almost every t and if the non-interpenetration inequality (1) holds at the initial instant t 0 , then it holds throughout I.

The decisive move of the Contact Dynamics strategy consists in complementing the contact laws (2) so as to obtain relationships in the same variables, containing the same stipulations as before and additionally the following ones:

• in all cases U α ∈ K α , • if U α ∈ interior K α , then R α = 0.
In other words, these relationships secure the implications

g α (t, q) ≤ 0 ⇒ n α . U α ≥ 0 (6) n α . U α > 0 ⇒ R α = 0. (7)
We propose to say that a package of information concerning the possible contact labelled α, if it possesses these properties, is a contact law of prospective type. The underlying idea is that such a law does not properly govern the values of U α and R α at the actual instant, but their limits on the right of this instant, assumed to exist. In fact, if n α . U α . > 0, the concerned instant is followed by a contactless time-interval. Since R α must vanish over this interval, the same should be true for its right-limit.

Once all the contact laws (2) have been complemented in that way, the Lemma shows that condition (1) becomes redundant, provided it is satisfied at the initial instant. It then becomes easy to imagine time-stepping procedures for the approximation of the evolution problem.

The very structure of the Coulomb model of friction induces to organize the step computation as an implicit scheme, at least with regard to the velocity function u, i.e. it is the unknown value of this function at the step end which, as an approximation, is entered into the law. A contact law of the prospective type is adequate here, since step end values mimick right-side limits. The possibility of contacts to get loose at every instant is then managed automatically, without needing any analysis of complementarity conditions.

Each time-step makes a nonlinear problem whose computational cost grows rapidly with the number of contacts in presence. As a compensation for this cost, the method accepts much larger steplengths than those required by the regularization techniques referred to in 1.3. Usually an iteration procedure à la Gauss-Seidel is applied. It amounts to treat cyclically a sequence of dynamical problems, each of them involving a single contact. In dense collections of bodies with slowly evolving contact list, computation efficiency may be greatly enhanced by taking as starting guess in iterations the contact forces calculated at the antecedent time step, whenever the same contact was already in effect.

The Contact Dynamics strategy may also be applied with other solvers used at each time step [START_REF] Johansson | Study of frictional impact using a nonsmooth equation solver[END_REF]. A single computation step is needed to check whether a given position of the system is that of a possible equilibrium. One just have to launch iterations with zero initial velocity; if the step-end velocity is found zero, the values obtained for the contact forces agree with equilibrium.

Nonsmoothness in time and collisions

The sudden occurrence at some instant t c of a new contact, i.e. a collision, is expected to generate a velocity jump, the dynamics of which cannot be governed by the differential equation [START_REF] Brogliato | Nonsmooth Mechanics[END_REF]. Even in the absence of events of this sort, dry friction at a contact point has long been known, in some cases, to forbid the existence of any smooth solution beyond a catastrophic instant, on the left of which the contact forces, as well as the derivative u , may become unbounded [START_REF] Génot | New results on Painlevé paradoxes[END_REF]. This is a dynamic analogue to the locking situations familiarly met in the statics of frictional systems. Around year 1900, such a lack of solution for an apparently well set problem seemed inadmissible to P. Painlevé and induced him to question the very concept of a contact force. Hence the misleading denomination of 'Painlevé's paradox' which has been given to the observation. Today, one is accustomed to see models reaching the limits of their validity domain. Such is the case for smooth dynamics at a catastrophic instant of any sort, leading to enlarge the function space where u is looked for [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of Painlevé's problem[END_REF].

A readily available mathematically framework is that of R n -valued functions of locally bounded variation on the time-interval I; notation lbv(I, R n ) [START_REF] Moreau | Bounded variation in time[END_REF]. With every u in this space, an R n -valued measure on I is classically associated that we shall denote by du: this is the differential measure (or Stieltjes measure) of u. In the smooth case, this measure admits the derivative u as density function relative to dt (the Lebesgue measure on I). A natural generalization of (3) therefore is

A(t, q) du = F (t, q, u) dt + α ds α , (8)
where the R n -valued measures ds α constitute the 'components', relative to the parametrization in use, of the respective contact impulses dS α which are E 3 -valued measures. This is a measure differential equation. Concerning the theoretical foundation of such an extension of Classical Dynamics, one may refer to [START_REF] Moreau | An expression of classical dynamics[END_REF]. On a time interval of smooth motion dS α = R α dt, while at the instant t c of a collision, each measure dS α is expected to present an atom whose value is nothing else than the percussion vector of the traditional theory of shocks, say P α . The latter is a priori unknown and the question arises whether it may be involved in some phenomenological relationship playing at instant t c a role similar to that of (2) on intervals of smooth motion. The laws of restitution of Newton, Poisson or others amount to pragmatic tricks of this sort, but it is now clear [START_REF] Stoianovici | A critical study of the concepts of rigid body collision theory[END_REF][12] that the empirical coefficients involved in such laws cannot in general be identified a priori. Apart the simple case of the collision of two, otherwise free, nearly spherical objects, the outcome of a collision can be predicted with some precision only through a detailed analysis of the contact interaction which actually involves a certain amount of deformation of the concerned bodies. When so studied, the process takes place on a nonzero time interval, contact forces and accelerations exist in the proper sense, leaving smooth mechanics valid. Another classical approach consists of a multiple scaling technique, using some microtime variable in the analysis of the process.

A critical issue, when restitution laws are used, is the energy balance of the instant process. A theoretical framework securing the dissipativity required by the Second Law of Thermodynamics may be found in [START_REF] Frémond | Rigid body collisions[END_REF].

The CD handling of nonsmooth evolutions

Formally, at each time-step of a CD computation scheme, say [t i , t f ] ('i' as in initial , 'f' as in final), t f = t i + h, the algorithm is ready to face collisions. Throughout this interval, the elements A, G α , W α are approximated by constants, namely the value they take at the 'midpoint' (t m = t i + h/2, q m = q i +h u i /2) and the element F by the value it takes at (t m , q m , u i ). It is also at (t m , q m ) that inequalities [START_REF] Abadie | Dynamic simulation of rigid bodies : modelling of frictional contact[END_REF] are checked in order to determine the set J m of the labels α to be treated as effective. By integrating both members of (3) over the interval, or both members of (8) as well , one obtains

u f = u i + h A -1 F + A -1 α p α , (9)
where the quantities p α are the covariant components of the contact impulses, i.e. the integrals of the ds α . The relationship that (2), ( 4), ( 5) establish between r α and u, does not commute in general with time integration because contact laws are nonlinear. When the motion one intends to approximate is smooth, the variation of r α and u over the time-step may be small enough for the decision of connecting p α with the final velocity u f by the said relationship to merely reflect the choice made in 2.1 of a time-stepping strategy of the implicit type. But if some collisions occur in [t i , t f ], i.e. J m includes labels it did not contain at the antecedent step, large variations of u are expected. The designation of any value of u to be connected with p α , constitutes an assertion about the collisional process, the phenomenological quality of which has to be checked in regard to physical circumstances.

To u f corresponds the final value U αf of the local velocity. The above choice amounts to connect with it, through the admitted contact law, the total local impulse on [t i , t f ]. If this interval contains a single collisional instant t c , the initial and final values of U α may be viewed as approximants of the left and right limits U - α (t c ) and U + α (t c ). Therefore, the collisions met in the course of the above computation are treated according to the following rule: for each contact α present at instant t c , the contact percussion P α is related, through the admitted contact law, to the right limit U + α (t c ). Contact laws here are assumed of the prospective type. In view of ( 6) and ( 7), this yields the implication P α = 0 ⇒ n α . U α = 0, expressing that the collision is soft, i.e. completely inelastic. But the possibility n α . U α > 0 with P α = 0 is left open. This model of soft collision in multicontact systems thus improves on the mere assumption of zero Newton restitution coefficient, because all the contacts present at instant t c are involved together in the dynamics.

An extension of the preceding enables CD algorithms to produce bouncy collisions: it consists in connecting, through the admitted contact law, the total local impulse with some weighted mean of U αi and U αf [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][25] [START_REF] Moreau | An introduction to unilateral dynamics[END_REF]. The choice of the weighting coefficients (possibly distinct for the normal and tangential components of velocities) directly determines some restitution rules, improving as before on those of Newton.

Anyway, such laws of restitution do not elude the criticism made in 2.2. Their advantage is only that their implementation in a CD algorithm adds nothing to the computation cost.

Numerical simulation in Granular Mechanics

What can be expected from computation

Granular mechanics has long been a topic in Civil Engineering as a part of Soil Mechanics. For a certain number of years, it has also attracted the attention of another public, coming from the Condensed Matter and Disordered Systems community. It is expected that the forms of statistical reasoning used in the latter domain could provide some insight into the many intriguing features of the granular behaviour. On this purpose, the model consisting of rigid grains which interact only through contact and Coulomb friction, without any adhesive effect, makes a sufficiently rich concept for investigating some fundamental questions.

In contrast dry sand, to which the above model is relevant, is not the most frequently met soil in Civil Engineering. In the numerical treatment of most problems, the soil material is viewed as a continuous medium, for which some empirical constitutive law should be available. Finite Element softwares are then applied to the corresponding boundary value problem. The difficulty lies in the determination of the constitutive law which anyway has to be tested against all available experimental information about the considered soil. Engineering needs may be met hrough this approach, but no insight into the intimate granular behaviour can be gained.

While the abbreviation FEM (for Finite Element Method) is used to refer to the above strategy, the representation of a granular sample as a collection of solids is called a Distinct Element Method (abbr. DEM, sometimes also read as Discrete Element Method). Numerical simulation may be conducted thanks to the methods referred to in the foregoing. Clearly the boundary problems of Civil Engineering cannot directly be handled that way since the number of grains involved in real soils exceeds the possibility of computers. The availability of scaling rules which would allow one to approximate these boundary problems through DEM with oversized grains is still an object of investigation [START_REF] Nouguier | Force fluctuations in a pushed granular material[END_REF].

Force chains

It is known from experiments with assemblies of photoelastic grains that the transmission of forces in the bulk of a granular material is concentrated in chains of grains more loaded than the surrounding ones.

A two-dimensional model of railway track is shown on Figure 1. In view of the elastic deformation of the rail, the force exerted by a passing wheel is distributed between the closest supporting concrete blocks, in relation with their distances. The blocks, as well as the ballast grains, are drawn with levels of gray corresponding to their respective loads, i.e. the sum of the normal components of the experienced contact forces. This makes visible that forces are transmitted to the ground across the ballast layer along some preferred chains of grains.

Deposition anisotropy

The response of a grain collection to some mechanical action strongly depends on the way it has been prepared or, more generally, on its past evolution. Some message from history is therefore written in the current state of the granulate. As a part of this message the statistical distribution of the directions of the normals at the grain to grain contacts is commonly viewed as a descriptor of a possible microstructural anisotropy. Here is an intriguing example.

A usual way of preparing granulates consists of producing a random rain of grains under gravity. If this rain is received on a fixed horizontal ground, a certain statistical anisotropy of the collected granular mass is expected, since the vertical direction plays a special role in the deposition process. The surprise was to discover experimentally that the most frequent directions of contact normals are not vertical nor horizontal but at some angle on both sides of the vertical. The physical experiment was performed with the model of two-dimensional granulate called a Schneebeli material, i.e. a collection of cylindrical rods of equal lengths stacked parallel and observed laterally; it was well reproduced by CD numerical simulations [START_REF] Daudon | A micromechanical comparison between experimental results and numerical simulation of a biaxial 2D granular material[END_REF]. Instead of circular objects, the two examples presented on Figure 2 concern convex irregular polygonal grains of random shapes, with two different distributions of sizes. In each sample, for all contact points detected in the rain-deposited layer, the normal directions -in an extended sense, since contact may affect a vertex -are recorded. Histograms of the distributions are constructed by dividing the 360 • range about the origin into 24 equal sectors. Each sector is filled in gray up to a radius proportional to the number of contact points whose normals have the corresponding directions. For legibility or for quantitative treatment, one may choose to fit a continuous curve to the distribution, the graph of a probability density expressed as a function of the angle θ by a linear combination of cosine and sine of multiples of 2θ. The same preference for two oblique directions as with circular grains is found.

Three-dimensional simulations yield similar results: in the granular layer created by raining spherical grains of dispersed sizes, contact directions of largest frequency make an angle with the vertical.

The stress tensor

In the Civil Engineering situations referred to in 3.1, the granular material is treated in the framework of classical Continuum Mechanics, with the Cauchy stress field as central concept. This macroscopic approach does not prevent one, at the stage of assessing the constitutive laws wich govern the medium behaviour, to refer to smaller scale features. Connecting with micromechanical quantities some tensor which could be identified as the stress is not trivial [START_REF] Goddard | Continuum modelling of granular assemblies[END_REF]. Since the very time of Cauchy, this question has mainly been addressed for assemblies of points subject to some interaction potentials. Granular media offer a very different sort of microstructure. On Figure 3 there is shown a part of a two-dimensional collection of irregular polygonal grains (grain to grain friction: 0.4) submitted to a 'biaxial test', i.e. the granulate is squeezed between frictionless rectilinear vertical boundaries in imposed motion while some containing pressure is exerted by horizontal ones. The representation of contact forces as line segments makes the force chains visible. The programme allows one to superimpose a circular probing area over which the average stress is computed according to the definition introduced in [START_REF] Moreau | Numerical investigation of shear zones in granular materials[END_REF], equivalent to some other ones found in literature. If this tensor is used in the standard way to express the normal and tangential components of the tension upon a cut with direction ϑ, the results are sinusoidal functions of ϑ whose graphs are drawn in gray. On the other hand, 31 diameters of the probing area have been selected. By computing, for each of them, the resultant of the contact forces exerted by the grains on one side upon the grains on the other side, one estimates the corresponding tensions. When the plots of the normal and tangential components of these vectors are compared with the preceding sinusoids, one concludes that the computed tensor, in spite of the small size of the sample and the inhomogenity due to force chains, plays with satisfactory precision the role expected from the Cauchy stress in continuous media. The cross drawn in the circular area represents the principal stresses.

Banks created by avalanches of dry grains

Figure 4 has been created in the course of the following two-dimensional simulation. Random polygonal grains are deposited at the top of a rigid incline, made of a collection of fixed grains with similar characteristics, so as to simulate roughness. Thanks to a landing plate placed at the bottom of the incline, a layer of some thickness accumulates, sporadically perturbed by avalanches. The angle that the free surface makes with horizontal fluctuates between two values: the angle of repose θ rep (approx. 25 • ), observed after an overall avalanche has remodelled the layer, and some larger value θ lim (approx. 29 • ), the limit angle, attained as the result of progressive accretion. At θ lim , the layer becomes instable. Similarly to that was done in 3.4, the programme allows one to choose a line segment and to calculate the resultant contact force transmitted across such an imagined cut.

For a line segment parallel to the free surface (the difference between θ rep and θ lim is neglected here), one finds a vertical force. This feature looks natural since the bank is made of layers successively deposited by avalanches. When such a layer stops, the supporting forces it experiences from the material beneath should equilibrate its weight, and this verticality of transmitted forces is expected to persist after other layers possessing the same property have been superimposed. This confirms the assumption commonly made in Civil Engineering that a freshly deposited bank is everywhere in a state of incipient failure.

For a vertical line segment, one finds a force parallel to the free surface. This observation, which seems novel, is connected in an essential way with the preceding one through the symmetry of the Cauchy stress tensor. It provides a convincing insight into the following subject which, in recent years, started a lot of speculation and controversy.

Let a conical pile be created by pouring grains from a source onto a rough rigid horizontal ground. Some experimentalists have been surprised to observe that the distribution of the pressure exerted upon the ground was not proportional to the height of the material above and even that a local minimum of ground pressure may be present at the vertical of the apex.

Actually, the fact that the free surface is not horizontal makes clear that the statics of the granular assembly has nothing to do with hydrostatics. In a three-dimensional CD numerical experiment [START_REF] Moreau | An introduction to unilateral dynamics[END_REF] not shown here, we have reproduced the construction of a conical pile of 14.000 spherical grains with dispersed sizes. Computed results reveal first that a test surface portion placed inside the peripheral bank of the pile and parallel to the free surface transmits a vertical force, indicative of a state of incipient failure. Secondly, a test surface in the shape of a cylinder coaxial with the conical pile is found to transmit average forces parallel to the free surface. As before, this couple of observations reflects the symmetry of the stress tensor. If such a cylinder is used to delimit some 'central core' of the pile, this core consequently has part of its weight supported by an arching action from the surrounding bank, which explains the experimental findings.

Figure 1 :

 1 Figure 1: A two-dimensional model of railway ballast.

Figure 2 :

 2 Figure 2: Histograms of contact directions.

Figure 3 :

 3 Figure 3: Quality of the stress tensor estimate.

Figure 4 :

 4 Figure 4: Two cuts in a bank.


-Salim Rao -

ESTIMATION OF THE PROPORTIONAL HAZARD PREMIUM UNDER RANDOM CENSORING

Keywords: Heavy tails, Proportional hazard premium, Random censoring

In this paper, we investigate the estimation of the proportional hazard premium for randomly right-censored risks. The asymptotic normality of the proposed estimator is established under mild conditions.

Introduction

Let X 1 , ..., X n be n ≥ 1 independent copies of a non-negative random variable (rv) X, with continuous cumulative distribution function (cdf) F. An independent sequence of independent rv's Y 1 , ..., Y n with continuous cdf G censor them to the right, so that at each stage j we only can observe Z j := min(X j , Y j ) and the variable δ j := 1 {X j ≤ Y j } (with 1 {•} denoting the indicator function) informing whether or not there has been censorship. This model is very useful in a variety of areas where random censoring is very likely to occur such as in biostatistics, medical research, reliability analysis, actuarial science,... In insurance, the worst scenarios are those caused by extreme events such as natural catastrophes, human-made disasters and financial crashes. These events increase the bill of insurance and reinsurance companies. A typical requirement for actuaries is the determination of adequate premiums for such risks. Usually, the insurer's claims data do not correspond to the underlying losses, because they are censored from above, since the insurer stipulates an upper limit to the amount to be paid out and the reinsurer covers the excess over this fixed threshold. This kind of reinsurance is called excess-of-loss reinsurance and the upper limit has distinct designations that are specific to each insurance type.

For instance, in life insurance, it is called the cedent's company retention level while in non-life insurance, it is called the deductible, where the losses should be treated separately.

Model and estimators

Let us assume that both F and G are heavy-tailed, that is there exist two constants γ 1 > 0 and γ 2 > 0, called tail indices or extreme value indices, such that

F (z) ∼ z -1/γ 1 ℓ 1 (z) and G(z) ∼ z -1/γ 2 ℓ 2 (z), as z → ∞, (2.1)
where ℓ 1 and ℓ 2 are slowly varying functions at infinity, i.e. lim z→∞ ℓ i (xz)/ℓ i (z) = 1 for every x > 0, i = 1, 2. Throughout the paper, we use the notation S(x) := S(∞) -S(x),

for any function S (x) of x > 0. If relations (2.1) hold, then we have, for any x > 0

lim z→∞ F (xz) F (z) = x -1/γ 1 and lim z→∞ G(xz) G(z) = x -1/γ 2 , (2.2)
and we say that F and G are regularly varying at infinity as well, with respective tail indices -1/γ 1 and -1/γ 2 , which we denote by F ∈ RV -1/γ 1 and G ∈ RV -1/γ 2 . Note that, in virtue of the independence of X and Y, the cdf of the observed Z ′ s, that we denote by H, is also heavy-tailed and we have -i+1:n] being the proportion of upper non-censored observations.

H ∈ RV -1/γ with γ := γ 1 γ 2 /(γ 1 + γ 2 ).
(i.e. δ [i:n] = δ j if Z i:n = Z j ), then Hill's estimator of γ 1 adapted to censored data is defined as γ (H,c) 1 := γ H / p, where γ H := k -1 k i=1 log (Z n-i+1:n /Z n-k:n ) represents Hill's estimator of γ, with k = k n being an integer sequence satisfying 1 < k < n, k → ∞ and k/n → 0 as n → ∞, (2.3) and p := k -1 k i=1 δ [n
In the excess-of-loss reinsurance treaty, the ceding company covers claims which do not exceed a (high) number R ≥ 0 (called retention level), while the reinsurer pays the part

(X i -R) + := max (0, X i -R) of each claim beyond R. Applying Wang's premium cal-
culation principle with a distortion function equal to x 1/ρ , one defines what is called the proportional hazard premium (PHP), where ρ ≥ 1 represents the distortion parameter or the risk aversion index. Then, the PHP of loss for the layer from R to infinity is defined as follows:

Π ρ (R) := ∞ R F (x) 1/ρ dx,
which may be rewritten into

Π ρ (R) = R(F (R)) 1/ρ ∞ 1 F (Rx) F (R) 1/ρ dx. Karamata theorem entails Π ρ (R) ∼ ρ 1/γ 1 -ρ R F (R) 1/ρ , 0 < γ 1 < 1/ρ, for large R. Since F ∈ RV -1/γ 1 , then F (x) ∼ F (h) (x/h) -1/γ 1 as x → ∞, where h = h n := H ← (1 -k/n
) with H ← (y) := inf {x : H (x) ≥ y} , 0 < y < 1, denoting the quantile function pertaining to H. We refer to [13,36,14,19,30,10,31,20,5,8,18] for other works linking risk theory with extreme-value statistics. The above formula leads us to derive a Weissman-type estimator for the distribution tail F for censored data as follows:

F (x) = x Z n-k:n -1/ γ (H,c) 1 F n (Z n-k:n ) .
In the context of randomly right censored observations, the nonparametric maximum likelihood estimator of F is given as the product limit estimator

F n (x) := Z i:n ≤x 1 - δ [i:n] n -i + 1 = Z i:n ≤x n -i n -i + 1 δ [i:n] , for x < Z n:n , which gives F n (Z n-k:n ) = n-k i=1 1 - δ [i:n]
n-i+1 . Thus, the distribution tail estimator is of the form

F (x) := x Z n-k:n -1/ γ (H,c) 1 n-k i=1 1 - δ [i:n] n -i + 1 ,
and consequently, we define the PHP estimator as follows:

Π ρ (R) := ρR 1/ γ (H,c) 1 -ρ R Z n-k:n -1/ ρ γ (H,c) 1 n-k i=1 1 - δ [i:n] n -i + 1 1/ρ .
The outline of the paper is as follows. In Section 3, we state our main result that consists in the asymptotic normality of the newly proposed estimator Π ρ (R), which we prove in Section 5. Finally, technical results are postponed to the Appendix.

Asymptotic distribution

It is well-known that the asymptotic normality of extreme value theory based estimators is adequately achieved within the second-order framework. Thus, it seems quite natural to suppose that cdf's F and G satisfy the well-known second-order condition of regular variation. That is, we assume that there exist two constants τ j ≤ 0 (called second-order parameters) and two functions A j , j = 1, 2, tending to zero and not changing sign near infinity, such that for any x > 0

lim t→∞ F (tx)/F (t) -x -1/γ 1 A 1 (t) = x -1/γ 1 x τ 1 /γ 1 -1 γ 1 τ 1 , lim t→∞ G(tx)/G(t) -x -1/γ 2 A 2 (t) = x -1/γ 2 x τ 2 /γ 2 -1 γ 2 τ 2 .
(3.4) Theorem 3.1. Assume that the second-order conditions of regular variation (3.4) hold, with 0 < γ 1 < 1/ρ and let k = k n be an integer sequence satisfying, in addition to

(2.3) , √ kA 1 (h) → λ 1 . Assume further that R/h → 1. Then √ k Π ρ (R) -Π ρ (R) (R/h) -1/ργ 1 R F (h) 1/ρ D → N µ, σ 2 , as n → ∞,
where

µ := ρλ 1 (1 -pτ 1 ) (1 -ργ 1 ) 2 + λ 1 ρ (γ 1 + τ 1 + ρ -2) (2 -ρ -γ 1 )
, and

σ 2 := γ 2 1 (1 -ργ 1 ) 2 p (2 -p) + ρ (p -1) (1 -ργ 1 ) + ρ 2 (1 -2p) p (1 -ργ 1 ) 2 .

Further work

When the studied phenomena X depends on a covariate, the associated risk measure also depends on X and one has to deal with conditional extreme-value analysis. This branch of statistics has become very active these past ten years, the main contributions to this domain are listed below:

• Theoretical issues: [33,34,32,16,23,25,12,11,22,26,21,24,4,17,9,15,3,1,39,38] • Quantile regression: [28,35] • Application to finance: [6,2,27,7,29,37] Our goal will be to adapt our methodology to this framework.

Proof of the main result

Before starting the proof of the main theorem, let us give a brief introduction on some uniform empirical processes under random censoring. To this end, let us define

H (j) (v) := P (Z ≤ v, δ = j) , j = 0, 1; v ≥ 0,
which have a prominent role to play in the random censorship setting. Their empirical counterparts are defined by

H (j) n (v) := 1 n n i=1 1(Z i ≤ v, δ i = j), j = 0, 1; v ≥ 0.
In the sequel, we will use the following two empirical processes

√ n H (j) n (v) -H (j) (v) , j = 0, 1; v ≥ 0,
which may be represented, almost surely, by a uniform empirical process. Indeed, let us define, for each i = 1, ..., n with θ := H (1) (∞) , the following rv

U i := δ i H (1) (Z i ) + (1 -δ i )(θ + H (0) (Z i )).
The rv's U 1 , ..., U n are independent and identically distributed according to the (0, 1)uniform law. The empirical cdf and the uniform empirical process based upon U 1 , ..., U n are respectively denoted by

U n (s) : = 1 n n i=1 1(U i ≤ s) and α n (s) := √ n(U n (s) -s), 0 ≤ s ≤ 1.
One has almost surely

H (0) n (v) = U n (H (0) (v) + θ) -U n (θ), for 0 < H (0) (v) < 1 -θ,
and

H (1) n (v) = U n (H (1) (v)), for 0 < H (1) (v) < θ.
It is easy to verify that we almost surely have √ n H

(1)

n (v) -H (1) (v) = α n (θ) -α n θ -H (1) (v) , for 0 < H (1) (v) < θ, (5.5) and √ n H (0) n (v) -H (0) (v) = -α n 1 -H (0) (v) , for 0 < H (0) (v) < 1 -θ.
(5.6)

Our methodology strongly relies on the well-known Gaussian approximation given in 6.1.

For our needs, we use the following form:

sup 1/n≤s≤1 n ζ |α n (1 -s) -B n (1 -s)| s 1/2-ζ = O P (1).
(5.7)

For the increments α n (θ)α n (θs), we will need an approximation of the same type as (5.7). Following similar arguments, we may show that, for every 0 < θ < 1 and

0 ≤ ζ < 1/4, we have sup 1/n≤s≤θ n ζ |{α n (θ) -α n (θ -s)} -{B n (θ) -B n (θ -s)}| s 1/2-ζ = O P (1).
(5.8)

The following Gaussian processes will be crucial to our needs:

B n (v) := B n (θ) -B n θ -H (1) (v) , for 0 < H (1) (v) < θ,
(5.9) and

B * n (v) := B n (v) -B n 1 -H (0) (v) , for 0 < H (0) (v) < 1 -θ.
(5.10)

Proof of Theorem 3.1. In the sequel, for two sequences of rv's, we write

V (1) n = o P V (2) n and V (1) n ≈ V (2)
n , as n → ∞, to say that V

(1) n /V

(2) n → 0 in probability and

V (1) n = V (2) n (1 + o P (1)) respectively. With the premium Π ρ (R) = R(F (R)) 1/ρ ∞ 1 F (Rx) F (R) 1/ρ dx,
and its estimator

Π ρ (R) = ρR 1/ γ (H,c) 1 -ρ R Z n-k:n -1/ ρ γ (H,c) 1 F n (Z n-k:n ) 1/ρ , it is easy to verify that √ k Π ρ (R) -Π ρ (R) (R/h) -1/(ργ 1 ) R F (h) 1/ρ = 5 i=1 S ni ,
where

S n1 := ρ 1/ γ (H,c) 1 -ρ F (Z n-k:n ) F (h) 1/ρ F n (Z n-k:n ) F (Z n-k:n ) 1/ρ × √ k        (R/Z n-k:n ) -1/ γ (H,c) 1 (R/h) -1/γ 1   1/ρ -1      , S n2 := F (Z n-k:n ) F (h) 1/ρ F n (Z n-k:n ) F (Z n-k:n ) 1/ρ √ k ρ 1/ γ (H,c) 1 -ρ - ρ 1/γ 1 -ρ , S n3 := ρ 1/γ 1 -ρ F (Z n-k:n ) F (h) 1/ρ √ k F n (Z n-k:n ) F (Z n-k:n ) 1/ρ -1 , S n4 := ρ 1/γ 1 -ρ √ k F (Z n-k:n ) F (h) 1/ρ -1 , and
S n5 := √ k ρ 1/γ 1 -ρ - F (R) /F (h) 1/ρ (R/h) -1/(ργ 1 ) ∞ 1 F (Rx) F (R) 1/ρ
dx .

We will represent the first three terms S ni , i = 1, 2, 3, in terms of the Gaussian processes B n and B * n and we will show that S n4 P → 0 while S n5 converges to a deterministic limit.

For the first term S n1 , we have γ

F n (Z n-k:n ) /F (Z n-k:n ) P → 1. It follows that S n1 = S (1) n1 + S (2)
n1 , where

S

(1)

n1 := (1 + o P (1)) ργ 1 1 -ργ 1 × √ k    Z n-k:n h 1/ ρ γ (H,c) 1 -1      R h 1/γ 1 -1/ γ (H,c) 1   1/ρ , and
S
(2)

n1 := (1 + o P (1)) ργ 1 1 -ργ 1 √ k        R h 1/γ 1 -1/ γ (H,c) 1   1/ρ -1      .
For S

(1) n1 , we use the mean value theorem, the consistency of γ (H,c) 1

and the fact that

Z n-k:n /h P → 1, to have S (1) n1 = (1 + o P (1)) 1 1 -ργ 1 √ k Z n-k:n h -1 .
Next,

S

(1)

n1 = (1 + o P (1)) γ 1 -ργ 1 n k B * n (h) .
In view of the consistency and asymptotic normality of γ (H,c) 1

and the assumption R/h → 1, we show, by applying the mean value theorem twice, that S

(2) n1 = o P (1). Thus, we end up with

S n1 = (1 + o P (1)) γ 1 -ργ 1 n k B * n (h) + o P (1).
(5.11)

By similar arguments and using the mean value theorem once again, we easily show that

S n2 = (1 + o P (1)) ρ (1 -ργ 1 ) 2 √ k γ (H,c) 1 -γ 1 , S n3 = (1 + o P (1)) γ 1 1 -ργ 1 √ k F n (Z n-k:n ) F (Z n-k:n ) -1 ,
and

S n4 = (1 + o P (1)) γ 1 1 -ργ 1 √ k F (Z n-k:n ) F (h) -1 .
We get, after a change of variables, that

S n2 = (1 + o P (1)) ρ (1 -ργ 1 ) 2 1 p n k ∞ 1 v -1 B * n (hv) dv - γ 1 p n k B n (h) + (1 + o P (1)) ρ √ kA 1 (h) (1 -pτ 1 ) (1 -ργ 1 ) 2 .
(5.12)

From Proposition 6.2, we infer that

S n3 = (1 + o P (1)) γ 1 1 -ργ 1 n k B n (h) + k n ∆ n + o P (1) .
(5.13) Now, we decompose S n4 into the sum of two terms

S

(1)

n4 := (1 + o P (1)) γ 1 1 -ργ 1 √ k F (Z n-k:n ) F (h) - Z n-k:n h -1/γ 1 , and
S
(2)

n4 := (1 + o P (1)) γ 1 1 -ργ 1 √ k Z n-k:n h -1/γ 1 -1 .
The second-order condition (3.4) of F and the fact that Z n-k:n /h

P → 1 yield that S (1) n4 = o P √ kA 1 (h) = o P (1) .
For S

(2) n4 , we, once again, apply the mean value theorem (with Z n-k:n /h

P → 1) to get S (2) n4 = -(1 + o P (1)) γ 1 -ργ 1 n k B * n (h) .
Consequently, we have

S n4 = -(1 + o P (1)) γ 1 -ργ 1 n k B * n (h) + o P (1) .
(5.14)

For the last term S n5 , we start by decomposing it into the sum of

S

(1)

n5 := - ργ 1 1 -ργ 1 1 (R/h) -1/(ργ 1 ) √ k    F (R) F (h) 1/ρ - R h -1/γ 1 1/ρ    , and
S
(2)

n5 := - F (R) /F (h) (R/h) -1/γ 1 1/ρ √ k ∞ 1 F (Rx) F (R) 1/ρ -x -1/γ 1 1/ρ dx.
By similar arguments as those used for S

(1) n4 , we show that (here we use the assumption that R/h → 1)

S (1) n5 = o P √ kA 1 (h) = o P (1) .
For S

(2) n5 , we first apply the mean value theorem to have

S (2) n5 = - 1 ρ √ k ∞ 1 F (Rx) F (R) -x -1/γ 1 ζ 1/ρ-1 (x)dx,
where ζ lies between F (Rx)/F (R) and x -1/γ 1 . Then we use Potter's bounds to get

S

(2)

n5 = (1 + o (1)) √ kA 1 (h) ρ (γ 1 + τ 1 + ρ -2) (2 -ρ -γ 1 )
.

Therefore

S n5 = (1 + o (1)) √ kA 1 (h) ρ (γ 1 + τ 1 + ρ -2) (2 -ρ -γ 1 )
+ o P (1) .

(5.15)

Finally, by gathering results (5.11) , (5.12) , (5.13) , (5.14) and ( 5.15) , we obtain the following asymptotic representation to the premium estimator:

√ k Π ρ (R) -Π ρ (R) (R/h) -1/ργ 1 R F (h) 1/ρ = o P (1) + γ 1 1 -ργ 1 k n ∆ n + 1 1 -ργ 1 n k Γ n (5.16) + ρ √ kA 1 (h) (1 -pτ 1 ) (1 -ργ 1 ) 2 + √ kA 1 (h) ρ (γ 1 + τ 1 + ρ -2) (2 -ρ -γ 1 )
,

where ∆ n is as defined in 6.18 and

Γ n := γ 1 1 - ρ p (1 -ργ 1 ) B n (h) + ρ p (1 -ργ 1 ) ∞ 1 v -1 B * n (hv) dv.
From (5.16) , we deduce that

√ k Π ρ (R) -Π ρ (R) / (R/h) -1/ργ 1 R F (h) 1/ρ is asymp- totically Gaussian with mean ρ (1 -pτ 1 ) (1 -ργ 1 ) 2 + 1 ρ (γ 1 + τ 1 + ρ -2) (2 -ρ -γ 1 ) lim n→∞ √ kA 1 (h) = µ,
and variance

lim n→∞ E γ 1 1 -ργ 1 k n ∆ n + 1 1 -ργ 1 n k Γ n 2 .
Note that from the covariance structure, we have the following useful formulas:

           E [B n (u) B n (v)] = min H (1) (u) , H (1) (v) -H (1) (u) H (1) (v) , E [B * n (u) B * n (v)] = min H (u) , H (v) -H (u) H (v) , E [B n (u) B * n (v)] = min H (1) (u) , H (1) (v) -H (1) (u) H (v) .
(5.17)

After elementary but very tedious computations, using these formulas with l'Hôpital's rule, we get as n → ∞,

h 0 E [B n (u) B n (h)] H 2 (u) dH (u) → -p, h 0 E [B n (h) B * n (u)] H 2 (u) dH (1) (u) → -p 2 , h 0 ∞ 1 E [B n (v) B * n (hu)] uH 2 (v) dudH (v) → -pγ, h 0 ∞ 1 E [B * n (v) B * n (hu)] uH 2 (v) dudH (1) (v) → -pγ, k n h 0 h 0 E [B n (u) B n (v)] H 2 (u)H 2 (v) dH(u)dH (v) → 2p, k n h 0 h 0 E [B * n (u) B * n (v)] H 2 (u)H 2 (v) dH (1) (u)dH (1) (v) → 2p 2 , and k n h 0 h 0 E [B n (u) B * n (v)] H 2 (u)H 2 (v) dH(u)dH (1) (v) → 2p 2 ,
Collection the results above with some further calculations leads to σ 2 .

Appendix

Proposition 6.1. There exists a probability space (Ω, A, P) with independent (0, 1)-uniform rv's U 1 , U 2 , ... and a sequence of Brownian bridges {B i (s); 0 ≤ s ≤ 1} (i = 1, 2, ...) such that, for every 0 < λ < ∞, we have as n → ∞

sup λ/n≤s≤1 n ζ |α n (s) -B n (s)| s 1/2-ζ =      O P (log n) when ζ = 1 4 , O P (1) when 0 ≤ ζ < 1 4 , sup 0≤s≤1-λ/n n ζ |α n (s) -B n (s)| (1 -s) 1/2-ζ =      O P (log n) when ζ = 1 4 , O P (1) when 0 ≤ ζ < 1 4 , and
sup λ/n≤s≤1-λ/n n ζ |α n (s) -B n (s)| (s(1 -s)) 1/2-ζ =      O P (log n) when ζ = 1 4 , O P (1) when 0 ≤ ζ < 1 4 .
In the next basic proposition, we provide an asymptotic representation to the Kaplan-Meier product limit estimator in Z n-k:n . This result will be of prime importance in the study of the limiting behaviors of many statistics based on censored data exhibiting extreme values.

Proposition 6.2. Assume that all second-order conditions (3.4) hold. Let k = k n be an integer sequence satisfying, in addition to (2.3) √ kA j (h) = O(1), for j = 1, 2, as n → ∞.

Then there exists a sequence of Brownian bridges {B n (s); 0 ≤ s ≤ 1} such that

√ k F n (Z n-k:n ) F (Z n-k:n ) -1 = n k B n (h) + k n ∆ n + o P (1) ,
where (6.18) with B n (v) and B * n (v) respectively defined in (5.9) and (5.10) . Consequently,

∆ n := h 0 B n (v) H 2 (v) dH (v) - h 0 B * n (v) H 2 (v) dH (1) (v) ,
√ k F n (Z n-k:n ) F (Z n-k:n ) -1 d → N (0, p(1 -p)) , as n → ∞, (6.19)
Proof. We have for any

x ≤ Z n-k:n , F n (x) -F (x) F (x) = x 0 d H (1) n (v) -H (1) (v) H (v) - x 0 H n (v) -H (v) H 2 (v) dH (1) (v) + O P 1 k .
Upon integrating the first integral by parts, we get

F n (x) -F (x) F (x) (6.20) = -H (1) n (0) -H (1) (0) + H (1) n (x) -H (1) (x) H (x) + x 0 H (1) n (v) -H (1) (v) H 2 (v) dH (v) - x 0 H n (v) -H (v) H 2 (v) dH (1) (v) + O P 1 k . Recall that √ n H n (v) -H (v) = √ n H 1 n (v) -H 1 (v) + √ n H 0 n (v) -H 0 (v) ,
which by representations (5.5) and (5.6) becomes

√ n H n (v) -H (v) = α n (θ) -α n θ -H (1) (v) -α n 1 -H (0) (v) .
On the other hand, by the classical central limit theorem, we have H

(1)

n (0) -H (1) (0) = O P n -1/2
. Using these results in (6.20) and then multiplying by √ k, we get

√ k F n (x) -F (x) F (x) = O P k n + O P 1 √ k + k n α n (θ) -α n θ -H (1) (x) H (x) + k n x 0 α n (θ) -α n θ -H (1) (v) H 2 (v) dH (v) - k n x 0 α n (θ) -α n θ -H (1) (v) -α n 1 -H (0) (v) H 2 (v) dH (1) (v) .
The Gaussian approximations (5.7) and (5.8) , in x = Z n-k:n , and the facts that k/n and 1/ √ k tend to zero as n → ∞, lead to

√ k F n (Z n-k:n ) -F (Z n-k:n ) F (Z n-k:n ) = n k B n (Z n-k:n ) + k n Z n-k:n 0 B n (v) H 2 (v) dH (v) - k n Z n-k:n 0 B * n (v) H 2 (v) dH (1) (v) + o P (1) .
Applying Lemma 6.1 completes the proof. The asymptotic normality property is straightforward. For the variance computation, we use the covariance formulas (5.17) and the results at the end of Section 5.

Lemma 6.1. Assume that the second-order conditions of regular variation (3.4) and let k := k n be an integer sequence satisfying (2.3). Then

(i) k n Z n-k:n h B n (v) H 2 (v) dH (v) = o P (1) . (ii) k n Z n-k:n h B * n (v) H 2 (v) dH (1) (v) = o P (1) . (iii) n k {B n (Z n-k:n ) -B n (h)} = o P (1) (iv) n k {B * n (Z n-k:n ) -B * n (h)} = o P (1) .
Proof. We begin by proving the first assertion. For fixed 0 < η, ε < 1, we have

P k n Z n-k:n h B n (v) dH (v) H 2 (v) > η ≤ P Z n-k:n h -1 > ε + P k n (1+ε)h h B n (v) dH (v) H 2 (v) > η .
It is clear that the first term the right-hand side tends to zero as n → ∞. Then, it remains to show that the second one goes to zero as well. Indeed, observe that

E k n (1+ε)h h B n (v) dH (v) H 2 (v) ≤ - k n (1+ε)h h E |B n (v)| dH (v) H 2 (v) .
From the first result of (5.17) , we have

E |B n (v)| ≤ H 1 (v). Then E k n (1+ε)h h B n (v) dH (v) H 2 (v) ≤ - k n (1+ε)h h H 1 (v) dH (v) H 2 (v) ,
which, in turn, is less than or equal to

k n H (1) (h) 1 H ((1 + ε) h) - 1 H (h)
.

Since H (h) = k/n, then this may be rewritten into

H (1) (h) H (h) H (h) H ((1 + ε) h) -1 .
Since H

(1) (h) ∼ pH (h) and H ∈ RV (-1/γ) , then the previous quantity tends to p 1/2 (1 + ε) 1/γ -1 as n → ∞. Being arbitrary, ε may be chosen small enough so that this limit be zero. By similar arguments, we also show assertion (ii) , therefore we omit the details. The last two assertions are shown following the same technique, that we use to prove (iv). Notice that, from the definition of B * n (v) and the second covariance formula in (5.17) ,

{B * n (v) ; v ≥ 0} d = B n H (v) ; v ≥ 0 ,
where {B n (s) ; 0 ≤ s ≤ 1} is a sequence of standard Brownian bridges. Hence

n k {B * n (Z n-k:n ) -B * n (h)} d = n k B n H (Z n-k:n ) -B n H (h) .
Let {W n (t) ; 0 ≤ s ≤ 1} be a sequence of standard Wiener processes such that B n (t) =

W n (t) -tW n (1) . Then n/k {B * n (Z n-k:n ) -B * n (h)} equals in distribution to

n k W n H (Z n-k:n ) -W n H (h) -H (Z n-k:n ) -H (h) W n (1)
.

By using the facts that H (h) = k/n and H (Z n-k:n ) /H (h) ≈ 1, we get .

Since E W n ǫH (h) ≤ ǫH (h) and P {A c n (ǫ)} < ǫ, thus P (|ϑ n | > η) ≤ η -1 ǫ 1/2 + ǫ which tends to zero as ǫ ↓ 0, as sought.

  and Z n-k:n /h P → 1, which, in view of the regular variation of F , implies that F (Z n-k:n ) /F (h) P → 1. Moreover, from (6.19) we have

  o P (1) .Next we show thatϑ n := n k W n H (Z n-k:n ) -W n H (h) = o P (1) .Let η > 0 be a fixed real number and show that P (|ϑ n | > η) → 0, as n → ∞. Since Z n-k:n /h P → 1, then for an arbitrary ǫ > 0 and sufficiently large n, the probability ofA n (ǫ) := {|Z n-k:n /h -1| ≤ ǫ} is close to 1.Next, we will use the following useful inequality:P (|ϑ n | > η) ≤ P {|ϑ n | > η, A n (ǫ)}+P {A c n (ǫ)}, where A c n (ǫ) denotes the complement set of A n (ǫ) . It is easy to verify that ϑ n may be rewritten intoW n H (h) ξ n + H (h) -W n H (h) H (h) , where ξ n := H (Z n-k:n ) /H (h)-1. Since H isregularly varying, then we may show readily that, in the set A n (ǫ) , we have |ξ n | ≤ ǫ too, therefore P (|ϑ n | > η) ≤ I n + P {A c n (ǫ)} + P {A c n (ǫ)} , where I n := P sup 0≤t≤H(h)ξn


1/4

Reconstructing Old Chinese uvulars in th e Baxter-Sagart system (ver. 0.97) Laurent Sagart Centre de Recherches Linguistiques sur l'Asie Orientale, Paris ---1. Uvular stops were first reconstructed for OC by Pan [START_REF] Wuyun | [END_REF] as the source of MC "laryngeals":

 *q > ʔ (凩損)  *qh-> x (鉑損)  *ɢ-> ɦi (卩尿)
2. Advantage: phonetic series which mix these initials can be treated as showing normal phonetic alternations between stops (B-S 0.97 reconstructions from here onward):

 勣 *ɢʷ(r)a > hju > yu2 'go; at'  巘 *qqʷ(r)a > 'wae > wa1 'impure, vile'

 䯂 *qʷʰ(r)a > xju > xu1 'great' compare:  借 *[dz]a > dzjoX > ju3 'chew'  杓 *[ts]a-s > tsjoH > ju4 'marshy ground'  鍍 *[ts h ]A > tshjaeX > qie3 'moreover'
As observed by Pan, stops and fricatives of the same place of articulation (like t and s) do not normally mix in phonetic series.

uvulars and labio-uvulars.

It is necessary to reconstruct both plain uvulars: *q-, *q h -, *ɢ-

 忌 *qqat > 'at > è 'repress; cease'  劉 *q h at > xjot > xie4 'cease, rest'  褐 *[ɢɢ]at > hat > hè 'coarse cloth'
(note: in the B-S system, singleton consonants note type-B/division-III initials, while doubled consonants note type-A/non division-III initials )

and labio-uvulars: *q w -, *q wh -, *ɢ w -

2/4

 顪 *qʷʰat-s > xjwojH > huì 'beard'  穢 *qʷat-s > 'jwojH > huì 'bad weeds; filth' Shi Jing rhyme-words from the 歲 series rhyme as OC -at, not -ot; so that the labiality in MC words of that series listed above cannot be explained as due to a rounded vowel.

4. Middle Chinese reflexes of *ɢ w -, *ɢɢ w -and *ɢɢ-.

The Middle Chinese initial *hj-(卩尿) occurred primarily in words with a rounded vowel and words with medial -w-. Its OC precursor must have been inherently labial or labialized. Baxter and I agree with Pan that it is the MC reflex of *ɢ w -. It corresponds to g-in WT:

 迄 *ɢ w aŋ > hjwang > wáng 'king' (WT gong 'a superior one')  婿 *[ɢ w ]ə[t]-s > hjw+jH > wèi 'stomach' (WT grod-pa 'stomach')  勣 *ɢʷ(r)a > hju > yu2 'go; at' (WT Hgro 'go')  咎 *[ɢ]ʷat > hjwot > yuè 'transgress' (WT Hgrod-pa go, to travel')  叮 *[ɢ]ʷa[n] >
hjwon > yuán 'long robe' (WT gon 'garment') *ɢ w -is a type-B initial. Its type-A counterpart is *ɢɢ w -. Baxter and I also agree with Pan that the MC reflex of *ɢɢ w -is hw-(余):

 公 *ɢɢ w aŋ > hwang > huáng 'sovereign'  畫 *[ɢɢ]ʷrek-s > hweaH > huà 'draw'  緩 *[ɢɢ]ʷa[n]ʔ > hwanX > huǎn 'slack'
We think OC *ɢɢ-goes to MC h-(余) :

 后 *[ɢɢ](r)oʔ > huwX > hòu 'sovereign' (WT mgo 'head', WB u 'head')  褐 *[ɢɢ]at > hat > hè 'coarse cloth'
5. Distinguishing velar and uvular series. OC *ɢɢ w -and *ɢɢ-are difficult to distinguish from *gg -and *gg-, which also evolve to MC h-. We propose basing the decision on the phonetic series. We observed that when a particular rhyme category has two competing series with the same velar initials, it is usually the case that only one of those two also includes MC reflexes of OC uvulars: '-(glottal stop), x-or h-. Here are some paired phonetic series with vowel -o-:

rhyme just velars Mixed: velars + '-/x-/h- -o GSR 109 饅 GSR 108 T -o GSR 113 彊 GSR 112 后 -oŋ GSR 1172 ^ GSR1173 「 -ok GSR 1225 援 GSR 1226 独 3/4
The pure velar series we reconstruct with OC velars: 饅 kk-, 彊 gg-, ^ kk-, 援 kk-; and the mixed velar+"laryngeal" series we reconstruct with OC uvulars: T q(q)-, 后 ɢɢ-, 「 qq-, 独 qq( h )-.

We think 后 and 褐 have uvular initials because they occur in such mixed series.

6. Why are there velars in uvular series ?

We assume that OC uvular stops evolve to MC velars when preceded by minor syllables, consisting of a consonant plus ə, as in these examples of word families from one phonetic series: Examples of MC y-from OC *l-are found in phonetic series with MC d-and/or th-(but no t-), and correspond to TB l-:

 to rest  劉 *qʰat > xjot > xiē 'cease, rest'  *Cə.qʰrat > khjet > qì 'to rest'  empty  虛 *q h (r)a > xjo > xù 'empty'  虛 *Cə.q h (r)a > khjo > xù 'ruins'  shadow  凩 *qraŋʔ > 'jaengX > yǐng 'shadow'  竺 *Cə.qraŋʔ > kjaengX > jǐng 'bright; image'  old man  勲 *Cə-qqoŋ > kuwng >
 余 *la > yo > yú '1sg.' (prob. a polite form. Mikir la 'here, where the speaker is') kloong. These l-clusters, however, are not supported by comparative evidence (for 勲, compare WB hông 'old'). With our proposal, there is no need of l-clusters in these words.

 以 *mə-ləʔ > yiX > yǐ 'take, use' (Jingpo la 'take')  揚 *laŋ > yang > yáng 'lift, raise' (WT lang 'rise, get up')  枼 *l[a]p > yep > yè 'leaf' (Jingpo lap 'leaf')  易 *N(ə)-lek-s > yeH >

  gōng 'father; prince'  模 *qqoŋ > 'uwng > wēng 'old man' 7. What is the MC reflex of *ɢ-? We think it is y-(卩蒜). That is, MC y-has two distinct OC sources: (1) OC *l-, and (2) OC *ɢ-.

  yì 'easy' (WT legs 'good, happy, etc.')  用 *loŋ-s > yowngH > yòng 'use' (WT longs 'to use') Examples of MC y-from OC *ɢ-are found in phonetic series without MC d-and/or th-, but with velars and/or laryngeals initials. They often correspond to TB words in y-or g-y-:  羊 *ɢ(r)aŋ > yang > yáng 'sheep' (Matisoff PTB *yang 'sheep') 4/4  洋 *[ɢ](r)aŋ > yang > yáng 'a great expanse of water' (WT yangs-pa 'wide, broad, large')  祥 *s-[ɢ]aŋ > zjang > xiáng 'auspicious' (WT g-yang 'happiness, blessing, prosperity) Note that all three examples above are from the phonetic series of 羊 (GSR 732), from which Karlgren arbitrarly excluded  羌 *Cə.qʰaŋ > khjang > qiāng 'Western tribes'  姜 *Cə.qaŋ > kjang > jiāng 'a family name' More examples of MC y-from OC *ɢ-:  亦 *[ɢ](r)Ak > yek > yì 'armpit' (Benedict TB *(g-)yak 'armpit')  夜 *[ɢ]ak-s > yaeH > yè 'night' (Benedict TB *ya 'night') Although this phonetic series (GSR 800) has no velars or laryngeals, the word 亦 'armpit' has a related words with a MC velar initial: 胳 *kkak (or *Cə.qqak ?) > kak > gē 'armpit'. 8. word-family contacts between MC y and laryngeals/velars. Our proposal that OC *ɢ-> MC y-makes sense of phonetic series and word-family contacts between MC laryngeals/velars and y-:  桾 *Cǝ.q(r)a > kjoX > jŭ 'lift, raise'  冉 *ɢ(r)a > yoX > yŭ 'give; for; and'  集 *qek > 'jiek > yì 'increase'  泡 *[ɢ]ik > yit > yì 'overflow'  国 *[ɢ](r)a[n]ʔ > yenX > yǎn 'overflow'  僻 *Cə.qʰ(r)a[n] > khjen > qiān 'exceed, err' If the only source of MC y-was *l-, the voiceless members in the above examples would need l-clusters, for instance Sagart (1999) had 桾 b k-la, Zhengzhang (2003) has qleg, 勲

Harbin, 27-29 September, 2007


In Memoriam: Richard A. Toupin,

1926-2017

I had the pleasure to interact personally with Richard Toupin over three consecutive days in Montreal. We spoke for more than twenty hours during coffee breaks, lunches, dinners and walks in downtown Montreal, on the occasion of the 4th Canadian Conference on Nonlinear Solid Mechanics (CanCNSM2013) at McGill Uni versity, 23-26 July, 2013. He was attending a symposium in his honor organized by David Steigmann and myself. We announced the symposium with the following words: "Dr. Richard Toupin is one ofthe 20th century's most eminent contributors to theoretical mechanics, applied mathematics and physics. His deep contributions to the Lagrangian school of Continuum Mechanics stand as milestones in the history of science. His mastery of math ematical methods allowed him to advance the works of the founders of his discipline: Gabrio Pio la and Eugène and François Cosserat. He gave us a framework in which to understand microstructured continua as well as electromechanical phenomena."

I had the honor ofToupin's presence at the conference and was delighted to meet this originai, good-natured person and to enjoy his many anecdotes about the history of our subject, the beauty of which was so greatly enhanced by his works. However, I had already known his deep and somehow mysterious mind via the papers which he had written and which I had the good fortune to read during my PhD studi es and also later in my career. Richard Toupin had seemed to me a voice out of the choir: a harmonious voi ce accompanied by those of great contemporaries like Raymond Mindlin, Ronald Rivlin, Alleo Pipkin, Paul Germain, Leonid Sedov and Ekhart Kroner. Richard Toupin is also rightly famous for his co-authorship, with C. Truesdell, of the monumental treatise: "The Classical Field Theories", published by Springer-Verlag in 1960. This masterwork remains the single most important source for theoretical mechanics to the present day.

I recall here some aspects of Richard Toupin's personal history, as he explained it to me at the Montreal meeting. He was of course of French background, his ancestors having lived in Quebec. Richard was very proud of his French origins. He and I joked that perhaps his inclination towards variational principles was due to his French heritage and to the fact that Lagrange and D'Alembert were working in Paris. I managed to obtain from Richard a detailed CV: Born August 20, 1926 in Miami, he got his BS at the University of South Carolina in 1946; his MS at the University of Hawaii 1949 with a dissertation having the title: "On the Motion of Singularities in the General Theory of Relativity"; and his PhD at Syracuse University 1961 with a dissertation entitled: "A Dynamical Theory of Elastic Dielectrics". Toupin's subsequent papers on this subject paved the way for the modern theory of electro-magneto-mechanical interactions in continuous media. 1973 [with G. Fan] for work on magnetic ink jet printing. His groundbreaking paper "Elastic Materials with Couple Stresses" was the most cited of the genre in the decade 1972-1982. Indeed this seminal work laid the foundation for research on scale effects in elasticity, currently a very active field. He was a Founding Member of the Society for Natural Philosophy, 1963-86, andSecretary of the Society, 1965-67. He was member of the editorial boards of: The American Mathematical Society, Archive for Rational Mechanics and Analysis, and International Journal of Engineering Science.

While I was reading his CV, Richard became amused by my surprise that he never had been a professor. I asked him if that was a choice or a circumstance. How could it be possible that he did not get any university position with his extraordinary curriculum and scientific standing? He did not answer directly except to say that it was not a choice but rather a circumstance of life. I understood later that around 1961 he had turned down an offer from Ronald Rivlin to become a professor in the Applied Mathematics Department at Brown, as he preferred to go to IBM Research instead.

While telling me that he had studied the textbook of theoretical physics by Landau, and that he had always believed that computers were bound to dramatically change our lives, Richard started telling me how much he enjoyed playing tennis and that his life had never been only science or career. Gently, Richard then returned to his "self-education": He chose Landau as his guide at first, upon reading the introductory chapter in Volume 2 of "The Classical Theory of Fields", where Lagrangian methods are described as the most powerful tool for formulating new physical theories. Richard remarked that Landau's title was basically replicated in his treatise with Truesdell, adding: "Of course if you already know the result, then you can also get it without resorting to variational principles". Being proud of his mastery of French, he recalled the famous quote of Lagrange: "On ne trouvera point de Figures dans cet Ouvrage. Les méthodes que j'y expose ne demandent ni constructions, ni raisonnemens géométriques ou méchaniques, mais seulement des opérations algébriques, assujetties à une marche réguliere et uniforme. Ceux qui aiment l'Analyse, verront avec plaisir la Méchanique en divenir une nouvelle branche, et me sauront gré d'en avoir étendu ansi le domaine". Richard completed his digression about his self-education by emphasizing the great influence on him of the works by Levi-Civita. He referred to himself as a champion of Levi-Civita tensor calculus and expressed the view that all the efforts made to change Levi-Civita's index notation were bound to be a total failure, until another bright mind like Levi-Civita himself finds a drastic further innovation. While concluding that he did not see anybody who was capable of improving on this notation, and that all the efforts which he had seen merely concealed the true content of the theory, he

  During the period 1946-47 he was U. S. Navy Navigator aboard the U.S.S Andromeda; and during 1949-50 Instructor of Physics and Mathematics, University of Hawaii, Honolulu. He was active in research in several public and state research laboratories, as documented in the following list of his positions: 1950-52 Research Staff Member, Applied Mathematics Branch,U. S. Naval Research Laboratory, Washington, D.C.; 1952-54 Research Assistant, Graduate Student Syracuse University, Syracuse, N. Y.; 1954-6: Research Staff Member, Applied Mathematics Branch, U. S. Naval Research Laboratory, Washington, DC.; 1956-62: Head, Theoretical Mechanics Section, U. S. Naval Research Laboratory, Washington, DC.; 1962-66: Research Staff Member, Mathematical Sciences Department, IBM Research Center, Yorktown Heights, N. Y.; 1966-68: Visiting Research Staff Member, IBM Research Laboratory, Zurich, Switzerland; 1968-74: Manager, Theoretical Mechanics, Mathematical Sciences Department, IBM Research Center, Yorktown Heights, N. Y.; 1974-81: Director, Mathematical Sciences Department, IBM Research Center, Yorktown Heights, N. Y.; 1982-83: Consultant, Theoretical and Applied Mechanics, IBM General Products Division Development Laboratory, San Jose, CA.; 1983-86: Manager, Scientific and Engineering Computation, IBM Science Center, Heidelberg.; 1986-88: FRG Visiting Scholar, Mathematics and Applied Mechanics, Stanford University, and Consultant in Applied Mechanics, IBM Almaden Research Center. Richard Toupin was the recipient of several prizes: The Hulbert Award, 1961; U. S. Naval Research Laboratory Outstanding Scientist of the Year; New York Academy of Sciences Annual Invited Evening Lecture on Mechanics 1965: "St. -Venant and a Matter of Principle"; and the IBM Outstanding Contribution Award,

started talking about the "Lecture Notes in Physics" by R. Feynman. He did not stop talking about the chapter on variational principles until dinner was finally served.

After the end of his time at the Theoretical Mechanics Section of the U. S. Naval Research Laboratory in Washington, D.C., Richard Toupin started what he described as a new life. Having read Trusedell's paper "The Computer: Ruin of Science and Threat to Mankind" in "An Idiot's Fugitive Essays on Science" (Springer-Verlag, 1984, pp. 594-631), I had naively attributed to Toupin the same attitude against computers. Richard, instead, was very proud of this second life: He modestly outlined his role in some of the most important advancements in computer technology. His CV documents the patents and invention disclosures which he had obtained While telling me the story of his works in transforming his equations into real objects, Richard started talking about his achievements in the research of innovative technological solutions: indeed the hard disks we are using in our computers are all the result of Richard's understanding of frictional and electrical phenomena. Probably Richard was too modest in underestimating his contributions to this field, but perhaps this can be understood as reflecting the need of IBM to guard some secrets about its new devices. Richard was very proud of his last papers: 1. Toupin's papers, lectures and presentations prove that he was ahead of his time. The subjects treated in them will become topical after some decades, and indeed many recent papers rediscover periodically the results presented there. A list of his papers and lectures may be found at: http://www.fdellisola.it/wpcontent/uploads/2013/08/CVToupin.pdf At the end of the Montreal conference, when we knew that our conversations were about to end, Richard invited me to visit him at his home in Florida. I regret that I did not manage to organize this visit in time. Richard's health deteriorated in 2014 and our email exchanges ended abruptly. However, before saying goodbye he added that he had had a happy life, full of fun, with a nice group of friends and that he did not regret any choice he had made. He did not believe himself to be a genius, he had not invented anything in particular, but he had worked seriously in service to science and knowledge by honoring the human mind. I did not have the chance to tell him what I write here: Richard Toupin was a great scientist. His family will miss his humanity; humanity will miss his mind.