Communication: Structure characterization of hard sphere packings in amorphous and crystalline states

The hard-sphere particle interaction limit is a tremendously versatile physical model, being widely used for structural studies of liquids,¹ glasses,^2,3 colloids,⁴ granular materials,⁵ and many others.^6,7 Its relevance in such a variety of physical systems suggests that many macroscopic properties arise by the fundamental fact of impenetrability of the systems constituents.⁸ The ultimate goal then becomes to establish relations between physical properties and the geometry of the arrangement of hard bodies in two or three-dimensions. Many decades of research, heavily driven by numerical experiments,² have led to various geometrical structure variables, with different levels of success in either uniquely characterizing each state or in deriving macroscopic physical properties from them.⁹ 

Our study is motivated by the structural phase transitions observed in molecular fluids and also replicated in hard-sphere systems under compression.¹⁰ We also consider analogous phenomena observed in granular materials, where the hard-sphere approximation is commonly used to successfully model complex rheological behaviours.⁵ As the volume fraction is increased, hard-spheres enter an entropy minimization driven phase where glass formation competes with the nucleation and growth of the crystalline phase.¹¹ Hard-sphere models are known to successfully reproduce the main structural properties of these states for various physical systems, either for crystallization,^12,13 or the amorphous solid phase transition.¹⁴ One of the main reasons for using the hard-sphere model over classical condensed matter systems is that simulations – as also colloidal suspensions experiments, which are a very good approximations of the hard-sphere model^13,15 – have particle-size spatial resolution, and thus the statics and dynamics can be studied from a microscopic perspective. Furthermore, due to the lack of long-range or non-binary interactions, and the simple geometry of the constituents, hard-sphere models are theoretically tractable.^12,16,17

Inspired by the highly ordered and easily describable crystalline phase, many researchers have searched for inherent geometrical relations in disordered (amorphous) packings. The straightforward approach is to analyze the static structure factor, which is a direct measure of the local microstructure of particles.¹⁸ The pair correlation function g(r) is also a popular quantity for the analysis of non-crystalline materials.^19,20 The problem with such quantities is that the detailed three-dimensional information is lost as a result of statistical averages, as also by considering only pairs of particles. In particular, they do not provide much information about the topology of the local structures in the particle-size scale, which are believed to distinguish different kinds of amorphous arrangements. It therefore becomes highly significant to exploit some other methods of three-dimensional characterization of these structures. Many attempts have already been made to quantify local or long-range ordering, providing further characterization of disordered packings.^21–31 

In the following we analyze local arrangements of particles recognizing the importance of caging and voids in the overall structure and properties of the arrangement. Here we extend the previous studies^22,23,26 by considering all particle triples in the particles neighbourhood that do not include the central particle. The distribution of voids allows us to clearly distinguish between the different structural phases, as also between different kinds of crystals and the relative number of each specific ordering, in systems presenting partial or many types of crystallization. Our analysis shows that in amorphous states there is a preferred local structure of four co-planar particles.

We consider systems of N non-overlapping spheres arranged in a three-dimensional cubic space of volume V. The spheres are located at positions x_i and have radii r. Our main parameter is the sphere volume fraction, defined as

To determine the neighbors of a particle, we first compute the weighted Delaunay triangulation of the set of points corresponding to the centers of the particles, {x_i}.³² Neighboring particles are then defined as those particles connected by the edges of the triangulation. For each particle we consider every possible combination of three neighbors, that is, all possible triangles that can be formed by the centers of any three of its neighbors. We refer to these triangles as neighbor-triangles. Notice that the neighbor-triangles do not contain the central particle. We then proceed to quantify the overall mobility of the particle by defining all channels through which the particle can pass. A channel is defined as the circle-area in the plane of a neighbor-triangle through which the central particle could move. This is computed by considering the Apollonius circle, i.e., the circle which is simultaneously tangent to all other three circles defined by the projection of the three spheres in the neighbor-triangle plane. There are at most eight possible Apollonius circles for each case, which are obtained analytically by solving a system of three quadratic equations.³³ From the set of eight possible solutions we choose the one which corresponds to the circle that does not contain any particle center of the neighbor-triangle, as it is the only one that corresponds to our definition of channel.⁴³ The radius of the respective channel is then defined as the radius of this circle, R_j, as shown in Fig. 1.

Having obtained R_j for all neighbor-triangles of every particle, we then compute the normalized probability distribution function of (scaled) channel sizes, f(R_j/r). The ratio R_j/r is calculated for all neighbor triples j with channel size R_j of every particle. Note that R_j/r has a direct physical interpretation, as less than unity corresponds to a closed channel, while greater than unity corresponds to an open channel, through which eventually the central particle could escape. Furthermore, the function f is well defined for spheres with any size distribution, since the radius of the central particle is scaled out. We analyze both the individual structure of f as also its evolution with volume fraction, for various particle systems.

In order to refine our definition of channels, we consider f(R_j/r) for very low volume fractions, where no structure is expected (see Fig. 2). When considering all triangles the distribution presents a recognizable wide tail structure, but after excluding from the distribution the channel sizes that correspond to non-acute neighbor-triangles, i.e., those where one of the angles is greater than 90°, the distribution becomes Gaussian, with high accuracy over three orders of magnitude.⁴⁴ The exclusion of non-acute triangles makes physical sense considering that channels defined by them cannot block the central particle, thus conflicting with our initial definition of a channel. For the rest of the analysis, non-acute triangles are never considered.

We use an event-driven molecular dynamics algorithm, as it is fundamentally suited for the simulation of hard spheres systems. The number of particles is by default N = 16³ = 4096, unless stated otherwise. Given the large amount of possible neighbor-triangles for each particle, the statistical significance rapidly increases with the number of particles in the system. We observed that 4096 particles were adequate, as increasing the number of particles did not produce any noticeable change in any of the results. Periodic boundary conditions are imposed to mimic an infinite system, i.e., a statistically homogeneous medium.

Starting from zero volume fraction, we compress the system towards a jammed state using a modification of the Lubachevsky-Stillinger algorithm,^34,35 which allows the radius of the particles to grow linearly in time with a dimensionless rate Γ.⁴⁵ The kinetic energy, E, is kept constant using a re-scaling thermostat procedure.^36,37

If the growing is sufficiently slow, Γ < 0.0007,³⁸ the monodisperse system stays in a gas-fluid state in approximate equilibrium during the densification phase, and exhibits a fluid-solid transition (crystallization) for volume fractions between ν_f ≈ 0.492 (freezing point) and ν_m ≈ 0.543 (crystal melting point). For infinitely slow compressions, it is expected that the system finally reaches a stable solid (crystalline) phase with close-packing fraction ν_cp ≈ 0.7405, corresponding to face-centered close packing. In our simulations, due to finite compression rates, we reach a crystalline phase with defects and different local arrangements, and packing fractions up to ν ≈ 0.73. This corresponds to a thermodynamically stable branch in the hard sphere phase diagram.¹⁰ On the other hand, for fast compression rates the system enters a metastable state for ν > ν_m, which extrapolates continuously from the fluid branch and is conjectured to end at some random close packing state, around ν_rcp ≈ 0.64, the interpretation of which is beyond the scope of this study, as its value depends on the details of the procedure.³⁹ 

We now observe the evolution of f(R_j/r) with ν for fast and slow compression rates.

As the volume fraction increases, the distribution of channel radii fundamentally changes, see Fig. 3(a). The distribution changes to non-Gaussian for fluid densities above ν ≈ 0.15–0.25. We speculate that this change corresponds to the percolation gas-to-fluid transition observed by Woodcock⁴⁰ at similar packing fractions, although we did not investigate this in detail. As the volume fraction increases, two smooth humps continuously grow, that at higher ν > 0.5 evolves into two well defined peaks, centered above R_j/r ≈ 0.15 and near R_j/r = 1. These values can be understood in terms of the geometry of the local arrangements: R_j/r ≈ 0.15 ideally corresponds to the channel size expected for three touching equal spheres, and thus the appearance and growth of this peak shows the appearance of triples in contact as well as the relative importance of density fluctuations. It is also the absolutely smallest possible channel size for equally sized spheres. The peak at unity, on the other hand, is obtained for three particles lying on the same plane with the central particle and practically touching it, i.e., when the channel essentially coincides with the central particle; we confirmed that the majority of particles corresponding to the peak at unity are indeed practically touching the central particle.

As expected for very slow compression, Γ = 16 × 10⁻⁶, the system exhibits (partial) crystallization near the melting point ν_m ≈ 0.54; crystallization at the freezing point is kinetically suppressed.⁴¹ The distribution is able to capture the crystallization transition by the development of two new peaks, at R_j/r ≈ 0.4 and R_j/r ≈ 1.4. This transition is shown with more detail in the inset of Fig. 3(a). The new peak at R_j/r ≈ 0.4 corresponds to a square crystalline arrangement. On the other hand, the peak at R_j/r ≈ 1.4 groups several distinct arrangements which can be distinguished as ν further increases, as shown below. This was confirmed for many runs with different initial particle velocities and positions.

Finally, as the maximum volume fraction is reached (i.e., as the pressure diverges), the distribution is mostly dominated by steep peaks, see Fig. 3(b). In analogy with Bragg peaks from common diffraction techniques, these peaks can be traced to the crystal structures present in the particles’ arrangement. Figure 4 shows f(R_j/r) for perfect FCC, HCP, BCC, and SC crystals of about 10³ particles, as a reference. They all present significant differences, which allows to distinguish between the types of crystals. A detailed discussion of these peaks is beyond the scope of this study.

Let us now focus on the case of fast compression, Γ = 16 × 10⁻³, for which crystallization is not happening. The channel size distributions for different volume fractions are shown in Figure 5. For volume fractions below ν ≈ 0.2, there is no appreciable difference with the slow compression case described in Subsection IV A. For values well above ν ≈ 0.2, the same peaks at R_j/r ≈ 0.15 and unity can be observed, but no other peaks are developed. Contrary to the crystallisation path at the melting point, f(R_j/r) does not show sharp signs of the glass transition at ν ≈ 0.57. The lack of defined peaks in the distribution, as in the previous slow compression cases, is a clear signal that the system remains, to very high degree, amorphous. Nevertheless, there exists partial order, as suggested by the high values of the peak at unity. While the peak at 0.15 can be interpreted in the same way as in the slow compression rate case, the peak at unity, on the other hand, is not as easily interpretable. While present in both crystalline and glassy configurations, its relative importance and shape are considerably different. The unity value corresponds to those configurations where three neighboring particles all touch the central particle and lie in the same plane with it. There are of course other cases where the value of R_j/r could be one. By looking at the distribution of distance of each neighbor-triangle particle to the center particle, we confirm that by far the most common case is when the three particles are indeed touching the central one. This suggests the existence of “hidden” local order in random close packings that cannot be easily measured by order parameters because such local planes, corresponding to R_j/r ≈ 1, are not oriented with respect to each other as in a periodic crystal structure. We have no explanation for this preference of the system structure, suggesting a direction of future research. Furthermore, we have confirmed that the same qualitative features are observed in polydisperse systems.⁴² 

Let us now take a look at the structure of f for values close to the minimal channel size, shown in the inset of Figure 5. The high values between the two ideal cases, corresponding to the three touching spheres (R_j/r ≈ 0.15) and a square arrangement (R_j/r ≈ 0.41), signal a significant presence of “intermediate” configurations. The drop at R_j/r ≈ 0.41 is due to exclusion of non-acute neighbor-triangles from the statistics. A similar distribution of channel sizes, up to R_j/r ≈ 0.41, was obtained in the studies of interstitial holes in random close packings of spheres.^22,23,26 They concluded that for monodisperse packings the spread in the channel sizes between 0.15 and 0.41 cannot be reduced to lead to a single distribution of channel sizes allied with mechanical stability.

The analysis of channel size distributions was shown to be able to distinguishing between gaseous, liquid, partially and fully crystallized, and glassy (random) jammed states. Unlike the usually computed pair-distribution functions or structure factors, the channel size distribution is highly sensitive to changes in volume fraction, and presents unique features for each phase. States of partial crystallization can be recognized and characterized by the development and position of specific peaks, which can be traced to specific crystalline configurations, and could be used to quantify the degree of crystallization of the system. On the other hand, we confirm that random glassy configurations of isotropically (rapidly) grown hard particle systems present a common structural feature, as shown by looking at the channel size distributions. The overpopulation of many three neighboring particles lying in the same plane as the central particle, almost touching it, could be considered a first microscopic trace of crystals in a plane. As different planes are not oriented with dominant relative angles as in a crystal, there is no appreciable global three dimensional ordering. Finally, we remark that our analysis should be easily extendable to other particle systems even in different dynamic regimes. Further work on the behaviour of the distribution near the state transitions could lead to deeper relations with previously known distribution functions or thermodynamic variables.

We would like to thank J. Finney and L. Woodcock for helpful discussions. This research is supported by the Dutch Technology Foundation STW, which is the applied science division of NWO, and the Technology Programme of the Ministry of Economic Affairs, Project No. STW-MUST 10120.