Wang-Landau algorithm for entropic sampling of arch-based microstates in the volume ensemble of static granular packings

We implement the Wang-Landau algorithm to sample with equal probabilities the static configurations of a model granular system. The"non-interacting rigid arch model"used is based on the description of static configurations by means of splitting the assembly of grains into sets of stable arches. This technique allows us to build the entropy as a function of the volume of the packing for large systems. We make a special note of the details that have to be considered when defining the microstates and proposing the moves for the correct sampling in these unusual models. We compare our results with previous exact calculations of the model made at moderate system sizes. The technique opens a new opportunity to calculate the entropy of more complex granular models.


I. Introduction
In the study of static packings there exists still a lack of predictive capabilities of the available theories. Assemblies of objects that pack (such as grains) can generally sample such packed configurations only by the external excitation of the system. These packings can be built by repeating a given packing protocol (e.g., homogeneous compression or deposition under an external field against a confining boundary) on an initial random configuration. Also, a Markovian or non-Markovian series can be constructed by exciting the system from the previous packing configuration. To what extent the series of packings obtained (using either type of protocol) can be modeled without information on the dynamics that drives the system to the packed configuration is still uncertain. The main reason for this is that the few statistical approaches that attempt to do this are strongly hinder by the poor current ability to generate such packed structures without using a dynamics to build the packings.
One might expect that the packing fraction and its fluctuations, among other properties, could be obtained from basic statistics without resourcing to a full molecular dynamic type of simulations (also known as "discrete element method", DEM). Although these types of simulations are powerful enough to predict the behaviour of most systems that pack, it is desirable to find a description that could neglect the detailed dynamics between consecutive packed configurations.
The use of the tools provided by the ensemble theory of statistical mechanics in problems of granular matter is still limited. Although many studies perform statistical analysis of configurations of a granular sample obtained during careful preparations in the laboratory or in molecular dynamictype simulations, very rarely sampling in a particular statistical ensemble is carried out either via an 070001-1 analytic or a computational calculation of a model system. This has prevented a direct assessment as to whether ensemble theory is appropriate to describe the behaviour of these peculiar systems.
One pioneering contribution to the topic is the idea that granular systems at mechanical equilibrium could be treated as ensemble members, putting forward the conjecture that the mean values of measurable quantities could be calculated using statistical mechanics for these ensembles [1]. In this scheme, each of the thermodynamic variables finds a counterpart, the volume taking the place of the energy, and the compactivity that of the temperature, amongst other transformations. However beautiful this description may seem, the computational challenges to generate ensemble samples in this context are extraordinaries [2][3][4][5].
One complication that prevents to a large extent the use of the machinery of statistical mechanics in this case is the fact that configurations, unlike in traditional liquid theories, have to be checked for the constraint of mechanical equilibrium. In a previous work [6], we have made a proposal on how to deal with this, at least in a first approximation, by describing the excitations of static granular systems under gravity in terms of its arches. Since arches are sets of grains that stabilize each other, these are the basic units of mechanically stable structures in the packing. Any static configuration can be described in terms of the arches formed by its grains, their arch shape, position, orientations, etc. We have considered, as an example, a model of a two-dimensional (2D) granular system composed of disks where arches are assumed to take a single possible structure and the arch-arch interactions (due to the interlocking of arches) is neglected. We calculated the exact entropy of this model (the noninteracting rigid arch model, NIRA) by constructing all possible configurations for moderate system sizes. Of course, generating each state is a cumbersome task if the system size is increased or if the number of degrees of freedoms (DoF) is increased by using more realistic models. Therefore, an alternative approach based on sampling the phase space with the desired probability is necessary.
In the present work, we will calculate the entropy of the NIRA model in the microcanonical ensemble [7] using entropic sampling through the Wang-Landau (WL) algorithm [8][9][10][11]. This approach allows us to obtain the entire entropy function for all possible volumes of the system in one single simulation for larger systems and potentially for more complex models. All derived properties, such as compactivity and volume fluctuations, can then be calculated through numerical differentiation. We pay particular attention to the different descriptions that can be realized for the NIRA model. Some of these representations do not provide a direct way of sampling the configuration space uniformly.
This work is organized as follows: in section II., we will review the WL algorithm. In section III., we will review the NIRA model and discuss different ways of representing it, along with the issues related to uniform sampling of the configurations. We then present a representation that allows very fast calculations of the entropy and we compare the results with the exact counting of all configurations for systems of moderate size. Finally, we discuss future directions to refine the arch-based ensemble volume function towards capturing detailed features of more realistic systems.

II. Wang-Landau algorithm
The WL method has revolutionized computational statistical mechanics [8][9][10][11]. WL is a pure statistical method that can retrieve the density of states (DoS) (hence the entropy) over a bounded region of the energy spectrum from the sole knowledge of the energy function. The spectacular computational performance achieved by this method stems from the fact that it presents no limitations for the system to tunnel between potential barriers, in stark contrast with classical Monte Carlo methods that underperform when they encounter deep valleys in the energy landscape.
WL finds the entropy of the system by means of a Markov chain in the energy landscape which is conveniently biased towards the less probable energies in a strongly history dependent manner. This is achieved by using the multicanonical approach [12] in which each possible configuration is sampled with a probability given by the inverse of the density of states for its given energy. Specifically, WL aims to obtain a flat histogram of visited energies E by forcing the system to go through all configurations with a probability which is inverse to the previous occurrence of that energy in the Markov chain. The method is ergodic and asymptotically fulfils the detailed balance condition [13]. There exists extensive literature on the WL algorithm but here, we only summarize its most relevant steps. In the following sections, we will refer to a configuration of the system as a fixed set of values of all its DoF.
In WL, one defines two histograms that are continuously updated as the Markov chain proceeds. These histograms are the entropy S(E), which is the output of the algorithm, and a control histogram H(E). After initializing H(E) = 0, S(E) = 1 and a starting configuration with energy E 0 , the rules to update these histograms are: I Propose a new configuration and calculate its energy E 1 . The new configuration is generally derived from the previous configuration by a change in the value of one of its DoF.
II Accept the new configuration according to a probability given by: III Update the two histograms in the correct energy bin E (E 1 if the new configuration is accepted and E 0 if otherwise), accordingly: H(E) = H(E) + 1 for the control histogram, and S(E) = S(E) + f for the entropy. Here, f is a correction that controls the precision of the algorithm, which will be decreased (see next step), usually starting at f = 1.
IV If the control histogram H(E) is flat enough according to some arbitrary criterion, decrease f (for instance by making f = f /2) and reset all entries of H(E) to zero.
V If f > ǫ (with ǫ a prescribed tolerance), return to step I, otherwise stop.
After each reduction of the correction term f , the entropy histogram is built with a finer grained precision. However, to speed up the initial estimates of S(E), f is set to a high initial value. Different approaches are followed to accelerate the final stages of refinement by decreasing f with alternative criteria [14].
In Monte Carlo approaches, the detailed balance condition ensures that the Markov chain has a limiting distribution [15]. Detailed balance can be stated as follows: where p µ is the probability distribution of configuration µ and P (µ → ν) the transition probability from configuration µ to configuration ν which can be written as: with S P (µ → ν) the selection probability, which is the probability that the algorithm generates a trial configuration ν starting from configuration µ; and A P (µ → ν) the acceptance probability, i.e., the probability that the algorithm will accept the trial configuration ν. Hence, since the target distribution in the entropic sampling is the inverse of the DoS, i.e., (1) and (2) implies In WL, detailed balance is not fulfilled in general. During the construction of entropy, the acceptance probability given in step II above (i.e., A P (µ → ν) = min [1, exp(S(E µ ) − S(E ν ))]) evolves and it is only when the entropy gets sufficiently refined that detailed balance is met.
Notice that the form chosen for A P requires that the selection probability be symmetric (i.e., S P (µ → ν) = S P (ν → µ)). Although this condition is simple to comply with in models like Ising, for arch-based descriptions of static granular packs this is non-trivial. Therefore, one must be careful to represent a system in such a way that trial moves between different configurations have the same backward and forward selection probability.
After reviewing the main characteristics of the arch-based ensemble in the next section, we will show that some natural representations of the microstates lead to non-symmetric selection probability schemes. We present, however, a way of representing the configurations that does allow for the direct use of WL. In all the WL simulations, we have used 300 bins for the histograms. The tolerance for the f correction was set to ǫ = 2 −15 . The histogram H(E) is considered flat (see step IV) whenever (H max − H min )/H max < 0.2, with H max and H min the maximum and minimum height of the histogram.

III. Arch-based microstates
In Ref. [6], we have introduced a way of describing the microstate of a static granular system under gravity by considering the arches that the grains form instead of the more traditional approach of using the particle positions. Arches are defined as sets of mutually stable grains. All other particles being fixed, the removal of any of the grains in the arch would induce the destabilization of the rest of the set. Any assembly of grains, static under gravity, can be split into a number of arches which are mutually exclusive [16,17].
The major difficulty in sampling static granular configurations is the fact that these are sparse (with zero measure) in the overwhelming number of possible particle positions. Moreover, there are not recipes to generate a static configuration from another static configuration by simply moving a grain from its position.
Since each arch is stable on its own right, the arch-based description warrants that any configuration proposed fulfils a basic stability constraint; i.e., that each set of grains identified as an arch has internal contacts that keep it stable. The problem is now moved to generating all possible combination of arches, including how many of them are of a given size, shape, orientation and position and also generating all arrangements of these that can be stable resting on each other. Of course, all of these DoF can be represented with different levels of approximation.
In Ref. [6], we have described the five general steps necessary to carry out an Edwards entropy calculation (i.e., the number of states associated to each given volume of the packing) within an archbased scheme. These are: 1. Define the microstate of the system in terms of arches.
2. Define the external constraints imposed to arches.
3. Define a volume function that yields the total volume of the microstate in terms of the arches.
4. Define an algorithm to generate all microstates defined in step 1 that comply with the external constraints of step 2. Or sample microstates with equal probabilities.
5. Calculate the volume of each microstate generated in step 4 using the function in step 3 and build the DoS.
Step 4 may constitute a significant limitation to the real possibility of calculating the density of states for systems of reasonable size. Generating all configurations is certainly impractical for most models (especially if they have continuous DoF). This paper demonstrates how to sample configurations by using WL (rather than generating all of them) to accomplish step 4.
We will focus on a model we have already solved exactly by counting every single possible configurations so as to have a reference system to compare with the WL algorithm results. This is the "non-interacting rigid arch model" (NIRA). In this model, only the number of arches n i of each size i (in number of grains) that form part of the packing is used to describe the microstate. All arches consisting of the same number of grains are considered to occupy the same volume (hence one single possible shape is assumed for each arch size; this is implied in the word "rigid" used to name the model). The total volume of the system is assumed to be the sum of the volume of the individual arches and the arch-arch interlocking DoF is not taken into account (hence the word "non-interacting"). To represent a two-dimensional pack of equal-sized disks, we have taken the volume v i of an arch of i grains to be the area under the regular polygon that inscribes all disks in a "regular" arch. The special cases of "arches" of one particle and twoparticle arches are considered separately [6]. The total volume V [{n i }] of the system is Where N is the total number of disks of diameter d in the system. It is important to mention that, typically, the maximum size that an arch can take is physically bounded (e.g., due to the size of the container that holds the granular sample). Hence, we will put a cutoff C to the largest arch allowed in the system. The cutoff C is an external constraint (that is imposed in step 2, above). Importantly, this cutoff imposes a limit to the correlations in the system which leads to an extensive entropy (see Ref. [6]) In counting microstates, one has to bear in mind that arches of the same size are indistinguishable in the NIRA model, whereas arches of a different size can be distinguished. Hence, if there are n i arches of i grains in a configuration (with i = 1, ...C, being C the size cutoff), then the number of permutation of arches that yield distinguishable microstates with these arches is where N A = i n i is the total number of arches of the configuration (including those "arches" of size 1). Despite all the simplifications, the model applied to a 2D system of equal-sized disks yields qualitative agreement with DEM simulations of tapped disks [6]. The NIRA model is in many respects similar to an ideal gas of excitations or quasi-particles (the arches) with a single DoF (their size). Figure 1(a) shows the entropy S calculated by counting all possible microstates for 500 disks using different C for the largest arch allowed [6]. As it is expected, if C increases, looser configurations are possible and hence states with higher volumes become more significant. However, the entropy for low volumes rapidly converges to a well defined curve. Part (b) of Fig. 1 shows the compactivity χ defined as χ −1 = ∂S/∂V . This is the analogue to the temperature in thermal systems [1]. The entropy presents a maximum as observed by others [2,3]. States for volumes beyond this maximum correspond to negative compactivities [see Fig. 1(b)]. This is caused by the inversion population of these volume bounded systems. Some authors suggest these negative χ macrostates may be inaccessible, however, this does not need to be the case; some preparation protocols may indeed lead to very low packing fractions [3,18]. An interesting prediction of the NIRA model is that systems constrained by different C achieve the same χ at different specific volume V /N (i.e., at different packing fractions). As a consequence, two samples of grains "equilibrated" to the same compactivity will show distinct packing fractions if the maximum arch size possible in each sample is different. Different values of C in practice may be achieved by using narrow containers or by changing the static friction coefficient of the grains. There have been some progress in the study of the equilibration of vibrated granular samples in "contact" [19]. However, there are still no attempts to couple static granular packs under gravity. Further developments in this direction may help validating this prediction of the NIRA model.
The configurations of the NIRA model are compatible with different representations. In the following subsections we will discuss some of these representations and their suitability for the implementation of the WL algorithm.

i. The arch size distribution representation
In our previous paper [6], we have used a vector {n i } that represents the number of arches consisting of i grains in the configuration, i.e., {n i } = (n 1 , n 2 , ..., n C ), with C the largest arch allowed in the system and n 1 the number of grains not forming arches. As an example, a possible configuration in a system of N = 10 grains and a cutoff arch length 070001-5 of C = 6 represented in this way could be: 1, 0, 0, 1, 0). (6) In this example, there are three grains not forming arches, two grains forming an arch of size two, and five grains forming another arch of size 5. In Ref. [6], we have swept all possible configurations and multiplied each by its analytical degeneration due to the different permutations of arches with repetitions given by Eq. (5).
It is difficult to propose an algorithm to move between configurations represented in this way and yet comply with the symmetry of the selection probability required by the WL algorithm. For example, consider a move that consists in removing a grain from one arch of size k and adding it to another arch of size k ′ . Such move would require subtracting 1 from the coordinate k of {n i } and adding 1 to the coordinate k − 1, since this arch will now be smaller by one grain. Additionally, the coordinate k ′ of {n i } needs to be reduced in 1 and the coordinate k ′ + 1 increased in 1, since this arch will now be part of the set of arches larger by one grain. There are different ways of selecting a grain to be moved and to select its arch of destination.
Let us consider, for instance, that all grains and destination arches are chosen with same probability. Now consider a move in the Markov chain that takes configuration (6) [i.e., µ = (3, 1, 0, 0, 1, 0)] into configuration ν = (2, 1, 0, 0, 0, 1). This corresponds to taking one grain that was not forming an arch and inserting it in the five-particle arch to make it a six-particle arch. The probability of selecting a particle from an "arch of size one" in this case is 3/10. The probability of choosing the arch of size five as the destination is 1/5 (there are four other arches in configuration µ plus the possibility of leaving the grain on its own without forming an arch with others). Hence the selection probability is S P (µ → ν) = 3/50. A similar analysis shows that to return to the original configuration S P (ν → µ) = 6/10 × 1/4 = 3/20 (there are 6 grains out of ten that can be taken from the six-grain arch and there are three possible other destination arches plus the case with the grain not forming an arch in the new configuration). Clearly, this selection probabilities are non-symmetric as they should be to apply the algorithm of section II.
A possible workaround to the previous representation is to multiply each coordinate of the vec-tor {n i } by the corresponding arch size. Therefore, each coordinate now indicates how many grains are involved in all arches of the given size. In this new representation, the configuration of Eq. (6) (10 particles with a cutoff C = 6) is written as In this case, a trial move may consist in randomly transferring a grain from one coordinate k to another coordinate k ′ . This is done simply by subtracting 1 to {n ′ k } and adding 1 to {n ′ k ′ }. In this representation, the selection probability of moving a particle from one arch of size k to another of size k ′ is S P (µ → ν) = 1/N ×1/(C −1) irrespective of k and k ′ . Hence, the selection probability is symmetric and suitable to implement the entropic sampling using WL.
Unfortunately, the representation of Eq. (7) does not tell apart two microstates that differ only in the permutation of two arches of different sizes. Therefore, the corresponding degeneracy given by Eq. (5) cannot be accounted for 1 . More importantly, the change of representations from Eq. (6) to Eq. (7) is not an exact mapping because these newly proposed moves allow for "fractions of arches" to exist since we do not request that the coordinates of {n ′ i } be a multiple of i after each move. For instance, the configuration (0, 0, 0, 0, 0, 10) is allowed in representation (7) but does not exist in representation (6). In the new representation such configuration implies that there is one arch of size 6 plus a fraction of an arch of size six. We can still treat these "fracrtions of arches" by assigning to them a fraction of the arch volume. However, there is a large number of new unrealistic configurations (there are many ways of chosing a numbers that are not conmensurate with the arch size) in this representation that will bias the result for the entropy. Figure 2 shows the entropy for the NIRA model using representation (7) for different maximum arch sizes C compared with the exact result obtained by counting all configurations and all permutations [6]. As we can see, not including all distinguishable permutations and including the new "fractional arches" gives a wrong entropy function. ii. The arch listing representation As previously discussed, other ways of representing the system should be carefully chosen in order to ensure that there exist moves that sample the different configurations uniformly. An alternative natural representation consists of using a vector, {m i }, with N coordinates, where each coordinate m i can take any value from 0 to C, the cutoff for the arch size, provided that i m i = N . The content of each coordinate indicates that the configuration has an arch of that size. The configuration of Eq. (6) in these representations can be expressed, for example, as (5, 0, 0, 1, 0, 1, 2, 0, 1, 0). (8) There are, of course, multiple permutations in Eq. (8) that lead to the same distribution of arch sizes compatible with Eq. (6). Trial moves in the system represented in this way can be done by subtracting one from a non-zero m i and adding one to any other coordinate that has a value smaller than C. This is equivalent to reducing the size of one arch in one particle and either creating a new arch of size one (if the new coordinate had a zero value) or increasing the size of another arch. Unfortunately, this algorithm has a non-symmetric selection probability.
where N A is the number of arches (i.e., the number of non-zero coordinates in the configuration) and  Fig. 2) using the binary arch representation to carry out the entropic sampling through the WL algorithm for 200 grains. The solid lines correspond to the exact counting of microstates from Ref. [6].
N C is the number of arches with the maximum allowed size C. Therefore, S P will depend on the total number of arches and the number of arches of size C.
Besides the non-symmetric S P , the system represented in this way and sampled with these trial moves clearly overestimates the number of states of a given volume. This is because, apart from the N A !/(n 1 !...n C !) permutations of the distinguishable arches [see Eq. One can, in principle, resolve these issues. The selection probability may be turned into a symmetric one by adapting the acceptance probability in Eq. (3). The degeneracy due to the presence of zeros in the representation (8) can be handled by switching to a representation without including these, and allowing for a vector of variable length. However, although less intuitive, there is a much simpler, suitable representation that we discuss in the next section.

iii. The binary arch representation
Finally, we present a representation that comply with the symmetric selection probability and simultaneously account for the permutation of distinguishable arches.
In this case, the system is chosen to be represented by a vector of N coordinates with binary values (zeros and ones). These coordinates are not associated to specific particles 1, 2, 3, etc. Rather, as we move along the vector from left to right, we can think the ones as representing a first grain in an arch (whatever its identity) and the following zeroes as the remaining particles. For instance, the configuration in equation 6 in this new representation can be given by The trial moves consist in picking a coordinate and changing the state of that coordinate (if 1 change to 0 and vice versa). This results in a symmetric selection probability of S P (µ → ν) = S P (ν → µ) = 1/N . In each move, the constraint imposed by the cutoff C must be checked and the trial configuration must be rejected whenever the constraint is not complied with.
In Fig. 3, the result for this representation and sampling strategy is plotted along with the exact result showing a remarkable agreement.

IV. Conclusions
We have been able to compute the entropy of a system of non-interacting rigid arches using a WL algorithm in the volume ensemble in different representations. We have exposed the difficulties in dealing with different representations of the configurations of arches and the mechanisms used to propose trial moves for the WL algorithm. These difficulties appear during the choice of a simple sampling scheme that ensure a symmetric selection probability of configurations. Additionally, the degeneracy due to distinguishable permutations of arches pose a further complication in the use of WL. The most suitable representation that we found for a noninteracting system of rigid arches resulted in a binary vector.
We believe that entropic sampling of arches through the WL algorithm has a great potential for testing the granular statistical mechanics hypothesis (such as equiprobability and ergodicity). Having a sampling algorithm like WL adapted for these types of models is crucial to continue the road map towards the refinement of an arch-based framework for static granular packs. In particular, the non-interacting condition is clearly a crude approximation and should be lifted, along with the introduction of a more accurate volume function.