{'fraction_non_alphanumeric': 0.030707425193941634, 'fraction_numerical': 0.007596047284817141, 'mean_word_length': 4.5267321679213985, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stress phase space for static granular matter

This paper proposes a phase space to compare the static packings of a granular system compatible to a macrostate that is set by the external stress. The nature of this phase space is analyzed, showing that the consideration of the allowed and forbidden regions and the internal degrees of freedom of every configuration (i.e. geometrical pattern) could be a relevant factor for the establishment of its probability and, therefore, of the expected properties of the sample. This is due to the fact that many combinations of forces acting on a particle can keep it in static equilibrium. Every set of forces can be considered equivalent to a microscopic stress field, but the kind of interaction and the geometrical restrictions mean that not all stress states can be represented by any set, whereas others can be represented by many sets. Consequently the points of the phase space are degenerate, and the density of states of each configuration strongly determines the most probable statistical distribution. It is shown how these functions just depend on the deviatoric stress. A first analysis of two-dimensional (2D) arrangements is included to clarify this assertion.


Introduction
Granular materials behave differently from any of the other familiar forms of matter [1,2,3]. In particular, a fundamental issue concerns the packing of these materials [4,1,2,3,5]: a granular system may be in a number of different microscopic states at fixed macroscopic constraints (densities, pressures, etc.) and many unusual properties are linked to its non-trivial packing.
Edwards was the first to propose that a statistical mechanics approach might be feasible to describe dense granular media [6]. Assuming that granular systems have entropy, it was claimed that the volume plays the role of energy (V-ensemble) [6,7,8,9,10]. Other formalisms are based on the energy of the whole system, such as [11,12,13]. Following this work, an approach based on considering the elastic potential energy of particles has recently been presented, referred to here as the elastic energy approach [14]. It is set up by considering the stress state of the system. In this respect, other previous approaches also consider the stress of the system together with the volume, Full Canonical Ensemble (or V-F ensemble) [15,16,17], or even alone (F-ensemble) [18,19]. This can be done via the forcemoment tensor and this approach is referred to here as the force-moment approach. It uses the concept of angoricity. On the other hand, the Force Network Ensemble [20,21,22,23] is not based on the stresses but directly on the forces between particles of a given contact network.
This paper deals with the establishment of a stress phase space suitable for the elastic energy and the force-moment approaches. This paper deals with the establishment of a stress phase space suitable for the elastic energy and the force-moment approaches. Both of them are based on canonical ensembles that are set up by using some functionals of the stresses. As they consider different macroscopic constraints, they lead to different statistical distributions. In practice, the main difference between them is whether the exponent of the Boltzmann factor decays as a linear or a quadratic function of the stresses. However, this paper also shows why it is not only the Boltzmann factor that determines the most probable distribution but also the degeneracy of the points of the phase space. This is precisely due to the existence of different configurations with respective allowed and forbidden regions in the phase space.

Revision of elastic energy and force-moment approaches
The elastic energy approach [14] was an attempt of expressing the Edwards' first formalism [6] in terms of elastic potential energy. As a consequence, a theoretical framework based on stresses and geometries was developed to explain the expected features of disordered arrangements of jammed granular media. The thermodynamical formalism was set up starting from the assumption that the most probable distribution of stresses of a compressed granular media can be identified with a specific value of the average elastic potential energy. When a driving process is able to switch from a given packing to any other possible packing of the ensemble, the ergodic hypothesis would be guaranteed not by averaging samples in time but in cycles. However, the physical restriction that the whole arrangement should be a solution of the elastic problem was not rigorously imposed in [14].
The force-moment approach [15,18,19] properly imposes that the assembly must fit a solution of the elastic problem and can be applied to arrangements with different values of equivalent stiffness and, therefore, of energy. This approach is based on the assumption that if particles are small in comparison with the scale of variation of the macroscopic stress field τ ij 1 , the latter must coincide with the average equivalent stress field over the whole systemσ ij . In consequence it is possible to set up an ensemble based on the additive force-moment tensor Σ ij =σ ij V , with V being the volume of the whole system. Then, the exponent of the Boltzmann factor used in the canonical ensemble is a linear function of the stresses, whereas the external control parameter, the tensorial temperature, is fixed by angoricity.
The elastic energy approach showed that consideration of the mechanical stability of geometrically compatible local packings (configurations) can be useful to determine why some of them are more likely than others and why some of them are barely present. The same methodology, revised and extended to general cases, is followed in this paper, since it is useful for a better definition of the random close packing [4] state or of the maximally random jammed state [24]: if the most probable statistical distribution of configurations were known, it would be possible to determine which the expected packing ratio of the sample is. However it is only possible if the constrains of the ensemble associated to a specific physical process are clearly established in a suitable phase space.

Equilibrium and phase space
In this section the τ (or T) phase space is proposed for those ensembles in which the stress (or the force-moment) state of the granular system is externally set (letters τ and T have been chosen because the Latin word for stress is tensio). Jammed and compressed packings of soft particles in static equilibrium with a given external stress field are supposed to correspond to a solution of a macroscopic elastic problem. A space described by stress (or force-moment) coordinates is referred to as phase space in the sense of that it is the space of all possible states of these systems (i.e. the possible solutions of the macroscopic elastic problem). Each solution or packing can be described by its contacts and forces networks. There are many packings of a given granular system that are in mechanical equilibrium and satisfy some macroscopic constraints (e.g. average elastic potential energy or force-moment). Each of these static packings is a microstate whereas the macrostate is precisely defined by reference to the macroscopic variables. The ensemble is therefore considered as a set of copies of a granular system that correspond to all possible packings for a given macrostate. Statistical equilibrium is identified with the most probable distribution in the τ (or T) phase space that satisfies the constraints of the ensemble. It is obtained by maximizing the Boltzmann's entropy S ∝ ln Ω 2 . The problem is not mechanical but statistical 1 τ ij is the stress field which fits a solution of the macroscopic elastic problem with average equivalent values of constitutive relationships.
2 Ω is the number of compatible microstates. It depends on the number of particles and the macroscopic variables one, because the role of the dynamics during the driving process is trivial for this definition: it is being compared, among the static solutions of a particular physical problem, the class that it is compatible with the macroscopic knowledge of a system. This is a Boltzmann approach to the equilibrium [25]. Note that this approach is only valid for those physical processes that make the packing switch from one contact network to another and when the size of the particles is small in comparison with the scale of variation of the expected macroscopic stress field (for instance, it cannot be applied to heaps, silos or systems with few particles).
At this point, it is interesting to compare these approaches with the Force Network Ensemble (FNE) theory [20,21,22,23], which uses a phase space defined by forces. It samples all possible force configurations for a given contact network with an equal probability. It can be applied to both hyperstatic and isostatic contact networks. However the approaches presented here sample all the equivalent heterogeneous stress (or force-moment) fields compatible to the constraints that are given by the ensembles, irrespective of the contact network and irrespective of the boundary forces 3 . It is also assumed that in the thermodynamic limit it is guaranteed that all the sampled heterogeneous fields actually correspond to several real contacts networks.
The equivalent stress of each domain (Voronoi cell) in a packing depends on the forces between the particles, which are strongly correlated each other (especially when the packing is close to the isostatic state or J point [26,27]). Equivalently, it can be asserted that the average stress of each domain depends on how it interacts with others: in series (equal stress, or Reuss case), in parallel (equal strain, or Voigt case), or mixed. In short, the equivalent stress state of a particle k in a given packing π depends on the geometry of the Voronoi tessellation, on the elastic properties of the particles and on the stress state of others (mainly on the stress state of its closer neighbors), so thatσ k,π ij =σ k,π ij(σ l,π ij ) (where l = k represents each of the z k,π particles interacting with the particle k in the packing π). If the packing is being changed during the driving process (tapping, anisotropic compression, etc) then these relationships also change, in such a way that ∀π ′ = πσ k,π ′ ij(σ l ij ) =σ k,π ij(σ l ij ) ). Well, this work starts from the hypothesis that, as the relationships between the particles are short-ranged (only close domains interact with each other, ∀k, π z k,π << N ) and as they are constantly changed during the driving process, then it is possible to be unconcerned about where (in the real space) each particle is or which its neighbors are. It is equivalent to assume that when all possible packings are sampled, the average relationship between the stresses can be neglected <σ k ij > =<σ k ij > (<σ l ij >) orσ k ij =σ k ij(σ l ij ) (< y > representes the average value of the function y (σ k ij ) over all the different packings compatible with some constraints). It means that in the thermodynamic limit, and as result of a complex mechanical interaction, a particle can be located at any accessible point in the stress (or force-moment) phase space. The probability of being located at each point depends on both internal and external factors (i.e., microscopic features of the particles and external macroscopic constraints). In consequence, within a statistical mechanics framework it is possible to obtain which the most probable distribution of points in the phase space is.
In addition, in the case of the ensembles reported in section 2 above, maximizing the Boltzmann's entropy leads to Maxwell-Boltzmann statistics. In both cases the exponent of the Boltzmann factor is represented by an additive and completely separable function of the phase space coordinates, , because both the total elastic energy and the total force moment are equal to the sum of the separate contribution of each particle. This, together with the assumption that in the thermodynamic limit and when all possible packings are sampled the average relationships between the stresses of the particles are negligible (σ k ij =σ k ij(σ l ij ) ; ∀k = l), makes it possible to use the τ (or T) phase space of one particle. Then the whole partition function is equal to the partition function of one particle to the power N and the "phase space" of a single particle can be used.
Now it is possible to analyze how the nature of the phase space of a single particle is. It has two important features: 1) its accessible region (given by the integration limits of the coordinates) and 2) the degeneracy of the points or density of states, which comes from the capability of geometrical patterns, or configurations, to cover regions of the phase space. Both these issues are analyzed below together with an explanation for why using a change of variables is interesting.

Coordinates (p, q, ω) and (P, Q, Ω)
If the ensemble is based on the total elastic potential energy of the granular system, a natural workshop is provided by the components of the stress tensor σ ij [14]. In the static case, for elastic continuum media the energy density is given by where S ijkl is the compliance tensor of the material.
Within granular media the particles are in equilibrium under a finite number of forces. If they are elastic, they store elastic potential energy, which for a particle k is given by where σ k ij is the actual stress field within the volume of the particle V k P , S k ijmn is the compliance tensor of the material,σ k ij is the stress field averaged not only over the volume of the particle but also over the respective Voronoi cell V k , andS k ijmn is the equivalent compliance tensor. f kl is the module of the compressive force exterted by particle l on particle k 4 and K D is the stiffness of the particles (i.e., the ratio between the applied force and the reduction of the diameter).
In order to translate a set of forces into an equivalent average stress field the average stress tensor (see [28]) can be used:σ where r kl is the position of the point of application of the force f kl .
According to the possible arrangements, some sets of forces can be equivalent to the same stress field. Therefore, the phase space determined byσ ij is a good workshop to compare microstates, and according to the value of the forces it is possible to calculate the elastic potential energy. This is the so-called τ phase space.
Nevertheless, the nature of granular media makes it preferable to use other stress state variables. For the case of 2D media, considering the invariants of the secondorder stress tensorσ ij , the variables (p, q, ω), (σ I , σ II , ω), or (σ, m, ω) are more suitable. The coordinates p, q, ω are based on Mohr's circle [29] and are defined as p and q 5 (or q ′ = q/p), the centrum and radius of Mohr's circle, relate to the invariants I 1 =σ ii and I 2 = ǫ ijσxiσyj of the Cauchy stress tensor (being ǫ ij the Levi-Civita symbol); therefore, they are also independent of the coordinate system. However, ω depends precisely on the orientation of the principal directions of the tensor (given by its eigenvectors). Other alternatives are to use the principal stresses σ I and σ II (also invariants) together with ω, or the minor principal stress σ and the anisotropy factor m (the ratio between both stresses m = σ I /σ II ).
Using these coordinates rather than the components of the Cauchy stress tensor is interesting because some physical restrictions and some magnitudes are conveniently expressed as functions of them; for instance, the elastic potential energy of an isotropic continuum medium does not depend on ω, but it often does when the medium is not isotropic (e.g., a particulate medium). In this paper, (p, q ′ , ω) are used for the stress phase space. In short, they express the stress level, the relative deviatoric component of the stress state, and the orientation, respectively.
The volume of the accessible region of the τ phase space of one particle V τ can be expressed as The integration limits 6 are those which guarantee that, for any state, the respective principal stresses are positive. It means that the stress field is compressive in any spatial direction (as it is expected for cohesionless materials). The principal stresses are the eigenvalues of the Cauchy stress tensor and, according to the characteristic equation, both of them are positive whenσ 2 xy ≤σ xxσyy .
On the other hand, the force-moment approach is based on the force-moment tensor Σ k ij =σ k ij V k = l r kl i f kl j . For a system of N particles, the total Σ ij is given by the sum of the force-moment tensor of every particle Σ ij = k Σ k ij . Therefore, the coordinates of the phase space can be directly Σ ij . It is being referred herein to this space as the T phase space.
However, according to the invariants of Σ ij , it is also possible to define the coordinates P, Q, Ω as P =Î 1 /2, In this paper, (P, Q ′ , Ω) are proposed for the forcemoment phase space. The integration limits and the volume element are similar to those used above: P ∈ [0, ∞), 1], Ω ∈ [0, 2π) and dV T = 2P 2 Q ′ dP dQ ′ dΩ. 6 Actually, ω should be evaluated within the interval [0, π), because the stress states of ω and ω + π are the same. Nevertheless, this doubled evaluation results in a more intuitive representation. Note that, in this way, a determined system of forces, equivalent to a stress field, would be represented by two symmetrical points in the phase space.

Elastic potential energy ensemble
Following the approach of elastic potential energy [14], the Boltzmann factor is can be expressed as: β relates to a quadratic form of the components of the external pressure tensor and it is determined by some values of the expected equivalent stiffness of a packing in statistical equilibrium. This approach does not mean that the energy is constant during the driving process, but that, once the kinetic energy has been dissipated and the packing is again in mechanical equilibrium with the external pressure, the elastic energy of the ensemble is the same that it was before the driving. As the energie of a each particle can be expressed as , the probability of the configurations actually depends on the functionsM α(q ′2 ,ω) 8 . The energy increases with the stress level and with the anisotropy factor m, so the statistical weight decays.

F-ensemble
Following the force-moment approach, the Boltzmann factor would depend on the force-moment tensor [15,16,18,19] as where χ ij is the inverse angoricity, a second-order tensor. If a coordinate system x, y according its principal directions, i.e., those in which χ xy = 0 9 is used, then the Boltzmann factor does not depend on Σ xy .
On the other hand, the physical meaning of the angoricity can be deduced by obtaining average values. For instance, by integrating over the volume of the T phase space of a single particle (analogous of Eq. 7 for Σ ij variables, instead ofσ ij ), it is obtained that the expected value of Σ xx is given by: The 3/2 and 1/2 powers on Σ xx come from the integration limits of Σ xy (Σ xy ∈ (− Σ xx Σ yy , + Σ xx Σ yy )). Acting in the same way, < Σ yy >= (3/2)(1/χ yy ) and < Σ xy >= 0. Using Σ ij = Vσ ij andσ ij = τ ij it is shown that χ ij is actually 3/2 of the inverse of the average external force-moment tensor V τ ij per particle, i.e., 7 In that work the function was expressed in terms of other variables Eα = E I(K D ,σ 2 ) M α(m,ω) 8 The reference system is established according to the orientation of the configuration, so that ω refers here to the angle of the principal stresses with respect to the configuration. 9 The inverse angoricity and the external stress are tensorial fields which can be related (Eq. 11). They share the principal directions, which are given by the eigenvectors of either tensor.
where V is the volume of the whole system. However, the exponent of the Boltzmann factor can also be expressed as where χ P = (χ xx + χ yy ) and χ Q ′ Ω = (χ yy − χ xx )/(χ yy + χ xx ) ∈ [0, 1). The advantage of writing the Boltzmann factor in this way is that, in the statistical weight of each point of phase space, it separates the contribution due to the force-stress level P from the contribution due to the relative deviatoric component Q ′ and the orientation Ω, depending on the external control parameters χ P and χ Q ′ Ω .
Accordingly, and as some physical restrictions are independent of the stress level, it is also interesting to consider the distribution of points in Q ′ , Ω planes, for any value of P or for all of them. The integration of the points over all the values of P and its representation over a polar Q ′ , Ω plane is referred to as the comparison plane (i.e. a plot which displays Q ′ sin Ω as ordinate plotted against Q ′ cos Ω as abscissa) 10 .
If every point in the T phase space of a single particle is supposed to have the same multiplicity, the partition function (when negative stresses are not allowed) is given by Integrating over all the values of P leads to the expected density in the comparison plane ρ (Q ′ ,Ω) , which just depends on χ Q ′ ,Ω .
For the particular case of isotropic compression, χ xx = χ yy = χ P /2, so that the exponent of the Boltzmann factor just depends on P , i.e., χ ij Σ ij = P χ P . This means that, if the granular system is in equilibrium with an external isotropic stress field, the probability of the particles in the ensemble will not depend on either Q ′ or Ω, so that the expected statistical distribution of particles according to their respective values of P will be satisfying ∞ 0 ρ T(P ) P 2 dP = 1.0. Note that the volume element is dV TP = P 2 dP and it affects the expected number of particles n P in each P ± ∆P/2 interval: Moreover, integrating over P , in the isotropic compression 10 In the τ phase space, the comparison plane plots the polar distribution of q ′ and ω, integrated for all p values. case every point in the comparison plane is supposed to be equally likely to be found On the other hand, for anisotropic stress fields, the higher the anisotropy ratio, the stronger the dependence of the Boltzmann factor on Q ′ and Ω. This is shown in Fig. 1, where it is plotted for different values of χ Q ′ Ω and P . As Ω is fixed according to the principal directions of the external stress field (via the angoricity), the expected local stress fields are not isotropically oriented but according to some prevailing directions.

Configurations.
Allowed and forbidden regions 6.1 The degeneracy of the points in the phase space Equilibrium configurations are defined herein according to the number of forces that make a particle be in static equilibrium. The points of application of the forces can be taken as vertices of a representative polygon (so it changes when the relative angles of the forces change). Then it is possible to establish a representative direction of the configuration in base of this polygon (e.g. any of the diagonals of a regular hexagon). The orientation of the configuration is established according to its representative direction and the external frame of reference. On the other hand, the equivalent stress field can be obtained from the value of the forces and their relative angles. It is a tensorial field and its principal directions define the orientation of the equivalent stress field, which can be considered in relation to either the representative direction of the configuration or the external frame of reference.
Irrespective the ensemble, it is highly relevant that some points of the stress phase space can be filled by several configurations while other points are forbidden for some of them. Moreover, some configurations have internal degrees of fredom because they can fill some of these points being oriented in several directions, having several representative polygons or even achieving the force balance with several values of the forces (it is only possible in hyperstatic configurations). In consequence the phase space must be characterized by a density of states function which comes from the sum of the allowed regions of all configurations in all possible orientations with respect to the principal directions of the external pressure. Then, Eq. 13 becomes where g (Q ′ ) is the total density of states function. In the next section below it is explained why this function is supposed to be independent of both Ω and P . Since this function is actually equal to the sum of the respective contribution of every configuration g α(Q ′ ) , it is also possible to express the partition function as The functions g α(Q ′ ) are equal to 0 within the respective forbidden regions but can take a normalized value within allowed regions (due to their internal degrees of freedom). Expressing the partition function in this way reflects exactly not only the stress state in which particles of every microstate are located but also the configuration to which they belong. Then, some expected properties, especially the associated packing ratio, can be directly measured if the distribution function is separable in configurations.

Physical reasons for allowed and forbidden regions
The functions g α(Q ′ ) depends on the kind of interaction between the particles and on the geometric compatibilities. Another important advantage of using the coordinates (P, Q ′ , Ω) is that the physical restrictions of the equilibrium are independent of the force-moment level. In consequence the ability of a configuration to fill a point only depends on Q ′ and Ω but it does not on P . Thus the comparison plane can be used to analyze the size and shape of the respective allowed region.
To illustrate this assertion, the simplest case of a 2D frictionless monodisperse system is analyzed. In this medium, only five configurations are possible, defined respectively by the number of forces L: 6, 5, 4, 3, and 2. The forces f kl acting on the particle are ordered according to the angle θ kl with respect to the representative direction of the configuration. There are three kind of physical restrictions: • Restrictions due to the kind of interaction In cohesionless materials, all forces must be positive. If the interacting forces are expressed as f kl = f kl r k −r l |r k −r l | (being r k , r l the positions of the particles k, l), then this restriction means that: ∀k, l f kl ≥ 0. (19) • Restrictions due to mechanical equilibrium In 2D there are just two restrictions per particle k: With i = x, y. As there is no friction, all forces are central, and there are no moments.
All these restrictions are independent of the force scale because if the forces are increased by X times, Eqs. 19, 20, and 21 are still satisfied. It makes the allowed regions to extend over the same Q ′ , Ω values in any constant-P plane. So, in the case of 2D systems it is not necessary to analyze the whole volume of the phase space, but just the comparison plane.
Then, the respective g α(Q ′ ) of a configuration is equal to the normalized integration over any orientation. These functions include the features of the particles (size distribution, shape, roughness, strength, etc.) because they determine which configurations are possible and which g α(Q ′ ) they have.
7 Qualitative analysis of a 2D monodisperse system

Configurations and allowed regions
Some of the possible equilibrium configurations of a 2D monodisperse granular system have been analyzed. This was done by numerically generating combinations of forces and relative angles which maintain the equilibrium and which correspond to the same value of P (P = 1.0) 11 . Then, points were plotted on the comparison plane.
The surveying of combinations of forces and angles for a configuration was done by taking possible values of some of them. The number of independent variables to choose depends on the configuration. In this work, some representative combinations of angles were fixed whereas the domain of the independent forces was divided into a uniform grid, according to a fixed interval. For each point of the grid, the other forces were obtained by applying the restrictions of force balance and of P = 1.0. The resulting groups which included negative forces were automatically discarded. The interval of the grid was taken with values of different orders of magnitude and it was observed that all the points were always located within the same region. When the elastic energy approach was set up [14] the boundaries of the six-force configuration were analytically obtained (although other variables were used in that paper) and that procedure gave the same result. Some configurations analyzed in this way are presented.
In the six-force configuration, the relative angles are fixed because inequalities of angles (Eqs. 21) become equalities. Thus the orientation of the configuration can be assigned to any of the three diameters of opposite forces. Once the orientation is fixed, it is possible to choose four independent forces and analyze which combinations of positive values are possible. Plotting all these points over the comparison plane reveals the allowed area of this configuration (Fig. 2). It is a key factor that there are three directions (with respect to the orientation of the configuration) in which the equivalent stress field can take any value of Q ′ between 0.0 and 1.0. This is the configuration which covers the largest area (in 2D systems).
In the five-force configuration it is possible to choose three forces and five angles, although angles are strongly bounded due to the inequalities (Eqs. 21). It is not clear if every combination of angles can be considered as an independent configuration, because the multiplicity of some (Q ′ , Ω) points would be too high. In any case, some representative cases are plotted in Figs. 3, 4, 5, and 6. It is important to realize that there are only points with high values of Q ′ for a few representative polygons. Indeed, the highest value of Q ′ = 1.0 is obtained precisely in the quasihexagonal arrangement (which coincides with a six-force configuration in which one force is null).
In the four-force configuration it is possible to chose two forces and four angles in order to sweep the phase space. Angles are somewhat less bound than in the fiveforce case. Some combinations of them lead to high values of Q ′ but only over a narrow area (indeed, only some lines), so that this configuration is in general quite unstable. Again, it is not clear whether the angles can change without breaking, which would extend the stability area of the configuration. This idea was used in a previous paper [14], where it was supposed that a transition between rhombic configurations happens until the system reaches a totally forbidden area. The result would be a four-pointed star area, covering some high values of Q ′ .

Numerical simulations
Two different and representative numerical molecular dynamics simulations are presented as qualitative analysis of the distribution in the T phase space. The purpose of these simulations is to check whether, at different stages of the compression process, the distribution of the points on the comparison plane is affected by the allowed regions of configurations or not. However the macroscopic constraints of the ensembles explained above are not totally imposed, because it is much more easy to perform simulations in which the displacement of the walls or the external pressure are controlled, rather than the force-moment or the final elastic energy of the packing. Nevertheless, as it is not being measured neither the entropy nor the ensemble average of any quantity, these simulations are valid for the objectives of this paper. In addition, after a cyclic compression schedule, if the variation of the volume of the sample is small, then the constraints of the force-moment ensemble are closely matched.
The LAMMPS code, a parallel particle simulator devel-oped at Sandia National Laboratory [30] was used (including the GRANULAR package). The systems consisted of 900 (30 × 30) and 4, 900 (70 × 70) frictionless disks 12 , which interacted via a viscous and linear (Hookean) contact law (i.e. f kl = (K H δ − γ m 2δ )n kl , being δ the overlapping between particles k and l and n kl = r k −r l |r k −r l | ). A normal damping coefficient γ ≃ 2 √ 2 K H /m was used in order to dissipate the kinetic energy and compare static solutions. Particles were initially arranged according to a regular square lattice. Some irregular walls were also generated, two of which were fixed while the others were subject to controlled motion or applied force. The compression schedule is described by a first isotropic compression stage and then some anisotropic compression cycles. The sole purpose of the irregularity of the walls was to avoid the whole crystallization of the sample. These walls were generated by locating particles in consecutive contact and in such a way that they separate from the average position of the wall according to a normal distribution. Both the duration of the isotropic compression stage and the period of the cycles of compression ∆t were chosen in such a way that ∆t/t c ∼ 10 4 , being t c the characteristic time t c = m/K H . In the case of 900 particles, an external pressure of p/K H ∼ 10 −3 was applied (p is given in units of force per unit of length) whereas in the sample of 4, 900 particles a displacement of the walls of ∆L ≃ 0.07L 0 was imposed (being L 0 ∼ 70 × D and D the diameter). Fig. 7 and Fig. 9 show the aspect of the granular packings, and Fig. 8 and Fig. 10 show how the respective population density on the comparison plane of the two simulations was. In the first case, this measure was taken after 4 anisotropic compression cycles whereas in the second case it was taken at the end of the isotropic compression stage. The first conclusion is that, although it is an isotropic compression state, the points are not uniformly distributed within the circle as expected for an isotropic compression state, providing evidence that the Boltzmann factor is not the only parameter which determines it (at least if the force-moment approach is followed, because results show that, the lower the value of Q ′ , the higher the density of points). Therefore, the density of states g (Q ′ ) strongly determines the expected distribution of particles in force-moment levels.
The number of points obtained in these simulations is not enough to ensure thermodynamic limit conditions. Therefore, obtaining histograms (which would be necessary to know how the g α(Q ′ ) functions are and how the distribution depends on the protocol) it is not still possible. Nevertheless, this plot shows how the allowed regions have a strong influence on the statistical distribution. It is quite interesting that many of the particles are precisely located at the boundaries of the allowed and forbidden regions of the main six-force arrangement.

Probability of configurations and transitions
According to this approach, not only the statistical weight but also the density of states determines which the statistical equilibrium of a system is, but it is necessary to emphasize that this state corresponds to the most probable stress distribution when all possible contact and forces networks are sampled. However this is not necessarily guaranteed by any protocol. In consequence the expected properties of a sample can be dependent on the protocol (e.g. see [31,32,33,13]). It can be explained in base of the constraints that the protocol is actually imposing to the static solutions and its capacity of to make particles switch from one configuration to another.
Although further research is needed on both issues (which are out of the scope of this work), the consideration in detail of the nature of the τ (or T) phase space may be useful to do that. The allowed regions of configurations and the density of states functions must affect not only the probability of finding a particle in a specific geometrical pattern and in a specific stress state but also the stability of configurations in a driving process. This could explain why, for instance, in 2D monodisperse media, protocols often cause most of the particles to arrange according to the six-force configuration, although other configurations may also be expected (because they also have considerable allowed areas).
When the external conditions of the whole system are modified during a driving cycle, another solution of the elastic problem can be achieved. As the forces on each particle can change, it is possible to draw its path in the τ (or T) phase space. If this path remains within the allowed region of its previous configuration, the particle can be kept arranged according to the same geometrical pattern with other values of the forces. The particle could also change into other configurations without breaking if they overlap their allowed regions. Finally, if the stress path reaches a region where the configuration is not possible, it will break locally 13 . Once this happens, the particle could fall into any configuration, according to the probability of each.

Conclusions
It is possible to compare different static packings (contact and forces networks) of granular systems compatible to the same external constraints, in a suitable τ (or T) phase space which considers the equivalent stress fieldσ k ij (or the force-moment tensor Σ k ij ) of the particles. When all possible packings are sampled, it is possible to use the phase space of one particle to establish which the most probable distribution of stresses (or force-moments) is. The nature of the phase space it is governed by its 13 Sometimes, the hyperstaticity of the whole sample can make a particle keep in the boundary of the allowed region while others change their stress state to find a collective solution. It occurs until a collective failure is produced. accessible region and the degeneracy of the points. Coordinates (p, q ′ , ω) (or (P, Q ′ , Ω)) based on the invariants of the stress tensor are more appropriate to analyze these issues. According to the scale invariance of the mechanical restrictions, for 2D systems it is not necessary to analyze the whole one-particle phase space but just the P distribution and the area in the comparison plane [polar representation of (Q ′ , Ω) for all P ]. The size and shape of the allowed regions and the degeneracy of states of each configuration are due geometrical incompatibilities, mechanical constraints, and the kind of interaction between the particles. It can be expressed as a density of states function g (Q ′ ) . More research is needed to analyze this function and how it can be separated according to the respective contributions of the different configurations g α(Q ′ ) . The statistical weight of the points in the phase space is fixed by the ensemble. Nevertheless, the most probable distribution is strongly affected by the density of states. The intrinsic properties of the medium, i.e., the features of the particles (size distribution, shape, roughness, etc.), determine these functions, whereas the external stress state fixes the angoricity or the elastic energy (the parameter controlling the statistical distribution). As granular systems do not tend to equilibrium by themselves, it is necessary to drive them. The constraints fixed by a driving process together with its capacity of sampling microstates determines the expected distribution. The shape and size of the allowed regions of configurations could determine the probability for a particle to switch between configurations.
10 Acknowledgements I thank Ricardo Brito (Complutense University of Madrid) and Rafael Jiménez (Technical University of Madrid) for stimulating discussions, technical support and sincere encouragement. Figure 1: Representation of the Boltzmann factor (using Eqs. 9 and 12) over two constant-P planes according to the force-moment approach. χ P = 1.0. The higher the P , the smaller the value of the Boltzmann factor, whereas χ Q ′ Ω sets the dependence on Q ′ and Ω. (a) P = 1.0 and χ Q ′ Ω = 0.33, (b) P = 1.0 and χ Q ′ Ω = 1.0, (c) P = 2.0 and χ Q ′ Ω = 0.33, and (d) P = 2.0 and χ Q ′ Ω = 1.0. If χ Q ′ Ω is 0.0 (meaning that the external stress is isotropic), every point within any constant-P plane is equally likely to be filled and uniform gray circles are expected. Figure 2: Allowed area of a six-force configuration. The only possible combination of angles corresponds to a regular hexagon. This is the configuration which covers the largest area in the constant-P plane. It is about 55% of the expected maximum area (the circle, which corresponds to an ideal, infinite-strength continuum medium). The highest value of Q ′ which can be resisted in three directions is 1.0, and some wide areas of high values of Q ′ are also covered Figure 3: Allowed area of a five-force configuration with pentagonal orientation of forces. The area is quite isotropic, and its size is about 32% of the expected maximum area. The highest value of Q ′ which can be resisted in five directions is less than 0.6 Figure 4: Allowed area of a five-force configuration with a slightly modified pentagonal orientation of forces. The area is slightly more anisotropic than in the regular pentagonal case. The size is about 31% of the expected maximum volume, and the maximum value of Q ′ is higher than the regular case in just one direction Figure 5: Allowed area of a mixed hexagonal-square configuration. This five-force configuration is usually found in the grain boundary of crystalline regions. The size is about 24% of the expected maximum area, and the highest value of Q ′ is again 1.0, being possible in only one direction Figure 6: Allowed area of a quasihexagonal five-force configuration, being equivalent to a six-force configuration in which one force is zero. The size is about 19% of the expected maximum volume, and the highest value of Q is P . However, it is possible in just two directions. As the shape is narrow and strongly anisotropic, this configuration is supposed to be quite unstable Figure 7: Arrangement of 900 elastic disks after the isotropic compression stage and some anisotropic compression cycles. The sample was initially ordered according to a square lattice, but the boundaries were irregular. Figure 8: Distribution of points in the comparison plane of the simulation of isotropic compression of 900 disks after the isotropic compression stage and 4 anisotropic compression cycles. The number of particles is far from the thermodynamic limit conditions and most of the points are located within the respective allowed region of the corresponding six-force configuration associated to the main crystalline domain (Fig. 7). Figure 9: Arrangement of 4, 900 elastic disks after the isotropic compression stage. The sample was initially ordered according to a square lattice, but the boundaries were irregular. In this case, two main crystalline domains (grains) are found in the arrangement. Figure 10: Distribution of points in the comparison plane of the 4, 900-disk compression simulation. It is still clear that the respective allowed regions of the two grains affect the statistical distribution of the points: two overlapping six-pointed stars can be discerned.