\section{\label{sec:level1}Introduction}
This paper discusses the structural characteristics of random packings of planar (2D) granular solids. The science of random granular packing has a very long history and statistical methods have been introduced into this field in the 60's \cite{Bernal1960, Mogami1965}. In 1989, Edwards and co-workers \cite{EdOa89a,EdOa89b} proposed a novel statistical mechanical approach, based on the entropy of configurations that assemblies of grains can take. In this approach, energy is replaced by volume, the Hamiltonian is replaced by a volume function and the temperature, characterising ensemble fluctuations, is replaced by compactivity. This latter quantity has been elusive to measure and quantify from bulk quantities. 
Much work followed, aiming both to use this formalism for derivation of macroscopic quantities, as in thermodynamics, and to test it experimentally and numerically  \cite{Knight-etal1995, Richard-etal2005, Song-etal2008}.


Fundamental to this approach is an identification of basic volume elements, which can be used as 'quasi-particles' for the volume function. Blumenfeld and collaborators have proposed that these quasi-particles are {\it quadrons} \cite{BaBl02,BlEd03,BlEd06} - structural elements that tessellate the granular space and can be used to quantify the local structure \cite{BaBl02,Fretal08,Fretal09}. Ball and Blumenfeld \cite{BaBl02} used the quadron description to construct a stress theory of isostatic granular materials.  
In contrast to other definitions of volume elements, such as Volonoi-based tesserations \cite{Song-etal2008}, quadrons preserve the connectivity information, since their construction is based on the force-carrying inter-granular contacts. They also allow for an unambiguous tensorial description of every structural element. A further advantage of this description is that, at least in 2D, there are as many quadrons as there are degrees of freedom in the statistical mechanical structural phase space \cite{BlEd03,BlEd06,Bl08a}. This makes it useful for treating these volumes as degrees of freedom in their own right.

It has been shown \cite{Fretal08} that the quadron volume probability density function (PDF), $P(V_q)$, is best understood in terms of its conditional PDF's $P(V_q\mid e)$, where $e$ is the order of the cell associated with the quadron (see Figure \ref{fig:quadron}). Based on geometrical considerations, and supported by numerical results, Frenkel et al. argued that the PDFs $P(V_q\mid e)$ should be independent of inter-granular friction and of the packing process (i.e. the history). 
A detailed analysis of this issue and an understanding of the shape of these PDFs from first principles is yet to emerge.

Here, we present a set of numerical experiments on assemblies of polydisperse discs, isotropically packed by the same process to different packing fractions $\phi$ and mean coordination numbers $\bar{z}$.
We investigate in detail the geometrical characterisation of the quadrons and the dependences  of the PDFs $P(V_q)$ and $P(V_q\mid e)$ on $\phi$ and $\bar{z}$.  We examine the decomposition proposed by Frenkel et al. \cite{Fretal08}, focusing on the significance of cell shapes, and we show that the universality suggested by them must be augmented with a mechanics-based consideration. 
This is supported by a free cell model, whose results are compared with the numerical experiments. 
We conclude with a discussion of the results and the new insight that they provide into structures of random granular packs in mechanical equilibrium. 

\section{Probability density function of quadron volumes in statistical mechanics}

As in thermodynamics, one can define the entropy, $S$, of a micro-canonical ensemble of granular assemblies in static equilibrium as the logarithm of  the number of ways, $M$, that a collection of $N$ grains fit into a given volume $V$,

\begin{eqnarray}
S(V,N)=\log M
\label{S}
\end{eqnarray}
where $M$ can be written as
\begin{eqnarray}
M=\sum_i \delta (V-V_i) G_i \to \int \delta (V-W)\ \Theta\ d^{N_{sdf}}\{\eta \}
\label{M}
\end{eqnarray}
Here $\delta$ is the delta function, $G_i$ is the number of configurations (microstates) in mechanical equlibrium having a total volume $V_i$,
$W$ is the volume function - a function that gives the volume in terms of $N_{sdf}$ internal structural degrees of freedom (SDF), $\{\eta \}$, for any grain configuration. $\Theta$ is a product of delta functions, constraining the possible configurations to those specified in the ensemble statistics, e.g. to assemblies in mechanical equilibrium, with a specified mean coordination number $\bar{z}$. From this micro-canonical description, one can readily pass to the canonical one \cite{Edwards_Rafi2007book}.
Assuming a Gibbs ensemble \cite{EdOa89a,EdOa89b}, the occurrence probability of a configuration of volume $V_i$ is 

\begin{equation}
P_i = \frac{1}{Z_v} G_i \exp \left( - \frac{ V_i }{X} \right)
\label{Pi}
\end{equation}
where $X = \left( \partial S / \partial V \right)^{-1}$ is a Lagrange multiplier, called compactivity, that characterizes the statistical fluctuations and is analogous to the temperature. 
$Z_v$ is a volumetric partition function

\begin{equation}
Z_v = \sum_i e^{-V_i/X} \to \int e^{-W\left(\left\{ \eta \right\}\right)/X} \Theta(\left\{ \eta \right\}) d^{N_{sdf}} \left\{ \eta \right\}
\label{Z}
\end{equation}

To use the  statistical mechanical formalism, one needs a set of elementary entities, whose summed volumes gives the volume of the system. The quadrons, discussed originally in \cite{BaBl02} and adapted to statistical mechanics in \cite{BlEd03,BlEd06}, are structural elements that fit this role nicely. They can be regarded as the 'quasi-particles' of this formalism\cite{BlEd03,BlEd06}. One also needs to identify the SDF. In terms of these, the volume function is  $\sum_q V_q\left(\left\{ \eta \right\}\right)$. It has been shown \cite{BlEd03,BlEd06}  that $N_{sdf}$ is equal to the number of  quadrons in 2D, but is larger in 3D\cite{Bl12}.

An explicit illustration of the quadron in 2D is shown in Figure \ref{fig:quadron}. First, one defines the centroids of grains $g$ and of cells $c$ as the mean position vectors of the contact points around them, respectively. Second, the contact points around every grain are connected to make polygons, whose edges are vectors, $\vec{r}^{gc}$, that circulate the grain in the clockwise direction. Third, one extends vectors $\vec{R}^{gc}$ from the centroids of grains $g$ to the centroids of the cells $c$ that surround them. A quadron is the quadrilateral whose diagonals are the vectors  $\vec{r}^{gc}$ and  $\vec{R}^{gc}$. Since every pair $gc$ corresponds exactly to one quadron, we shall denote a quadron by the index $q$ for  convenience. 
The quadron's shape is  quantified by the following structure tensor

\begin{equation}
C^q = \vec{r}^q\otimes\vec{R}^q
\label{QuadTensor}
\end{equation}
and its volume is
\begin{eqnarray}
V_q= \frac{1}{2}\mid\vec{r}^q \times \vec{R}^q\mid = \frac{1}{2} (r_x^q R_y^q - r_y^q R_x^q) 
\label{QuadVolume}
\end{eqnarray}

\begin{figure}[!bht]
\begin{minipage}[t]{0.3\textwidth}
\includegraphics[width=1.0\textwidth]{Fig01_quadron.eps}
\caption{Illustration of quadrons (shaded) and cells (yellow) in 2D.}
\label{fig:quadron}
\end{minipage}
\end{figure} 

It has been shown that, in 2D, there are as many quadrons as there are degrees of freedom \cite{BlEd03,BlEd06,Bl08a}, $N_q = N_{sdf}$. This makes it possible to use the quadron volumes to span the $N_q$-dimensional phase space, $\vec{V}$. 
We define the `density of states'  $g\left(\vec{V} \right)d^{N_q} \vec{V}=g\left(\left\{ V_q \right\}\right)\prod_q{d V_q}$ as the number of states whose volumes are between $\vec{V}$ and $\vec{V}+d\vec{V}$. 

The volume partition function can be then written as

\begin{equation}
Z_v = \int e^{-\sum_{j=1}^{N_q} V_q/X} g(\left\{ V_q \right\}) \prod_{j=1}^{N_q} dV_{q,j}
\label{PartFun}
\end{equation}
where the integration is understood as taken only over the systems that satisfy the constraints in $\Theta(\left\{ \eta \right\})$ and $W=\sum_{q=1}^{N_q} V_q $.

Let us now consider a model system - an ideal quadron gas (IQG) - where the quadrons are uncorrelated. The general partition function (\ref{PartFun}) becomes
\begin{equation}
Z^{iqg}_v = \left[ \int e^{-V_q/X} g_1(V_q) dV_q \right]^{N_q}
\label{PartFunIQG} 
\end{equation}
where $g_1(V)$ is a one-quadron density of states \cite{BlEd06}. 
The propability density function (PDF) of finding a quadron of volume between $V_q$ and $V_q+dV_q$ is
\begin{equation}
P(V_q) = \frac{1}{Z^{iqg}_v} e^{- V_q /X} g_1(V_q)
\label{PDFQuad}
\end{equation}

In the following, we investigate this PDF in our numerical experiments and relate it to the cell structures. The combination of the exponential term and the density of states, whose exact form is unknown, makes exact calculations impossible without additional information. Nevertheless, it is possible to gain insight from such a study, as we show in the following.


\section{Numerical experiments}

Our numerical experiments have been carried out using the Discrete Element Method (DEM) (\cite{Cundall1979,Mat-Chang-2011}). 
The method consists of using an incremental time marching scheme, wherein the motions of 2D discs (grains) are each computed using Newton's second law. We postulate a repelling harmonic interaction potential, characterized by normal and tangential spring constants, $k_n$ and $k_s$, activated upon contact and overlap between discs. We set $k_s/k_n=1/4$ in this study.

Our aim is to study polydisperse systems and we have chosen the distribution of the radii of the discs to be log-normal distribution due to its wide use in civil engineering and soil sciences \cite{Mitchell-Soga-2005},
\begin{eqnarray}
P (D) =\frac{1}{\sqrt{2 \pi} \sigma D} \exp \left( \frac{(\ln D - \ln D_0)^2}{2 \sigma^2}   \right)
\end{eqnarray}
where we assign $D_0=1.$ and $\sigma=0.2$, which give for the mean grain size and mode, respectively, $\bar{D}=1.02$ and $D_{mode}=0.961$ (see Figure \ref{fig:size_distrib}). 

\begin{figure}[!bht]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig02_size_distrib.eps}
\caption{Grain diameter probability density function.}
\label{fig:size_distrib}
\end{minipage}
\end{figure} 

The packing protocol of our systems is as follows. 
First, we construct three random packs in a double periodic domain, engineered to be on the verge of jamming. The packs consist of about $21400 \pm 1000$ discs and are made at packing fractions $\phi=0.76, 0.82$ and $0.84$. The small variation in the numbers of particles is due to the different densities of the three initial configurations.
These configurations are then used as initial states for the packing procedure; loose initial state (LIS), intermediated initial states (IIS), and dense initial state (DIS), respectively. 

Once an initial state is set, we assign all the particles the same friction coefficient, $\mu$, and apply to the system a slow isotropic stress $\sigma_c$ by changing the periodic length on both sides. 
We limit the applied stress to a level that corresponds to an average overlap discs of $\delta=\sigma_c / k_n =10^{-5}$.
No gravitational force is applied and the compression continues until the fluctuations of both grain positions (per mean grain diameter) and inter-granular forces (per mean average contact force) are below very small thresholds - $10^{-9}$ and $10^{-6}$, respectively. This procedure is carried out for each initial state at five different values of the inter-granular friction coefficients: $\mu=0.01, 0.1, 0.2, 0.5$ and $10$, giving altogether 15 assemblies.

Using this procedure, we have computed the packing fractions, the mean coordination numbers, and studied the structures in these systems. 
For the determination of $\bar{z}$, we disregard `rattlers', i.e. grains with one or no force-carrying contact. The results are shown in Figure \ref{fig:packing_fraction}. 
The upper and lower bounds of the coordination number $\bar{z}_{max}=4$  and $\bar{z}_{min}=3$ correspond to the isostatic states for discs with friction coefficients of $\mu=0$ and $\infty$, respectively.
The three systems with $\mu=0.01$ converge into states with ($\bar{z}$, $\phi$)=(3.941, 0.840),(3.944, 0.842) and (3.943, 0.843), respectively, which are very close to the ideal frictionless jammed state. However, with increasing friction coefficient, the difference between the values of $\bar{z}$ and $\phi$ in the final states of systems started from different initial states increases (see Figure  \ref{fig:packing_fraction}).

\begin{figure}["here"]
\begin{minipage}[t]{0.45\textwidth}
\includegraphics[width=1.0\textwidth]{Fig03_zg-phi.eps}
\caption{Coordination number vs. packing fraction for 15 systems, generated from three different intial states, LIS, IIS, and DIS. Note the convergence for $\mu\to 0$.}
\label{fig:packing_fraction}
\end{minipage}
\end{figure} 

Turning to the analysis of cell structures in these systems, consider the three examples shown in Figures \ref{fig:cell_structure1} - \ref{fig:cell_structure3}, constructed from LIS, with $\mu=10$ (\ref{fig:cell_structure1}) and $0.01$ (\ref{fig:cell_structure2}), and from DIS with $\mu=10$ (\ref{fig:cell_structure3}). 
While traces of the initial condition can be detected in both the first two systems, for higher friction one clearly ends up with typically larger cells and a correspondingly lower $\bar{z}$.  
Note that the packing fraction in Figures \ref{fig:cell_structure1} and \ref{fig:cell_structure3} are very similar (see also Figure \ref{fig:packing_fraction}), but the cells in  \ref{fig:cell_structure1} are noticeably bigger typically. 
Also note the considerable number of rattlers in the large cells, an issue that will be discussed in detail below. 
To quantify the structural differences, we next study the quadron volume distribution. 

\begin{figure}[!bht]
\begin{minipage}[t]{0.3\textwidth}
\includegraphics[width=1.0\textwidth]{Fig04_cell_loose_mu10.eps}
\caption{Example of an assembly with $\mu=10$, generated from LIS.}
\label{fig:cell_structure1}
\end{minipage}
\end{figure} 

\begin{figure}[!bht]
\begin{minipage}[t]{0.3\textwidth}
\includegraphics[width=1.0\textwidth]{Fig05_cell_loose_mu01.eps}
\caption{Example of an assembly with $\mu=0.01$, generated from LIS.}
\label{fig:cell_structure2}
\end{minipage}
\end{figure} 

\begin{figure}[!bht]
\begin{minipage}[t]{0.3\textwidth}
\includegraphics[width=1.0\textwidth]{Fig06_cell_dense_mu10.eps}
\caption{Example of an assembly with $\mu=10$, generated from DIS.}
\label{fig:cell_structure3}
\end{minipage}
\end{figure} 

 
\section{Structural characteristics}

\subsection{\label{sec4-1} Quadron volume statistics}

The grain size distribution is the same for all packs, making it convenient to normalize the quadron volumes by the mean grain volume $\bar{V}_g$, $v\equiv V_q/ \bar{V}_g$. The overall PDFs of quadron volumes are plotted in Figure \ref{fig:PDF of quadron} for all friction coefficients and for LIS and DIS. 
The corresponding means are plotted against the mean coordination numbers, for all initial states, in Figure \ref{fig:mean quadron volume vs zg}. All the points fall nicely on one curve {\it regardless of initial state}. 

\begin{figure}["here"]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig07_vqvg_all.eps}
\caption{PDF of the quadron volumes, normalized by the mean grain volume, $v=V_q/\bar{V}_g$ for all the 10 systems generated from LIS and DIS.}
\label{fig:PDF of quadron}
\end{minipage}
\end{figure} 

\begin{figure}["here"]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig08_vqvg-zg.eps}
\caption{Mean quadron volume $\bar{v}$ vs. mean coordination number $\bar{z}$ for all initial states.}
\label{fig:mean quadron volume vs zg}
\end{minipage}
\end{figure} 

On close inspection, this very weak dependence, if any, on initial state can also be observed in Figure \ref{fig:PDF of quadron}. 
This implies that the dependence of the relation between $\bar{v}$ and $\bar{z}$ on the initial state, seen in Figure  \ref{fig:packing_fraction}, can be made to disappear under the right parameterization. A clue to such a parameterization is the fact that the quadron volumes are computed {\it ignoring the rattlers}. 
Following this idea, we can derive the relation between $\bar{v}$ and $\bar{z}$ as follows

\begin{eqnarray}
\phi'=\frac{V_s'}{V}=\frac{\sum^{N_g'} V_g}{\sum^{N_q} V_q}
=\frac{N_g' \bar{V}_g}{N_q \bar{V}_q}
=\frac{1}{\bar{v} \bar{z}}
\label{eq:model03}
\end{eqnarray}

\begin{eqnarray}
\bar{v} = \frac{1}{ \phi' \bar{z} }= \frac{1}{(1-\eta_{r}) \phi \bar{z}}
\label{eq:mod_phi}
\end{eqnarray}
where $\eta_{r}=N_{r}/N$ is the rattlers fraction and $\phi'$ is the packing fraction after removing all the rattlers. 
This elegant result suggests that there is indeed a unique relationship between $\bar{z}$ and the packing fraction, but only if the latter is corrected for the rattler fraction. 
This is verified by a direct plot of these two quantities, Figure \ref{fig:mean zg vs. mod phi}.
It is also important to note that this relation is independent of initial state as Figure \ref{fig:mean zg vs. mod phi} shows.

\begin{figure}["here"]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig09_zg-phi_mod.eps}
\caption{Mean coordination number $\bar{z}$ vs. the rattlers-free packing fraction $\phi'$.}
\label{fig:mean zg vs. mod phi}
\end{minipage}
\end{figure} 

Comparing Figures \ref{fig:mean zg vs. mod phi} and \ref{fig:packing_fraction}, we see that the rattler-free packing fractions are much smaller than the packing fractions commonly reported in the literature. 
The removal of rattlers also makes sense on mechanical grounds, since it does not affect the force transmission in our systems (nor in any system experiencing no body forces).

We conclude that the seeming sensitivity of $\bar{z}(\phi)$ to the initial state (e.g. Figure \ref{fig:packing_fraction}), commonly seen in the literature,  stems directly from the variation in $\eta_{r}$. Plotting $\eta_{r}$ as a function of $\bar{z}$ in our systems (Figure \ref{fig:r_rattlers vs. mean zg}), we see that it is the small differences between the curves that gives rise to the observed differences in the conventional plots.

\begin{figure}["here"]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig10_Rrattlers-zg.eps}
\caption{ The rattlers fraction $\eta_{r}$ vs. the mean coordination number $\bar{z}$.}
\label{fig:r_rattlers vs. mean zg}
\end{minipage}
\end{figure} 

\begin{figure}["here"]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig11_vqvq_pdf.eps}
\caption{The PDF of the normalized quadron volume $u=V_q/\bar{V}_q$ for  for all the 10 systems generated from LIS and DIS.}
\label{fig:vq_PDF_normalized_by_mean_vq}
\end{minipage}
\end{figure} 

It is important to comment here that this result is not really universal. It is correct for the specific packing protocol used here. While we believe that, for any given protocol, the plots of $\bar{z}(\phi')$ would collapse onto one curve, there is no reason that this curve should be universal.  In other words, we expect different protocols to display different $\bar{z}(\phi')$ curves, probably depending on the rate of rattlers generation. 
However, these results do provide a universal protocol-independent insight - the relation between the $\bar{z}$ and the packing fraction is directly linked to the mechanical stability of the structure and the way that forces are transmitted. We will discuss this insight more in the concluding section.  

Turning to consider the PDFs of the quadron volumes more closely, it is natural to expect the mean quadron volume to increase with $\mu$, simply because the cells get bigger. To get insight into the shapes of the PDFs, we scale the quadron volumes by their means, $u\equiv  V_q / \bar{V}_q=v/\bar{v}$. This simple scaling is sufficient to collapse all the PDFs almost perfectly onto one curve for all the systems, independent of friction and initial state (Figure \ref{fig:vq_PDF_normalized_by_mean_vq}).  

Moreover, the collapsed PDF has an exponential tail, which  is a signature of the Boltzmann-like factor of equation \ref{PDFQuad}. 
Our best fit to the collapsed curve has the following form
\begin{eqnarray}
P(u) = f(u) e^{\alpha u}
\end{eqnarray}
where $\alpha$ is a slope of the tail and $f(u)$ is a rational function. Comparing this expression with equation \ref{PDFQuad}, we identify
\begin{eqnarray}
X=\bar{V}_q / \alpha
\end{eqnarray}
where we used $u=v/\bar{v}=V_q/\bar{V}_q$.
This implies that the compactivity $X$ is only a function of $\bar{V}_q$, in agreement with the equipartition principle derived recently in \cite{Bl12}.  
Intriguingly, combining this observation (or the equipartition principle of \cite{Bl12}) with equation (\ref{eq:mod_phi}), we find a expression relating directly the mean coordination number $\bar{z}$, the compactivity $X$ and the rattler-free packing fraction $\phi'$,

\begin{equation}
\frac{1}{ \phi' \bar{z} X }= \frac{\alpha}{\bar{V}_g}
\label{eq:EOS}
\end{equation}
where the right hand side is a constant that depends only on the grain size distribution.

Thus, the quadron description makes it possible to collapse the statistics of all the systems, given the correct normalisation. As such, it gives better insight into the general characteristics of granular packs - characteristics that are independent of both the inter-granular friction coefficient and the initial state. 
In the next section, we explore the reasons for the apparent unified nature of the quadron volumes statistics.
 

\subsection{\label{sec4-2} Cell order statistics}

To understand better the structural characteristics, we use a recently-proposed decomposition of the quadron volumes into conditional distributions \cite{Fretal08}

\begin{eqnarray}
P\left(v\right)=\sum_e Q(e) P\left(v\mid e\right)
\label{eq_decomposition}
\end{eqnarray}
\\
Here $Q(e)$ is the occurrence probability of cells of order $e$ and $P(v|e)$ is the conditional PDF of the normalized quadron volume, given that it belongs to a cell of order $e$.
$Q(e)$ is essential to the understanding of random granular packing \cite{BlToMaSoon} and it is this PDF that we wish to focus on next. 
As expected from Euler topological relation (below) and as can be observed from Figures \ref{fig:cell_structure1}-\ref{fig:cell_structure3}, a lower value of $\bar{z}$ must be accompanied with a higher mean value of $\bar{e}$ and corresponds to a lower packing fraction $\phi'$. 

\begin{figure}["here"]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig12_q_e.eps}
\caption{The cell order probability $Q(e)$ for  for the 10 systems generated from LIS and DIS.}
\label{fig:q(e)_dense}
\end{minipage}
\end{figure} 

In Figure \ref{fig:q(e)_dense} we plot $Q(e)$ for all the systems generated from LIS and DIS. The plots make evident two points. One is that the higher the inter-granular friction the larger the fraction of high-order cells. The other is  that $Q(e)$ is hardly dependent on the initial state for any $\mu$.

There is a direct relation between the mean cell order, $\bar{e}$, and the mean coordination number $\bar{z}$ and it can be derived from Euler's topological relation for a planar graph,

\begin{equation}
N_V - N_E +N_C = 1
\label{Euler}
\end{equation}
In this relation, $N_V$, $N_E$ and $N_C$ are, respectively, the numbers of the graph's vertices, edges and cells. In this relation we disregard the one large cell making the outside of the system, whose inclusion gives a topological characteristic 2, rather than 1, on the right hand side. This topological characteristic is negligible for sufficiently large systems. 
There are several ways to choose the graph in granular assemblies. One is to take the grain centres as vertices and the lines connecting touching grains as edges. This construction is convenient for discs, but can be awkward for grains of arbitrary shapes. Another, more generally applicable, construction is to take the contact points as vertices and the lines connecting contact points around grains as edges. This graph is consistent with the definition of the quadrons (Figure \ref{fig:quadron}). 
Both graphs give the same result for packs of $N(\gg1)$ grains

\begin{eqnarray}
\bar{e} = \frac{2 \bar{z}}{ \bar{z} -2} + O\left(\frac{1}{\sqrt{N}}\right)
\label{eq:mean_e}
\end{eqnarray}
where the rightmost term is a (negligible) boundary correction. Figure \ref{fig:Euler's relation} shows that the numerical systems satisfy this relation very well, showing that the sizes of the numerical pack are sufficiently large to ignore finite size effects.

\begin{figure}["here"]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig13_zg-e_relation.eps}
\caption{The mean cell order $\bar{e}$ vs. the mean coordination number $\bar{z}$.}
\label{fig:Euler's relation}
\end{minipage}
\end{figure} 

$Q(e)$ is sensitive to the value of $\bar{z}$, which is a function of $\mu$ (see Figure \ref{fig:q(e)_dense}), while $P\left(v\mid e\right)$ has been shown by Frenkel et al \cite{Fretal08}  (see also below) to be hardly dependent on $\mu$. This suggests that there may be a parameterization of $Q(e)$ that collapses all the curves corresponding to the collapse of $P\left(u^q \right)$ in Figure \ref{fig:vq_PDF_normalized_by_mean_vq}.

We find that all the curves collapse if we plot $\bar{e} e Q(e)$ as a function of $e'=(e-\bar{e})/\bar{e}^2$.
We note that the integrals of the curves in Figure \ref{fig:normalized_qe} are exactly 1. The collapsed curve appears to be fitted nicely by a truncated Gaussian

\begin{eqnarray}
\bar{e} e Q(e) = 
\frac{\sqrt{2/\pi}}{\sigma_r (1+erf(e'_{min}/\sqrt{2}\sigma_r))} 
exp\left(-e'^2/2\sigma_r^2\right)
\label{eq:gaussian}
\end{eqnarray}
with $erf(x)$ being the error function, $e'_{min}=0.076$ and $\sigma_r=0.082$. However, we note small deviation from the Gaussian form at the large $e$-tail, for which we currently have no explanation.

\begin{figure}[!bht]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig14_normalized_qe_gauss.eps}
\caption{The normalized cell probability $\bar{e}e Q(e)$. It is fitted well by a Gaussian form, except for small, but consistent, deviations in the large-$e$ tail (highlighted by a dash lined ellipse).}
\label{fig:normalized_qe}
\end{minipage}
\end{figure} 


\subsection{\label{sec4-3} Conditional quadron volume distributions}

To gain further insight into the universal properties of the structure, let us consider in more detail the conditional PDFs $P\left(v\mid e\right)$. 
These PDFs were studied by Frenkel et al. \cite{Fretal08}, who argued, on the basis of geometrical considerations, that they should be independent of inter-granular friction. Their argument was based on the observation that, given a collection of $N$ arbitrary grains, the number of ways to arrange $e$ grains into a cell of order $e$  depends only on the grains shapes and not on the inter-granular  friction. 
While our results seem to provide a support to this argument, a closer look shows a systematic $\mu$-dependence of $P\left(v\mid e\right)$ for $e> 6$ (Figures \ref{fig:cPDF_e5} and \ref{fig:cPDF_e6}). We can also see this effect in the behaviour of the mean of the conditional quadron volume as a function of $e$, $\bar{v}(e)=\sum_{q\in e}v P\left(v\mid e\right)$, shown in Figure \ref{fig:Mean Vq by cell order}. We include in the figure calculated values of quadron volumes of regular polygonal cells (RPC) of order $e$, 

\begin{eqnarray}
v^{RPC}(e)=\frac{V_q^{RPC}(e)}{\bar{V}_g}
=\frac{1}{\pi} cot \frac{\pi}{e}
\label{eq:regular_polygon}
\end{eqnarray}
which is clearly an upper bound for $\bar{v}(e)$. 
It is constructive to consider the ratio of the observed mean quadron volume to the regular polygon value, $\gamma(e)\equiv\bar{V}_q(e)/\bar{V}_q^{RPC}(e)$, shown in Figure \ref{fig:vq_vrp}. We observe that $\gamma(e)$, which is always below unity, has a minimum between $e=5$ and $e=6$ and that it decreases with $\mu$. These can be understood as follows. Carrying out a cell shape analysis (to be reported elsewhere), where cells were approximated to lowest order as ellipses, there is a strong correlation between $\gamma(e)$ and the mean aspect ratio of the ellipses for any given $e$.  
The implication is that the increase of $\gamma(e)$ with $\mu$ is due to the decrease of mechanical stability of elongated cells as $\mu$ decreases.
The behaviour for any particular $\mu$ depends then on two effects. As $e$ increases, more cell configurations can be realised, allowing for more elongated cells and this accounts for the initial decrease of $\gamma(e)$ with $e$.
However, the more elongated the cell, the less stable it is, leading to fewer elongated cells as $e$ increases. 
The competition between these two effects gives rise to the observed minimum.
To test this understanding, we use next an analytical model.

\begin{figure}[!bht]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig15_conditional_vq1.eps}
\caption{The conditional PDF of the quadron volumes for $e \le 5$ collapse nicely for all the systems.}
\label{fig:cPDF_e5}
\end{minipage}
\end{figure} 
\begin{figure}[!bht]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig16_conditional_vq2.eps}
\caption{The conditional PDF of the quadron volumes for $e = 6, 8$ collapse, but not as sharply as those for $e \le 5$ (Figure \ref{fig:cPDF_e5}).}
\label{fig:cPDF_e6}
\end{minipage}
\end{figure} 
\begin{figure}[!bht]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig17_mean_vq-mean_vg_for_cell.eps}
\caption{The variation of the mean conditional quadron volume, $\bar{v}$, with cell order $e$. The quadron volume of the regular polygon cell  (dashed line) form an upper bound.}
\label{fig:Mean Vq by cell order}
\end{minipage}
\end{figure} 

\begin{figure}[!bht]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig18_mean_vq_per_cell.eps}
\caption{The variation of the ratio of the mean quadron volume to regular polygon volume, $\gamma$ ($<1$), as a function of cell order $e$.}
\label{fig:vq_vrp}
\end{minipage}
\end{figure} 


\section{A free cell model}

We believe that the geometry-based argument of Frenkel et al \cite{Fretal08}, mentioned above, should be augmented with stability considerations. Indeed, the number of arrangements of $e$ grains, chosen from a given collection of $N$ grains, is independent of $\mu$ if we care only about the geometry. However, not all these shapes are stable under all possible combinations of compressive forces and it is this $\mu$-dependent stability that constraints the possible shapes. We demonstrate this effect, using a 'free cell model' (FCM). 
The model consists of a single cell, surrounded by $e$ discs, chosen from the same distribution of sizes as in the numerical simulation. We construct a large number of cells for different values of $\mu$ and analyse their shape distributions. 
The cell construction is as follows.
First, a random set of $e$ discs is generated from the size distribution. 
We place the first disc at the origin and and another one touching it along the $x$-axis.
We then generate a random angle ($\theta_2$ in Figure \ref{fig:free_cell}) and place the third disc in contact with the second, making sure that it does not touch the first disc. 
This process is repeated, making sure that each added disc makes contact only with the previous one,  until we reach the penultimate, $(e-1)$th, disc. 
We then determine the size and position of the last disc. If it fits in the remaining opening  then the cell is kept, otherwise it is discarded and we attempt to generate a new cell.
Figure \ref{fig:six_cells} shows examples of free cells, generated by this procedure.
The FCM allows us to generate quickly and efficiently a large number of cells and analyse them statistically.
By comparing this geometry-based model with the DEM results, we expect to identify the effects of mechanical stability for high order cells.

A comparison between the resulting PDFs and those from DEM simulation for all $e=3-8$ are shown in Figure \ref{fig:pdf_vq_freecell}. 
The PDFs are similar, but with consistent differences. 
The mean $v$ for different $e$ (Figure \ref{fig:vq_vrp})  are lower than of the numerical ones for $e>6$. 
This systematic difference is the direct consequence of the cell stability constraint, which has not been considered by Frenkel et al. The DEM cells are closer to regular polygons because they must be stable under compressive forces. The more frictional the grains, the more elongated cells that they can support. This is the reason for the minimum in $\gamma(e)$  decreasing and shifting to higher $e$ with increasing $\mu$.
 
\begin{figure}[!bht]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig23_free_cell.eps}
\caption{The construction of a cell configuration studied by the free cell model.}
\label{fig:free_cell}
\end{minipage}
\end{figure} 

\begin{figure}[!bht]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=0.8\textwidth]{Fig25_six_cells_s.eps}
\caption{ Examples of free cell configurations for $e=6$.}
\label{fig:six_cells}
\end{minipage}
\end{figure}

\begin{figure}[!bht]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig28_vq_free_cell.eps}
\caption{A comparison of the PDF of quadron volumes, generated both with the free cells model and with the DEM numerical experiments.}
\label{fig:pdf_vq_freecell}
\end{minipage}
\end{figure} 

\section{Conclusions and Discussions}


To conclude, we have studied in detail the effects of inter-granular friction and initial conditions on the structural characteristics of 2D granular assemblies. We have shown that structure can be analysed quantitatively, using the quadron description\cite{BaBl02,BlEd03}. In particular, we have established that a number of distributions collapse onto single curves, pointing to a `universal'-like structural characteristics that are independent not only of initial state, but also of the friction coefficient. These are significant results, whose implications are discussed at the end of this section.

The raw quadron volume distribution, which is key for the statistical mechanics of granular matter\cite{BlEd03,BlEd06,HiBl12}, changes systematically towards larger values with increasing friction coefficient $\mu$ and hence with decreasing mean coordination number $\bar{z}$. However, when normalised by the mean quadron volume, it collapses nicely onto a single curve. Following the insight proposed by Frenkel et al.\cite{Fretal08}, we have traced this to a similar collapse of the conditional PDFs $P(V_q/\bar{V}_g\mid e)$ for every cell order $e$, except for small, but insightful, deviations for large $e$, which are discussed below.

The common wisdom in the soil mechanics literature is that the relation between the the mean coordination number, $\bar{z}$, and the packing fraction $\phi$ depends on the initial state, as shown in Figure \ref{fig:packing_fraction}.
However, prompted by the observation that the collapse appears once the rattlers are disregarded in the quadron partition of the granular space, we have examined this relation as a function of rattlers-free packing fraction $\phi'$. We discover that, with this correction, the dependence on the initial states disappears and all the relations collapse onto a single curve, $\bar{z}(\phi')$. 
This is consistent with our observation that the plot of $\bar{V}_q$ as a function of $\bar{z}$ is also independent of the initial state.
Since rattlers do not support forces in our systems, these results carry a fundamental significance - they suggest that the packing fraction, the mean coordination number, the quadron volume distribution and generally the structure are linked directly to the manner in which stresses are transmitted.
While this gives insight into the structure of the system, it should be remembered that it is difficult to estimate the packing fraction of rattlers in real experiments and it is unclear at this stage how to use this insight for engineering applications.

Considerations of mechanical stability also suggest that the argument proposed by Frenkel at el.\cite{Fretal08}, leading to independence of the conditional quadron volume distributions of inter-granular friction, needs to be modified. Their argument assumes that cell shapes are independent of mechanics and therefore that arranging an $e$-sided cell depends only on the grain shape distribution. This, in turn, eliminates friction as a factor in the occurence probability of cell structures. However, as friction decreases,  elongated cells are less stable mechanically, reducing the number of available cell shapes. Indeed, we can observe a small increase in the occurrence probability of large quadron volumes in $P(V_q/\bar{V}_g\mid e)$ for $e>6$ as $\mu$ increases (Figure \ref{fig:cPDF_e6}). 

This effect is seen most clearly when plotting the ratio of the mean quadron volume to that of a regular polygon as a function of $e$ (Figure \ref{fig:vq_vrp}). The initial drop reflects the departure of the cell from a regular polygon due to the increased number of possible geometrical configurations, while the subsequent increase shows the effect of mechanical stability in limiting such a departure. The more pronounced increase for lower $\mu$ is evidence for our argument above.

To  illustrate the effect of mechanical stability, we have constructed a free cell model (FCM) that takes into account only geometric considerations and disregards mechanical stability. Indeed, this model shows the same initial decrease of the above ratio with $e$ and none of the subsequent increase (Figure \ref{fig:vq_vrp}). 
Thus, neglecting the mechanical stability effect, the FCM provides a low bound for the mean quadron volume, as illustrated in Figure \ref{fig:VqvsZg}.

\begin{figure}["here"]
\begin{minipage}[t]{0.4\textwidth}
\includegraphics[width=1.0\textwidth]{Fig30_discussion.eps}
\caption{Simple prediction of the mean quadron volume in terms of the mean coordination number.}
\label{fig:VqvsZg}
\end{minipage}
\end{figure} 

We have further found that the raw cell order PDF, $Q(e)$, changes systematically with $\mu$, but hardly at all with the initial state. 
However, all the PDFs also collapse to a single curve under the  transformations $Q(e)\to\bar{e}eQ(e)$ and $e\to\left(e-\bar{e}\right)/\bar{e}^2$.  
We note that the collapse cannot be perfect at the large $e$ tail due to the aforementioned constraint of mechanical stability. This friction-specific behaviour, which is evidenced by Figures \ref{fig:cPDF_e6} and \ref{fig:vq_vrp}, is hardly noticed at all in Figure \ref{fig:normalized_qe}, demonstrating that the effect of mechanical stability is small compared to the geometrical effect discussed by Frenkel et al. \cite{Fretal08}.

Our results are significant for several reasons. 
First, they are a step towards a systematic understanding of structural organisation of granular matter in response to a specific packing procedure. In particular, they suggest that the structural characteristics of granular matter are best understood via the statistics of cell configurations. 
Second, they make it possible to progress on the statistical mechanics of granular matter, where knowledge of the quadron volumes distribution is key. 
Third, the results shed light on the common wisdom in the soil mechanics and civil engineering communities that the initial state affects the final relation between $\bar{z}$ and $\phi$. 

Our findings show that this relation is directly linked to mechanical stability of the packing backbone - for any specific pack generation process, the initial state affects the fraction of rattlers, which do not participate in the stress transmission. Once the rattlers are disregarded, there emerges a unique relation between $\bar{z}$ and the corrected rattlers-free packing fraction $\phi'$.
Fourth, our analysis has a significant ramification for the famous packing problem. Different packing procedures lead to different structures. For example, procedures not constrained by mechanical stability would give rise to different $Q(e)$'s than those that are. Specifically, the former would tend to have lower packing fractions due to a higher occurrence probability of large cells, which are mechanically unstable. This reinforces the need for a more accurate definition of the packing problem, which we discuss elsewhere\cite{BlToMaSoon}.

Taken together, these results show that the packing process has a fundamental effect on the structural characteristics while the initial conditions and the inter-granular friction are details that only modify this effect and can be scaled away. The collapsed distributions show how to predict such modifications.
It should be emphasised that the collapses, which we have found, do not imply universality in the traditional sense because the specific forms of the collapsed curves depend on the packing procedure. 
For example, it has been demonstrated that, using a different packing protocol, it is possible to generate packs dominated by cells of order 3 \cite{HiBl12}, which cannot be scaled to match the ones we obtained here.  There is in principle an infinite number of possible procedures to generate packs, depending on an astronomically large number of parameters. It would be impossible to map all the possible processes to a manageable parameter space. Each such procedure gives rise to its own characteristic $Q(e)$ and it would be impossible to collapse all of these $Q(e)$ onto one curve. 

To fully understand how structural characteristics depend on the packing procedure, one needs a good model of both the geometric effects and the limitations that mechanical stability poses on the cell shape distribution. We are currently developing such a model and will report results in a later publication.
Intriguingly, the relevance of both volumetric and stress effects is reminiscent of the inter-dependence demonstrated recently for the statistical mechanical understanding of granular matter \cite{Bl12} and it would be interesting to find out whether the two issues are related or not.  





\nocite{*}

 