F A Gilabert

J.-N Roux

A Castellanos

Computer simulation of model cohesive powders: influence of assembling procedure and contact laws on low consolidation states

Keywords: numbers: 81.05.Rm: Porous materials, granular materials, 83.10.Rs: Computer simulation of molecular and particle dynamics, 61.43.Hv: Fractals, macroscopic aggregates, 47.57.J-: Colloidal systems

Molecular dynamics simulations are used to investigate the structure and mechanical properties of a simple two-dimensional model of a cohesive granular material. Intergranular forces involve elasticity, Coulomb friction and a short range attraction akin to the van der Waals force in powders. The effects of rolling resistance (RR) at intergranular contacts are also studied. The microstructure of the cohesive packing under low pressure is shown to depend sensitively on the assembling procedure which is applied to the initially isolated particles of a granular gas. While a direct compression produces a final equilibrated configuration with a similar density to that of cohesionless systems, the formation of large aggregates prior to the application of an external pressure results in much looser stable packings. A crucial state variable is the ratio P * = P a/F0 of applied pressure P , acting on grains of diameter a, to maximum tensile contact force F0. At low P * the force-carrying structure and force distribution are sensitive to the level of velocity fluctuations in the early stages of cluster aggregation. The coordination number of packings with RR approaches 2 in the limit of low initial velocities or large rolling friction. In general the force network is composed of hyperstatic clusters, typically comprising 4 to a few tens of grains, in which forces reach values of the order of F0, joined by barely rigid arms, where contact forces are very small. Under growing P * , it quickly rearranges into force chain-like patterns that are more familiar in dense systems. Density correlations are interpreted in terms of a fractal structure, up to a characteristic correlation length ξ of the order of ten particle diameters for the studied solid fractions. The fractal dimension in systems with RR coincides, within measurement uncertainties, with the ballistic aggregation result, in spite of a possibly different connectivity, but is apparently higher without RR. Possible effects of micromechanical and assembling process parameters on mechanical strength of packings are evoked.

I. INTRODUCTION

Granular materials are currently being studied by many research groups [START_REF]Physics of Dry Granular Media[END_REF][START_REF] Kishino | Powders and Grains 2001[END_REF][START_REF]The Physics of Granular Media[END_REF][START_REF]Powders and Grains 2005[END_REF], motivated by fundamental issues (such as the relations between microstructure and global properties) as well as by practical needs in civil engineering and in the food and drug industries. The relation of their mechanical behavior in quasistatic conditions to the packing geometry, which depends itself on the assembling procedure, tends to escape intuition and familiar modelling schemes.

The configuration of the contact networks is hardly accessible to experiments, even though particle positions are sometimes measured [START_REF] Richard | [END_REF]6,7,8], and some experimental quantitative studies on intergranular contacts carried out in favorable cases (such as millimeter-sized beads joined by capillary menisci [6,9,10,11]). Intergranular † LMSGC is a joint laboratory depending on Laboratoire Central des Ponts et Chaussées, École Nationale des Ponts et Chaussées and Centre National de la Recherche Scientifique. * Electronic address: gilav@us.es forces are also, most often, inaccessible to measurements. Consequently, computer simulation methods of the "discrete element" type, as introduced 30 years ago [12], have proved a valuable tool to investigate the internal states of granular systems. Simulation methods like molecular dynamics [13] or "contact dynamics" [14,15] have been gaining an increasingly large constituency of users and a wide range of applications, as witnessed e.g., by recent conference proceedings [START_REF]Powders and Grains 2005[END_REF].

Dry assemblies of grains interacting via contact elasticity and friction, such as sands or glass beads, might form stable packings of varying solid fraction (typically between 58% and 64% for monosized spheres if they do not crystallize), which deform plastically in response to changes in stress direction, rather than stress intensity. Their elastic or elastoplastic properties have been studied by discrete simulation (see, e.g., [16,17]), and, in agreement with laboratory experiments and macroscopic modelling [START_REF] Wood | Soil Behaviour and Critical State Soil Mechanics[END_REF], found to depend sensitively on the initial density. Numerical simulation also stressed the importance of additional variables such as coordination number [START_REF] Agnolin | [END_REF] and fabric [START_REF] Bathurst | [END_REF]21], and it has often been applied to the study of quasi-static stress-strain behavior of gran-ular assemblies (refs. [16,17,21,22] are a few examples among a large literature).

Cohesive grains exhibit much larger variations in their equilibrium densities, and they are sensitive to stress intensity as well as direction : on increasing the confining pressure, the specific volume of a clay can irreversibly decrease by a factor of 4 [START_REF] Mitchell | Fundamentals of soil behavior[END_REF]. Likewise, series of experiments carried out in the Seville group on model powders [START_REF] Castellanos | [END_REF] (xerographic materials) in which the strength of van der Waals attraction is controlled by additives covering part of the grain surfaces, reveal a similar variation of porosity with confining pressure. It is notable that such packings of particles of rotund shape and nearly the same size can stay in mechanical equilibrium at much lower solid fractions (down to 25-30%) than cohesionless granular systems.

Despite this wider variety of equilibrium structures and mechanical behaviors, cohesive granular materials have much less frequently been investigated by numerical simulation than cohesionless ones. Some of the recent numerical studies, such as those of refs. [25,26] have investigated the packing structures of spherical beads deposited under gravity, depending on micromechanical parameters, including adhesion strength. Another set of publications report on simulations of the dynamical collapse and compaction, both in two [27,28,29] and three [30] dimensions, the main results being the relations between density and pressure increments, and their dependence on micromechanical parameters. Some works focussed on the fracture of bound particle assemblies in static [31,32] or dynamic [33] conditions, others on wet bead packs in which cohesion stems from liquid bridges joining neighboring particles, investigating the structure of poured samples [34] or the shear strength [11] of such materials. These two latter types of studies deal with relatively dense materials, as does the numerical biaxial compression test of [35]. Flow of cohesive materials has also been addressed in recent publications [36,37,38].

Yet, numerical studies of the mechanics of loose, solidlike cohesive granulates are quite scarce. This contrasts with the abundant literature on the geometry of model loose particle packings and colloidal aggregates, which tend to form fractal structures. Refs. [39,40] are useful overviews of aggregation processes and the geometric properties of the resulting clusters, as obtained by numerical simulation. In such processes, particle aggregates are usually regarded as irreversibly bound, rigid solids, while the interaction between separate clusters reduces to a "sticking rule", so that both intra-and inter-aggregate mechanical modelling is bypassed. Interestingly, one simulation study [41] shows that structures resulting from geometric deposition algorithms are not always stable once a mechanical model is introduced.

It seems that numerical simulations of both geometric and mechanical properties of loose granular assemblies forming solid aggregates are still lacking.

The present paper addresses part of this issue. It re-ports on numerical simulation studies of cohesive granular materials, with the following specificities:

• the assembling process is simulated with the same mechanical model as applied to solid-like configurations, and its influence on the packing microstructure is assessed ;

• special attention is paid to loose particle packings in equilibrium under vanishing or low applied pressure;

• both geometric and mechanical properties are investigated ;

• isotropic and homogeneous systems are studied, as representative samples for bulk material properties.

We consider a simple model system in two dimensions, introduced in section II, along with the numerical simulation procedure. Despite its simplicity we shall see that this model yields results that are amenable to comparisons with experimental situations.

Section. III is devoted to the important issue of the procedure to prepare samples, and its influence, as well as that of micromechanical features such as rolling resistance (RR), on final density and coordination number in solid packings in equilibrium. In Section IV we investigate the force distributions and force patterns of the equilibrated loose configurations under vanishing or low applied pressure. Some specific aspects of the forcecarrying structures in low density assemblies will be studied and related to the assembling process. In Section V, we characterize the geometry and density correlations in loose samples, resorting to the fractal model traditionally employed for colloidal aggregates. Finally we conclude in Section VI with a few remarks about future improvements and further developments of this work, some of which will be presented in a forthcoming publication [START_REF] Gilabert | Computer simulation of model cohesive powders: plastic compaction under isotropic loads[END_REF].

II. MODEL MATERIAL A. System definition, equations of motion

We consider a two-dimensional model material: an assembly of N disks with diameters (d i ) 1≤i≤N uniformly distributed between a/2 and a. The maximum diameter, a, will be used as unit of length. The mass of grain i is

m i = d 2
i /a 2 and its moment of inertia I i = m i d 2 i /8, i.e. disks are regarded as homogeneous bodies and the mass of a disk of maximum diameter a is the unit of mass.

The disks are enclosed in a rectangular cell the edges of which are parallel to the axes of coordinates x 1 and x 2 , with respective lengths L 1 and L 2 . Periodic boundary conditions are used, thereby avoiding wall effects. Neighboring grains, say i and j, might interact if they are brought into contact or very close to each other, hence a force F ij and a moment Γ ij exerted by i onto j at the contact point. Simulations do not model material deformation in a contact region, but consider overlapping particles, and the contact point is defined as the center of the intersecting surface of the two disks. In the case of an interaction without contact, the force will be normal to the surfaces at the points of nearest approach, and therefore carried by the line of centers. Let r i denote the position of the center of disk i. r ij = r j -r i is the vector joining the centers of i and j, and h ij = | r ij | -(d i + d j )/2 their overlap distance. The degrees of freedom, in addition to the positions r i , are the angles of rotation θ i , velocities v i , angular velocities ω i = θi of the grains (1 ≤ i ≤ N ), the dimensions (L α ) α=1,2 of the cell containing the grains and their time derivatives, through the strain rates: ǫα = -Lα /L 0 α , in which L 0 α denotes the initial size for the corresponding compression process. The time evolution of those degrees of freedom is governed by the following equations.

m i d 2 r i dt 2 = N j=1 F ij (1)
I i dω i dt = N j=1 Γ ij (2)
M d 2 ǫ α dt 2 = σ I αα -σ M αα σ M αα = 1 A N i=1   m i v 2 i,α + j =i F (α) ij r (α) ij   (3)
In Eqns. ( 1) and [START_REF] Kishino | Powders and Grains 2001[END_REF], only those disks j interacting with i, i. e. in contact or very close, will contribute to the sums on the right-hand side. In Eqn. (3), σ I αα is the externally imposed stress component, σ M αα is the measured stress component, resulting from ballistic momentum transport and from the set of intergranular forces F ij , A = L 1 L 2 denotes the cell surface area, and M is a generalized inertia parameter.

Stresses σ 11 and σ 22 , rather than strains or cell dimensions, are controlled in our simulation procedure. Note that compressions are counted positively for both stresses and strains. Eqn. (3) entails that the sample will expand (respectively, shrink) along direction α if the corresponding stress σ M αα is larger (resp., smaller) than the requested value σ I αα , which should be reached once the system equilibrates. This barostatic method is adapted from the ones initially proposed by Parrinello and Rahman [43,44,45] for Hamiltonian, molecular systems.

The choice of the "generalized mass" M is rather arbitrary, yet innocuous provided calculations are restricted to small strain rates. In practice we strive to approach mechanical equilibrium states with good accuracy, and choose M in order to achieve this goal within affordable computation times. We usually attribute to M a value equal to a fraction of the sum of grain masses (3/10 in most calculations), divided by a linear size L of the cell. This choice is dimensionally correct and corresponds to the appropriate time scale for strain fluctuations in the case of a thermodynamic system.

B. Interaction law

The contact law in a granular material is the relationship between the relative motion of two contacting grains and the contact force. As we deal with particles that may attract one another at short distance without touching, the law governing intergranular forces and moments is best referred to simply as the interaction law.

Although the interaction we adopted is based on the classical linear "spring-dashpot" model with Coulomb friction for contact elasticity, viscous dissipation and sliding, as used in many discrete simulations of granular media [13,37,38,46], some of its features (short-range attraction and rolling resistance) are less common ; moreover, one can think of different implementations of the Coulomb condition, depending on which parts of normal and tangential force components are taken into account. Therefore, for the sake of clarity and completeness, we give a full, self-contained presentation of the interaction law below.

We express intergranular forces in a mobile system of coordinates with axes oriented along the normal unit vector nij (along r ij ) and the tangential unit vector tij (n ij , tij is a direct base in the plane), and use the convention that repulsive forces are positive.

The intergranular force F ij , exerted by grain i onto its neighbor is split into its normal and tangential components, F ij = N ij nij + T ij tij thus defining scalars N ij and T ij . N ij comprises a static term depending on the distance between disk centers, combining contact elasticity and distant, van der Waals type attraction, as shown on Fig. 1a, and a velocity-dependent viscous term N v ij . T ij (Fig. 1b) is due to the tangential elasticity in the contact, and is limited by the Coulomb condition. If disks i and j are not in contact, both the tangential component of force F ij and the viscous part of the normal component vanish, while i and j still attract each other if the gap (h ij ≥ 0) between their surfaces is smaller than the attraction range D 0 (0 ≤ h ij ≤ D 0 ):

F ij = N a ij nij with N a ij = -F 0 (1 - h ij D 0 ) nij ( 4
)
This expression is a linear approximation of a realistic van der Waals force law (see Fig. 1a), and contains two essential parameters, maximum attractive force F 0 , range, D 0 . Typically, F 0 is of the order of γl, γ being a superficial energy, l the typical size of asperities [47] and D 0 is in the nanometer range. In the case of contacting disks (h ij < 0), the attractive term N a ij is kept constant, equal to -F 0 , while strains in the contact region result in normal (N e ij ) and tangential (T ij ) elastic forces. It is also assumed that a viscous normal term N v ij opposes relative normal displacements. One thus writes:

F ij = (N e ij + N v ij -F 0 ) nij + T ij tij (5)
The different terms introduced in Eqn. ( 5) are defined according to the following models. First,

N e ij = -K N h ij
is the linear elastic unilateral repulsion, due to the normal deflection -h ij in the contact as the disks are pressed against each other. K N is the normal stiffness coefficient, related to the elastic moduli of the material the grains are made of. The viscous normal force opposes the normal relative receding velocity δv N ij = nij • ( v j -v i ) as long as the contact persists. The relative normal motion of two disks i and j in contact is that of an oscillator with viscous damping, and η ij is the damping coefficient. We choose its value as a constant fraction ζ of the critical damping coefficient,

η ij = ζ 4K N m i m j m i + m j . ( 6
)
This is equivalent to the choice of a constant restitution coefficient in normal collisions if F 0 = 0. In the presence of attractive forces the apparent restitution coefficient in a collision will depend on the initial relative velocity, and will be equal to zero for small values, when the receding velocity after the collision will not be able to overcome the attraction and separate the particles. The minimum receding velocity for two particles of unit mass (i.e., of maximum diameter a) to separate is V * √ 2, with

V * = F 0 D 0 . (7)
The elastic tangential force in contact i, j is linearly related to the elastic part δu T ij of the total relative tangential displacement ∆u T ij , as

T ij = K T δu T ij ,
and is subject to the Coulomb inequality. K T is the tangential stiffness coefficient. ∆u T ij can be updated for all closed contacts according to

d∆u T ij dt = ( v ij • tij )
and vanishes as soon as the contact opens. Its elastic part satisfies

dδu T ij dt = H µN e ij K T -|δu T ij | ( v ij • tij )
in which H denotes the Heaviside function. This last equation introduces the friction coefficient µ. It is important to note that the Coulomb inequality,

|T ij | ≤ µN e ij , (8)
applies to the sole repulsive elastic component of the normal force (see Fig. 1b). We chose not to implement any tangential viscous force. The moment that disk i exerts onto its contacting neighbor j, of radius R j , in its center, is denoted by Γ ij in Eqn. [START_REF] Kishino | Powders and Grains 2001[END_REF]. It is first due to the tangential contact force, then to a possible moment Γ r ij of the force density distribution within the contact region. One thus writes:

Γ ij = -T ij R j + Γ r ij . (9)
Γ r ij is most often neglected on dealing with smooth, convex particle shapes, because the contact region is very small on the scale of the particle radius.

To model rolling resistance (RR), like in [48], we introduce a rotational stiffness parameter K r and a rotational friction parameter µ r in contacts, so that rolling elasticity and rolling friction are modelled just like sliding elasticity and friction. One thus writes

Γ r ij = K r δθ ij ,
while enforcing the inequality

K r |δθ ij | ≤ µ r N e ij . (10)
This involves the definition of δθ ij as the elastic part of the total relative rotation ∆θ ij . The total relative rotation angle satisfies

d∆θ ij dt = ω j -ω i ,
while the equation for δθ ij is

dδθ ij dt = H µ r N e ij K r -|δθ ij | (ω j -ω i ).
Parameters K r and µ r are often related to the size of a contact region [27]. K r is dimensionally the product of a stiffness by the square of a length, which is of the order of the contact size. In the following we set K r to 10 -4 a 2 K N , while µ r , which has the dimension of a length, is chosen equal to 10 -2 µa. The motivation for the introduction of RR into our model is twofold. First, cohesive particles are usually small (typically less than 30µm in size) and irregular in shape. Contacts between grains are likely to involve several asperities, and hence some lateral extension, of the order of the distance between asperities, however small the normal deflection -h. Then, it will be observed that even quite a small rotational friction has a notable influence on the microstructure of cohesive packings.

C. Control parameters and dimensional analysis

In this section we present the dimensionless parameters which express the relative importance of different physical phenomena. Such parameters enable qualitative comparisons with real materials, bearing in mind that the present model is admittedly an idealization of real powders and that our simulations do not aim at quantitative accuracy.

Dimensionless numbers related to contact behavior are the reduced interaction range D 0 /a, the friction coefficient µ, the viscous damping parameter ζ, and the stiffness parameter κ.

Under the attractive force -F 0 , the elastic deflection of one contact is

h 0 = F 0 /K N (11)
The stiffness parameter κ def = aK N /F 0 characterizes the amount of elastic deflection h 0 under contact force F 0 , relative to grain size a (h 0 /a = κ -1 ). A suitable analogous definition for Hertzian spheres in three dimensions would be κ = (Ea 2 /F 0 ) 2/3 .

The dimensionless number h 0 /D 0 is the ratio of elastic to adhesive stiffnesses, and its physical meaning is similar to that of the Tabor parameter λ = (1/D 0 )(γ 2 a/E 2 ) 1/3 [START_REF] Maugis | Contact, adhesion and rupture of elastic solids[END_REF] for a Hertzian contact between spheres of diameter a when the material Young modulus is E and the interfacial energy is γ (more precisely, the equilibrium normal deflection h 0 , due to adhesion, in the contact between an isolated pair of grains, satisfies λ ∼ (h 0 /D 0 ) 1/3 in this case).

The viscous damping parameter, ζ, corresponds to a normal restitution coefficient e N = exp[-πζ/ 1 -ζ 2 ] in the absence of cohesion (F 0 = 0).

In our calculations we set ζ = 0.8, corresponding to a high viscous dissipation in collisions, or a very low restitution coefficient in binary collisions. Models with a constant ζ were adopted in other published simulation works [46,[START_REF] Somfai | [END_REF], although little is known about dissipation in collisions. ζ is known to influence the packing structures obtained in the initial assembling stage [46,51], but we did not investigate its effects in the present study. The simulations reported in [25,26] use the viscous force model introduced in ref. [52], with a choice of parameters corresponding to strongly overdamped dynamics (i.e., analogous to ζ ≫ 1 in our case).

In addition to those control parameters determined by the contact behavior, other dimensionless numbers are introduced by the loading or the process being applied to the material. The effect of the external pressure, compared to the adhesion strength, is characterized by a dimensionless reduced pressure P * :

P * def = P a/F 0 . (12)
In the present paper, we focus on the assembling process and the low P * range. As we shall see below (Sec. III) low density, tenuous structures are then stabilized by adhesion, and the relevant force scale is F 0 . However, as briefly reported in [53], such structures tend to collapse upon increasing P * . These phenomena will be the subject of another paper [START_REF] Gilabert | Computer simulation of model cohesive powders: plastic compaction under isotropic loads[END_REF]. Wolf et al. [29] introduced a dimensionless stress proportional to P * , and observed, in numerical simulations, stepwise increases in pressure to produce large dynamical collapse effects around P * = 1. The importance of P * was also stressed in simulations of cohesive granular flow, in which the effects of cohesion on rheological laws were expressed in terms of a cohesion number defined as 1/P * [37]. In three dimensions, P * should be defined as a 2 P/F 0 . For large reduced pressures, externally imposed forces dominate the adhesion strength, and one should observe behaviors similar to those of confined cohesionless granular materials. For P * > 1, the relevant force scale is aP . The influence of κ, which should then be defined as κ def = K N /P , so that the typical contact deflection h satisfies h/a ∝ κ -1 , was studied in simulations of grains without adhesion [START_REF] Combe | Deformation characteristics of geomaterials[END_REF]. Whatever the reference force used to define it, the limit of rigid grains is κ → +∞. With relatively soft grains (say, κ below 10 3 ), a significant number of additional contacts appear in dense configurations, due to the closing of gaps between near neighbors. Such a κ parameter defined with reference to pressure, in the case of contacts ruled by Hertz's law between spherical grains made of a material with Young modulus E, should be chosen as κ = (E/P ) 2/3 , in order to maintain h/a ∼ κ -1 .

In order to stay within the limit of rigid grains both for small and large P * , we choose quite a large value of κ = K N a/F 0 : κ = 10 4 or κ = 10 5 .

Table I KN a F0 and D0/a to 10 -2 or 10 -1 . In the absence of cohesion, or for values of P ≥ F0/a, κ is defined as KN /P eters, adopted as a plausible choice for realistic orders of magnitudes, some calculations were also performed with deliberately extreme choices, such as very large RR (µ r = 0.5a) or absence of friction (µ = 0 and µ r = 0), in order to better explore some connections between micromechanics and macroscopic properties. The corresponding results will be described in Section IV.

The definition of dimensionless parameters, suitably generalized to three-dimensional situations as P * = a 2 P/F 0 and κ ≃ (Ea 2 /F 0 ) 2/3 for spherical particles of diameter a, enables one to discuss qualitative features and orders of magnitude in the model system defined with the parameters of Table I with comparisons to some cohesive packings studied in the laboratory.

When adhesive forces are due to liquid menisci joining neighboring particles, we should take F 0 ∼ γa, where γ is the surface tension. P * = 1 corresponds then to confining pressure P in the range of 10-100 Pa for millimeter-sized particles, taking standard values for γ. Those are rather low pressures in practice, which are comparable, e.g., to the ones caused by the weight of a typical laboratory sand sample. Thus wet granular materials are commonly under reduced pressures P * of order 1 or larger, and are not observed with much lower solid fractions than dry ones [6,9,10,11].

The cohesive powders studied in refs. [START_REF] Castellanos | [END_REF][START_REF] Watson | [END_REF]56,57] are xerographic toners with typical particle diameter a ∼ 10 µm. F 0 , the van der Waals attractive force, is a few tens of nN, and the range D 0 is several nanometers [58]. Therefore, a reduced pressure P * = 0.01 would correspond to about 1 Pa in the experimental situation [59]. This is an initial state of very low consolidation stress, which is present in a powder under gravity, provided a controlled gas flow, going upwards through the powder, counterbalances part of its weight [59]. As to contact stiffnesses, our values of h 0 /a would correspond to E ∼ 0.1 GPa (for K N = 10 4 F 0 /a) or 3.2 GPa (for K N = 10 5 F 0 /a), while the ratio D 0 /a would imply an interaction range of 10 nm. This gives us the correct orders of magnitudes for the toner particles, those being made of a relatively soft solid (polymer, such as polystyrene) with E ∼ 3-6 GPa. Xerographic toner particles appear to undergo plastic deformation in the contacts [56,60,61,62]. Plastic deflections of contacts are accounted for in the model of ref. [35], applied to the simulation of a biaxial compression of a dense powder. In our study, for simplicity's sake, and because we expect macroscopic plasticity of loose samples to be essentially related to the collapse of tenuous structures, we ignored this feature.

D. Equilibrated states

Although numerical simulations of the quasistatic response of granular materials requires by definition that configurations of mechanical equilibrium should be reached, equilibrium criteria are sometimes left unspecified, or quite vaguely stated in the literature. Yet, in order to report results on important, often studied quantities like the coordination number or the force distribution, it is essential to know which pairs of grains are in contact and which are not. Due to the frequent occurrence of small contact force values, this requires forces to balance with sufficient accuracy. We found that the following criteria allowed us to identify the force-carrying structure clearly enough. We use the typical intergranular force value F 1 = max(F 0 , P a) to set the tolerance levels. A configuration is deemed equilibrated when the following conditions are fulfilled:

• the net force on each disk is less than 10 -4 F 1 , and the total moment is lower than 10 -4 F 1 a;

• the difference between imposed and measured pressure is less than 10 -4 F 1 /a;

• the kinetic energy per grain is less than 5 •10 -8 F 1 a.

We observed that once samples were equilibrated according to those criteria, then the Coulomb criterion (8), as well as the rolling friction condition (10) were satisfied as strict inequalities in all contacts. No contact is ready to yield in sliding, and with RR no contact is ready to yield in rolling either.

III. ASSEMBLING PROCEDURE

It has been noted in experiments [START_REF] Mitchell | Fundamentals of soil behavior[END_REF] and simulations [16,[START_REF] Agnolin | [END_REF]46] that the internal structure and resulting behavior of solid-like granular materials is sensitive to the sample preparation procedure, even in the cohesionless case.

In the case of powders, it has been observed that the sedimentation in dry nitrogen (to minimize the capillary effects of the humidity on the interparticle adhesion) of a previously fluidized bed produced reproducible states of low solid fractions (down to 10 -15%) [63,64]. This initial state under such a low consolidation, as we commented in II C, plays a decisive role on the evolution of the dynamics of powder packing. That is, appreciable differences in initial states will lead to considerable ones in final packings [57]. This is mainly due to the role of aggregation, which we shall analyze in the second part of this section.

The motivation of this section is to investigate the dependence on packing procedure in a cohesive granular system, the first step being to obtain stable equilibrated configurations with low densities. For comparison, some simulation results are presented for the same model material with no cohesion.

Specimens were prepared in two different ways, respectively denoted as method 1 and method 2, and the resulting states are classified as type 1 or type 2 configurations accordingly.

Due to our choice of boundary conditions, our samples will be completely homogeneous, under a uniform (isotropic) state of stress. This choice is justified by the complexity of seemingly more "realistic" processes, such as gravity deposition, due to the influence of many material (such as viscous dissipation, as recalled in Section II C) and process parameters. Both pouring rate and height of free fall should be kept constant during such a pluviation process in order to obtain a homogeneous packing [51,65] with cohesionless grains. Cohesive ones, because of the irreversible compaction they undergo on increasing the pressure, would end up with a density increasing with depth. Hence our choice to ignore gravity in our simulations. Our final configurations should be regarded as representative of the local state of a larger system, corresponding to a local value of the confining stress.

A. Method 1

In simulations of cohesionless granular materials, a common procedure [16,17,66] to prepare solid samples consists in compressing an initially loose configuration (a "granular gas"), without intergranular contacts, until a state of mechanical equilibrium is reached in which interparticle forces balance the external pressure (further compaction being prevented by the jamming of the particle assembly). We first adopted this traditional method, hereafter referred to as method 1, to assemble cohesive particles.

In this procedure, disks are initially placed in random non-overlapping positions in the cell, with zero velocity.

We denote such an initial situation as the I-state. Then the external pressure is applied, causing the cell to shrink homogeneously. Thus contacts gradually appear and the configuration rearranges until the system equilibrates at a higher density.

Examples of equilibrated configurations are shown on Fig. 2, with and without cohesion. This state is characterized by its solid fraction (Φ = A -1 i πd 2 i /4) and its coordination number z, defined as the average number of interactions (contacts and distant attractions) for a particle in the packing, when the applied pressure is significantly smaller than F 0 /a (P * ≪ 1), as in the case of small powder samples assembled under gravity. With the values indicated above (at the end of Section II C) for toner particles, F 0 /a 2 (the relevant pressure scale in 3D) is of the order of 100 Pa, which corresponds to a normal consolidation stress in a cohesive powder with 34% solid fraction [59].

In the absence of cohesion, the value of the applied pressure does not affect the properties of the packing (apart from setting the scale of intergranular forces) provided the typical contact deflection, aP/K N , is small enough (rigid particle limit). We set this ratio to the value of F 0 /(aK N ) in the cohesive case, i.e., equal to κ -1 (see table I), so that typical contact forces are of the same order of magnitude (due either to P or predominantly to F 0 ) in both cases.

Effects of the initial solid fraction in the I-state, and of cohesion, friction and rolling resistance parameters on Φ and z were measured in 3 sets of samples, with Φ I = 0.13, 0.36 and 0.52. Each set consisted of configurations with the initial disorder (particle radii and initial positions) abiding by the same probability distribution, and the same number of particles (N = 1400). The values of the friction coefficient, µ, used in these tests were 0.15 and 0.5. The values of Φ and z in these samples are listed in tables II and III. Each one is an average on the different samples, and the indicated uncertainly is equal to the standard deviation.

Tables II andIII show that the introduction of cohesion reduces the solid fraction at equilibrium, but this is a limited effect (less than 10% density reduction), which is quite insufficient to account for experimental observations. Unlike powders or clays, 2D particle packings with Φ ≥ 0.7 cannot undergo very large plastic density increases.

Theses tables also show that the increase of the friction coefficient and/or the inclusion of rolling resistance in the model tend to hinder motions and stabilize looser, less coordinated configurations, which results in a decrease of Φ and z.

However, the observed differences are rather small, especially in cohesionless systems. To evaluate the influence of RR with a given value of µ, we define ∆X

(µ) = 1 -X (µ) RR /X (µ)
no RR , as the relative average decrease of the quantity X due to the existence of RR. For example, for µ = 0.15 results differ by a mere ∆Φ (0.15) = 0.24% and ∆z (0.15) = 1.3%, and for µ = 0.5, these variations are with the usual convention that the width of the lines joining the centers of interacting pairs of disks is proportional to the normal force, on scale aP (left) and F0 (right). Red, green, and blue lines distinguish compressive, tensile, and distant interactions in the cohesive case, while rattlers appear in grey in the cohesionless sample.

Non-Cohesive samples no RR RR ΦI µ = 0.15 µ = 0.5 µ = 0.15 µ = 0.5 Solid fraction 0.130 ± 0.001 0.8262 ± 0.0007 0.811 ± 0.001 0.8238 ± 0.0014 0.803 ± 0.002 0.3631 ± 0.0006 0.8256 ± 0.0005 0.811 ± 0.001 0.8231 ± 0.0013 0.805 ± 0.002 0.5244 ± 0.0012 0.8236 ± 0.0007 0.8092 ± 0.0005 0.8215 ± 0.0005 0. ∆Φ (0.5) = 0.84% and ∆z (0.5) = 4.6%. Comparing the effect of RR on Φ with µ = 0.5 for µ = 0.15 and µ = 0.5, one has ∆Φ (0.5) /∆Φ (0.15) = 3.5. Likewise, for coordination numbers z, one observes ∆z (0.5) /∆z (0.15) ≃ 3.53. This shows a clear correlation of variations introduced by friction and RR. The data in the non-cohesive case also exhibit very little dependence on initial density Φ I .

Results on cohesive systems show similar variations with the parameters of the contact model (friction and RR), but depend somewhat more sensitively on Φ I .

More refined information on the contact network is provided by the distribution of local coordination numbers, i.e. the proportions x k of particles interacting with k neighbors, for 0 ≤ l ≤ 6 (higher values were not observed). This distribution is depicted on Figure 3, for both cohesive and non-cohesive samples. These results gather information from all the statistically equivalent simulated samples, and slight corrections were applied in order to ignore the contacts with "rattler" particles in the non-cohesive case. Such particles are those that Cohesive samples no RR RR ΦI µ = 0.15 µ = 0.5 µ = 0.15 µ = 0.5 Solid fraction 0.130 ± 0.001 0.7635 ± 0.0023 0.751 ± 0.001 0.757 ± 0.002 0.709 ± 0.001 0.3631 ± 0.0006 0.727 ± 0.001 0.7232 ± 0.0012 0.710 ± 0.002 0.688 ± 0.001 0.5244 ± 0.0012 0.737 ± 0.002 0.733 ± 0.002 0.7248 ± 0.0002 0.733 ± 0.002 Coordination number 0.130 ± 0.001 3.563 ± 0.005 3.163 ± 0.004 3.189 ± 0.008 3.059 ± 0.003 0.3631 ± 0.0006 3.345 ± 0.009 3.103 ± 0.006 3.253 ± 0.003 2.971 ± 0.006 0.5244 ± 0.0012 3.189 ± 0.008 3.059 ± 0.003 3.096 ± 0.002 2.851 ± 0.001 TABLE III: Solid fractions and coordination numbers obtained at the preparation of the specimens in equilibrated samples under P * = 0.01, F0/(KN a) = 10 -5 for cohesive particles, using method 1. are free to move within the cage of their near neighbors, and transmit no force once the system is equilibrated. If they happen to be in contact with the backbone (i.e., the force-carrying structure), then the forces carried by such contacts should be below the tolerance set on the equilibrium requirement, and can safely be ignored. This is how the population of rattlers is identified. We observe that it can involve up to 18% of the total number of grains in the absence of cohesion (see Fig 2(a)).

This contrasts with the cohesive case, for which nearly all the grains are captured by the force-carrying structure because of attractive forces, and the rattlers are virtually absent. The particles with one contact equilibrate when the deflection of that contact is h 0 , as defined in (11). With RR, such a particle is entirely fixed. Without RR, it is only free to roll without sliding on the perimeter of its interacting partner, because such a contact is able to transmit a tangential force smaller than or equal to µK N h 0 = µF 0 .

Without cohesion, the coordination of the forcecarrying structure can be characterized with a coordi-nation number z * , different from z:

z * = z 1 -x 0 , ( 13
)
z * is the average number of contacts bearing nonnegligible forces per particle on the backbone. Without cohesion, the backbone (or set of non-rattler grains) is the rigid part of the packing. With cohesion and RR, the whole interacting contact network is to be considered in order to study the rigidity properties of the system, and there are nearly no particles to eliminate. With cohesion and no RR, we observe in the samples obtained by the presently employed procedure (method 1) that the network of interparticle contacts or interactions is also rigid, apart from the free rolling of isolated grains with only one contact. (The rigidity properties of equilibrated samples are discussed below in Section IV and Appendix A).

Cohesive samples in equilibrium also comprise a small number of pairs of particles interacting without contact, i.e. separated by a gap smaller than the range of attraction, D 0 . These are only a small fraction, below 1%, of interacting pairs. Such pairs do not contribute to dissipation, since the frictional and viscous force components are only present in true contacts between neighboring grains. We observed that the time necessary to equilibrate the sample tend to increase when such distant interacting pairs are more numerous.

In addition to the elimination of free rattlers, the most notable effect of cohesion on local coordination numbers (Fig. 3) is to increase the proportion of disks with 2 contacts. Without cohesion, the Coulomb condition restricts the angle between the directions of the 2 contacts to values between π -2ϕ and π, where ϕ is the friction angle (tan ϕ = µ). Thus if µ is small, a disk with two contacts should have its center close to the line of centers of its two partners. The increase of the population of 2-coordinated disks as µ is raised from 0.15 to 0.5 (see Fig. 3) in cohesionless systems corresponds to a less severe geometric restriction on contact angles. With cohesion, contacts may transmit a tangential force reaching µF 0 while the normal force is equal to zero. Consequently, a disk might be in equilibrium with two contact points in arbitrary positions on its perimeter. As there is no geometric constraint on the angle between the two contact directions, 2-coordinated disks are easier to stabilize, and their proportion raises from about 5% without cohesion to above 15% with cohesion in the case µ = 0.15. A population of disks with one contact (therefore carrying a vanishing normal force, with deflection -h = h 0 ) is also present. Those particles are fixed by a small rolling resistance, but are free to roll on their interacting neighbor without RR. Such a rolling motion is not damped in our model. Therefore, on waiting long enough, they should eventually stop after a collision, in a stable position with 2 contacts. Such a collision is bound to happen because the contact network is completely connected. However, we stop our calculations when the kinetic energy is below a set tolerance (see Section II D), and we do not wait until all freely rolling disks reach their final position. Hence the remaining population of disks with one contact in samples without RR.

The final configuration, with this preparation method, depends somewhat on the rate of compaction in the assembling stage. The latter is related to the choice of the dynamical parameter M , the "mass" with which the changes in cell dimensions are computed with Eqn. [START_REF]The Physics of Granular Media[END_REF]. The slight influence of the initial solid fraction, Φ I , also relates to such dynamical effects: a lower value of Φ I entails larger colliding velocities, which favors larger final solid fractions.

Although some of the aspects of the model (in particular the homogeneous shrinking imposed through the periodic cell dimensions in a dynamical regime) do not correspond to experimental conditions, configurations of type 1 should be regarded as typical results of fast assembling processes, in which the particles are requested to balance the external pressure before stable loose structures can be built. When the toner particles mentioned at the end of Section II C are first fluidized, and then settle under their own weight, a rough estimate of the settling time, assuming particles are settling individually in air, and fall over distances of order 1 cm, is ∼ 1 s. Fig. 4, with the value T 0 ∼ 10 -5 s corresponding to such particles, shows that the duration of the "method 1" compression process is a few milliseconds. In practice, due to the presence of the surrounding fluid, the packing of a powder in a loose state by settling and compaction of an initially fluidized state is therefore considerably slower than this numerical process.

In the next section, we consequently turn to the opposite limit, in which the external confining pressure is felt only after large, tenuous contact networks are formed.

B. Method 2 1. Numerical procedure.
The second method to prepare numerical samples allows for aggregate formation before imposing an external pressure. Along with method 1, its different stages are schematically presented on Fig. 4.

The aggregation phenomenon plays an important role in the experimental preparation procedure of refs. [67,68] in which powder particles in a fluidized bed collide and stick to each other. Then they settle under their weight when the upwards air flow is abruptly shut off. The numerical method was designed to reproduce, in some idealized way, the final state of a set of colliding particles in the absence of external force fields. In the initial disordered low-density configuration (the same I-state as in method 1), particles are now attributed random velocities drawn according to a Maxwell distribution, with mean quadratic velocity V 0 .

We performed systematic sets of simulations of disk packings with V 0 = 9.48V * (see Eqn. (7)). V 0 is thus large enough for the initial kinetic energy to overcome potential energy barriers in the process of aggregation. (The dependence of the final packing structure on this initial velocity of agitation, or "granular temperature", in the assembling stage in systems with small RR will be studied in Section IV B 6).

Once launched with such random velocities the particles are left to interact and stick to one another within a cell of constant size, forming larger and larger aggregates, as appears on the image marked "B2" on Fig. 4. Eventually, all particles are connected to one another by adhesive contacts, and reach an equilibrium position. At this stage, the two degrees of freedom of the cell are set free, and the stress-controlled calculation proceeds with σ 11 = σ 22 = 0 (or P * = 0) until an equilibrium state is reached. This relaxation step does not lead to any rearrangement of the contact structure, it only entails a very small increase of the solid fraction (hence the values slightly larger than Φ I given below). The final equilibrium structure exhibits large density inhomogeneities, as apparent on Fig. 4, which are characteristic of aggregation processes [39], and will be quantitatively studied in Section V. Point C2 corresponds to the stage when all disks are assembled in a unique aggregate, then equilibrated at P * = 0 (both aggregation and equilibration stages take place between A and C2). The time unit is T0 = ma/F0. Note the duration of the preparation process with method 2, and the difference in final equilibrated states compared to method 1.

Unlike cohesionless systems, which are devoid of any "natural" state of stress, clusters of cohesive particles can exist in a well defined state of mechanical equilibrium in the absence of any external force. Once the state at zero pressure is obtained, we subsequently apply the same load P * = 0.01 as in method 1, which results in further compression and notable changes in the packing structure: Φ increases from values close to Φ I up to the 0.45-0.55 range (see fig. 4). Nevertheless, the final solid fraction under P * = 0.01 is considerably lower than the one obtained with method 1.

It should be noted on Fig. 4, which summarizes the assembling procedures, that the aggregation stage makes method 2 computationally quite costly because of the time necessary for clusters to merge, and especially for the stabilization of loose samples in equilibrium config-urations (lower contact numbers implying lower rates of energy loss as well as larger and slower fluctuations of soft, tenuous structures). In an attempt to limit the influence of compaction dynamics, which results in denser samples when the lower density of the initial state allows the compaction process to accelerate more (as noted in Sec. III A), we tested the effect of limiting the maximum strain rate ǫmax . Without any limitation, we obtained a maximum value ǫ ≃ 0. gave a final state close to the original one. The others two values produced similar results, with a relative decrease in density of about 10% compared to the original procedure. We chose to enforce condition ǫmax = 0.05 T -1 0 , to save computational time. This value has been applied to prepare all samples studied in the following.

Fig. 4 shows that method 2 succeeds in stabilizing open structures. Final solid fractions agree with the typical values observed in powders if one uses the correspondence between 2D and 3D packing fractions suggested by Campbell in [69]:

Φ 3D = 4 3 √ π Φ 3/2 2D ≃ 0.752Φ 3/2 2D . (14)
Numerical samples under P * = 0.01, with solid fractions around 45%, would correspond to a powder consolidated in the laboratory under 1 Pa with a solid fraction of about 23%. This is in satisfactory agreement with the experimental results of Ref. [59].

We therefore regard method 2 as an appropriate way to reach an essential objective of this work, since stable loose structures are obtained.

Although we perform simulations of a mechanical model, the final configurations exhibit at first sight (Fig. 4) similar features as those obtained with geometric algorithms implemented in numerical studies of colloid aggregation models [28,40]. We are not aware of similar results in the literature, at least with equilibrium requirements comparable to those of Sec. II D.

Tenuous, fractal-like contact networks contain denser regions and large cavities. Such heterogeneities produce long range density correlations, to be analyzed in Sec. V. Without tensile contact forces, the walls of the cavities, comprising particles that are pushed towards the hole by the resultant of contact normal forces, would tend to buckle in.

We regard method 2 as yielding typical results for assembling processes in which particles form tenuous aggregates before they are packed in a structure that is able to sustain a confining stress. In the sequel, we focus on the tenuous structures obtained with method 2.

2. Global characterization of loose packings at P * = 0 and P * = 0.01.

We simulated four samples with 1400 disks and three of 5600 for Φ I = 0.36, rather than lower initial densities, in order to achieve statistical significance at affordable computational costs, and to check for possible size effects. This set of samples will be denoted as series A.

Samples with Φ = 0.13 (series A0), which require the initial cell to shrink more before a stable network can resist the pressure, request longer calculations. Although some samples were prepared at P * = 0, we do not use them any more in the following, except for the values showed in Table IV.

To accelerate the numerical assembling procedure, we also created samples with K N = 10 4 F 0 /a, using an intermediate value of Φ I = 0.26, and softer contacts, such that κ = 10 2 in the initial aggregation stage (recall the time step is proportional to A/ √ K N ). Once equilibrium was reached with P * = 0, we slowly changed the stiffness parameter from κ = 10 2 to κ = 10 4 , and recorded the final equilibrated configuration. This procedure is about ten times as fast as the normal one, and generates similar structures and coordination numbers as series A prepared with the same I-state density. We shall refer to this set as series B.

In table IV we list the corresponding results for solid fractions and coordination numbers. In such data we did not find a significant difference between the two different sample sizes, and therefore we did not distinguish between sizes in the presentation of statistical results. The tenuous networks obtained with method 2 collapse on changing the pressure: Table V gives the new values of Φ and z after the compaction caused by the pressure increase from P * = 0 to P * = 0.01.

Structural changes between P * = 0 and P * = 0.01 are shown on Fig. 5, which illustrates by means of four selected snapshots the mechanism of the closing of pores in a 1400 disks sample of series A. The first image corresponds to equilibrium at P * = 0, and the fourth one to equilibrium under P * = 0.01. The two others show intermediate, out of equilibrium configurations during the collapse. One may appreciate how the denser regions grow and merge while pores shrink. Fig. 5 also makes it quite evident that the size of 1400 disk samples is not very much larger than the scale ξ of density heterogeneities (typical diameter of large pores or dense regions, which will be studied in Sec. V). These systems will exhibit large fluctuations in their mechanical properties: the rectangular shape of the final configuration displayed on Fig. 5 shows that the disorder is large enough for the mechanical response of the system to become anisotropic. Isotropy should be recovered in the limit of large sample sizes, L ≫ ξ.

Finally, Fig. 6 displays the histogram of local coordination numbers (percentage of particles interacting with k others, 0 ≤ k ≤ 6), for the same samples as those of Tables IV and V (µ = 0.5, Φ I = 0.36). It is remarkable that this distribution, in spite of the large difference in sample geometries, remains rather close to the one observed in the denser packings made with method 1 (compare P * = 0.01, µ = 0.5 case), just like global coordination numbers take very similar values in samples prepared with both methods (see Tables III andV), in spite of the very different solid fractions.

An essential conclusion of the present study is therefore, for one given material, the absence of a general relation between the density of a cohesive packing and its coordination number, in spite of previous claims [26]. Both quantities are determined, rather, by the conjunction of micromechanical laws and sample preparation history.

Effects of micromechanical parameters

Adhesion should enhance the role of sliding friction and rolling friction, because the limiting values for tangential contact forces and rolling moments are both proportional to the elastic repulsive part of the normal force, N e (

|T ij | < µN e ij , | Γ r ij | < µ r N e ij )
. Consequently, contacts with the equilibrium value h 0 of the elastic deflection for an isolated pair of grains transmit no normal force, but are able to sustain tangential force components as large as µF 0 and rolling moments as large as µ r F 0 . Those values might turn out to be large in comparison to the typical level of intergranular forces under low external pressure (P * ≪ 1). Therefore, even very low values of µ and µ r should affect the final structure of equilibrated packings considerably more than in the cohesionless case. This is indeed the case for the coordination numbers observed in our simulations (see tables II and III) which dropped more significantly, upon introducing the small level of RR we have been using, in cohesive systems than in cohesionless ones.

On Fig. 7 we show the configurations at P * = 0 (a) and P * = 0.01 (b) of the same sample assembled using method 2 with RR (left) and without RR (right). The denser regions in the inhomogeneous packings are joined by slender "arms" (see Figs. 5 and7). Such arms can in principle reduce to a chain of particles in the presence of rolling resistance. Such chains are otherwise destabilized by a rolling mechanism, hence the difference in the thickness of the arms with or without RR (see the blown-up detail in Fig. 7-a), the lower coordination numbers of configurations assembled with RR (Tables IV andV). This might also explain the greater fragility of equilibrium configurations with RR, in which a larger compaction step (see Table V) is necessary, on applying P * = 0.01, before a new stable structure is reached.

Another important parameter is the initial velocity of agitation, V 0 . Its influence has been assessed on one 1400 disks sample, with Φ I = 0.36. The changes of coordination number with V 0 at P * = 0 are presented on Fig. 8.

Low velocity values produce more tenuous aggregates (z ∼ 2), since even a small level of RR is able to slow down local rearrangements and stabilize tree-like structures (i.e., devoid of flops) immediately after the collisions between particles or small clusters.

A large kinetic energy cannot be absorbed by the RR, and as a result disks are able to rotate, which leads to better connected structures (z ∼ 3). In a sense, a large V 0 kills the effects of RR, and packings are similar to those made without RR in such cases.

We therefore conclude that the connectivity of loose samples with RR assembled by aggregation depends on the initial magnitude of velocity fluctuations and on the level of rolling friction.

As figure 8 shows, the same trend was found on reducing contact stiffness parameter κ, as a larger translational and rotational compliance creates more contacts.

V 0 is analogous to the particle fluctuating velocity in experiments on gas-fluidized beds of xerographic toners under gravity [67]. Such velocities are larger than the gas velocity by two orders of magnitude. Typically, one has v gas ∼ 1 -4 mm/s, while V * , deduced from the contact parameters with relation ( 7) is about 1 cm/s. Such a value is therefore comparable to the particle fluctuation velocity.

Of course, such a comparison is only indicative, because the influence of V 0 on packing structures depends on µ r , and is also very likely to be affected to some extent by the viscous dissipation model we have adopted. Both rolling resistance and viscous forces are micromechanical features for which no accurate physical identification is available. Yet, it seems plausible that powder packings, because of their initial agitated states, stabilize in better connected states than predicted by geometric aggregation models.

We now turn our attention in the next section to the forces in the contact networks, in particular the loose ones formed with method 2.

IV. MECHANICAL CHARACTERIZATION OF CONTACT NETWORKS

Many numerical studies, in the past 15 years, have addressed the issue of contact network geometry and force distribution in cohesionless systems [70]. The image of force chains, i.e. a pattern in which larger intergranular forces tend to line up on the scale of several grains, was evidenced in experiments [71,72] and simulations [73,74], and the p.d.f of contact force values has often been measured and studied. An interpretation of the mechanical role of "force chains" [75] is that they carry the essential part of deviatoric stress, while the contacts carrying the lower forces are less sensitive to stress orientation and laterally stabilize the strong force chains against buckling.

The main features of the distribution of forces and their spatial correlations have been reproduced by approximate models [76] based on local equilibrium rules on each grain, supplemented by inequality constraints. One important such constraint is released in cohesive systems, in which normal force components can have either sign. It is therefore worth investigating how the usual features of force-carrying structures in equilibrated granular packings are affected by the presence of negative normal forces. One may also wonder to what extent the considerable difference in the density fields will affect the force patterns, given that the coordination of the force networks, as observed previously, does not seem to be very sensitive to density levels and density fluctuations.

A. Force scale and force distribution

The first obvious distinction between cohesive and noncohesive systems is the appearance of a new force scale F 0 , in addition to the one provided by the confining pressure, i.e., aP , the ratio of those characteristic forces defining the reduced pressure, P * . It is especially interesting to investigate the values and spatial organization of forces in systems with P * ≪ 1, as little information is to be found in the literature on this issue: numerical studies of loose cohesive systems [29] tend to focus on density and geometry of packings as a function of applied stresses. Some information on force networks is provided in a recent publication [11] on bead assemblies with capillary cohesion, but the confining stress is considerably higher is that study (P * of order 1) than in the present one.

In the absence of cohesion, the distribution of force values is usually presented in a form normalized by its average, which itself scales with the applied pressure. This scaling can be made more quantitative on using a general relation between pressure P and the average normal contact force F N def = N ij and particle diameter d, which is known in the literature on powders as the Rumpf formula. We write it here in a form involving the spatial dimension D, which is valid both for D = 2 and D = 3:

P = 1 π zΦ d D-1 F N . (15)
In ( 15), d stands for the typical grain diameter. This relation can be made more accurate if one notes that it stems from the standard formula for stresses in an equilibrium configuration (see the r.h.s. term in Eqn. 3). To derive the formula, defining P = 1 D D α=1 σ αα , the average pressure, one assumes h ij ≪ R i + R j and then neglects correlations between particle radii and forces, assuming

N ij (R i + R j ) ≃ F N d . (16)
Then, with a simple transformation of the sum, one obtains

P = 1 π d d D zΦF N . (17)
With D = 2 and our diameter distribution (for which d = 3a/4 and d 2 = 7a 2 /12) this yields

F N = 7πa 9 P zΦ . (18)
We found relation [START_REF] Wood | Soil Behaviour and Critical State Soil Mechanics[END_REF] to be remarkably accurate in all our simulations, with or without cohesion, with configurations obtained by either method 1 or method 2, thereby checking that the correlations between particle sizes and contact forces could safely be neglected on writing (16). Without cohesion, Eqn. ( 18) yields the correct scale for forces, i.e. the frequency of occurrence of intergranular forces larger than a few times F N is very small. With cohesion, when P * = 0 or P * ≪ 1, contact forces of order F 0 are quite common, as shown on Fig. 9, on which normal force distributions are represented. Hence Eqn. ( 18) cannot be used to predict "typical" contact forces. The presence of forces of order F 0 explains the sensitivity of type 1 and type 2 samples with P * ≪ 1 to friction coefficient and rolling resistance: densities and coordination numbers (tables III, IV and V), in cohesive systems prepared under P * = 0 or P * = 0.01 with µ = 0.15 and with µ = 0.5, or with and without RR, differ significantly. Otherwise, if contact forces were of order of the average F N , the value of which is correctly predicted by [START_REF] Wood | Soil Behaviour and Critical State Soil Mechanics[END_REF], thresholds µF 0 or even µ r F 0 would be very large compared to typical forces and moments, and become irrelevant.

(It should be recalled that Rumpf's name is often associated (as in Ref. [11]) to a means to predict the macroscopic tensile strength of a powder. As the essential ingredient of the Rumpf approach [77] is Eqn. 15, we refer here to that relation (like in [START_REF] Castellanos | [END_REF]), as the Rumpf formula).

Normal force distributions in cohesionless, cohesive type 1 and cohesive type 2 samples, the latter being obtained with Φ I = 0.36 (series A), are shown on Fig. 9. Those distribution functions are roughly symmetric about 0, decay approximately exponentially at intermediate values, and vanish at -F 0 , and F 0 . In type 2 samples without RR, for P * = 0, there is a finite proportion of contacts carrying vanishing forces, about one fourth in the A series (Φ = 0.36). In addition to this Dirac mass, there might be a power-law divergence near 0, with an exponent our level of statistics is not sufficient to resolve accurately (about 0.6 to 0.8 in the range of forces between 10 -3 F 0 and 10 -2 F 0 ). This proportion of zero forces is smaller, down to 9%, with RR, and drops as P * reaches 0.01, to 7% and 3%, respectively, without and with RR. It is worth pointing out that the corresponding contacts carry zero total forces, i.e. vanishing normal components (-h = h 0 , see (11)) and no tangential elastic displacement either. In principle we cannot distinguish them from forces below the numerical tolerance defined in Sec. II D. However, as we shall argue below in Section IV B, under P * = 0 one could expect all contact forces to vanish, and non-zero forces are related to the small, but finite degree of force indeterminacy.

Before turning our attention to such features and to the spatial organization of forces, let us briefly discuss the differences between sets of (type 2) samples A and B. B samples, which are obtained with the "accelerated" procedure and Φ I = 0.26 (see Sec. III), exhibit, due to their specific history, larger forces at P * = 0, with as many as 10% of the contacts transmitting normal forces N such that |N | > F 0 /10, while this proportion lies below 2% in A samples. On the other hand, B samples are looser, with more open contact networks under P = 0 and a larger proportion of contacts (about one third in configurations without RR) carrying vanishing forces. In the following we shall use them to illustrate qualitative tendencies in very loose samples.

When the pressure is increased to P * = 0.01, differences in force distributions between A and B samples, despite their different solid fractions (see table V), have considerably decreased, as shown on Fig. 10. The influence of such differences in the aggregation stage as those between our samples A and B are therefore expected to fade out after the systems are compressed to higher pressures and densities.

B. Packing structure and force patterns

The spatial organization of forces in type 2 samples, which we now discuss, is related to the distribution of force values, and should determine the ability of given configurations and contact networks to support stress increments. We first discuss systems without RR, then with the small RR values we adopted in most cases (see Table I). We emphasize the role of force indeterminacy and assembling history (the collisions by which cohesive clusters were built) in the final force patterns in equilibrium under vanishing or low applied stress. Extreme cases of systems with large RR on the one hand, or without friction on the other hand, are useful reference situations, which we briefly examine and discuss. We conclude this part with a discussion of the main physical implications of the relationships between force patterns, assembling process, geometry and micromechanical parameters

Qualitative aspects of force networks with no RR

It is instructive to represent the forces carried by the contact network with a visualization of positive (repulsive) and negative (attractive) normal forces, as was done on Fig. 2(b), showing the force network in one type 1 sample. Figs 11 and12 respectively correspond to equilibrated samples prepared with method 2 under P * = 0 (immediately after the aggregation stage) and under P * = 0.01, without RR. They are represented here with (approximately) the same scale. Both belong to the (Φ I = 0.26) B series. Line widths, which are proportional to the intensity of the total interaction force, i.e. to ||F|| = N n + T t , witness the presence, in spite of the low pressure, of many forces of order F 0 (which cor-respond on the figures to line thicknesses comparable to particle radii). Stressed clusters, in loose type 2 samples under P * = 0, are separated by large parts of the interacting network in which contacts carry vanishing forces: the corresponding normal deflection is h 0 (Eqn. (11)) and there has been no elastic relative tangential displacement. Attractions (green) and compressions (red) have to compensate for Eqn. [START_REF] Wood | Soil Behaviour and Critical State Soil Mechanics[END_REF] to hold true. This compensation appears to operate on a smaller scale in type 2 samples, because internal forces were previously balanced within isolated particle clusters. Such a local balance of forces is quite conspicuous at P * = 0 (Fig. 11), in which internal stresses in small clusters often take the form of a peripheral tension compensating a radial compression, or the other way round. This contrasts with samples prepared with method 1 (Fig. 2(b)), in which the spatial distribution of forces is more similar to the familiar "force chain" pattern of cohesionless systems, although there are of course compressive and tensile "chains". Unstressed regions are rather scarce in type 1 samples, although some areas with smaller forces are still present. The structure of type 2 samples under P * = 0.01 (Fig. 12) is somewhat intermediate: isolated stressed clusters are still present, but elongated, force-chain-like structures emerge.

To characterize such force patterns in a slightly more quantitative way, one can evaluate a threshold force F perc such that contacts carrying a force F with ||F|| ≥ F perc percolate through the sample. Such a criterion was used to identify a "strong" subnetwork of force chains in [75]. One observes F perc of the order of the tolerance 10 -4 F 0 in Φ I = 0.26 samples with no RR and P * = 0, which shows that stressed regions are isolated "islands" within the network. F perc raises to slightly less than 0.1F 0 under P * = 0.01.

Configurations of series A, assembled with Φ = 0.36, possess the same qualitative features, although quantitatively slightly weaker, due to their higher density. For instance, local stressed regions are somewhat less isolated, with a threshold force F perc between 10 -3 F 0 and 8 × 10 -3 F 0 at P * = 0.

Force indeterminacy (without RR)

The presence of large interaction forces of order F 0 in equilibrated samples is not obviously necessary, and is related to the assembling process. Let us imagine particles are brought very slowly, one by one, within interaction range of the previous network, thus gradually building a unique cluster in equilibrium in the absence of external stress. One could expect, rather, each new contact to stabilize with N = T = 0 and h = -h 0 . The existence of non-zero interaction forces in equilibrium is related to the hyperstaticity or force indeterminacy of the contact network. On writing all equilibrium equations for grains and collective degrees of freedom (i.e., setting acceleration terms to zero in Eqns. 1, 2 and 3) and regarding all contact forces as unknowns, the degree of force indetermi-FIG. 11: (Color online) Sample of type 2 (N=1400), in equilibrium under P * = 0 after aggregation stage, with solid fraction Φ = 0.26 (series B). Same conventions as on Fig. 2(b), except for the blue color corresponding to contacts carrying a total force below tolerance 10 -4 F0 (deflection h0 and no mobilization of tangential force). Note the large number of such interactions and the local compensation of attractions and repulsions in small prestressed clusters. To help visualize unstressed regions, disks only interacting at contacts bearing forces below tolerance are filled in light grey. nacy h (or degree of hyperstaticity) is the number of remaining independent unknowns, which cannot be determined by the equilibrium requirement. If h = 0, knowing that some equilibrium forces exist (since an equilibrium state has been found), then one would necessarily have all interaction forces equal to zero under P * = 0 (since this is one obvious possible solution). The notion of force indeterminacy has been recently discussed by different groups in the context of granular materials, essentially because of the special case of rigid frictionless grains, for which the contact network is generically such that forces are uniquely determined [73,[START_REF] Roux | Proceedings of the Saint-Venant Symposium on Multiple Scale Analysis and Coupled Physical Systems[END_REF][START_REF] Moukarzel | [END_REF]80,81]. The degree of force indeterminacy is linked to the number of degrees of freedom, equal to 3N (or 3N + 2 if the cell sizes can change), to the number of contacts N c , the number of distant interactions N d and the number of independent mechanisms or floppy modes k (also called degree of hypostaticity [81]) by the following relation (written here for a fixed cell)

3N + h = 2N c + N d + k (no RR). ( 19
)
A proof of this simple result (which is classical in structural engineering), and the relation of numbers h and k to the rigidity and stiffness matrices of the contact network, are recalled in Appendix A. Mechanisms are those sets of velocities (or small displacements, dealt with as infinitesimal) which entail no relative velocities (or small relative displacements) in contacts. For distant interactions, only normal relative velocities are relevant, hence their particular treatment in [START_REF] Agnolin | [END_REF]. In Appendix A we explain how we determine whether a given configuration is rigid, i.e., devoid of mechanisms (apart from the two global translational motions of the whole set of grains, rendered possible by the periodic boundary conditions).

It is customary to relate the level h of force indeterminacy to the coordination number z in granular materials. However, this is not possible in general, which motivated our recalling [START_REF] Agnolin | [END_REF] in its complete form. ( 19) can be rewritten, neglecting the very scarce distant interactions, as

h = N (z -3) + k.
Hence, in the absence of floppy modes, h = N (z -3). However, there are still a few floppy modes on structures like those of Fig. 11 at the end of the aggregation stage, and this relation, which predicts a small degree of hyperstaticity relative to the number of grains (see Tables IV andV), is only approximate. Some mechanisms are due to the (exceptional) 1-coordinated disks and others, less trivial, are associated with larger parts of the structure which are connected to the rest of the packing via one single 2-coordinated disk. This floppiness is obviously related to the assembling process: before any external pressure is applied, nothing really requests the aggre-gates to possess a rigid backbone. The free motion of mechanisms in assembling method 2 is largely responsible for the very long equilibration time (see Fig. 4): such motions entail no restoring force and no dissipation of kinetic energy. Floppy modes in the final state obtained with our criteria (Section II D) being scarce (typically, a few such mechanisms per 1400-disk sample), we conjecture that they would disappear entirely on adopting stricter equilibrium requirements in terms of kinetic energy. If a mechanism survives, it should generically be in motion with a non-vanishing velocity, as a residual effect of the initially agitated state. As the connected aggregate partly folds onto itself, such motions should eventually create new adhesive contacts, thereby reducing k, until the network becomes rigid. Once some rigid aggregate is formed in the assembling process, it will keep the same shape and structure, unless the collisions and perturbations it subsequently undergoes cause it to break, because of the limited tensile strength of contacts or because of the Coulomb inequality. This is the reason why the initial mean quadratic velocity of isolated grains in method 2 should be compared to V * , as given by (7).

It is easy to see that the closing of one contact can convert an aggregate from floppy to hyperstatic, the simplest example thereof being the "double triangle" structure of Fig. 13. By [START_REF] Agnolin | [END_REF], this small structure, which is rigid (k = 3 counting the free motions of an isolated object in 2D) has a degree of force indeterminacy h = 1. This is δ FIG. 13: Three (grey) disks initially forming an isostatic structure, when a fourth one (coming from the left) adheres to one of them (dotted position), it can roll (this is a floppy mode) until another (fifth) contact is formed, stabilizing it in the final position drawn with continuous line. The final "double triangle" structure is hyperstatic. Final forces (see main text) were computed for different initial velocities of the mobile disk and for different values of impact parameter δ.

how the self-stressed clusters of Fig. 11 are formed. Such structures have a strong influence on force values and force distribution. In particular, we show now that they entail specific correlations between normal and tangential force components in contacts.

Local patterns and specific force orientations

Fig. 14 shows all contact force values as points in the N, T plane for a 5600 disk sample without RR equilibrated under P * = 0. Fig. 14 displays a striking Xshaped distribution in the N, T plane, corresponding to a ratio T /N of ± √ 3 This cross pattern fades away in systems which have rearranged to support P * = 0.01, although the corresponding specific T /N ratios are still overrepresented, as shown on Fig. 15. The cross pattern of Fig. 14 corresponds to angles θ = 60 • and θ = 120 • on Fig. 15, and the second graph shows that θ = 120 • still corresponds to a peak in the distribution once the sample has been compressed (and rearranged) to P * = 0.01. As other characteristic features of force patterns in loose type 2 samples, this correlation between tangential and normal force components is stronger in the more tenuous networks of series B samples, for which the data are also represented on Fig. 15. The difference between both sample series tends to disappear on compressing to P * = 0.01 (Fig. 15(b)).

The prevalence of ratio | T N | ≃ √ 3 is in fact easy to understand. Many disks are in equilibrium with two contact forces, with two other disks which are themselves contacting each other, as on Fig. 16. (Occasionally, a third contact might be present, bearing a much smaller force, which we neglect in the present argument). In such a situation, without RR, the three equations expressing the balance of forces and moments on disk D involve four unknown force components. Labels corresponding to contacts with disks marked 1 and 2 like on Fig. 16, one obtains, on counting positively repulsive nor-FIG. 14: Values of normal and tangential contact forces in a 5600 disk, type 2 sample, in equilibrium under P * = 0, with Φ = 0.36 (A configuration). In addition to the remarkable cross-shaped pattern, marked with dashed lines of slopes ± √ 3, note the large number of very small forces, the numerous points with |T | ≪ |N | and the relevance of the value of the friction coefficient (µ = 0.5 here), as a small number of forces approach the Coulomb cone. mal forces on disk D and tangential forces with a positive moment:

N 1 = N 2 T 1 = -T 2 T 1 N 1 = T 2 N 2 = R(R 1 + R 2 + R) R 1 R 2 (20)
The ratio in [START_REF] Bathurst | [END_REF] varies for the radius distribution we are using in the present study but its most frequent value, corresponding to 17 shows the same graph as that of Fig. 14, in the case of a loose packing of disks with the same radii. In agreement with formula (20), the "X" shape is sharply defined.

R 1 = R 2 = R, is √ 3. Fig.
To understand the frequency of occurrence of very small T /|N | values, let us now consider again the smallest cluster with force indeterminacy, without RR, which comprises four disks and 5 contacts, as schematized on Fig. 18. Fig. 18 shows graphically that the balance of contact forces implies that the tangential force within the contact corresponding to the common base of the two triangles should be very small, thereby explaining the "dense line" along the N axis on Fig. forming a triangle equilibrate with zero contact forces, since there is no force indeterminacy. On simulating the collision of a fourth disk with such a triangle (as already sketched on Fig. 13) all four particles having the same radius a/2, we could observe final equilibrium situations with contact forces depending on the impact velocity, provided of course a hyperstatic structure like that of Fig. 18, with 5 contacts, was assembled. Tensile forces equal to -0.133 × F 0 in contacts C 1 and C 2 of Fig. 18 were created for a contact with an initial velocity due to the sole acceleration of the distant attractive force over distance D 0 = 10 -3 a (within a range of impact parameter δ defined on Fig. 13). Larger, repulsive forces were 16: The bottom disk, marked D, of radius R, is in contact with two other disks 1 and 2, themselves touching, whose radii are R1 and R2. At equilibrium, contact forces on disk D should be carried by the dotted line joining its 2 contact points, which determines the ratio of tangential to normal force components. observed for higher initial approaching velocities. Selfbalanced forces of order F 0 therefore naturally appear in the assembling process.

R 1 R 2 R 2 1 D FIG.

Systems with small RR

With the small level of rolling resistance we have chosen, µ r = 0.005a (see Table I), the general features of the force patterns in systems without RR are only slightly altered, as apparent on Fig. 19, which shows the interaction forces in an equilibrated sample of series B with RR under P * = 0. Like in type 2 systems devoid of RR under P * = 0, forces of order F 0 only exist in isolated regions. Note, however, that the small forces outside these regions with self-balanced stresses do not vanish, but are of order µ r F 0 /a, a feature which is further commented below.

In principle, the discussion of force indeterminacy and rigidity is quite different with RR. Contacts now carrying

D D C 1 C 2 C 1 C 2 1 2
FIG. 18: Hyperstatic 4-disk cluster, with 5 contacts. The force at the contact point between 1 and 2 should be carried by the continuous line joining this point to the intersection of the dotted lines. Those lines are respectively defined, as on Fig. 16, by the two contact points C1 and C2 of lower disk D with disks 1 and 2, and the two contact points C ′ 1 and C2 ′ of upper disk D ′ with disks 1 and 2. Note that the continuous line is close to the line of centers, hence a small value of ratio |T /N | in the contact between 1 and 2.

a moment, the analog of relation 19 becomes

3N + h = 3N c + N d + k ( RR) (21)
With RR (and as we could check with the method of Appendix A), all connected clusters are rigid. One may therefore use directly k = 0 (ignoring the two global translations) in (21). All independent loops contribute 2 to the degree of force indeterminacy, and the coordination number corresponding to isostaticity (no floppiness, no hyperstaticity) is equal to 2. To find self-balanced forces in a loop, note that each one of the two contacts of any particle in the loop will carry opposite forces (whence two independent force components); the resulting torques are then to be compensated by the rolling moments at the contacts, to be determined with a number of equations (one per particle in the loop) equal to the number of unknowns (one per contact in the loop). However, those moments are severely limited by inequality (10).

The constant force F transmitted around the loop should then be of the same order of magnitude as µ r F 0 /a , hence small. If we had used the same threshold for blue contacts on Fig. 19 than on Fig. 11, then all contacts within a loop, because they carry forces above the tolerance level 10 -4 F 0 , would have appeared as red (compressive) or green (tensile). Resetting the threshold to 10 -2 F 0 , of the order of µ r F 0 /a, thus enabled us to distinguish the hyperstatic clusters analogous to the previous case without RR from the new source of hyperstaticity, the effects of which are limited by the smallness of the RR parameter µ r . We could check that the force threshold N perc for percolation, as defined above in paragraph IV B 1, is close to 0.01F 0 in that case. If clusters made with RR, which are (infinitesimally) rigid could not be broken, no loop should appear because two independent clusters do not generically collide simultaneously in several points. The existence of loops in the final structure therefore witnesses the fragility of tenuous structures which form with a small level of RR (which is further confirmed by the large scale changes observed between P * = 0 and P * = 0.01).

Other features of force distributions and force patterns in systems without RR, such as the correlations between normal and tangential contact force components, can still be observed with the small rolling friction level µ r = 5.10 -3 a. The graphs of Figs. 14 and15, if drawn for configurations prepared in the same way with a small RR, are very similar. The small RR level used in simulations therefore only introduces small quantitative differences in that respect, at least for the parameters of the assembling procedure defined in Section III B. In the next paragraphs, we investigate, first, as an instructive limiting case, the effects of large RR, and then the situations in samples with low RR assembled with different initial random velocities (as on Fig. 8). and19, obtained with a large rolling resistance: µ r = 0.5a in a 5600 disk sample. The resulting structure has very few, large loops, hence an extremely small degree of hyperstaticity, and most contacts carry but very small forces. The characteristic prestressed clusters of Figs. 11, 12 and 19 have disappeared. Such packings with large RR therefore approach the limit in which a simple geometrical rule is adopted to aggregate particles: in the present case one recovers the results of the ballistic aggregation algorithm, stipulating that particles or clusters move on straight-line trajectories and join to form larger, rigid objects as soon as they touch. This results in isostatic, loop-free structures with coordination number 2. The resulting contact network has no force indeterminacy, and is consequently not prestressed. Our introduction, in the previous simulations, of a finite rolling resistance (and a finite friction coefficient) changes those structures in two respects: first, they form better connected structures, with loops ; second, they carry significant self-balanced forces, of the order of the maximum tensile force in a contact. Those effects are however dependent on the initial conditions for aggregation, as we now report.

Effects of initial velocities in aggregation process

As shown on Fig. 8, the initial mean quadratic velocity V 0 in the aggregation stage of assembling method 2 determines the final coordination number of systems with RR. Isostatic, loop-free networks are formed with the small RR level (µ r /a = 0.005) used in our systematic simulation series provided V 0 is small enough. The resulting force network, as displayed as an inset on Fig. 20, approaches a tree-like, loop-free structure, in which all contact forces vanish under P = 0.

Table VI shows the dependence of coordination numbers and force values on initial velocity parameter V 0 /V * . One distinguishes three populations of contacts or interactions: those with repulsive, negative and vanishing normal forces (i.e., below tolerance level 10 -4 F 0 ), and, likewise, between the average number of contacts per grains of each kind, respectively contributing z + , z -and z 0 to the coordination number z. N + (respectively, N -) is the average value of repulsive (attractive) normal forces, N

+ (resp., N (2)
-) the quadratic average. These results illustrate the dependence of the force distribution on the initial velocity parameter. Force indeterminacy and significant non-vanishing forces appear as V 0 /V * reaches V0/V * z z+ z-z0 10 2 N+ N

(2) + values of a few units, with V 0 /V * = 9.5 corresponding to the simulation series labelled A and to the force distribution shown on Fig. 9. This set of results therefore bridges the gap between our mechanical studies of cohesive particle aggregation, with the parameters given in Table I and the preparation method of Section III B, involving parameter V 0 , and the results of geometric algorithms, which are more traditional in the field of colloid aggregation [40].

10 2 N-N (2)
Geometric changes due to the breaking and rearrangements of clusters as they aggregate lead to better connected and presumably less fragile structures, which carry forces of the order of the maximum tensile force.

The special case of frictionless disks

As a complementary study of the opposite extreme case to that of large RR, we ran some exploratory simulations of frictionless, cohesive grains (also devoid of RR). In the limit of rigid disks (κ → ∞), one knows then that such assemblies are devoid of force indeterminacy: h = 0 [73,[START_REF] Moukarzel | [END_REF]80,81]. As a consequence, once large clusters are formed under no external pressure, all contacts should bear normal forces equal to zero. Such a situation is depicted on Fig. 21. The aggregate represented on Fig. 21 is obviously floppy. The analog of relation [START_REF] Agnolin | [END_REF] is

2N + h = N c + N d + k (µ = 0, no RR) ( 22
)
(It is customary, on counting degrees of freedom for frictionless disks or spheres, to discard rotations, which are all irrelevant, thereby reducing the number of degrees of freedom to 2N in the l.h.s. of ( 22) ; an alternative is to regard each rotational degree of freedom as an independent mechanism). Formula ( 22) in the frictionless case yields for h = 0 a number of floppy modes equal to 2N -N c -N d = (4 -z)N/2. The configuration of fig. 21 has a coordination number z = 3.14, hence a number of mechanisms larger than 40% of the number of particles. Such aggregates are therefore very floppy, although particles are firmly tied to their contacting neighbors. Large parts of the particle cluster of Fig. 21 are connected to the rest of the structure by only one or two contacts, thereby allowing large scale motions maintaining all contacts. Not surprisingly, the application of a small pressure P * = 0.01 to the system of Fig. 21 produces a very large compression step, resulting in the configuration shown on Fig. 22. The coordination number is now 4.01 (corresponding to a very small degree of hyperstaticity, due to finite contact stiffnesses aK N /F 0 = 10 4 , as well as to distant interactions), and the force-carrying network has a rigid subset.

We conclude that assemblies of frictionless, cohesive particles are rather singular, and do not seem capable of forming stable loose structures under a non-vanishing confining pressure. It could of course be conjectured, like in the frictional case, that floppy networks as shown on Fig. 21, with some residual motion, would gradually form better coordinated structures and eventually become rigid, but such an evolution is too slow to be efficiently followed in our simulations.

C. Discussion

The study of force values, force distributions and spatial force patterns we have been presenting here opens quite a few perspectives that are worth pursuing in more detailed and quantitative form. In particular, we have left the investigation of elastic moduli and vibrational eigenmodes of the tenuous structures formed with method 2 for future work.

However, two qualitative conclusions can be drawn, which might have broad physical relevance.

First, essentially by direct inspection of force patterns, we observed that, in loose configurations under relatively low pressure if compared to the tensile strength of bonds (as expressed by P * ≪ 1), local arrangements of grains tend to form isolated self-stressed clusters where forces are of the order of the maximum tensile force in a contact, F 0 . Those clusters comprise any number of grains between a few units to a few tens, keep the memory of the assembling process, and strongly influence the force distributions. These features are more apparent at lower densities. The degree of force indeterminacy h might be a useful indicator, but is not sufficient in itself, as it is related to the coordination number, which is very similar in type 1 (dense) and type 2 (loose) systems, and, moreover, does not account for the role of inequalities ( 8) and (10). As a general rule, loose cohesive systems tend to have a wider force distribution when h is larger, whereas the opposite behavior was observed for confined cohesionless granular materials [74]. Dense hyperstatic clusters in loose packings are connected by regions which bear very small forces. On increasing the applied pressure by small amounts, important changes occur, in which these prestressed regions merge together and large forces tend to organize in locally preferred directions, as in "force chains". Such structures are likely to play an important role in the mechanics of loose cohesive granular assemblies. Our mechanical study stresses the different effects of the two physical origins of forces -interparticle attraction and applied pressure -which tend to create different geometries, force patterns and force distributions.

Second, the structure of the loose packings and the forces they carry are strongly influenced by the assembling conditions. The relative duration of compression and aggregation processes might produce results as different as type 1 or type 2 configurations. The velocity of agitation in the initial assembling stage affects the final coordination number, as shown in Sec. IV B 6. Such parameters affect the force patterns as well, and those are also modified if contacts are initially modelled as soft (κ = 10 2 ), as in the procedure leading to configurations B.

We expect that mechanical strength properties will also be sensitive to the aggregation process.

V. GEOMETRIC CHARACTERIZATION

A. Introduction

Aggregation processes are well-known to produce fractal structures, which have been studied for many years, in particular with numerical simulations (see ref. [40] for a review). Universal fractal regimes due to various types of aggregation processes (ballistic, diffusionlimited, reaction-limited) are most conveniently observed in very low density samples. Indeed, an object of fractal dimension d F extending over distance L in d dimensions (d > d F ) will have an apparent volume fraction Φ ∼ L dF -d , which vanishes as L → ∞. Starting from N isolated particles in a finite volume with periodic boundary conditions, an aggregation process cannot produce a fractal geometry over arbitrarily large length scales. In practice, for low enough values of Φ, the aggregation process will begin just like in the Φ → 0 limit, when correlations between particles can be neglected. Later on, the crowding and interpenetration of clusters will prevent the fractal behavior to extend to larger scales [82], and a classical geometric model to describe this situation is a dense packing of fractal domains (sometimes called "blobs") of typical diameter ξ. ξ is the upper limit of the fractal regime, and is related to Φ (see the discussion of eqn. 1 in [40]) as

ξ ∼ Φ -1 d-d F , (23)
a relation which should be independent of the total sample diameter L, provided L ≫ ξ. This "fractal blob" model is reminiscent of semi-dilute polymer solutions [START_REF] De Gennes | Scaling Concepts in Polymer Physics[END_REF] and has been employed in many different physical situations, e.g., silica aerogels [START_REF] Hasmy | [END_REF]. It has been shown to describe experimental results on the packing of cohesive powders [START_REF] Castellanos | [END_REF]59]. If such a geometric description applies to our loose systems, then ξ should be of the order of the typical size of large density inhomogeneities (dense regions or holes) in the samples depicted on Fig. 5.

B. Definitions

Self-similarity is conveniently detected on studying the density autocorrelation function (DACF), as follows. Let χ(r) denote the indicator function of solid particles, taking values 1 if point r is within a solid disk, and zero otherwise. Then we define the DACF as:

C(r) = χ(R)χ(R + r) R = 1 A χ(R)χ(R + r))dR, (24)
with an average over the origin position R over the whole sample surface, of area A. On computing C(r) periodic boundary conditions should be accounted for, so that position R + r stays within the simulation cell. Isotropy ensures that C(r) is only dependent on distance r = ||r|| in the large sample limit (or on taking its ensemble average). C(r), by construction, takes the value Φ (the solid fraction) for r = 0, and tends to Φ 2 as r → ∞.

In practice it is convenient to calculate, rather than C(r) -Φ 2 , its Fourier transform, a function of the magnitude k = ||k|| of wavevector k by isotropy, which we denote as I(k). I(k) is simply related to the Fourier transform χ of the field χ(r) by:

I(k) = | χ(k)| 2 A . (25)
The notation I(k) is of course reminiscent of the scattering intensity per unit volume for wavevector k (as used in e.g. small-angle X-ray or neutron scattering experiments), which is equal to I(k), up to a "contrast factor", replaced by 1 in (25).

A fractal structure with dimension d F in 2D should have a power-law decreasing scattering intensity over some range of k:

I(k) ∝ k -dF ( 2π ξ ≪ k ≪ 2π a ). (26)
An exponential cut-off of the decreasing power law behavior of C(r) around r ∼ ξ is sometimes used [85,86]:

C(r) -Φ 2 = Φ r ℓ dF -2 e -r/ξ (27)
where the length ℓ, introduced to make C(r) appropriately dimensionless, is a constant of the order of the average particle radius. Then the corresponding form of I(k), is given in terms of Gauss's hypergeometric function 2 F 1 [a, b; c; x] [87]:

I(k) = Cst+ Φ2πℓ 2 Γ(d F ) ξ dF ℓ dF 2 F 1 1 + d F 2 , d F 2 ; 1; -ξ 2 k 2 (28
) In 2D, as soon as the particles form one continuous aggregate, the empty space is split into a set of disconnected holes or pores. The distribution of sizes and shapes of such holes is another way to characterize the system geometry.

C. Procedure

To compute I(k) from the configurations obtained in simulations, we first discretized the density field χ(r), i.e., we considered its values on the points of a regular mesh, with spacings ∆x and ∆y along the edges of the rectangular cell of the order of a/100. χ(k) was then evaluated using a two-dimensional FFT algorithm, from which I(k) was deduced by formula (25) and orientationally averaged on binning values of wavevectors k according to k = ||k||.

The field χ(r) = 0 defines a set of holes. We characterize a hole labelled as H by the value of its equivalent radius R H . R H is defined as the radius of a disk with the same radius of gyration as the hole. Specifically, if N H is the number of mesh nodes in the hole, which are labelled as i, 1 ≤ i ≤ N H , and have coordinates x i , y i , on denoting as (x c H ,y c H ) the coordinates of the mass center of the hole, one has:

R H = 2 N H NH i=1 [(x i -x c H ) 2 + (y i -y c H ) 2 ] (29)
Holes have complicated shapes and may be characterized by other quantities such as eccentricity or higher geometrical moments, but such refinements lie outside the scope of this paper.

For each sample, we record the first weighted moment (or mass average) of the distribution of hole equivalent radii, R w , defined as

R w = n H=1 N H R H n H=1 N H ( 30
)
where n is the total number of holes in the sample. In loose cohesive samples, we obtained a rapid power-law decay for the shape of this distribution. Definition (30), rather than a simple number average, ensures that the very small cavities (formed by three or four disks in contact) do not dominate in the evaluation of the average and R w indeed characterizes the large pores in the loose packings. However, this definition can only be applied when holes do not percolate through the aggregate. Thus, we have restricted the calculation of R w to samples with a non-vanishing confining pressure, P * > 0, in which case we regard it as an independent measurement of length scale ξ.

D. Results

Functions I(k) are shown on Fig. 23, along with their fits by Eqn. 28, for P * = 0 or 0.01 with and without RR, for the configurations of series A (parameters of Table I and Φ I = 0.36). The FFT calculations have been averaged over different characteristics density maps, and the bars denote the standard errors. To carry out these tories and stick to one another, forming rigid objects, as soon as they touch. The tenuous, loop-free structure of such objects, as previously commented, is retrieved in our simulations on using a large rolling friction coefficient or a small level of initial velocity fluctuations. If measured on such samples as those of Fig. 20, the same result was obtained, as expected: d F = 1.56 ± 0.04. With small RR or larger initial velocities, our simulations produce structures with, apparently the same fractal dimension, but a larger coordination number. Another observation from Fig. 23 is the presence of a slight bump (maximum) in I(k), for 2π k ≃ 10a, (which is not present in the fitted function 28). Such a feature is analogous to the peak in the structure factor of dense particle assemblies, and is likely related to the packing of aggregates. As the aggregates are requested to be mechanically rigid they tend to be multiply connected and, at least in 2D, nearly impenetrable: the maximum in the structure factor is a signature of steric exclusion.

The evaluation of the first weighted moment R w of the distribution of hole equivalent radii, for P * = 0.01 yields R w /a = 6.6 ± 0.2, and R w /a = 5.47 ± 0.14 with RR. As expected, these results are similar to the values of ξ given in Table VII.

VI. CONCLUSIONS AND FINAL REMARKS

Cohesive packings can form equilibrium structures at very small solid fractions in qualitative agreements with experiments in fine and ultrafine powders. Equilibrated configurations are sensitive to the level of applied pressure, relative to contact tensile strength F 0 , as expressed by the dimensionless number P * . They crucially depend on the assembling procedure, as even a low pressure (P * ≪ 1) can lead to rather compact states if applied to a initial "cold" (i.e., with vanishing or low velocities) granular gas of isolated particles, as in method 1. If, on the other hand, particles are given some random motion and have time to stick to one another before having to sustain some stress, as in method 2, tenuous particle networks and open structures are obtained. The initial random motion, which is ballistic in our simulations, could be diffusive in practical situations in which fine particles are dispersed in a fluid. Some random relative motion of different particles is also present, due to velocity fluctuations, in situations of global flow or suspension sedimentation.

Under low pressure, such loose packings carry self-balanced forces of the order of the maximum tensile force F 0 in hyperstatic, well-connected lumps joined by thinner arms where many contacts carry vanishing or very small forces. Such structures are sensitive to the magnitude of initial velocities with which particles collide on forming aggregates. In general force networks differ from the usual "force chain" patterns of cohesionless systems, and are associated with different force distributions. The force balance is strongly influenced by the structure of small aggregates that are first created on assembling the system. They evolve very fast as the system rearranges when P * grows even by small amounts (from P * = 0 to P * = 0.01).

Due to the limited strength of contacts with respect to tangential relative displacement and rolling, forcecarrying structures therefore differ from the ones obtained with geometry-based algorithms in which any particles or clusters that join form one unique rigid, unbreakable object. The result of such algorithms is however retrieved, in the presence of rolling resistance, if large strength properties are attributed to contacts (to the RR parameter µ r in particular) or if initial velocities of colliding grains are kept low enough. In such limits isostatic, loop-free clusters are formed with coordination number 2.

Micromechanical parameters do otherwise influence the structure of packings and the initial (self-balanced) forces they carry, especially those without rolling resistance.

The study of density correlations show that loose configurations can be regarded as dense packings of selfsimilar blobs of typical size ξ (about 10 times as large as the average diameter in our case), as in fractal clusters produced by colloid aggregation models. The estimated value of the fractal dimension, with RR, is compatible with the 2D result for ballistic aggregation, even when the connectivity (coordination number) is different. We thus expect different structures of the same density and fractal dimension to possess, due to the difference in loop numbers and self-stresses, different mechanical properties.

The fractal dimension appears to be larger in systems without RR. Thus systems without RR seem to exhibit systematic qualitative peculiarities, and since a small level of rolling resistance is likely to exist in all realistic models, this feature should preferably be included in numerical studies.

The effect of a growing pressure, as well as pressure cycles, on the packing density and internal state will be investigated in a forthcoming publication [START_REF] Gilabert | Computer simulation of model cohesive powders: plastic compaction under isotropic loads[END_REF]. Other immediately related perspectives are the study of macroscopic tensile and shear strength in relation to geometric characterizations and self-balanced forces.

Degrees of force indeterminacy h, of velocity indeterminacy k and their relations are properties of the rigidity matrix G, which is defined as follows. First, let us denote as U a displacement vector for all degrees of freedom in the system, U = ((ũ i , ∆θ i ) 1≤i≤n , (ǫ α ) 1≤α≤2 ) .

(A1) in which one conveniently separates out in the displacement u i of grain i the part due to the global strain, thus writing u i = -ǫ • r i + ũi . U has dimension 3N + 2 for N disks and 2 strain increments. Then for each one of the N c contacts, say between i and j, the relative displacement of the contact point (with notations R i,j for the radii, and tij for unit tangential vectors as in Sec. II B)

δu ij = ũi + ∆θ i × R i tij -ũj + ∆θ j × R j tij + ǫ • r ij , (A2)
can be regarded as providing 2 coordinates to one 2N cdimensional vector of relative displacements, δu. As (A2) expresses a linear dependence of δu on U, one has defined a 2N c × (3N + 2) matrix, which is the rigidity matrix G:

δu = G • U (A3)
All coordinates of u and δu are to be thought of as small (infinitesimal) increments, for which the system geometry is fixed. The degree of displacement (or velocity) indeterminacy k is by definition the dimension of the null space of G. The relevant definition of relative displacements includes all relative motions that are associated to forces or moments. In the presence of RR, one should include all relative rotations δθ i -δθ j into the components of δu, the dimension of which thus raises to 3N c . On the other hand, in the absence of friction the tangential relative displacement of the contact point becomes irrelevant, and δu should only include normal relative displacements. In general all distant attractions between close neighbors should be dealt with similarly, because only normal forces are transmitted between such pairs. For future use we just denote as M the appropriate dimension of the relative displacement vector. On writing δu it is most convenient to use a local basis for each contact, with normal and tangential directions as coordinate axes. Increments of contact forces, and possibly moments (with RR), are related via the contact law to δu. Together they define a contact force vector f , the dimension of which is equal to that of δu. f , in a system with RR, also includes rolling moments at contacts.

Externally applied forces and torques onto the grains, as well as stresses, define together a vector of external forces F ext :

F ext = ((F i , Γ i ) 1≤i≤N , (Aσ αα ) 1≤α≤2 ) .
(A4)

A denotes the surface area of the sample, so that the work of the load for small displacements is just F ext • U.

The equilibrium relations, stating that contact forces f balance the load F ext , just read (as one easily checks):

F ext = T G • f , (A5)
with the transposed rigidity matrix, T G. That matrices appearing in relations (A3) and (A5) are transposed to each other is just a statement of the theorem of virtual work : the work of external forces in any displacement vector is F ext • U = f • δu, provided F ext is related to f by (A5) and δu is related to U by (A3). By definition, the degree of force indeterminacy h is the dimension of the null space of T G.

The rank of matrix G is r = N fk, with N f the number of degrees of freedom (the dimension of displacement or external load vectors). This rank r is also the dimension of the range of the matrix, which is the orthogonal subspace, within the M -dimensional space of relative displacements, to the null space of its transpose T G in the dual space of contact forces. Hence r = Mh. We have obtained

N f + h = M + k,
which yields, according to the appropriate definition of relevant relative motions, relations [START_REF] Agnolin | [END_REF], (21), and (22).

Assuming elastic behavior in the contact (i.e., strict inequalities in ( 8) and (10), which, as noted in Section II D, is the general case at equilibrium), in a quasistatic experiment contact force increments ∆f relate to relative displacement increments ∆δu with a contact stiffness matrix K:

∆f = K • ∆δu.
K is a square, diagonal matrix, containing coefficients K N , K T and (with RR) K r for each contact. Thus K only contains positive elements, except for the (very scarce) distant interactions, which contribute the negative normal stiffness -F 0 /D 0 in our model. If ∆f balances some load increment ∆F ext , while δu corresponds to the N fdimensional displacement vector U, one then has:

∆F ext = K • U,
where one has introduced the stiffness matrix K:

K = T G • K • G. ( A6
)
(K is traditionally called dynamical matrix in the context of solid-state physics and interactions of atoms or ions in a crystal [START_REF] Ashcroft | Solid State Physics[END_REF]). Unlike G and T G, K is always a square, symmetric matrix. It has to be positive definite in order for the equilibrium state to be stable, because it expresses the elastic energy associated with small displacements. (In fact, the full stiffness matrix also contains a small non-symmetric correction to (A6) [START_REF] Bagi | On the Concept of Jammed Configurations from a Structural Mechanics Perspective[END_REF] due to the effect of contact forces prior to the application of the load increment, which we ignore here.)

By construction, the null space of G is contained in the null space of K, and coincides with it in the absence of distant attractions, because K is then a positive matrix. In practice, the positiveness of K can be investigated with the Cholesky algorithm. We applied this method (in a form suitable for sparse matrices, stored in a "skyline" form) to the stiffness matrix of the contact networks of the simulated equilibrium configurations. This is how, on finding that K was positive definite, we could conclude that the contact structure was devoid of mechanisms (or floppy modes, eigenmodes of K with eigenvalue zero) in all cases with P * = 0.01. On the contrary, stiffness matrices associated with contact structures without RR at P * = 0 usually possess some mechanisms, although we

FIG. 1 :

 1 FIG. 1: Graphical representation of the model for the adhesive elastic contact force as a function of the distance between the surfaces of particles i and j, hij. (a) The elastic normal force consists of a repulsive Hookean part N e ij plus a linearized attractive part N a ij . (b) The elastic tangential force is limited by the Coulomb cone (adhesion shifting its tip to -F0 on the normal force axis).

  FIG.2:(Color online) Aspect of force-carrying structures in cohesionless and cohesive samples. Contact forces are displayed with the usual convention that the width of the lines joining the centers of interacting pairs of disks is proportional to the normal force, on scale aP (left) and F0 (right). Red, green, and blue lines distinguish compressive, tensile, and distant interactions in the cohesive case, while rattlers appear in grey in the cohesionless sample.

FIG. 3 :

 3 FIG.3: Distribution of local coordination numbers (percentage of total particle number), without (left graph) and with (right graph) cohesion.

FIG. 4 :

 4 FIG.4: Solid fraction versus time for both preparation procedures, showing some aspects of the configurations at different stages. Point A is the initial state (or ΦI). Aspects of configurations are shown for intermediate states B1 and B2, and for final equilibrated states C1 and D2 (at P * = 0.01). Point C2 corresponds to the stage when all disks are assembled in a unique aggregate, then equilibrated at P * = 0 (both aggregation and equilibration stages take place between A and C2). The time unit is T0 = ma/F0. Note the duration of the preparation process with method 2, and the difference in final equilibrated states compared to method 1.

1 0

 1 15 T -1 0 . Using the samples with Φ I = 0.13 (the lowest value used in this work) with K N = 10 5 F 0 /a, three different values for ǫmax were tested: 0.10 T -1 0 , 0.05 T -and 0.015 T -1 0 . The condition ǫ ≥ 0.10 T -1 0

FIG. 5 :FIG. 6 :

 56 FIG.5: Configuration of 1400 disk sample of series A without rolling resistance. Note the gradual closing of pores as the external pressure is increased from P * = 0 (first image) to P * = 0.01 (last image) going through two intermediate stages.

FIG. 7 :

 7 FIG. 7: Typical configurations of 1400 disk samples of series A with (left) and without (right) rolling resistance, at P * = 0 (a) and P * = 0.01 (b). Note the difference in local structure of thin "beams" joining dense regions with or without RR.

3 FIG. 8 :

 38 FIG.8: Final coordination number z versus initial quadratic average velocity in agitation stage of method 2, normalized by characteristic velocity V * . The arrow points to the value most often used in our calculations.

5 - 1 5 - 1 5 FIG. 9 :

 515159 FIG. 9: Distribution of normal forces for series A samples. The non-cohesive case (9(a)) is normalized by the average repulsive elastic part. The cohesive cases (9(b) and 9(c)) are normalized by F0 (note that the average of the elastic part of N is Ne ≃ F0 in cohesive cases with P * ≪ 1).

FIG. 10 :

 10 FIG. 10: Comparisons between probability distribution functions of normal force values in samples of type A (histogram, in black) and B (shaded histogram, grey) without RR under P * = 0.01.

FIG. 12 :

 12 FIG.12:(Color online) Sample of Fig.11, with same scale and color conventions, in equilibrium under P * = 0.01. The solid fraction increased to Φ = 0.39. The threshold force (used to distinguish blue lines and grey disks) was set to 0.01F0.

14 .FIG. 15 :

 1415 FIG.15: Histograms of angle θ, between normal vector n and total contact force F. Conventionally, θ = 0 • for a repulsive normal force and T = 0, and θ = 180 • for a tensile normal force and T = 0. Shaded histograms (grey) correspond to B configurations (ΦI = 0.26), bold-line non-shaded ones (black) to A ones (ΦI = 0.36)

FIG. 17 :

 17 FIG. 17: Values of normal and tangential contact forces in a 5600 disk, type 2 sample of monodisperse disks in equilibrium under P * = 0. Note the sharp"X" shape on blown-up detail of small forces.

5 .

 5 Effect of a large rolling resistance.

Fig 20

 20 Fig 20 shows the analog of Figs.11 and 19, obtained with a large rolling resistance: µ r = 0.5a in a 5600 disk sample. The resulting structure has very few, large loops, hence an extremely small degree of hyperstaticity, and most contacts carry but very small forces. The characteristic prestressed clusters of Figs. 11, 12 and 19 have disappeared. Such packings with large RR therefore approach the limit in which a simple geometrical rule is adopted to aggregate particles: in the present case one recovers the results of the ballistic aggregation algorithm, stipulating that particles or clusters move on straight-line trajectories and join to form larger, rigid objects as soon as they touch. This results in isostatic, loop-free structures with coordination number 2. The resulting contact network has no force indeterminacy, and is consequently not prestressed. Our introduction, in the previous simulations, of a finite rolling resistance (and a finite friction coefficient) changes those structures in two respects: first, they form better connected structures, with loops ; second, they carry significant self-balanced forces, of the order of the maximum tensile force in a contact. Those effects are however dependent on the initial conditions for aggregation, as we now report.

FIG. 19 :

 19 FIG. 19: (Color online) Same as Fig. 11, in a (B series) sample with RR, Φ = 0.26. Threshold force 10 -2 F0.

FIG. 20 :

 20 FIG. 20: (Color online) Same as Fig. 19, in a sample with large RR, µr = 0.5a, N = 5600, and Φ = 0.26. Inset: force network in N = 1400 sample obtained with low initial mean quadratic velocity V0 and small RR (corresponding to the bottom left point on Fig 8).

FIG. 21 :

 21 FIG. 21: (Color online) Same as Fig. 11, in a sample with N = 1400, Φ = 0.26, no RR and no friction (µ = 0).

FIG. 22 :

 22 FIG. 22: (Color online) Sample of Fig. 21 under P * = 0.01. Φ increased to 0.72.

FIG. 23 :

 23 FIG.23:(Color online) Scattering functions I(k) of samples with and without RR for P * = 0 and P * = 0.01, averaged over 4 samples of 1400 disks and 2 with 5600 disks. Fits of data points with Eqn. 28 are drawn with continuous lines. Both with and without RR, I(k) is larger for P * = 0, corresponding to larger density fluctuations.

TABLE IV :

 IV ± 0.0003 (series A0) 0.1303 ± 0.0003 3.197 ± 0.002 0.1304 ± 0.0003 2.656 ± 0.007 0.2649 ± 0.0006 (series B) 0.265 ± 0.001 3.123 ± 0.004 0.2652 ± 0.023 2.963 ± 0.006 0.361 ± 0.007 (series A) 0.3616 ± 0.0003 3.1407 ± 0.0016 0.361 ± 0.009 2.660 ± 0.004 Values of Φ and z on equilibrating configurations at P * = 0 with µ = 0.5. Samples with ΦI = 0.26 correspond to κ = 10 4 , the others to κ = 10 5 .

		no RR		RR
	ΦI (z = 0)	Φ	z	Φ	z
	0.1301 no RR		RR
	ΦI (z = 0)	Φ	z	Φ	z
	0.2649 ± 0.0006 (series B) 0.448 ± 0.006 3.235 ± 0.003 0.42 ± 0.01 3.085 ± 0.005 0.361 ± 0.007 (series A) 0.472 ± 0.008 3.175 ± 0.003 0.524 ± 0.008 2.973 ± 0.004

TABLE V :

 V Values of Φ and z in equilibrated configurations at P * = 0.01. These results are averaged over the whole set of samples prepared with µ = 0.5, for ΦI = 0.26 (series B) on the one hand, and for ΦI = 0.36 (series A) on the other hand. Series A0, prepared with ΦI = 0.13, yielded very similar results but due to computational costs the number of samples was too small to record data in statistical form.

  fits, we have applied the Levenberg-Marquardt method for nonlinear least-squares fittings[START_REF] Moré | Lecture Notes in Mathematics[END_REF]. This fitting procedure yields values of d F and ξ listed in TableVII. As expected, the fractal dimension is conserved in the compaction between P * = 0 and 0.01, but the fractal range shrinks. The marked difference in d F caused by the introduction of a small level of RR is remarkable. = 0.01 dF 1.925 ± 0.024 1.93 ± 0.04 1.53 ± 0.04 1.51 ± 0.04 ξ/a 8.29 ± 0.15 6.07 ± 0.2 9.3 ± 0.4 5.06 ± 0.21

	no RR		RR
	P * = 0	P * = 0.01	P * = 0	P
				While self-
				similar clusters are very nearly dense (d F approaching
				2) without RR, more open fractal structures are stabi-
				lized on small scales with µ r = 0.005a. This value of the
				fractal dimension obtained with RR appears to coincide,
				within the error bar, with the value d F = 1.55 ± 0.02 ob-tained for the ballistic cluster aggregation model [39, 89],
				assuming particles or clusters move on rectilinear trajec-

*

TABLE VII :

 VII Fractal dimension and fractal blob size obtained on fitting the data of I(k) to Eqn.(28).

argued that their number k must be small. Aknowledgements: This work has been supported by the Ministerio de Educación y Ciencia of the Spanish Government under contract number BFM2003-1739. J.-N. Roux wishes to thank Dietrich Wolf for useful contacts and discussions.