Emilien Azéma
email: azema@lmgc.univ-montp2.fr

Farhang Radjaï
email: radjai@lmgc.univ-montp2.fr

Stress-strain behavior and geometrical properties of packings of elongated particles

Keywords: numbers: 45.70.-n, 83.80.Fg, 61.43.-j

We present a numerical analysis of the effect of particle elongation on the quasistatic behavior of sheared granular media by means of the Contact Dynamics method. The particle shapes are rounded-cap rectangles characterized by their elongation. The macroscopic and microstructural properties of several packings subjected to biaxial compression are analyzed as a function of particle elongation. We find that the shear strength is an increasing linear function of elongation. Performing an additive decomposition of the stress tensor based on a harmonic approximation of the angular dependence of branch vectors, contact normals and forces, we show that the increasing mobilization of friction force and the associated anisotropy are key effects of particle elongation. These effects are correlated with partial nematic ordering of the particles which tend to be oriented perpendicular to the major principal stress direction and form side-to-side contacts. However, the force transmission is found to be mainly guided by cap-to-side contacts, which represent the largest fraction of contacts for the most elongated particles. Another interesting finding is that, in contrast to shear strength, the solid fraction first increases with particle elongation, but declines as the particles become more elongated. It is also remarkable that the coordination number does not follow this trend so that the packings of more elongated particles are looser but more strongly connected.

I. INTRODUCTION

Since a few years, the research for a better understanding of the complex rheology of granular media is enriched by an increasing focus on nonspherical particles [1][2][3][4][5][6][7][8][9]. The wide-spread use of spherical or disk-like particles has been motivated by the fact that the rheology of granular media is basically governed by the collective contact interactions of the particles so that the particle shape can be viewed as a secondary effect. In practice, both in experiments and discrete element simulations, the spherical or circular particles, such as glass beads and disks, are easier to handle and the results are generally more directly amenable to theoretical analysis. However, owing to the fast progress in experimental and numerical techniques during the last decade, there is now a wide scope for the investigation of materials composed of more complex particle shapes. In this respect, the model granular media with spherical particles provide a reference material for understanding the rheology when the particle shapes deviate from a spherical or circular shape [10,11].

A wide variety of particle shapes can be found in nature and industry: elongated and platy shapes, e.g. in biomaterials and pharmaceutical applications, angular and facetted shapes, e.g. in geomaterials, and nonconvex shapes, e.g. in sintered powders. The behavior under various types of loading is strongly influenced by particle shape. Rounded particles enhance flowability whereas angular shape is susceptible to improve shear strength. In many applications, the particle shapes need to be optimized in order to increase performance [12][13][14][15][16][17]. These trends are generally explained in qualitative terms and linked with the jamming of the particles.

The effect of particle shape is mediated by the specific granular texture (or fabric) induced by each particle shape. For example, it is found that hard ellipses can be jammed even though they are underconstrained [18][19][20][21][22][23]. In general, the anisometric or elongated particle shapes, such as spheroids and sphero-cylinders, tend to develop orientational order affecting force transmission and frictional behavior [3,[24][25][26]. This "nematic" ordering occurs while, in contrast to liquid crystals, the particles interact only via contact and friction [27].

In a sheared granular material, the local equilibrium structures are generically anisotropic in terms of contact directions and forces [5,[28][29][30][31][32][33][START_REF] Radjaï | 6th International Conference on the Micromechanics of Granular Media[END_REF]. It was recently shown that the fabric anisotropy in a sheared granular assembly crucially depends on particle shape [10,11]. In the case of polygonal and polyhedral particles, due to large contact area of side-to-side contacts, the fabric anisotropy appears to be marginal compared to force anisotropy [10,11]. Those contacts play a major role in force transmission by accommodating long force chains that are basically unstable in a packing composed of spheres.

The force and fabric anisotropies are at the origin of the enhanced shear strength of materials composed of nonspherical particles [11,24,32,[START_REF] Azéma | [END_REF]. The particle shape affects the compactness and dilatancy of granular materials. A nontrivial effect, evidenced recently by experiments and numerical simulations for spheroids, is the finding that the solid fraction is not a monotonous function of the aspect ratio [19][20][21][22]36]. The solid fraction increases linearly to a maximum and then declines in inverse proportion to the aspect ratio [37]. In powder processing, the particle shape appears also to be an important parameter controlling the flowability, discharge rates and compaction of powders [38,39].

In this paper, we use contact dynamics simulations to investigate the rheology of large packings of elongated particles with increasing aspect ratio. The particles are rectangles with rounded caps to which we will refer as Rounded-Cap Rectanglular (RCR) particles. These particles may be considered as 2D analog of spherocynlinders. The RCR shape can be characterized by a single aspect ratio α or, as we shall see, by an elongation parameter η varying from 0 to 1 as the particle shape varies continuously from a circle to an thin line. We are interested both in the properties of the static packings of RCR particles prepared by isotropic compaction without friction and in the stress-strain behavior under biaxial compression with finite friction between particles.

The macroscopic behavior is studied in terms of the internal angle of friction and solid fraction for different values of η. We find a nonmonotonous variation of the solid fraction and a nearly linear increase of the internal angle of friction with η. In order to understand the origins of this behavior and the role of particle shape, we perform a detailed analysis of the microstructure and stress transmission. We consider the organization of the particles and contacts in the simulated packings, as well as the stress transmission by means of a harmonic representation of the stress tensor in terms of force and fabric anisotropies. The microstructure is increasingly dominated by a short-range nematic ordering of particle orientations as η increases. We show that the internal angle of friction is influenced by this ordering via an increasing anisotropy of friction forces and contact orientations with the elongation parameter. For all values of the latter, the harmonic approximation provides an excellent fit to the shear stress.

An important feature of RCR particles is that, like polygonal particles, they have lineal edges and can thus form side-to-side contacts as well as side-to-cap and capto-cap contacts. Hence, in a packing of RCR particles, the texture can characterized by the networks of these various contact types, and the influence of the shape parameter on force transmission and shear strength may be analyzed in terms of these contacts and more specially the side-to-side contacts which are expected to play a stabilizing role in the packing.

In the following, we first introduce our numerical approach in Section II. Then, in Section III, the stress-strain behavior is presented for different values of η. The microstructure is analyzed in Section IV in terms of connectivity, orientations of the particles and the contact network. We also introduce the harmonic approximation of the stress tensor allowing us to track the origins of the internal angle of friction via force and fabric anisotropies. In Section V, we present an additive decomposition of the connectivity, anisotropies and forces as a function of different contact types. The force distributions are presented in Section VI. In section VII, we analyze the structure of force networks with cap-to-cap, cap-to-side and side-to-side contacts. We conclude with a summary and discussion of the most salient results of this work.

II. SYSTEM DESCRIPTION AND NUMERICAL PROCEDURES

In this section, we briefly introduce the contact dynamics (CD) method used for the simulations and the numerical procedures for sample preparation.

A. Contact Dynamic method

The simulations were carried out by means of the contact dynamics (CD) method [2,[START_REF] Jean | Proceedings of Contact Mechanics International Symposium (Presses Polytechniques et Universitaires Romandes[END_REF][START_REF] Moreau | [END_REF][42][START_REF] Moreau | Novel approaches in civil engineering[END_REF][START_REF] Radjai | 14e Congres Francais de Mécanique[END_REF][START_REF] Dubois | Actes du sixième colloque national en calcul des structures -CSMA-AFM-LMS[END_REF][START_REF] Renouf | [END_REF][START_REF] Radjai | Physics of dry granular media[END_REF][START_REF] Radjaï | [END_REF]. The CD method is a discrete element approach based on a nonsmooth approach in which an integrated form of the equations of dynamics. The integration time interval corresponds to the time step and may involve discontinuous variation of the velocities due to collisions. The frictional and collisional interactions are described as complementarity relations between the relative velocities between particles and the corresponding momenta at the contact points without elastic or viscous regularization. Thus, the condition of geometrical contact between two particles is expressed by the following mutually exclusive alternatives:

f n 0 and u n = 0, f n = 0 and u n > 0.

(

) 1
where f n is the normal contact force and u n the relative normal velocity between two particles in contact is counted positive when they move away from each other.

In the same way, the Coulomb friction law involves three mutually exclusive conditions:

f t = -µf n and u t > 0, -µf n f t µf n and u t = 0, f t = µf n and u t < 0, (2)
where u t is the sliding velocity at the contact, µ is the friction coefficient and f t is the friction force. Remark that this relation cannot be reduced to a (mono)valued functional dependence between the two variables as assumed in the Molecular Dynamics (MD) method. The above formulation is implicit in the sense that the complementarity relations should be satisfied for the velocities at the end of each time step. An iterative algorithm based on a nonlinear Gauss-Seidel scheme is used to solve the system of equations and complementarity relations for contact forces and particles velocities. The uniqueness is not a priori guaranteed for perfectly rigid particles. However, by initializing each step of calculation with the forces calculated in the preceding step, the set of admissible solutions shrinks to fluctuations which are basically below the numerical resolution.

The CD method is particularly suitable for the simulation of rigid nearly undeformable particles. In this limit, the MD method requires steep interaction potentials and thus very small time steps. Nevertheless, several comparisons between the two methods suggest that both methods are equally valid and efficient for the simulation of granular materials [2,[START_REF] Radjai | Physics of dry granular media[END_REF]49].

B. Simulation of RCR particles

We model the RCR particle as a juxtaposition of two half-disks of radius R ′ with one rectangle of length L and width 2R ′ ; see Fig. 1. The shape of a RCR particle is a circle of radius R ′ for L = 0. The aspect ratio α = (L + 2R ′ )/(2R ′ ) is 1 in this limit and increases with L for a fixed value of R ′ . In this paper, we use an alternative parameter describing the deviation of the particle shape from a circle. Let R be the radius of the circle circumscribing the particle. We have R = L/2 + R ′ . The radius R ′ is also that of the inscribed circle. Hence, the deviation from a circular shape can be characterized by ∆R = R -R ′ = L/2. We use the dimensionless parameter η defined by

η = ∆R R = α -1 α . (3)
It varies from η = 0, for a circle, to 1 corresponding to a line. We will refer to η as the elongation parameter as in rock mechanics [START_REF] Folk | Petrology of Sedimentary Rocks[END_REF].

For the detection of the contacts between two RCR particles, we use the schema shown in Fig. 2. Three types of contact can be distinguished: cap-to-cap (cc), cap-toside (cs) and side-to-side (ss). The contacts between the particles are thus detected separately for the pairs of circles and rectangles. In general, in the CD method ss contact between two rectangles is treated as composed of two point contacts and the contact laws (1) and ( 2) are applied separately to each point. The choice of these points does not affect the resultant force and its point of application. Hence, for RCR particles, as shown in Fig. 2, ss contact is composed of four point contacts : two points due to the rectangle-rectange interface and two points due to the cc contacts. Thus, four forces are calculated by the CD algorithm but only their resultant and application point are material.

The detection of line contacts between rectangles was implemented through the so-called shadow overlap method devised initially by Moreau [15,[START_REF] Dubois | Actes du sixième colloque national en calcul des structures -CSMA-AFM-LMS[END_REF] for polygons. The reliability and robustness of this method have been tested in several years of previous applications to granular materials [3,10,11,15,[START_REF] Azéma | [END_REF][START_REF] Cholet | World Congress on Railway Research[END_REF][START_REF] Azéma | [END_REF]. This detection procedure is fairly rapid and allows us to simulate large samples composed of RCR particles. For our simulations, we used the LM GC90 which is a multipurpose software developed in Montpellier, capable of modeling a collection of deformable or undeformable particles of various FIG. 2: Representation of cap-to-cap, cap-to-side and sideto-side contact and they will be referred as cc contacts, cs contacts and ss contact, respectiveley.

shapes (spherical, polyhedral, or polygonal) by means of the contact dynamics (CD) method [START_REF] Dubois | Actes du sixième colloque national en calcul des structures -CSMA-AFM-LMS[END_REF].

C. Sample preparation

We prepared 8 different packings of 13000 RCR particles with η varying from 0 to 0.7 by steps of 0.1. The radius R of the circumscribing circle defines the size of a RCR particle. In order to avoid long-range ordering in the limit of small values of η, we introduce a size polydispersity by taking R in the range [R min , R max ] with R max = 3R min with a uniform distribution in particle volume fractions.

All samples are prepared according to the same protocol. A dense packing composed of disks (η = 0) is first constructed by means of a layer-by-layer deposition model based on simple geometrical rules [53][54][55][56]. The particles are deposited sequentially on a substrate. Each new particle is placed at the lowest possible position at the free surface as a function of its diameter. This procedure leads to a random close packing in which each particle is supported by two underlying particles and supports one or two other particles. For η > 0, the same packing is used with each disk serving as the circumscribing circle of a RCR particle. The RCR particle is inscribed with the given value of η and random orientation in the disk.

Following this geometrical process, the packing is compacted by isotropic compression inside a rectangular frame of dimensions l 0 × h 0 in which the left and bottom walls are fixed, and the right and top walls are subjected to a compressive stress σ 0 . The gravity g and friction coefficients µ between particles and with the walls are set to zero during the compression in order to avoid force gradients and obtain isotropic dense packings. Fig. 3 displays snapshots of the packings for several values of η at the end of isotropic compaction.

The isotropic samples are then subjected to vertical compression by downward displacement of the top wall at a constant velocity v y for a constant confining stress σ 0 acting on the lateral walls. The friction coefficient µ between particles is set to 0.5 and to zero with the walls. The simulations were run with a time step of 2.10 -4 s. The CPU time was 5.10 -4 s per particle and per time step on an AMD processor. Since we are interested in quasistatic behavior, the shear rate should be such that the kinetic energy supplied by shearing is negligible compared to the static pressure. This can be formulated in terms of an inertia parameter I defined by [57]:

I = ε m p , (4)
where ε = ẏ/y is the strain rate, m is the particle mass and p is the mean pressure. The quasistatic limit is characterized by the condition I ≪ 1. In our simulations, I was below 10 -3 .

III. STRENGTH AND DILATANCY

In this section, we consider the stress-strain and volume-change behavior as a function of the shape parameter η. We need to evaluate the stress tensor and solid fraction during deformation from the simulation data. For the stress tensor, we start with the tensorial moment M i of each particle i that is defined by [START_REF] Moreau | Friction, Arching, Contact Dynamics[END_REF][START_REF] Staron | [END_REF]:

M i αβ = c∈i f c α r c β , (5)
where f c α is the α component of the force exerted on particle i at the contact c, r c β is the β component of the position vector of the same contact c, and the summation runs over all contact neighbors of particle i (noted briefly by c ∈ i). The average stress tensor σ in the volume V of the granular assembly is given by the sum of the tensorial moments of individual particles divided by the volume [START_REF] Moreau | Friction, Arching, Contact Dynamics[END_REF][START_REF] Staron | [END_REF]:

σ = 1 V i∈V M i = 1 V c∈V f c α ℓ c β , (6)
where ℓ c is the branch vector joining the centers of the two touching particles at the contact point c. Remark that the first summation runs over all particles whereas the second summation involves the contacts, each contact appearing only once. Under biaxial conditions with vertical compression, we have σ 1 ≥ σ 2 , where the σ α are the stress principal values. The mean stress p and stress deviator q are defined by:

p = 1 2 (σ 1 + σ 2 ), (7)
q = 1 2 (σ 1 -σ 2 ). ( 8
)
For our system of perfectly rigid particles, the stress state is characterized by the mean stress p and normalized shear stress q/p. The strain parameters are the cumulative vertical, horizontal and shear strains ε 1 , ε 2 and ε q , respectively. By definition, we have

ε 1 = h h0 dh ′ h ′ = ln 1 + ∆h h 0 , ( 9
)
where h 0 is the initial height and ∆h = h 0 -h is the total downward displacement, and

ε 2 = l l0 dl ′ l ′ = ln 1 + ∆l l 0 , (10)
where l 0 is the initial box width and ∆l = l -l 0 is the total change of the box width. The cumulative shear strain is then defined by

ε q ≡ ε 1 -ε 2 . (11)
Finally, the cumulative volumetric strain ε p is given by

ε p = ε 1 + ε 2 = V V0 dV ′ V ′ = ln 1 + ∆ν ν (12)
where V 0 = l 0 h 0 is the initial volume and ∆ν = ν -ν 0 is the cumulative change of solid fraction.

Figure 4 shows the normalized shear stress q/p as a function of shear strain ε q for different values of η. The jump observed at ε q = 0 reflects both the rigidity of the particles and high initial solid fraction of the samples (see below). In all cases, the shear stress passes by a peak before relaxing to a stress plateau corresponding to the so-called "residual state" in soil mechanics [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF]. We remark that the residual shear stress increases with η.

The internal angle of friction ϕ * , representing the shear strength of the material, is defined from the mean value (q/p) * of the normalized shear stress in the residual state by [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF] sin

ϕ * = q p * . ( 13
)
Fig. 4 shows the variation of sin ϕ * as a function of α and η. We see that the shear strength is an increasing nonlinear function of the aspect ratio, but, interestingly, it varies linearly when plotted versus the elongation parameter. Hence, we have

sin ϕ * = sin ϕ * 0 + k η = sin ϕ * 0 + k 1 - 1 α (14)
This observation indicates that the evolution of shear strength reflects more directly shape elongation than aspect ratio. In the following, we will use η as shape parameter. Figure 6 (a) displays the cumulative volumetric strain ε p as a function of ε q for different values of η. Starting with an initially dense state, all packings dilate and hence the volume increases. For η ≤ 0.4, a plateau is reached beyond ε q = 0.3, corresponding to a state of isochoric deformation. For larger aspect ratios, the dilatation continues even at very large deformations. This is an indication of an inhomogeneous dilation due to the formation of shear bands in the bulk, which is enhanced by the elongated shape of the particles. Since different parts of the packing undergo differential volume change, longer shearing is required to reach a fully dilated state for the whole packing. The initiation and evolution of shear bands for different particle elongations will be reported in more detail elsewhere.

Figure 6 (b) displays the solid fraction ν as a function of η at different levels of shear deformation ε q . It is remarkable that, at all levels of deformation, the solid fraction increases with η, reaches a maximum at η ≃ 0.4 and then declines as η further increases. We note that solid fractions as large as 0.90 are reached for η = 0.4 in the initial state. A similar nonmonotonous behavior was observed for packings of ellipses or ellipsoidal particles [19,20,22]. This is somewhat a counterintuitive finding as the shear strength (a monotonous function of η) does not follow the trend of solid fraction (nonmonotonous). This behavior is clearly not related to shear localization since it is observed at all levels of deformation including the initial isotropic state ν 0 = ν(ε q = 0) where the packings are homogeneous.

A rapid fall-off of solid fraction for elongated parti-

ν ε q = 0 ε q = 0.1 ε q = 0.2 ε q = 0.3 ε q = 0.4 (b)
FIG. 6: Cumulative volumetric strain εp as a function of shear strain εq (a); Solid fraction as a function of particle shape parameter η (b) at different levels of shear strain.

cles in 3D was observed at very large aspect ratios and attributed to the excluded volume due to disorder, predicting a fall-off in inverse proportion to the aspect ratio [37,[START_REF] Blouwolff | [END_REF]62]. The initial rapid increase of solid fraction, as observed in Fig. 6, reveals that excluded-volume effects are not the prevailing mechanism at low aspect ratios. In this limit, slight deviations from spherical shape have strong space-filling effect on the packing although the excluded volume increases at the same time and becomes dominant at very large aspect ratios.

The volumetric deformation can also be expressed in terms of the dilatancy angle ψ defined by [START_REF] Wood | Soil behaviour and critical state soil mechanics[END_REF]:

sin ψ = ∂ε p ∂ε q . ( 15
)
The plot of ψ as a function of ϕ, the so-called "stressdilatancy diagram", is shown in Fig. 7 for different values of η. We observe a linear correlation between ϕ and ψ irrespective of the value of η. We have

ϕ ≃ ϕ * + ψ. ( 16
)
This is a particularly simple relation compared to several models proposed in soil mechanics [START_REF] Wood | Soil behaviour and critical state soil mechanics[END_REF][START_REF] Radjai | The Physics of Granular Media[END_REF]. It reflects the "non-associated" character of granular plasticity, an associated behavior implying simply ϕ = ψ, which is unrealistic for granular materials [START_REF] Wood | Soil behaviour and critical state soil mechanics[END_REF][START_REF] Taylor | Fundamentals of soils mechanics[END_REF][START_REF] Vermeer | Constitutive relation for soils[END_REF][START_REF] Vermeer | [END_REF]. According to relation (16), the dilatancy angle vanishes in the residual state where ϕ = ϕ * . Recent work on cohesive and granular packings of polygonal particles in 2D is in agreement with this correlation [START_REF] Azéma | [END_REF]68]. In the following, we focus on the microstructure and force transmission that provide a key to a better understanding of the physical mechanisms underlying the effect of particle shape on the shear strength.

IV. GRANULAR TEXTURES

In this section, we investigate the general organization (texture) of our packings of RCR particles in terms of particle orientations and contact network. This will allow us to quantify the effect of the elongation parameter and its connection with shear strength.

A. Particle orientations

The principal feature of elongated particles is their orientational degree of freedom. The particle orientation is represented by a unit vector m as shown in the inset to Fig. 9. In 2D, it is parametrized by a single angle ϑ. Let D(ϑ) be the set of particles with direction ϑ ∈ [ϑ -δϑ/2; ϑ + δϑ/2] for angle increments δϑ, and N p (ϑ) its cardinal. The probability P ϑ (ϑ) of the orientations of particles is given by

P ϑ (ϑ) = N p (ϑ) N p , (17)
where N p is the total number of particles. Figure 8 displays a polar representation of P ϑ (ϑ) for η = 0.7 at various stages of deformation. In the initial state, corresponding to an isotropic stress state, the distribution is anisotropic with privileged direction ϑ p close to π/2. This means that for elongated particles, the particle orientations are not fully correlated with the stress state so that the resulting particle orientation anisotropy depends on details of the assembling process that can not be controlled by simply subjecting the particles to isotropic stresses from the boundary.

The priviliged direction rotates as a result of vertical compression and becomes horizontal (parallel to the minor principal stress direction) in the residual state. The

ε q = 0 ε q = 0.1 ε q = 0.2 ε q = 0.3 ε q = 0.4
FIG. 8: Polar representation of the probability density P ϑ of particle orientations ϑ for η = 0.7 at various stages of deformation εq. The symbols are the simulations data. Solid lines represent harmonic fit according to equation (18).

distribution are nicely fitted by harmonic approximation corresponding to the lowest order terms of the Fourier expansion of P ϑ (ϑ) [START_REF] Rothenburg | Powder and Grains[END_REF][START_REF] Reddy | [END_REF]:

P ϑ (ϑ) = 1 2π (1 + a p cos(2(ϑ -θ σ )), (18)
where a p represents the anisotropy of the distribution and θ σ is the major principal stress direction. The choice of θ σ as reference direction is motivated by the observation that the privileged orientation of the particles tends to align itself with the minor principal stress direction.

Hence, negative values of a p mean that the particles are preferentially oriented perpendicular to the major principal stress. We plot the particle orientation anisotropy a p in Fig. 9 as a function of η at different stages of shear ε q . We see that the particle orientations are isotropic for η ≤ 0.4 at the initial state and they become increasingly anisotropic as η increases beyond 0.4. At nearly all stages of shear, a p is negative, and at most advanced stage, i.e. corresponding approximately to the residual state, it becomes nearly independent of η. The large absolute value of a p (∼ 0.35) suggests that many particles are aligned in horizontal layers just as in nematic order. One example of this nematic ordering is shown in Fig. 10 for two different values of η. This ordering may be attributed to the favored mechanical equilibrium of the particles under the action of vertical stress and enhanced by boundary alignment of the elongated particles [3,25,26,[START_REF] Zuriguel | Proc. R. Soc. A[END_REF]. This point will be analyzed more deeply below.

B. Particle connectivity

The primary statistical quantity describing the contact network is the coordination number z (average number -0,1 0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 -0,4 -0,2 0,0 0,2 0,4 a p q = 0 q = 0.1 q = 0.2 q = 0.3 q = 0.4 m FIG. 9: Particle orientation anisotropy ap as a function of shape parameter η at different stages of shear εq.

of contacts per particle). For our elongated particles, each side-to-side contact is counted as one contact even if side-to-side contacts are treated as four point contacts belonging to the contact segment (see section II). The floating particles with no force-bearing contact (i.e. with less than two active contacts) are thus removed from the statistics. The fraction of floating particles decreases linearly with η from 17% for η = 0 to 10% for η = 0.7.

Figure 11 displays the evolution of z as a function of η, in the initial and residual states. The initial-state value z = 4 corresponds to a frictionless packing of circular particles in the isostatic state with z = 2d where d is space dimension (indeed, the packings where prepared by setting the friction coefficient to zero). However, as η increases, z increases to a plateau value of ∼ 5.6. This is in agreement with recent work showing that large disordered jammed packings are isostatic only for disks or spheres [18][19][20][21][22][23][START_REF] Roux | [END_REF]. For nonspherical or noncircular particles, we have z d(d + 1).

Numerical results for frictional or frictionless systems of rigid or deformable disks and spheres [18,[73][74][75][76][77], as well as with ellipses and spheroids [19][20][21][22][23] confirm this point. In the residual state, the mean value of z is below that in the isotropic state, and it grows from 3 to 5 with η. It is interesting to note that z does not follow the solid fraction which, as we have seen before, is nonmonotonous as a function of η. This means that for large aspect ratios, the packings are loose but well connected.

The connectivity of the contact network can be characterized in more detail by the fraction P (c) of particles with exactly c contact neighbors. Fig. 12 shows P (c) in the residual state for different values of η. The distribution is increasingly broader as η becomes larger. The particles can have as many as 10 contact neighbors at η = 0.7. This is allowed both by the geometry and polydispersity of the particles as shown by a typical grey-level map of particle connectivities in Fig. 13. For η = 0, we observe a peak at c = 3. This peak slides gradually to c = 5 at η = 0.7 as observed also for z, which is, by definition, the mean value of c for force-bearing contacts:

z = c = ∞ c=2 c P (c).

C. Force and texture anisotropies

Equation (6) shows that the expression of stress tensor is an arithmetic mean involving the branch vectors ℓ and contact force vectors f . This means that for the analysis of stress transmission and shear strength from a particlescale viewpoint we need a statistical description of these quantities.

A common approach used by various authors is to express branch vectors and contact force orientations in terms of the contact direction, i.e. in the local contact frame (n, t), where n is the unit vector perpendicular to the contact plane, and t is an orthonormal unit vector oriented along the tangential force [4, 5, 24, 28-30, 32, 35, 78-80] ; see figure 14(a). The components of the branch vector and contact force are expressed in this frame as:

ℓ = ℓ n n + ℓ t t, f = f n n + f t t, (19)
where ℓ n and ℓ t are the normal and tangential compo- Remark that only for disks or spherical particles we have ℓ = ℓn where ℓ is the length of the branch vector. A consequence of noncircular or nonspherical particle shape is to dissociate the contact frame from the branch vector frame (n ′ , t ′ ), where n ′ is the unit vector along the branch ℓ and t ′ is the orthoradial unit vector [11,81] ; see figure 14(b). We express the components of the branch vector and contact force also in this frame:

ℓ = ℓ n ′ n ′ , f = f n ′ n ′ + f t ′ t ′ , (20)
where f n ′ and f t ′ are the radial and orthoradial components of the contact force, and ℓ n ′ = ℓ.

In two dimensions, let θ and θ ′ be the orientations of of n and n ′ , respectively. From the numerical data, we can evaluate the probability density functions P θ (θ) and P θ ′ (θ ′ ) of contact and branch vector orientations, respectively, as well as the angular averages of the force components f n (θ), f t (θ), f n ′ (θ ′ ), and f t ′ (θ ′ ) and branch vector components ℓ n (θ), ℓ t (θ), ℓ n ′ (θ ′ ). In the absence of an intrinsic polarity for n and n ′ , all these functions are π-periodic. The insets to Figs. 15 17 display polar representations of these functions for η = 0.7 at the end of shearing. All these angular functions are generically anisotropic. The peak values occur along the axis of compression (θ = π/2) for P θ , P θ ′ , f n and f n ′ , and along the axis of extension (θ = 0) for ℓ n and ℓ n ′ . The maxima for the tangential components occur in the direction of π/4 with respect to the axis of compression.

The simple shapes of the above functions suggest that they can be approximated by their Fourier expansions up to the second term [11,24,30,32,[START_REF] Radjai | Granular Materials: Fundamentals and Applications[END_REF]:

         P Θ (Θ) = 1 2π {1 + a * c cos 2(Θ -Θ * c )}, ℓ n * (Θ) = ℓ {1 + a ln * cos 2(Θ -Θ ln * )}, ℓ t * (Θ) = ℓ a lt * sin 2(Θ -Θ lt * ), f n * (Θ) = f {1 + a f n * cos 2(Θ -Θ f n * )} f t * (Θ) = f a f t * sin 2(Θ -Θ f t * ), (21)
where Θ stands either for θ or for θ ′ depending on the local frame used. The ℓ is mean length of branch vectors and f is the mean force. (a * c , a ln * , a lt * , a f n * , a f t * ) = (a c , a ln , a lt , a f n , a f t ) and (Θ * c , Θ ln * , Θ lt * , Θ f n * , Θ f t * ) = (θ c , θ ln , θ lt , θ f n , θ f t ) are the anisotropy parameters and the angle of privileged direction of each function in the frame (n, t).

In the same way, we have

(a * c , a ln * , a lt * , a f n * , a f t * ) = (a ′ c , a ln ′ , a lt ′ , a f n ′ , a f t ′ ) and (Θ * c , Θ ln * , Θ lt * , Θ f n * , Θ f t * ) = (θ ′ c , θ ln ′ , θ lt ′ , θ f n ′ , θ f t ′ ) in the frame (n ′ , t ′ ).
Note that, by construction, we have a lt ′ = 0 and θ lt ′ = 0. In the following, we will refer to a c as contact anisotropy, to a ′ c as branch vector orientation anisotropy, to (a ln * , a lt * ) as branch length anisotropies and to (a f n * , a f t * ) as normal and tangential or radial and orthoradial force anisotropies depending on the local frame [9,[START_REF] Azéma | [END_REF]. These harmonic approximations are well fit to our data as shown in the insets to Figs. 15, 16 and 17. In practice, it is convenient to estimate the above anisotropies from the following fabric and force tensors [31]:

                               F * αβ = 1 π π 0 n * α n * β P Θ (Θ)dΘ, χ ln * αβ = 1 ℓ π 0 ℓ n * (Θ)n * α n * β P Θ (Θ)dΘ, χ lt * αβ = 1 ℓ π 0 ℓ t * (Θ)n * α t * β P Θ (Θ)dΘ, χ f n * αβ = 1 f π 0 f n * (Θ)n * α n * β P Θ (Θ)dΘ, χ f t * αβ = 1 f π 0 f t * (Θ)n * α t * β P Θ (Θ)dΘ, (22)
where α and β design the components in the considered frame. Note that, by construction, we have χ lt ′ αβ = 0. From equations ( 21) and (22), assuming that in a sheared state Θ * c = Θ f n * = Θ f t * = Θ σ , Θ ln * = Θ lt * = 0 or θ σ , the following relations are easily obtained:

           a * c = 2(F * 1 -F * 2 )/(F * 1 + F * 2 ), a ln * = 2(χ ln * 1 -χ ln * 2 )/(χ ln * 1 + χ ln * 2 ) -a * c , a lt * = 2(χ l * 1 -χ l * 2 )/(χ l * 1 + χ l * 2 ) -a * c -a ln * , a f n * = 2(χ f n * 1 -χ f n * 2 )/(χ f n * 1 + χ f n * 2 ) -a * c , a f t * = 2(χ f * 1 -χ f * 2 )/(χ f * 1 + χ f * 2 ) -a * c -a f n * , (23)
where χ l * = χ ln * + χ lt * , χ f * = χ f n * + χ f t * and the indices 1 and 2 refer to the principal values of each tensor. By construction, we have (

F * 1 +F * 2 ) = 1, (χ l * 1 +χ l * 2 ) = ℓ and (χ f * 1 + χ f * 2 ) = f
. Note that a * c , a f n * and a f t * are always positive whereas a ln * and a lt * are negative. We have Θ ln * = 0 and Θ lt * = 0,.

Figure 15 displays the variation of contact anisotropy a c and branch vector orientation anisotropy a ′ c , both averaged in the residual state, as a function of η. We observe two distinct behaviors: a c increases quickly from 0.3 to 0.7 with η whereas, after a slight increase, a ′ c declines to nearly 0 for η = 0.7. It is often admitted that the shear strength in granular media is a consequence of the buildup of an anisotropic geometrical structure due to friction and steric exclusions between particles [4,5,[START_REF] Satake | Proceedings of the IUTAM symposium on deformation and failure of granular materials[END_REF][START_REF] Radjai | Granular Materials: Fundamentals and Applications[END_REF][START_REF] Oda | [END_REF][START_REF] Cambou | Powders and Grains 93[END_REF]. But, here we have two different structural anisotropies a c and a ′ c that vary in opposite directions as η is increases. Hence, when the granular structure is complex as in our packings of nonspherical particle c as a function of shape parameter η averaged in residual state. The error bars represent the standard deviation in the residual state. The inset shows the angular probability densities P θ (θ) in black and P θ ′ (θ ′ ) in red for η = 0.7 calculated from the simulation data (points) together with the harmonic approximation (lines).

-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

η -0.3 -0.2 -0.1 0.0 0.1
a ln a lt a ln' FIG. 16: Normal and tangential branch length anisotropies a ln and a lt and branch length anisotropy a ln ′ as a function of shape parameter η in the residual state. The error bars represent the standard deviation in the residual state. The inset shows the angular average functions ℓ ln (θ), ℓ lt (θ) and ℓ ln ′ (θ) in black, red and green, respectively, for η = 0.7 calculated from the simulation data (points) and approximated by harmonic fits (lines).

shapes, the choice of the statistical representation of the granular structure has to be specified [11,[START_REF] Azéma | [END_REF]. This point will be addressed in more detail in section VII.

The branch vector length anisotropies a ln , a lt and a ln ′ , averaged in the residual state, are plotted in Fig. 16 as a function of η. These parameters are negligibly small at small values of η, i.e. for nearly circular particles, and decline to negative values as η is increased. This means that the particles tend to form longer branch vectors with their neighbors in the direction of extension, suggesting that the particles touch preferentially along their minor axes when the contact orientation is close to the compression axis, and along their major axis when the contact orientation is close to the extension axis; see section VII. It is also remarkable that a ln ≃ a ln particularly for η 0.4.

The normal and tangential force anisotropies a f n and a f t are plotted in Fig. 17 as a function of η, together with the radial and orthoradial force anisotropies a f n ′ and a f t ′ , averaged in the residual state. In contrast to contact anisotropy, we see that a f n and a f n ′ grow together slowly until η = 0.4, then a f n remains nearly constant whereas a f n ′ increases. On the other hand, the anisotropy a f t ′ of orthoradial forces grows much faster with η than the anisotropy a f t of tangential forces. Remarkably, from η > 0.4 the orthoradial force anisotropy is higher than radial force anisotropy (a f t ′ > a f n ′ ) whereas, even if the tangential force anisotropy increases with η, it is still below the normal force anisotropy (a f t < a f n ) and remains always below the contact anisotropy (a f n < a c ). In other words, the force anisotropy described in terms of branch vectors reflects more sensitively the effect of particle shape elongation than in terms of contact normal vectors. We will see below that this behavior is related to the mobilization of friction (section V) and contact types (section VII).

V. GEOMETRICAL AND MECHANICAL ORIGINS OF SHEAR STRENGTH

The stress tensor as formulated in Eq.( 6) is a function of discrete microscopic parameters attached to the contact network. For sufficiently large systems, the dependence of volume averages on individual discrete parameters vanishes [24,30,[START_REF] Landau | Statistical physics : course of theoretical physics[END_REF] and the discrete sums can be replaced by integrals. According to Eq. ( 6), we have

V σ αβ = c∈i f c α r c β = N c f α ℓ β , (24)
where N c is the total number of contacts. By writing the average on the right hand side in integral form, we get

σ αβ = n c Ω f α ℓ β P ℓf df dℓ, (25)
where P ℓf is the joint probability density of forces and branch vectors, n c is the number density of contacts and Ω is the integration domain in the space (ℓ, f ). The integral appearing in equation ( 25) can be reduced by integrating first with respect to the forces and branch vector lengths. Considering the components of the forces and branch vectors in one of the two local frames (n, t) or (n ′ , t ′ ) and neglecting the branch vector-force correlations, we get [9,11,30,[START_REF] Radjaï | 6th International Conference on the Micromechanics of Granular Media[END_REF]:

σ αβ = n c π 0 { ℓ n * (Θ) n * α (Θ) + ℓ t * (Θ) t * β (Θ)} { f n * (Θ) n * α (Θ) + f t * (Θ) t * β (Θ)}P (Θ) dΘ. ( 26
)
The expression of the stress tensor by equation ( 26) makes appear explicitly the average angular functions representing the fabric and force states. Using the harmonic approximation (21), equation ( 26) can be integrated with respect to space direction Θ and the stress invariants p and q extracted. Assuming that the stress tensor is coaxial with the fabric and force tensors (22), we get the following simple relations:

q p ≃    1 2 (a c + a ln + a lt + a f n + a f t ) in (n, t) 1 2 (a ′ c + a ln ′ + a f n ′ + a f t ′ ) in (n ′ , t ′ ). ( 27
)
The assumption of coaxiality is natural since, even if the preferential orientations of the forces and branch vectors are not fully correlated, we observe that shearing tends to align the contacts and forces with the principal directions of the stress tensor. Figure 18 displays the residual-state value of the normalized shear stress (q/p) * as a function of η calculated both directly from the simulation data and from equation (27) separately for the two local frames by using the values of various anisotropies estimated from the simulation data. As we see, for both local frames, equation ( 27) provides an excellent fit to the data for all values of η. Note, however, that the second expression in equation ( 27) is more simple than the first expression (4 anisotropic parameters vs 5 anisotropic parameters) and the resulting fit appears to be more accurate.

The two equations ( 27) are interesting as they reveal distinct origins of shear strength in terms of force and texture anisotropies with two different decompositions. Various anisotropies do not contribute equally to the shear strength. The dominant anisotropies are texture anisotropies for projection on contact frame since a c + a ln + a lt > a ′ c + a ln ′ and force anisotropies for projection on the branch vector frame since a f n ′ + a f t ′ > -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 η 0.1 0.2 0.3 0.4 0.5 0.6 0.7 (q/p) * 0.5(a c +a ln +a lt +a fn +a ft ) 0.5(a' c +a l' +a fn' +a ft' ) FIG. 18: Normalized shear stress (q/p) * in the residual state as a function of η together with two analytical expressions given by equation (27). The error bars represent the standard deviation in the residual state.

a f n + a f t . The fact that the texture anisotropy prevails in the contact frame may be attributed to its strong correlation with particle orientations. Geometrically, for a particle oriented along a direction ϑ, more contacts can be formed with the flat side of the particle with normals oriented along ϑ + π/2 than with its rounded caps. This is consistent with the observation that the particle orientations are strongly anisotropic with an anisotropy a p of negative sign (preferred direction along the extension axis) and the contact normal anisotropy a c is positive (along the compression axis) and increases with aspect ratio; see section IV.

In Fig. 19 two maps of radial forces are shown for packings with η = 0.2 and η = 0.7, respectively. In the presence of long parallel sides, the strong force chains are more tortuous. Hence, the stability of such structures requires strong activation of tangential forces. Indeed, we remark that the orthoradial force anisotropy is above the radial force anisotropy (a f t ′ > a f n ′ ) for the most elongated particles in contrast to the tangential force anisotropy which is below the normal force anisotropy (a f t < a f n ). As a result of the increasing activation of tangential forces, the fraction of sliding contacts (i.e. contacts where |f t | = µ|f n |) grows with η as shown in Fig. 20. The contributions of side-to-side and cap-to-side contacts to force anisotropy and friction mobilization, which are major effects of particle shape, will be analyzed in section VII.

VI. FORCE DISTRIBUTIONS

The force chains and spatially inhomogeneous stress distributions are well-known features of granular media. A well-known observation is that a large number of contacts transmit very weak forces, a signature of the arching effect, whereas a smaller fraction of contacts carry strong force chains [31]. Force transmission has been investigated by experiments and numer- ical simulations for disks, ellipses and polygonal particles in 2D as for spherical, cylindrical and polyhedral particles in 3D [11,25,[START_REF] Azéma | [END_REF][START_REF] Liu | [END_REF][87][88][89][90][91][92][93][94]. In close correlation with shear strength and solid fraction, the stress transmission is strongly influenced by particle shape. In particular, one expects that elongated particle shapes will influence mainly the distribution of weak forces by enhancing the arching effect.

The probability density function (pdf) of radial forces normalized by the mean radial force f n ′ is shown in Fig. 21 in log-linear and log-log scales at large strains (the data are cumulated from several snapshots in the residual state) for all values of η. As usually observed, in all packings the number of forces above the mean f n ′ falls off exponentially whereas the number of forces below the mean vary as a power-law:

P (f n ′ ) ∝ e -α n ′ (η)(1-f n ′ / f n ′ ) , f n ′ > f n ′ , f n ′ f n ′ β n ′ (η) , f n ′ < f n ′ , (28)
where α n ′ (η) and β n ′ (η) are the exponents which decrease with η from α n ′ (0) ≃ 1.69 to α n ′ (0.7) ≃ 0.88, and from β n ′ (0) ≃ 0.13 to β n ′ (0.7) ≃ -0.49. Figure 22 shows the pdf P (f t ′ ) of orthoradial forces normalized by by the mean orthoradial force f t ′ in each packing. These distributions are also characterized by an exponential falloff for the forces above the average force f t ′ and a power law for the forces below f t ′ :

P (f t ′ ) ∝ e -α t ′ (η)(1-|f t ′ |/ |f t ′ | ) , |f t ′ | > |f t ′ | , |f t ′ | |f t ′ | β t ′ (η) , |f t ′ | < |f t ′ | , (29)
with the corresponding exponents α t ′ (η) and β t ′ (η), which decrease from α t ′ (0) ≃ 1.09 to α t ′ (0.7) ≃ 0.73, and from β t ′ (0) ≃ -0.37 to β t ′ (0.7) ≃ -0.73. These distributions show clearly the larger inhomogeneity of stress transmission in a granular packing composed of elongated particles. We find that (the results not shown here), as with circular particles, the contacts can be classified into strong and weak networks. Evaluating q/p separately for each network, it is found that the shear stress is almost totally sustained by the strong contact network. Maps of strong and weak radial networks for radial forces are displayed in Fig. 19 for η = 0.2 and η = 0.7. The fraction of floating particles (less than two contacts) decreases linearly with η from 17% for η = 0 to 10% for η = 0.7. Hence, more particles are involved in the contact network for more elongated particles, but it is remarkable that the proportion of weak forces grows from 60% for η = 0 to 70% for η = 0.7. In other words, although the number of strong contacts decreases with η, stronger force chains occur with more elongated particles. Although we focused here on the networks of radial force components, we basically obtain the same conclusions for the normal force components.

We also remark that the packings are increasingly more inhomogenous in the sense that as particle elongation increases, the packing involves less strong force chains in number but with stronger forces. This decreasing force homogeneity in spite of increasing connectivity (fig. 11), means that force distributions are controlled by more subtle details of the microstructure than the density of contacts or solid fraction. As we shall see below, this is related to the role of various contact types in the contact network.

VII. EFFECT OF CONTACT TYPES

Remembering that RCR particles are clumps of two disks with one rectangle, in this section we revisit the results of the previous sections in the light of the organization of cap-to-cap (cc), cap-to-side (cs) and side-to-side (ss) contacts in each packing. The side-to-side or sideto-cap contacts do not transmit torques, but they are able to accommodate force lines that are usually unsustainable by cap-to-cap contacts. For this reason, it is worth while trying to isolate their respective roles with respect to the shear strength. The proportions of these contact types and their contributions to the structural anisotropy and force transmission are key quantities for understanding the effect of particle shape on the shear strength properties of granular media [11,[START_REF] Azéma | [END_REF].

In the residual state, the proportions of different contact types are nearly constant. Fig. 23 shows the proportions k cc , k cs and k ss of cc, cs and ss contacts averaged over the residual state as a function of η. We see that k cc declines with η from 1 (for disks) to 0.2 for η = 0.7. At the same time, k cs and k ss increase from 0 to 0.6 and to 0.2, respectively. Interestingly, k cs ≃ k cc for η ≃ 0.4, and k ss ≃ k cs for η = 0.7. In this way, as the particle elongation increases, the packing passes from a contact network dominated by cc contacts to a contact network dominated by the complex contacts cs and ss.

To identify the impact of each contact type on the shear strength, we proceed by additive decomposition of the stress tensor by considering the expression (25) of the stress tensor and grouping the contacts according their types:

σ = σ cc + σ cs + σ ss , (30)
where σ cc , σ cs and σ ss are obtained from the expression of the stress tensor Eq. ( 25) by restricting the summation to cc, cs and ss contacts, respectively. The corresponding stress deviators q cc , q cs and q ss are then calculated and normalized by the mean pressure p. Fig. 24 shows (q/p)* q cs /p q cc /p q ss /p FIG. 24: Shear strength (q/p) * for cs, ss and cc contacts as a function of η, together with the harmonic approximation fits in (n, t) frame (---) and (n ′ , t ′ ) frame (...).

q cc /p, q cs /p and q ss /p averaged in the residual state as a function of η. We see clearly that q cc /p follows a trend opposite to that of q cs /p. For η < 0.3, (q/p) * is dominated by cc contacts. For η ≃ 0.3, cc and cs contacts participate equally to the shear stress, and for η > 0.3, the cs contacts dominate (q/p) * . Remarkably, q ss /p ≃ 0 for η < 0.4. As we shall see below, the ss contacts participate to the strong force chains only in the case of the most elongated particles. In this way, the growth of the number of cs and ss contacts shown in Fig. 23 is clearly at the origin of a gradual consolidation of the packings as η increases.

In order to get further insight into the organization of different contact types, it is useful to consider partial connectivities P cs (c), P cc (c) and P ss (c) defined as the fraction of particles with exactly c contacts of cs type, cc type and ss type. These functions are displayed in Fig. 25 for all our packings. Note that, by definition P cc (c) ≡ P (c) for η = 0. We see that P cc gets narrower as η increases whereas P cs and P ss get broader. It should be noted that, even for the most elongated particles, a particle can have at most two ss. This means that, the elongated particles tend mainly to pile up like bricks. On the other hand, the peak of P cs slides to the larger values as η increases. For η = 0.7, most particles have three or four cs contacts (for nearly 40%).

We now consider the anisotropy of the branch vectors and contact forces supported by the three contact types at the contact and branch frames. Following the same procedure as for the stress tensor (see equation (30), we perform an additive decomposition of the fabric and force tensors:

       F * = F * cc + F * cs + F * ss , χ ln * = χ ln * cc + χ ln * cs + χ ln * ss , χ f n * = χ f n * cc + χ f n * cs + χ f n * ss , χ f t * = χ f t * cc + χ f t * cs + χ f t * ss , (31)
where the indices refer to the partial contributions of cc, cs and ss contacts. The corresponding anisotropies of each tensor can be extracted. In principle, the principal directions of these partial tensors do not coincide with those of the overall tensors at all stages of shearing. But, in practice, in the residual state, the principal directions coincide so that the global anisotropy of each tensor is the sum of its partial anisotropies:

q γ p ≃    1 2 (a cγ + a lnγ + a ltγ + a f nγ + a f tγ ) in (n, t) 1 2 (a ′ cγ + a ln ′ γ + a f n ′ γ + a f t ′ γ ) in (n ′ , t ′ ), (32)
where γ stands alternatively for {cc, cs, ss}. This decomposition is nicely verified by our numerical date as shown in Fig. 24.

Since the contact orientation anisotropy expressed in (n, t) frame and the force anisotropy expressed in (n ′ , t ′ ) frame provide respectively fine descriptions of the geometrical and force organizations (see section V), we restrict here our analysis to the contribution of various contact types to the contact orientation anisotropy a c and the radial force anisotropies a f n ′ and a f t ′ . Figure 26 shows the variation of the partial contact anisotropies a ccc , a ccs and a css due to cc, cs and ss contacts in the residual state as the function of η. The anisotropy a css of ss contacts increases slowly with η from 0 to 0.18. At the same time, a ccc decreases and at η = 0.7 we have a css = a ccc . Hence, although the ss contacts represent at η = 0.7 nearly 20% of contacts, their contribution to the contact anisotropy remains modest and of the same order as cc contacts. The variation of the contact anisotropy a c is thus largely governed by that of a ccs .

Figure 27 shows the partial radial force anisotropies a f n ′ cc , a f n ′ cs and a f n ′ ss , as well as the partial orthoradial force anisotropies a f t ′ cc , a f t ′ cs and a f t ′ ss in the residual state as the function of η. As for contact anisotropies, the cs contacts carry most of the radial and orthoradial force anisotropies. The ss contacts contribute modestly to the global force anisotropies only for η ≥ 0.4. The anisotropy a f n ′ cc declines with η, mainly due to their low number, and a f t ′ cc stays nearly constant.

A map of contact forces projected along the branch vectors is displayed in Fig. 28 in different colors according to the type of contact. For η = 0.7, we see that the network of very strong zigzag force chains is composed mostly of cs and ss-contacts and occasionally mediated by cc contacts. In contrast, for η = 0.2, the cc contacts appear clearly to be correlated in the form of long chains across the packing rarely mediated by cs-contacts. In all cases, the strong force chains are mostly parallel to the direction of compression.

In order to recognize quantitatively the roles of cc, cs and ss contacts with respect to the force network, we plot in Fig. 29 their respective proportions k cc , k cs and k ss alternatively for the strong and weak networks in the residual state as as a function η. Notice that the data are normalized for each network. The proportion of cc contacts declines rapidly in both networks as η increases whereas that of cs and ss contacts grow. We also remark that the cs contacts are slightly more numerous in the weak network than in the strong network. The proportions have nearly the same value in the two networks for cc and ss contacts.

VIII. CONCLUSIONS

In this paper, we investigated the effect of particle elongation on the quasistatic behavior of sheared granular materials by means of Contact Dynamics simulations. The particle shapes are rounded-cap rectangles (RCR) characterized by their elongation η defined as deviation from a reference circular shape, or alternatively by their aspect ratio. As the elongation increases from 0 to 1, the particle shape varies continuously from a disk to an increasingly thin rectangle with rounded caps. The macroscopic and microstructural properties of several packings of 13000 particles, subjected to biaxial compression, were analyzed as a function of η.

An interesting finding is that the shear strength is an increasing linear function of elongation, suggesting that the parameter η is a "good" shape parameter for our 2D granular packings. In order to understand the microscopic origins of this behavior, we performed an additive decomposition of the stress tensor based on a harmonic approximation of the angular variation of average local branch vectors, contact normals and forces. This approximation of the shear strength in terms of texture and force anisotropies turns out to be in excellent agreement with our numerical data in the investigated range of the elongation parameter (η ∈ [0, 0.7]). Given the evolution of various anisotropies with particle elongation, we find that both force and texture anisotropies contribute to the increase of shear strength, but the increasing mobilization of friction force and the associated anisotropy seem to be the key effect of particle elongation. In particular the proportion of sliding contacts increases strongly as the particles become more elongated. This effect is correlated with a local nematic ordering of the particles which tend to be oriented perpendicular to major principal stress direction. This ordering is enhanced beyond η = 0.4 but remains essentially of local nature. In this respect, the fraction of side-to-side contacts increases at large particle elongations. However, the force transmission is found to be mainly guided by cap-to-side contacts, which represent the largest fraction of contacts for the most elongated particles and carry a large part of the shear strength.

In contrast to shear strength, the solid fraction is not a monotonous function of particle elongation; It first increases with particle elongation, then declines as the particles become more elongated. In other words, small deviation from circular shape favors the space-filling aptitude of the particles but beyond a characteristic elongation the excluded-volume effects prevail and lead to increasingly larger pores which cannot be filled by the particles. It is remarkable that the coordination number does not follow the solid fraction but increases with particle elongation, so that the packings of the most elongated particles are loose but well connected.

Some features discussed in this paper can legitimately be attributed to the two-dimensional geometry of the particles. For example, rounded-cap-cylinders (spherocylinders), as three-dimensional analogs of RCR particles in 2D, do not undergo spontaneously a nematic ordering. However, we expect that 3D packings of roundedcap-cylinders behave in many ways as our 2D packings with increasing particle elongation. In particular, the excluded volume effect is reinforced by particle elongation and it leads to a similar nonmonotonous dependance on the elongation as in 2D [19,20,22]. In any case, it would be highly instructive to investigate 3D packings of elongated particles, along the same lines as in this paper. It is also obvious that more work required to assess the proper role of friction in 2D for RCR particles since friction mobilization seems to underlie to a large extent the shear strength, we are performing presently more simulations with lower values of friction coefficient. In the same way, we consider the effect of cohesion between particles with respect to the shear strength of packings of elongated particles.

We specially thank I. Zuriguel and F. Dubois for fruitful discussions.

This work was done as part of PPF CEGEO research project (www.granuloscience.org/CEGEO).

R∆RFIG. 1 :

 1 FIG.1: Shape of a Rounded-Cap Rectangle (RCR).

7 FIG. 3 :

 73 FIG. 3: Examples of the generated packings at the initial state.

7 FIG. 4 :FIG. 5 :

 745 FIG.4: Normalized shear stress q/p as a function of cumulative shear strain εq for different values of the shape parameter η.

FIG. 7 :

 7 FIG. 7: Dilatancy angle ψ versus internal angle of friction ϕ for different values of η.

FIG. 10 :

 10 FIG. 10: (Color online) Color level map of particle orientations for η = 0.2 (up) and η = 0.7 (down) in the residual state.

FIG. 11 : 7 FIG. 12 :

 11712 FIG. 11: Initial and residual coordination numbers as a function of shape parameter η. The error bars represent the standard deviation in the residual state.

FIG. 13 :FIG. 14 :

 1314 FIG. 13: (Color online) Color map of particle connectivities. Color intensity is proportional to coordination number.

FIG. 15 :

 15 FIG.15: Contact anisotropy ac and branch vector anisotropy a ′ c as a function of shape parameter η averaged in residual state. The error bars represent the standard deviation in the residual state. The inset shows the angular probability densities P θ (θ) in black and P θ ′ (θ ′ ) in red for η = 0.7 calculated from the simulation data (points) together with the harmonic approximation (lines).

FIG. 17 :

 17 FIG. 17: Normal and tangential force anisotropies an and at and radial and orthoradial force anisotropies a ′ n and a ′ t as a function of η in the residual state. The error bars represent the standard deviation in the residual state. The inset shows the angular average functions fn (θ) and f n ′ (θ) in black and green, respectively, (a) and ft (θ) and f t ′ (θ) in red and blue, respectively, (b) for η = 0.7 calculated from the simulation data (points) together with the harmonic approximation (lines). The error bars represent the standard deviation in the residual state.

FIG. 19 :FIG. 20 :

 1920 FIG.19: Map of radial forces for η = 0.2 (up) and η = 0.7 (down). Line thickness is proportional to the radial force. We represent the strong network in black and the weak network in red lines (see section VI). The floating particles excluded from the force network are in white.

7 FIG. 21 :

 721 FIG. 21: Probability distribution function of radial forces f n ′ normalized by the average radial force f n ′ in log-linear (up) and log-log (down) plots for different values of η.

7 FIG. 22 :

 722 FIG. 22: Probability distribution function of orthoradial forces f t ′ normalized by the average orthoradial force |f t ′ | in log-log for all values of η.

FIG. 23 :

 23 FIG.23: Proportions of side-to-side (ss), cap-to-side (cs) and cap-to-cap (cc) contacts as a function of η in the residual state. The error bars represent the standard deviation in the residual state.

FIG. 25 :

 25 FIG. 25: Partial connectivity diagrams for all packings in the residual state.

FIG. 26 :FIG. 27 :

 2627 FIG. 26: Partial contact orientation anisotropies accc, accs and acss of cc, cs and ss contacts as the function of η in the residual state. The error bars represent the standard deviation in the residual state.

FIG. 28 : 1 ηFIG. 29 :

 28129 FIG.28:(Color online) Snapshot of radiale forces for η = 0.2 (up) and η = 0.7 (down). Line thickness is proportional to the radial force. The cap-to-cap, cap-to-side and side-to-side contacts are in black, in red (dark gray) and in green (light gray).


Emmanuel Russ

Yannick Sire

A NOTE ON PRECISED HARDY INEQUALITIES ON CARNOT GROUPS AND RIEMANNIAN MANIFOLDS

Keywords: AMS numbers 2000: Primary: 46E35. Secondary: 35R11, 42B35, 58J35 Hardy inequalities, Lie groups, Riemannian manifolds, fractional powers

We prove non local Hardy inequalities on Carnot groups and Riemannian manifolds, relying on integral representations of fractional Sobolev norms.

Introduction

In the whole paper, when two quantities A(f ) and B(f ) depend on a function f ranging in some space V , the notation

A(f ) B(f ) ∀f ∈ V
means that there exists C > 0 such that A(f ) ≤ CB(f ) for all f ∈ V . Moreover, the notation

A(f ) ∼ B(f ) ∀f ∈ V
means that there exists C > 0 such that C -1 B(f ) ≤ A(f ) ≤ CB(f ) for all f ∈ V .

1.1. Setting of the problem. The simplest Hardy inequality on R n asserts that, if n ≥ 3, (1.1)

R n u 2 (x) |x| 2 dx R n |∇u(x)| 2 dx = u Ḣ1 (R n ) , ∀u ∈ D(R n ).
A non local version of (1.1) can be given, where the Ḣ1 norm in the right hand side is replaced by a Ḣs norm for 0 < s < n 2 (see [START_REF] Bahouri | Precised Hardy inequalities[END_REF]):

(1.2)

R n u 2 (x) |x| 2s dx u 2 Ḣs (R n ) , ∀u ∈ D(R n ).
When 0 < s < 1, it is well-known (see for instance [START_REF] Adams | Pure and Applied Mathematics[END_REF]) that u Ḣs (R n ) can be represented by means of an integral quantity involving first order differences of u, and (1.2) can therefore be rewritten as

(1.3) R n u 2 (x) |x| 2s dx R n ×R n |u(x) -u(y)| 2
|x -y| n+2s dxdy.

These Hardy inequalities (i.e. the local and the non local version) were transposed to the framework of the Heisenberg group in [START_REF] Bahouri | Precised Hardy inequalities[END_REF][START_REF] Bahouri | Trace and trace lifting theorems in weighted Sobolev spaces[END_REF]. More precisely, in the Heisenberg group H d (d ≥ 1), the following Hardy inequality was established in [START_REF] Bahouri | Trace and trace lifting theorems in weighted Sobolev spaces[END_REF]:

(1.4)

H d u 2 (x) ρ 2 (x) dx ∇ H u 2 2 , ∀u ∈ D(H d ),
where ρ(x) denotes the distance of x to the origin and ∇ H stands for the gradient associated to the vector fields Z 1 , ..., Z 2d (see [START_REF] Bahouri | Trace and trace lifting theorems in weighted Sobolev spaces[END_REF] and the notations therein). The non local version of (1.4), which was proven in [START_REF] Bahouri | Precised Hardy inequalities[END_REF] (where it was derived from precised inequalities involving Besov norms) says that, for 0 < s < d + 1, (1.5)

H d u 2 (x) ρ 2s (x) dx u 2 Ḣs , ∀u ∈ D(H d ).
When 0 < s < 1, an integral representation for the fractional Sobolev homogeneous norm was proven in [START_REF] Coulhon | Sobolev algebras on Lie groups and Riemannian manifolds[END_REF] (note that an analogous representation holds in any connected Lie group wth polynomial volume growth, and even in any unimodular Lie group if one works with the inhomogeneous version of this norm), so that (1.5) can be rewritten as

(1.6)

H d u 2 (x) ρ 2s (x) dx H d ×H d |u(x) -u(y)| 2 ρ(y -1 x) 2d+2+2s dxdy.
Hardy inequalities in local versions on more general Lie groups, namely Carnot groups, were obtained in [START_REF] Kombe | Rellich and Uncertainty principle inequalities on Carnot Groups[END_REF]. In the present paper, we establish non local versions of these Hardy inequalities on Carnot groups, in the spirit of (1.3) and (1.6). We also investigate the similar problem on Riemannian manifolds.

1.2. The case of Lie groups. We now describe the general framework for Lie groups. Let G be a unimodular connected Lie group endowed with the Haar measure. By "unimodular", we mean that the Haar measure is left and right-invariant. If we denote by G the Lie algebra of G, we consider a family

X = {X 1 , ..., X k }
of left-invariant vector fields on G satisfying the Hörmander condition, i.e. G is the lie algebra generated by the X ′ i s. A standard metric on G , called the Carnot-Caratheodory metric, is naturally associated with X and is defined as follows: let ℓ : [0, 1] → G be an absolutely continuous path. We say that ℓ is admissible if there exist measurable functions a 1 , ..., a k : [0, 1] → C such that, for almost every t ∈ [0, 1], one has

ℓ ′ (t) = k i=1 a i (t)X i (ℓ(t)).
If ℓ is admissible, its length is defined by

|ℓ| = 1 0 k i=1 |a i (t)| 2 dt 1 2
. For all x, y ∈ G, define d(x, y) as the infimum of the lengths of all admissible paths joining x to y (such a curve exists by the Hörmander condition). This distance is left-invariant. For short, we denote by |x| the distance between e, the neutral element of the group and x, so that the distance from x to y is equal to |y -1 x|.

For all r > 0, denote by B(x, r) the open ball in G with respect to the Carnot-Caratheodory distance and by V (r) the Haar measure of any ball.

We denote

-∆ G = - k i=1 X 2 i the sub-laplacian on G and ∇ G = (X 1 , ..., X k ) the associated gradient.
The Lie group G is called a Carnot group if and only if G is simply connected and the Lie algebra of G admits a stratification, i.e. there exist linear subspaces

V 1 , ..., V k of G such that G = V 1 ⊕ ... ⊕ V k which [V 1 , V i ] = V i+1 for i = 1, ..., k -1 and [V 1 , V k ] = 0. By [V 1 , V i ] we mean the subspace of G generated by the elements [X, Y ] where X ∈ V 1 and Y ∈ V i .
Recall that the class of Carnot groups is a strict subclass of nilpotent groups. Moreover, if G is a Carnot group, there exists n ∈ N, called the homogeneous dimension of G, such that, for all r > 0,

(1.7) V (r) ∼ r n
(see [START_REF] Folland | Hardy spaces on homogeneous groups[END_REF]). The Heisenberg group H d is a Carnot group and n = 2d + 2.

Let G be a Carnot group, denote by δ the Dirac distribution supported at the origin and let u be a solution of

-∆ G u = δ. Define N(x) = u(x)
1 2-n for x = 0 and N(0) = 0. The function N is an homogeneous norm on N by [START_REF] Folland | Subelliptic estimates and function spaces on nilpotent Lie groups[END_REF]. Kombe [Kom06] proved the following Hardy inequality on G: let α > 2n and u ∈ D(G \ {0}), then the following holds (1.8)

n + α -2 2 2 G u 2 (x) |∇ G N(x)| 2 |N(x)| 2 N α (x)dx ≤ G |∇ G u(x)| 2 N α (x)dx.
We prove the following non local version of (1.8):

Theorem 1.1. Let G be a Carnot group with homogeneneous dimension n ≥ 3. Then for all α > 2n and all s ∈ (0, 1), (1.9) where C > 0. Then, for all α < C -1 and all u ∈ D(M \ ρ -1 (0)),

G u 2 (x) |∇ G N| |N| s N α dx G×G |u(x) -u(y)| 2 |y -1 x| n+s N α (x) dx dy ∀u ∈ D(G\{0}).
(1.12)

C -1 -α 2 2 M u ρ 2 (x)ρ α (x)dx ≤ M |du(x)| 2 ρ α (x)dx.
Moreover, if the codimension of ρ -1 (0) is greater than 2-α, (1.12) holds for all function u ∈ D(M) (see [START_REF] Carron | Inégalités de Hardy sur les variétés riemanniennes non compactes[END_REF], Théorème 1.4 and Remarque 1.5) .

We provide here a non local version of (1.12). To state this result, we introduce some extra assumptions on M. The first one is a Faber-Krahn inequality on M. For any bounded open subset Ω ⊂ M, denote by λ D 1 (Ω) the principal eigenvalue of -∆ on Ω under the Dirichlet boundary condition. If p ≥ n, consider the following Faber-Krahn inequality: there exists C > 0 such that

(1.13) λ D 1 (Ω) ≥ Cµ(Ω)
2 p for all bounded subset Ω ⊂ M.

Let Λ p > 0 be the greatest constant for which (1.13) is satisfied. In other words,

Λ p = inf λ D 1 (Ω) µ(Ω) 2 p
, where the infimum is taken over all bounded subsets Ω ⊂ M. The Faber-Krahn inequality (1.13) is satisfied in particular when an isoperimetric inequality holds on M: namely there exists C > 0 and p ≥ n such that, for all bounded smooth subset Ω ⊂ M,

(1.14) σ(∂Ω) ≥ Cµ(Ω) 1-1 p ,
where σ(∂Ω) denotes the surface measure of ∂Ω. If M has nonnegative Ricci curvature, then (1.14) with p = n and (1.13) with p = n are equivalent. More generally, if M has Ricci curvature bounded from below by a constant, (1.13) with p > 2n implies (1.14) with p 2 ([Car96], Proposition 3.1, see also [START_REF] Coulhon | Sobolev inequalities on graphs and on manifolds[END_REF] when the injectivity radius of M is furthermore assumed to be bounded). Note that there exists a Riemannian manifold satisfying (1.13) for some p ≥ n but for which (1.14) does not hold for any p ≥ n ([Car96], Proposition 3.4).

It is a well-known fact that (1.13) implies a lower bound for the volume of geodesic balls in M. Namely ([Car96], Proposition 2.4), if (1.13) holds, then, for all x ∈ M and all r > 0,

(1.15) V (x, r) ≥ Λ p 2 p+2 p 2 r p .
We will also need another assumption on the volume growth of the balls in M. Say that M has the doubling property if and only if there exists C > 0 such that, for all x ∈ M and all r > 0, (D)

V (x, 2r) ≤ CV (x, r).

There is a wide class of manifolds on which (D) holds. First, it is true on Lie groups with polynomial volume growth (in particular on nilpotent Lie groups). Next, (D) is true if M has nonnegative Ricci curvature thanks to the Bishop comparison theorem (see [START_REF] Bishop | Geometry of manifolds[END_REF]). Recall also that (D) remains valid if M is quasi-isometric to a manifold with nonnegative Ricci curvature, or is a cocompact covering manifold whose deck transformation group has polynomial growth, [START_REF] Coulhon | Variétés riemanniennes isométriques à l'infini[END_REF]. Contrary to the doubling property, the nonnegativity of the Ricci curvature is not stable under quasi-isometry.

We prove the following theorem.

Theorem 1.2. Let M be a complete non compact Riemannian manifold. Assume that (1.13) holds and that M has the doubling property. Assume also that C > 0 and ρ : M → [0, +∞) are such that (1.10) and (1.11) hold. Then, if α < C -1 and ρ -1 (0) has codimension greater than 2α, one has, for all s ∈ (0, 1), (1.16)

M u 2 (x)ρ s(α-2) (x)dx M ×M |u(y) -u(x)| 2 d(x, y) p+s ρ α (x)dxdy ∀u ∈ D(M\ρ -1 (0)).

Proof of Theorem 1.1

In order to prove Theorem 1.1, we need to introduce an operator L N α on L 2 (G). Let L 2 (G, N α ) denote the L 2 space on G with respect to the measure N α dx and H 1 (G, N α ) the Sobolev space defined by

H 1 (G, N α ) := f ∈ L 2 (G, N α ); X i f ∈ L 2 (G, N α ) ∀1 ≤ i ≤ k . Define now the operator L N α on L 2 (G, N α ) by L N α u := -N -α k i=1 X i (N α X i u).
The domain of L N α is given by

D (L N α ) = u ∈ L 2 (G, N α ); N -α X i (N α X i u) ∈ L 2 (G, N α ) ∀1 ≤ i ≤ k . One has, for all u ∈ D (L N α ) and all v ∈ H 1 (G, N α ), G L N α u(x)v(x)N α (x)dx = G k i=1 X i u(x)X i v(x)N α (x)dx.
The operator L N α is therefore clearly symmetric and nonnegative on L 2 (G, N α ), and the spectral theorem allows to define the usual powers (L N α ) β for any β > 0 by means of spectral theory.

By the definition of L N α , (1.8) means, in terms of operators in L 2 (G, N α ), that, for some λ > 0, (2.17)

L N α ≥ λµ,
where µ is the multiplication operator by |∇ G N | |N | . Using a functional calculus argument (see [START_REF] Davies | One-parameter semigroups[END_REF], p. 110) one deduces from (2.17) that, for any s ∈ (0, 2),

(L N α ) s/2 ≥ λ s/2 µ s/2
which implies, thanks to the fact (L N α ) s/2 = ((L N α ) s/4 ) 2 and the symmetry of (

L N α ) s/4 on L 2 (G, N α ), that G |u(x)| 2 µ(x) s N α (x)dx ≤ C G (L N α ) s/4 u(x) 2 N α (x)dx = C (L N α ) s/4 u 2 L 2 (G,N α )
. The conclusion of Theorem 1.1 follows now from the estimate

(L N α ) s/4 u 2 L 2 (G) ≤ C G×G |u(x) -u(y)| 2 |y -1 x| n+2s N α (x) dx dy,
which is exactly the estimate for M = N α provided in Lemmata 3.2 and 3.3 in [START_REF] Russ | Non local poincaré inequalities on lie groups with polynomial volume growth[END_REF] (remember that (1.7) holds).

Proof of Theorem 1.2

The proof of Theorem 1.2 relies again on estimates for the powers of a suitable operator. Namely, define

L ρ α u = ρ -α div(ρ α ∇u),
where ∇ is the gradient induced by the Riemannian metric and div is the divergence operator on M. As before, L ρ α is a nonnegative symmetric operator on L 2 (M, ρ α dx) and L β ρ α is defined by spectral theory for all β > 0. If µ denotes the multiplication operator by ρ -2 , (1.12) means that L ρ α ≥ cµ in L 2 (M, ρ α dx). Spectral theory then yields L s/2 ρ α ≥ cµ s/2 , which means that

M u 2 (x)ρ s(α-2) (x)dx L s/4 ρ α u 2 L 2 (M,ρ α dx)
, and we are therefore left with the task of checking

L s/4 ρ α u 2 L 2 (M,ρ α dx) M ×M |u(y) -u(x)| 2 d(x, y) p+s ρ α (x)dxdy.
To that purpose, one first notices that

L s/4 ρ α u 2 L 2 (M,ρ α dx) +∞ 0 t -1-s/2 t L ρ α (I + t L ρ α ) -1 u 2 L 2 (M,ρ α dx) dt
and it is therefore enough to show that

+∞ 0 t -1-s/2 t L ρ α (I + t L ρ α ) -1 u 2 L 2 (M,ρ α dx) dt M ×M |u(y) -u(x)| 2 d(x, y) p+s ρ α (x)dxdy.
The proof follows the same lines as the one of Lemma 3.3 in [START_REF] Russ | Non local poincaré inequalities on lie groups with polynomial volume growth[END_REF] and we will therefore be sketchy, only indicating the main differences. Using (D), one can pick up a countable family x t j , j ∈ N, such that the balls B x t j , √ t are pairwise disjoint and

(3.18) M = j∈N B x t j , 2 √ t .
By (D), there exist constants C > 0 and κ > 0 such that for all θ > 1 and all x ∈ G, there are at most C θ 2κ indexes j such that

|x -1 x t j | ≤ θ √ t.
For fixed j, one has

t L ρ α (I + t L ρ α ) -1 u = t L ρ α (I + t L ρ α ) -1 g j,t
where, for all x ∈ M, g j,t (x) := u(x)m j,t and m j,t is defined by

m j,t := 1 V x t j , 2 √ t B(x t j ,2 √ t)
u(y)dy.

Note that, here, the mean value of u is computed with respect to the Riemannian measure on M. Since (3.18) holds, one clearly has

t L ρ α (I + t L ρ α ) -1 u 2 L 2 (M,ρ α dx) ≤ j∈N t L ρ α (I + t L ρ α ) -1 u 2 L 2 (B(x t j ,2 √ t),ρ α dx) = j∈N t L ρ α (I + t L ρ α ) -1 g j,t 2 L 2 (B(x t j ,2 √ t),ρ α dx) ,
and it is therefore enough to ckeck

j∈N t L ρ α (I + t L ρ α ) -1 g j,t 2 L 2 (B(x t j ,2 √ t),ρ α dx) (3.19) M ×M |u(y) -u(x)| 2 d(x, y) p+s ρ α (x)dxdy.
As in [START_REF] Russ | Non local poincaré inequalities on lie groups with polynomial volume growth[END_REF], this is a consequence of L 2 off-diagonal estimates for L ρ α and upper estimates for the functions g j,t . Let us recall the off-diagonal estimates for L ρ α for completeness:

Lemma 3.1. There exists C with the following property: for all closed disjoint subsets E, F ⊂ M with d(E, F ) =: d > 0, all function f ∈ L 2 (M, ρ α dx) supported in E and all t > 0,

(I + t L ρ α ) -1 f L 2 (F,ρ α dx) + t L ρ α (I + t L ρ α ) -1 f L 2 (F,ρ α dx) ≤ 8 e -C d √ t f L 2 (E,ρ α dx) .
The proof of Lemma 3.1 is analogous to the one of Lemma 3.1 in [START_REF] Russ | Non local poincaré inequalities on lie groups with polynomial volume growth[END_REF].

As far as estimates for g j,t are concerned, set, for all k ≥ 1, C j,t 0 = B x t j , 4 √ t and C j,t k = B x t j , 2 k+2 √ t \ B x t j , 2 k+1 √ t , and g j,t k := g j,t 1 C j,t k , k ≥ 0, where, for any subset A ⊂ M, 1 A is the usual characteristic function of A. We then have: Lemma 3.2. There exists C > 0 such that, for all t > 0 and all j ∈ N:

A. |u(x)u(y)| 2 ρ α (x)dx dy.

g j,
The proof of Lemma 3.2 is analogous to the one of Lemma 3.4 in [START_REF] Russ | Non local poincaré inequalities on lie groups with polynomial volume growth[END_REF], the only extra ingredient being the lower bound (1.15) applied with the balls B(x t j , 2 k+2 √ t). We then conclude the proof of (3.19) in the same way as for the conclusion of the proof of Lemma 3.3 in [START_REF] Russ | Non local poincaré inequalities on lie groups with polynomial volume growth[END_REF].

Yannick Sire-Université Paul Cézanne, LATP, Faculté des Sciences et Techniques, Case cour A Avenue Escadrille Normandie-Niemen, F-13397 Marseille, Cedex 20, France et CNRS, LATP, CMI, 39 rue F. Joliot-Curie, F-13453 Marseille Cedex 13, France.

1. 3 .

 3 The case of Riemannian manifolds. A general principle was developed in[START_REF] Carron | Inégalités de Hardy sur les variétés riemanniennes non compactes[END_REF] to derive Hardy inequalities on Riemannian manifolds. Let us recall here an example of such an inequality. Let M be a Riemannian manifold, denote by n its dimension, by dµ its Riemannian measure, by d the exterior differentiation and by ∆ the Laplace-Beltrami operator. For all x ∈ M and all r > 0, let B(x, r) be the open geodesic ball centered at x with radius r, and V (x, r) its measure. Assume that ρ : M → [0, +∞) satisfies (1.10) |dρ| ≤ 1 on M, and (1.11) ∆ρ ≤ -C ρ in the distribution sense,


Hermann Pfitzner

THE DUAL OF A NON-REFLEXIVE L-EMBEDDED BANACH SPACE CONTAINS l ∞ ISOMETRICALLY

Keywords: 1991 Mathematics Subject Classification. Primary 46B20; Secondary 46B03, 46B04, 46B26 embedded Banach spaces, isometric copies of c 0

A Banach space is said to be L-embedded if it is complemented in its bidual such that the norm between the two complementary subspaces is additive.

This note is an afterthought to a result of Dowling [START_REF] Dowling | Isometric copies of c 0 and ℓ ∞ in duals of Banach spaces[END_REF] according to which a dual Banach space contains an isometric copy of c 0 if it contains an asymptotic one. (For definitions see below.) It is known ( [START_REF] Pfitzner | L-summands in their biduals have Pe lczyński's property (V * )[END_REF] or [4, Th. IV.2.7]) that the dual of a non-reflexive L-embedded Banach space contains c 0 isomorphically. For a special class of L-embedded Banach spaces the construction of the c 0 -copy has been improved so to yield an asymptotic one [START_REF] Pfitzner | A note on asymptotically isometric copies of l 1 and c 0[END_REF]Prop. 6]) and it turns out that this improvement is possible in the general case which together with Dowling's result yields isometric copies of c 0 in the dual of an L-embedded Banach space.

All this is perhaps known -or at least not surprising -to experts in the field but the final result, i.e. an isometric copy of c 0 (and hence of l ∞ ) in the dual of an L-embedded space X, is optimal in the category of Banach spaces and therefore it seems worthwile proving it explicitely. As in [START_REF] Pfitzner | L-summands in their biduals have Pe lczyński's property (V * )[END_REF] we will prove a bit more by constructing the c 0 -copy within the context of Pe lczyński's property (V * ) i.e. the c 0 -basis will be constructed so to behave approximately like biorthogonal functionals on the basis of a given l 1 -basis in X, see (2) and (3) below where in particular the value cJ (x n ) in (2) is optimal. (For the definition and some basic results on Pe lczyński's property (V * ) see [START_REF] Harmand | M -ideals in Banach Spaces and Banach Algebras[END_REF].)

Preliminaries: A projection P on a Banach space Z is called an L-projection if P z + z -P z = z for all z ∈ Z. A Banach space X is called L-embedded (or an L-summand in its bidual) if it is the image of an Lprojection on its bidual. In this case we write X * * = X ⊕ 1 X s . Among classical Banach spaces, the Hardy space H 1 0 , L 1 -spaces and, more generally, the preduals of von Neumann algebras or of JBW * -triples serve as examples of L-embedded spaces. A sequence (x n ) in a Banach space X is said to span c 0 asymptotically isometrically (or just to span c 0 asymptotically) if there is a null sequence

(δ n ) in [0, 1[ such that sup(1 -δ n )|α n | ≤ α n x n ≤ sup(1 + δ n )|α n | for all (α n ) ∈ c 0 .
X is said to contain c 0 asymptotically if it contains such a sequence (x n ). Recall the routine fact that if (x * n ) in X * is equivalent to the canonical basis of c 0 then α n x * n makes sense for all (α n ) ∈ l ∞ in the w * -topology of X * and by lower w * -semicontinuity of the norm an estimate

α n x * n ≤ M sup |α n | that holds for all (α n ) ∈ c 0 extends to all (α n ) ∈ l ∞ .
The Banach spaces we consider in this note are real or complex, the set IN starts at 1. To a bounded sequence (x n ) in a Banach space X we associate its 'James constant'

c J (x n ) = sup c m where the c m = inf n≥m |αn|=1 n≥m α n x n
form an increasing sequence. If (x n ) is equivalent to the canonical basis of l 1 then c J (x n ) > 0 and more specifically, c J (x n ) > 0 if and only if there is an integer m such that (x n ) n≥m is equivalent to the canonical basis of l 1 . Roughly speaking, the number c J (x n ) may be thought of as the 'approximately best l 1 -basis constant' of (x n ); precisely speaking, there is a null sequence

(τ m ) in [0, 1[ (determined by c m = (1 -τ m )c J (x n )) such that ∞ n=m α n x n ≥ (1 -τ m )c J (x n ) ∞ n=m |α n |
for all (α n ) ∈ l 1 and c J (x n ) cannot be replaced by a strictly greater constant. If one passes to a subsequence

(x n k ) of (x n ) then c J (x n k ) ≥ c J (x n ) hence it makes sense to define cJ (x n ) = sup n k c J (x n k ).
The standard reference for L-embedded Banach spaces is the monograph [4, Ch. IV]. For general Banach space theory and undefined notation we refer to [START_REF] Diestel | Sequences and Series in Banach Spaces[END_REF], [START_REF] Johnson | Handbook of the Geometry of Banach Spaces[END_REF], or [START_REF] Lindenstrauss | Classical Banach Spaces I and II[END_REF].

The main result of this note is Theorem 1. Let X be an L-embedded Banach space and let (x n ) be equivalent to the canonical basis of l 1 . Then there is a sequence

(x * n ) in X * that generates l ∞ isometrically, more precisely α n x * n = sup |α n | for all (α n ) ∈ l ∞ (1)
and there is a strictly increasing sequence

(p n ) in IN such that lim |x * n (x pn )| = cJ (x m ) (2) x * n (x p l ) = 0 if l < n (3)
In particular, the dual of a non-reflexive L-embedded Banach space contains an isometric copy of l ∞ .

In order to prove the theorem we first state and prove Dowling's result in a way which fits our purpose.

Proposition 2. Let (ε n ) be a null sequence in [0, 1[, let (N n ) be a sequence of pairwise disjoint infinite subsets of IN and let (y * n ) in the dual of a Banach space Y span c 0 such that α n y * n ≤ sup(1 + ε n )|α n | and y * n → 1 (4) for all (α n ) ∈ c 0 .
Then the elements

x * n = k∈Nn y * k 1 + ε k (5)
generate l ∞ isometrically (as in [START_REF] Diestel | Sequences and Series in Banach Spaces[END_REF]).

Proof of the proposition: With (ε n ), (N n ) and (y * n ) as in the hypothesis of the statement define x * n by [START_REF] Johnson | Handbook of the Geometry of Banach Spaces[END_REF]. Then x * n ≤ 1 for all n ∈ IN by the first half of (4). For the inverse inequality we have that

x * n ≥ 2 y * m 1 + ε m - y * m 1 + ε m - k∈Nn,k =m y * k 1 + ε k ≥ 2 y * m 1 + ε m -1
holds for all m ∈ N n hence x * n ≥ 1 by the second half of (4) which proves

x * n = 1.
Similarly we show (1): First, "≤" of (1) follows from the first half of (4); second, by the just shown inequality we have

α n x * n ≥ 2|α m | -α m x * m - n =m α n x * n ≥ 2|α m | -sup |α n |
for all m ∈ IN hence "≥" of (1).

Proof of the theorem:

Let (δ n ) be a sequence in ]0, 1[ converging to 0. Suppose (x n ) is an l 1 -basis and write c = cJ (x n ) for short.

Observation: Given τ > 0 there is a subsequence (x n k ) of (x n ) such that c J (x n k ) > (1 -τ )c
and by James' l 1 -distortion theorem there are blocks of the x n k which span l 1 almost isometrically that is to say there are pairwise disjoint finite sets

A l ⊂ {n k | k ∈ IN}, a sequence of scalars (λ n ) such that k∈A l |λ k | = 1 and such that the sequence (z l ) defined by z l = zl / zl and zl = k∈A l λ k x k satisfies (1 -2 -m ) ∞ l=m |α l | ≤ ∞ l=m α l z l ≤ (1 + 2 -m ) ∞ l=m |α l | for all m ∈ IN; furthermore zl → c J (x n k ) whence the existence of l ′ such that |c -zl ′ | < τ .
By induction over n ∈ IN we will construct finite sequences (y

(n) * i
) n i=1 in X * , a sequence (ỹ n ) in X, pairwise disjoint finite sets C n ⊂ IN and a scalar sequence (µ n ) such that, with the notation y n = ỹn / ỹn ,

k∈Cn |µ k | = 1, ỹn = k∈Cn µ k x k , |c -ỹn | < δ n , (6) |y
(n) * i (y i )| > 1 -δ i ∀i ≤ n, (7)
y (n) * i (y l ) = 0 ∀l < i ≤ n, (8)
y (n) * i (x p ) = 0 ∀p ∈ C l , ∀l < i ≤ n, (9)
m i=1 α i y (n) * i ≤ max i≤m (1 + (1 -2 -n ) δ i ) |α i | ∀m ≤ n, α i scalars. ( 10
)
For n = 1 we use the observation with τ = δ 1 and choose l 1 such that | zl 1 -c| < δ 1 . Then we choose y

(z l 1 ) = z l 1 . It remains to set C 1 = A l 1 , µ k = λ k for k ∈ C 1 and ỹ1 = zl 1 .
For the induction step n → n + 1 we recall that (P * ) |X * is an isometric isomorphism from X * onto X ⊥ s , that X * * * = X ⊥ ⊕ ∞ X ⊥ s and that (P * x * ) |X = (x * ) |X for all x * ∈ X * . Let (z l ) be as in the observation above with τ = δ n+1 and let z s ∈ X * * \ X be a w * -accumulation point of the z l . Then z s ∈ X s and z s = 1 by the proof of [8, Lem 1] (or by some general folklore argument). Choose t ∈ ker P * ⊂ X * * * such that t = 1 and t(z s ) = z s . Put

E = lin({P * y (m) * i | i ≤ m ≤ n} ∪ {t}) ⊂ X * * * , F = lin({y i | i ≤ n} ∪ {z s } ∪ {x p | p ∈ l≤n C l }) ⊂ X * *
and choose η > 0 such that

(1 + η)(1 + (1 -2 -n ) δ i ) < 1 + (1 -2 -(n+1) ) δ i and η < (1 -2 -(n+1) ) δ n+1
for all i ≤ n. The principle of local reflexivity provides an operator R : E → X * such that

α i = 0 if m < i ≤ n + 1) by n+1 i=1 α i y (n+1) * i (11) ≤ (1 + η) n i=1 α i P * y (n) * i + α n+1 t = (1 + η) max n i=1 α i P * y (n) * i , α n+1 t = (1 + η) max n i=1 α i y (n) * i , α n+1 t (10) ≤ (1 + η) max max i≤n (1 + (1 -2 -n ) δ i )|α i |, |α n+1 | ≤ max i≤n+1 (1 + (1 -2 -(n+1) ) δ i )|α i |.
Since z s is a w * -cluster point of (z l ) we have

|y (n+1) * n+1 (z l )| > |z s (y (n+1) * n+1 )| -δ n+1 (12) = |t(z s )| -δ n+1 = 1 -δ n+1
for infinitely many l; furthermore, an l n+1 can be chosen among these l so to obtain [START_REF] Lindenstrauss | Classical Banach Spaces I and II[END_REF] holds and (7, n + 1) holds for i = n + 1. For i ≤ n, (7, n + 1) follows from

| zl n+1 -c| < δ n+1 . Set C n+1 = A l n+1 , ỹn+1 = zl n+1 , µ k = λ k for k ∈ C n+1 . Then
y (n+1) * i (y i ) = (P * y (n) * i )(y i ) = y (n) * i (y i ) (7) > 1 -δ i .
Condition (8, n + 1) holds for i = n + 1 by

y (n+1) * n+1 (y l ) = (Rt)(y l ) (12) = t(y l ) = 0 ∀ l < n + 1
and it holds for i < n + 1 by

y (n+1) * i (y l ) = (P * y (n) * i )(y l ) = y (n) * i (y l ) (8) = 0 ∀ l < i.
The proof of (9, n + 1) works like the one of (8, n + 1). This ends the induction.

Now we define y

* i = 1 1+δ i lim n∈U y (n) * i
for all i ∈ IN where U is a fixed nontrivial ultrafilter on IN and where the limit is understood in the w *topology of X * . Then by w * -lower semicontinuity of the norm and by [START_REF] Simons | On the Dunford-Pettis Property and Banach spaces that contain c 0[END_REF] [START_REF] Pfitzner | L-summands in their biduals have Pe lczyński's property (V * )[END_REF] and (y * i ) satisfies (4) for ε n = 0.

α i y * i ≤ sup(1 + δ i ) |α i | 1 + δ i = sup |α i | for all (α i ) ∈ l ∞ . In particular, y * i ≤ 1 hence y * i → 1 by
Let (N n ) be a sequence of pairwise disjoint infinite subsets of IN such that (i n ) increases strictly where i n = min N n . By the proposition the sequence defined by

x * n = i∈Nn y * i generates l ∞ isometrically and we have |x * n (y in )| (8)
= |y * in (y in )|

≥ 1δ in . By construction of the y i there is, for each n ∈ IN, an index

p n ∈ C in such that |x * n (x pn )| ≥ (1 -δ in ) ỹin (6)
≥ (1δ in )(cδ in ) which will yield "≥" of (2). In order to show "≤" of (2) suppose to the contrary that x * nm (x pn m ) > κ+c for appropriate subsequences, all m and κ > 0. According to an extraction lemma of Simons [START_REF] Simons | On the Dunford-Pettis Property and Banach spaces that contain c 0[END_REF] we may furthermore suppose that

j =m |x * n j (x pn m )| < κ/2 for all m. Then given (α m ) and θ m such that θ m α m = |α m | we obtain α m x pn m ≥ ( j θ j x * n j )( m α m x pn m ) ≥ (κ + c) m |α m |-m j =m |α m | |x * n j (x pn m )| ≥ ((κ/2)+c) m |α m |
which yields the contradiction c J (x pn m ) > c and thus shows "≤" and all of (2) whereas (3) follows from (9) via y * i (x p ) = 0 for p ∈ C l , l < i. The last assertion of the theorem is immediate from the fact that nonreflexive L-embedded spaces contain l 1 isomorphically [4, IV.2.3] Remarks: 1. It is not clear whether (3) can be obtained also for l > n. What can be said by Simons' extraction lemma (used in the proof) is that, under the assumptions of the theorem and given ε > 0, it is possible (after passing to appropriate subsequences) to obtain in addition to (3) that l-1 n=1 |x * n (x p l )| = n =l |x * n (x p l )| < ε for all l. In case cJ (x n ) = 1 = lim x n (which happens when the x n span l 1 almost isometrically) this can be improved to

n =l |x * n (x p l )| = ( n |x * n (x p l )|) -|x * l (x p l )| ≤ x p l -|x * l (x p l )| → 0. ( 13
)
One might also construct straightforward perturbations of the x * n in order to get (3) for l = n but then it is not clear whether these perturbations can be arranged to span c 0 isometrically, not just almost isometrically.

Since in general L-embedded spaces do not contain l 1 isometrically (see below, last remark) it is in general not possible to improve (2) and (3) so to obtain x * n (x p l ) = c(x m ) if l = n and = 0 if l = n. 2. As already alluded to in the introduction, the construction of c 0 in this paper bears much resemblance to the one of [START_REF] Pfitzner | L-summands in their biduals have Pe lczyński's property (V * )[END_REF]. A different way to construct c 0 is contained in [START_REF] Pfitzner | Separable L-embedded Banach spaces are unique preduals[END_REF] but it seems unlikely that this construction can be improved to yield an isometric c 0 -copy. 3. It follows from (2) that in L-embedded spaces the sup in the definition of cJ is attained by the James constant of an appropriate subsequence. For general Banach spaces this is not known although it can be shown by a routine diagonal argument that each bounded sequence (x n ) admits a c J -stable subsequence (x n k ) (meaning that cJ (x n k ) = c J (x n k )) whose James constant is arbitrarily near to c(x n ). 4. Each normalized sequence (x n ) in an L-embedded Banach space that spans l 1 almost isomorphically contains a subsequence each of whose w *accumulation points in the bidual attains its norm on the dual unit ball. To see this let (x * n ) and (x pn ) be the sequences given by the theorem and by Simons' extraction lemma (see (13) above), let x s be a w * -accumulation point of the x pn and let x * = x * n ; then x * = 1 and on the one hand x s = 1 by [START_REF] Pfitzner | A note on asymptotically isometric copies of l 1 and c 0[END_REF] and on the other hand x s (x * ) = lim x * (x pn )

(13) = lim x * n (x pn ) (2)
= c J (x n ) = 1. It would be interesting to know whether this remark holds for the whole sequence (x n ) instead of only a subsequence (x pn ). A kind of converse follows from [9, Rem. 2] for separable X: If x s ∈ X s attains its norm on the dual unit ball then it does so on the sum of a wuC-series. 5. Let us finally note that the presence of isometric c 0 -copies in X * does not necessarily entail the presence of isometric copies of l 1 in X even if X is the dual of an M-embedded Banach space. This follows from [4, Cor. III.2.12] which states that there is an L-embedded Banach space which is the dual of an M-embedded space (to wit the dual of c 0 with an equivalent norm) which is strictly convex and therefore does not contain l 1 isometrically although it contains, as do all non-reflexive L-embedded spaces, l 1 asymptotically ( [START_REF] Pfitzner | A note on asymptotically isometric copies of l 1 and c 0[END_REF], see [START_REF] Dowling | The optimality of James's distortion theorems[END_REF] for the definition of asymptotic copies and the difference to almost isometric ones).

( 1 -=

 1 η) e * * * ≤ Re * * * ≤ (1 + η) e * * * , (11) f * * (Re * * * ) = e * * * (f * * ), (12) for all e * * * ∈ E and f * * ∈ F . We define y (n+1) * i = R(P * y (n) * i ) for i ≤ n and y (n+1) * n+1 Rt and obtain (10, n + 1) (with