Structural analysis of disordered dimer packings

Abstract

Jammed disordered packings of non-spherical particles show significant variation in the packing density as a function of particle shape for a given packing protocol. Rotationally symmetric elongated shapes such as ellipsoids, spherocylinders, and dimers, e.g., pack significantly denser than spheres over a narrow range of aspect ratios, exhibiting a characteristic peak at aspect ratios of α_max ≈ 1.4–1.5. However, the structural features that underlie this non-monotonic behaviour in the packing density are unknown. Here, we study disordered packings of frictionless dimers in three dimensions generated by a gravitational pouring protocol in LAMMPS. Focusing on the characteristics of contacts as well as orientational and translational order metrics, we identify a number of structural features that accompany the formation of maximally dense packings as the dimer aspect ratio α is varied from the spherical limit. Our results highlight that dimer packings undergo significant structural changes as α increases up to α_max manifest in the reorganisation of the contact configurations between neighbouring dimers, increasing nematic order, and decreasing local translational order. Remarkably, for α > α_max our metrics remain largely unchanged, indicating that the peak in the packing density is related to the interplay of structural rearrangements for α < α_max and subsequent excluded volume effects with unchanged structure for α > α_max.

---

1 Introduction

Jammed disordered particle packings have been used as a model to understand the structures of liquid crystals, glasses, self-assembly of nanoparticles, biological systems and granular media.¹ While there has been considerable recent progress in our understanding of jammed sphere packings,² the effect of particle shape on the properties of jammed packings has been much less explored.³ Considering one of the simplest macroscopic observables of packings – the packing density – one finds that many non-spherical shapes pack denser than spheres, which achieve maximal packing densities of ϕ_j ≈ 0.64 for a wide range of packing protocols (although denser packings can also be achieved for specific protocols, see the discussion in ref. 4). For example, many polyhedra,^5–9 ellipsoids,^10–13 spherocylinders,^14–20 and dimers,^21,22 as well as irregular shapes such as those composed of a number of overlapping spheres^23,24 achieve packing densities ϕ_j ≥ 0.7, with the densest disordered packings so far found for tetrahedra at ϕ_j ≈ 0.78.⁵ Plotting the packing density as a function of a continuous shape descriptor, such as the aspect ratio α (for rotationally symmetric elongated shapes), exhibits a non-monotonic behaviour with a peak at α ≈ 1.4–1.5 for ellipsoids, spherocylinders, and dimers, with some variations due to the packing protocol. For larger aspect ratios, the packing density decreases, following, e.g., an approximate scaling behaviour ϕ_j ∼ 1/α for spherocylinders.²⁵

In this study, we revisit dimer packings simulated with the MD platform LAMMPS using a gravitational pouring protocol. Our goal is to identify structural features that characterize the peak in the packing density by focusing on details of the contact statistics as well as positional and orientational order metrics. In this context, it is important to emphasize the role of the protocol in the packing generation. The interplay between the packing density and the degree of order that arises by tuning the protocol parameters has been widely discussed for spheres, most notably in the critique of the well-posedness of the concept of “random close packing”.²⁶ For non-spherical particles, the protocol dependence is manifest in the relatively large variance of results reported for ϕ_j for the same shape, e.g., for spherocylinders.^14–20 Our viewpoint is thus to focus on packings generated by a specific protocol, namely the widely used pouring under gravity, and understand how shape variation changes the structural features of these packings.

Previous studies of ordering effects in random packings of elongated particles obtained inconsistent results, which might be due to different protocols and boundary conditions used. For example, simulations of prolate ellipsoids by pouring into a container under gravity found considerable nematic order, whereby the ellipsoids’ symmetry axes (the semi-major axes) tend to lie within the plane normal to the gravity direction.^10,11,27 This ordering effect has been explained as a result of the particles' tendency to minimize the gravitational potential energy.¹⁰ On the other hand, simulations that compress or inflate the non-spherical particles from an initial random state such as the Lubachevsky–Stillinger algorithm (applied to ellipsoids^12,28) or a mechanical contraction algorithm (applied to spherocylinders^29–31) do not find any significant order as is also observed with other geometric simulation methods.^32,33 While 3D experiments of ellipsoids¹³ and elongated colloids³⁴ did not observe any signatures of order, experiments of asymmetric dumbbells in 2D showed strong orientational correlations between neighbours due to mutual restrictions on positions.³⁵ The order characteristics of dimers in 3D have so far not been investigated to our knowledge.

The remainder of this article is organized as follows: in Section 2 we present the details of our simulation method with LAMMPS. In Section 3, we present results on our analysis of the packing fraction, contact number, and orientational/positional order metrics. Finally, we conclude in Section 4 with a discussion of our results.

2 Simulation method

Disordered packings of monodisperse frictionless symmetric dimers in three dimensions are generated with the molecular dynamics platform LAMMPS.^36,37 The dimers are obtained by overlapping two identical spheres with diameter d and mass m. We study dimers with aspect ratios α in the range 1.0005 ≤ α ≤ 2, where α is given as the ratio of the length over the width, see Fig. 1. In this packing protocol, N = 12 000–15 000 monodisperse dimers are poured under gravity into a three-dimensional box of side length ≈20d. The lateral (x̂–ŷ-plane) boundary conditions are chosen to be periodic and the box is bounded in the ẑ-direction by a rough surface at the bottom (implemented by the “fix wall/gran hertz/history” command). During a simulation run, a gravitational force acts on the dimers in the ẑ-direction. The pouring protocol makes use of LAMMPS’ “fix pour” command, which repeatedly inserts particles into the simulation box every few timesteps within a specified insertion region 30–40d above the bottom and releases them until N particles have been added overall. In the insertion region, particles are added with random positions and orientations and without any overlap. Particles are only inserted again after the previously inserted particles have fallen out of the insertion region under the gravitational force.

Fig. 1 Dimer shape defined by the aspect ratio: (a) α = 1.05, (b) α = 1.4, (c) α = 2.

LAMMPS treats a dimer defined by a fixed distance between its two constituent spheres as an independent rigid body (implemented by the “fix rigid/small” command). The total force and torque on each dimer rigid body are computed as the sum of the forces and torques on its constituent spheres in every time step. The coordinates, velocities, and orientations of the constituent spheres are then updated so that the dimer moves and rotates as a single entity.

LAMMPS can natively implement different models for calculating the contact forces between the spheres. In this study, a Hookean model is chosen because of its convenience to dissipate residual kinetic energy and hence to reach a static state quickly.³⁸ In the Hookean model, when two spheres i and j having positions r_i and r_j, respectively, are in contact, they experience a relative normal compression with overlap δ = d − r_ij, where r_ij = r_i − r_j and r_ij = |r_ij|. The resulting force is F_ij = Fⁿ_ij + F^t_ij, where F^n,t_ij are the normal and tangential contact forces, respectively, given as:³⁸

Here, n_ij = r_ij/r_ij, v_n,t are the normal and the tangential components of the relative velocity of the spheres i and j, and K_n,t and γ_n,t are the elastic and viscoelastic constants, respectively. The quantity Δs_t denotes the elastic tangential displacement between the spheres.³⁸ The total force F^tot_i on sphere i in a gravitational field g = −gẑ is then given as:where the sum runs over all j spheres in contact with sphere i.

Throughout the investigation we set our basic units as d = 1, m = π/6, and g = 1. Distances, times, velocities, forces and elastic constants are then measured in units of d, , , mg, mg/d, respectively. We generally use a normal spring constant K_n = 2 × 10⁵mg/d unless otherwise indicated. Additionally, we simulate also harder dimers with K_n = 2 × 10⁶mg/d and softer ones with K_n = 2 × 10⁴mg/d to examine the effect of particle hardness on the contact number of the dimers at small aspect ratios. We set γ_t = 0 and the remaining parameters used are given in Table 1. The choice of most of these values follows the discussion in ref. 39.

Table 1 Material parameter values and time step Δt used in the simulations

K n(mg/d)	K t/Kn	γ n (mg/d)	Δt
2 × 10^4	2/7	15	0.003
2 × 10^5	2/7	50	0.001
2 × 10^6	2/7	150	0.0003

---

Simulations are run until the system reaches a static equilibrium when the kinetic energy per particle is less than 10⁻⁸mgd for small K_n and up to three orders of magnitude less for large K_n. For example, when K_n = 2 × 10⁵mg/d the simulation takes 3–8 × 10⁶Δt to reach equilibrium, which depends on the chosen aspect ratio and also on the random initial configurations when particles are poured into the container. For further details of any of the LAMMPS commands used, we refer the reader to the LAMMPS documentation.³⁷

3 Structural analysis

3.1 Packing fraction

We calculate the packing fraction of the dimer packings for various aspect ratios. The packing density is determined for the bulk region shown in Fig. 2. The particles within 5–8d from the container floor are excluded from the bulk region since they can be highly crystallized. The thickness of this crystallized region depends on many factors such as the box dimension and the pouring height. Excluding the particles within 5–8d provides results that are largely unaffected by the crystallization. The particles within 5d from the upper-most particles have also been excluded from the bulk because their Voronoi volumes can not be decided accurately due to deficiencies in their neighbourhood.

Fig. 2 The bulk region shown in the x̂–ẑ-plane.

In order to determine the packing density in the bulk region, we calculate the Voronoi volume of each dimer in the bulk, which is defined as the space that is closer to the surface of a given dimer than to that of any other dimer. While a formal parametrization of the Voronoi volume of a dimer is analytically tractable,⁴⁰ a straightforward computational method makes use of LAMMPS’ built-in routine to determine the Voronoi volume of the individual spheres in the packing using a conventional Voronoi tessellation. The Voronoi volume W_i of a dimer is then found by summing the Voronoi volumes of its two constituent spheres. The bulk volume V_b occupied by N_b dimers in the bulk is calculated as . We then obtain the packing fraction as ϕ_j = N_bV_α/V_b where V_α is the volume of a dimer with aspect ratio α. The volume V_α is found by subtracting the overlap volume from the sum of its constituent sphere volumes. The overlap volume contains two equal spherical caps whose volume can be calculated exactly, see Appendix A. Note that a dimer is considered to be part of the bulk region only if the centres of both constituent spheres are within the bulk. All average quantities discussed in the following are calculated for dimers in the bulk only.

We plot the packing fraction ϕ_j of the dimers as a function of the aspect ratio α in Fig. 3. As can be seen from Fig. 3, the packing fraction ϕ_j has a non-monotonic relationship with α, i.e., it increases as α increases until reaching a peak at ϕ_j = 0.707 for α = α_max = 1.4, beyond that it decreases. These results are in agreement with previous studies^21,22 and also show reasonably good agreement with results from a mean-field calculation,⁴⁰ shown in Fig. 3. Systematic deviations between our simulations and the mean field theory are in particular visible in the behaviour for larger aspect ratios α > 1.5, which are likely due to the strong mean-field assumptions. In fact, the mean-field theory relies on a reduction of higher-order positional correlations to pair correlations and also neglects orientational correlations between particles. The latter become more significant for particles of larger aspect ratios, see Section 3.3.1.

Fig. 3 The packing fraction ϕ_j as a function of the dimer aspect ratio α. Simulation values of ϕ_j are shown averaged over 10 independent simulation runs for α ≥ 1.1 (dots), and for a single run for α < 1.1 (diamonds).

3.2 Contact and coordination numbers

We introduce the contact number z as the average number of contact points of a dimer and the coordination number z_c as the average number of neighbours of a dimer, whereby a neighbour is defined as another dimer with which at least one contact point is shared. While z = z_c for smooth convex shapes like spheres, ellipsoids, and spherocylinders, z ≥ z_c for concave shapes like dimers, since two particles can share more than one contact point. In general, two dimers A and B share a contact point if the separation vector of two spheres i and j, with sphere i in dimer A and sphere j in dimer B, satisfies r_ij ≤ d, which can be detected with high numerical precision. Two dimers can thus share up to four different contact points. Due to the soft interaction potential the contact “point” is strictly a small overlap region, which creates some complications at small dimer aspect ratios, see below.

In Fig. 4(a) we show the behaviour of z_c as a function of α and the associated distributions of z_c for a set of aspect ratios. We observe a smooth increase of z_c(α) for α > 1 with a maximum at z_c = 8.34 for α = 1.5 followed by a slight decay. The qualitative behaviour is in line with the results of ref. 22, where dimer packings were generated using an energy minimization protocol, although our values of z_c are consistently larger over the range of aspect ratios. The distributions P(z_c) are approximately symmetric and Gaussian (Fig. 4(a), inset).

Fig. 4 (a) The coordination number z_cvs. α and distributions P(z_c) for various aspect ratios (inset). (b) The contact number z vs. α and distributions P(z) (inset). The values of z_c and z are shown averaged over 10 independent simulation runs for α ≥ 1.1 and α = 1 (dots), and for a single run for 1 < α < 1.1 (diamonds).

On the other hand, the contact number z does not exhibit such a smooth increase, see Fig. 4(b). First establishing the baseline for sphere packings at α = 1 with our protocol, we find that z = 6.14 for spheres. This value is slightly above the isostatic value of z = 2d_f = 6, where d_f denotes the degrees of freedom of a particle, generally found for disordered sphere packings using a variety of packing protocols.⁴ We suspect that this difference is due to the gravitational packing protocol and the interaction potential with non-zero softness, see also the comparable values found in the studies of sphere packings^21,39 using a similar protocol. Deforming spheres into dimers, the smallest aspect ratio of dimers for which we are able to report the contact number reliably is α = 1.05, for which we find z = 10.39. For larger aspect ratios, z decreases slightly, but then remains unchanged at z = 10.28 for α > 1.2. The difference with the isostatic value z = 2d_f = 10 is approximately of the same magnitude as the difference for spheres using our packing protocol. By comparison, the studies in ref. 22, 41 and 42 find that dimers are almost exactly isostatic, which is thus in line with our findings. The observation of a constant z for all aspect ratios of dimers is an important difference with the behaviour of convex elongated shapes such as ellipsoids and spherocylinders, which are hypostatic (z < 2d_f) at small aspect ratios and show a smooth increase upon shape deformation from the sphere like the coordination number z_c here.^28,30

We highlight that for very small aspect ratios α ∈ (1,1.05) the calculation of z is unreliable, since our particle model leads to incorrect contact detections: the overlap regions due to the particle softness can extend far enough into the dimer as to create a contact with an interior sphere as illustrated in Fig. 5.

Fig. 5 Illustrations of “double” and “cusp” contacts shown in 2D as discussed in ref. 41. (a) Double contact: the yellow sphere is embedded into the red dimer so deeply that it contacts both red spheres. (b) Cusp contact: the yellow sphere contacts both red spheres by covering the cusp point (black point) of the red dimer.

Such problematic contact configurations for dimers were also identified in the recent work by Shiraishi et al.^22,41 and separated into “double” and “cusp” contacts, see Fig. 5. Shiraishi et al. investigated the contact number of dimer packings using a compression protocol with soft particle interactions for various packing densities ϕ. For large enough values of the excess packing density Δϕ = ϕ − ϕ_j, where ϕ_j denotes the packing density at jamming onset, “double” and “cusp” contacts were observed. In their analysis, these contacts could thus be avoided by setting an upper limit for Δϕ at each aspect ratio studied and they observed that this upper limit approaches zero as α → 1. In our case, the occurrence of these configurations depends on the stiffness value K_n as shown in Fig. 6, where it can be seen that the threshold aspect ratio, at which double and cusp contacts occur, is shifted to smaller aspect ratios for larger K_n. For any value of K_n, double and cusp contacts will occur at sufficiently small aspect ratios and thus the contact number very close to the sphere shape can not be reliably established. For K_n = 2 × 10⁵ we see that double and cusp contacts do not occur for α ≥ 1.05, which is the lower limit of α used in our contact number analysis.

Fig. 6 The fraction of double (solid lines) and cusp contacts (dashed lines) in the dimer packings for small α and three normal spring constants K_n.

In order to refine our analysis of the packing microstructure, we define five distinct contact configurations according to the number of contact points that are shared by two neighbouring dimers, see Table 2. Excluding the regime α ∈ [1,1.05), we determine how the fraction of each configuration type changes as a function of α, see Fig. 7. We see that even though the average number of contacts z is approximately constant over this range of α, the underlying contact configurations change significantly with α. Most notably, the two most common contact configurations, Type 1 and Type 2, increase and decrease, respectively, as α increases up to around α_max and remain approximately unchanged for α > α_max. The remaining contact configurations confirm this trend, showing the strongest variations in the regime α < α_max. Overall, we see that contact configurations, in which spheres of neighbouring dimers only have one contact point (Type 1 and Type 3) increase, while those with multiple contact points (Types 2, 4 and 5) decrease as the packing becomes denser up to the packing density peak at α_max. This trend is somewhat counter-intuitive, since the Types 2, 4 and 5 configurations correspond to more optimal local arrangements between two dimers, which locally reduce the packing density. Similar results for the fractions of these five configuration types have been found for packings of shapes composed of four overlapping spheres.⁴³

Table 2 Five distinct contact configurations of two dimers. We show illustrations for aspect ratios α = 1.2 and α = 2. The total number of contact points for each type is: one (Type 1), two (Types 2 and 3), three (Type 4), four (Type 5)

α	Type 1	Type 2	Type 3	Type 4	Type 5
1.2
2

---

Fig. 7 The fractions of the five contact configuration types of Table 2 for packings of dimers with different α.

Rather than excluding the aspect ratio regime where the problematic double and cusp contacts occur it might be tempting to re-assign such contacts and thus infer the properties of the small aspect ratio regime in an ad hoc way. For example, a double contact as in Fig. 5(a), which creates two overlaps of sphere pairs and is thus counted as two contact points, could be counted as only one, effectively ignoring the incorrect overlap with the interior sphere. This can be done likewise for other contact configurations, which require a careful consideration of the relative position and orientation of the overlapping dimer pair, see the full discussion in Appendix D. Re-assigning contacts in this way leads to a rapid but smooth decrease of z to the corresponding value of spheres z ≈ 6 as α → 1 (Fig. 21), but also exhibits seemingly unphysical behaviour, such as sharp peaks in the fractions of the Types 1–5 contact configurations around α ≈ 1.05, i.e., at the aspect ratio where double and cusp contacts start to occur (Fig. 22).

3.3 Order metrics

We employ several order metrics to measure global and local ordering in the dimer packings at various aspect ratios. The nematic orientational order parameter and the orientational pair correlation function are used to evaluate orientational ordering. Translational ordering is investigated with bond orientational order parameters, the radial distribution function and bond angle distributions. All calculations are made for the particles within the bulk volume so as to discard the crystallized region observed at the bottom of the container.

3.3.1 Metrics for orientational order. The nematic orientational order parameter S has traditionally been applied to identify different ordered phases of liquid crystals by characterising the average molecular orientation.⁴⁴S is defined as:where is the second Legendre polynomial and β_i the angle between the orientation of dimer i and the so-called director, which specifies the average orientation of the particles. The dimer orientation is described by the unit vector u⁽ⁱ⁾ measured along the dimer's long axis.

We apply this parameter to the dimer packings to quantify the global orientational order. When all u⁽ⁱ⁾ are randomly oriented, S = 0, while if all u⁽ⁱ⁾ are oriented in a plane normal to the director, S = −0.5, which corresponds to a perfect oblate phase. When all u⁽ⁱ⁾ are aligned with the director, we have perfect nematic order with S = 1.

In order to determine the director and S, we first evaluate the tensor Ω defined as:

Denoting by λ_max the eigenvalue of Ω with the largest absolute value, we identify the director as the eigenvector corresponding to λ_max. For all aspect ratios, we find that the director is aligned with the ẑ-axis (gravity direction). We then obtain S directly as:

S = λ_max. (5)

We also determine the orientational pair correlation function S₂ in order to quantify local ordered structures at a radial distance r from a reference particle. S₂ is calculated as:

where cos β_ij = u⁽ⁱ⁾·u^(j) and n_i(r) denotes the set of particles in a spherical shell of width Δ(r) = 0.025d at a distance r from the centre of dimer i in the bulk. The expression |n_i(r)| refers to the size (cardinality) of the set n_i(r). We note that the spherical shell considered in S₂ can extend into the boundary region beyond the bulk and thus include particles in partially crystallized regions, although the effect on the average should be small. In general, due to the non-periodic boundary conditions in the ẑ-direction our packings are not rotationally invariant and thus the restriction to a radial coordinate is only an approximation.

We present the dependence of S and S₂(r) on the aspect ratio α in Fig. 8. We see that S changes rapidly as α increases from the sphere value, reaching its minimum at around α_max and remaining approximately constant for α > α_max, in line with the behaviour of z_c and the different contact types. Interestingly, the behaviour of S(α) as α → 1 appears almost singular, but the range of values is not sufficient to identify a clear power-law. The minimum of S at ≈−0.16 indicates slight oblate ordering, where the dimers' long axes are oriented close to the horizontal plane normal to the direction of gravity. This ordering is thus in agreement with that observed in simulation studies of prolate ellipsoids using also pouring under gravity.^10,11,27 In order to compare the magnitude of the orientational ordering with these studies, we also calculated the order parameter χ used in ref. 10, 11 and 27, which is defined in eqn (17). We find a maximum of χ ≈ 0.32 for α = 1.4. By comparison, in ref. 10 the maximum is χ ≈ 0.4 for α ≈ 1.5, while ref. 11 and 27 find χ ≈ 0.25 and χ ≈ 0.5, respectively, for α ≈ 1.5. Note that in ref. 11 and 27, χ monotonically increases upon further elongation over the observed range of aspect ratios.

Fig. 8 The nematic orientational order parameter S vs. the aspect ratio α. Values of S are shown averaged over 10 independent simulation runs for α ≥ 1.1 (dots), and for a single run for α < 1.1 (diamonds). Inset: The orientational pair correlation function S₂vs. r/d for various aspect ratios.

The plot of S₂ in Fig. 8(inset) demonstrates how orientational correlations become more long-range for larger aspect ratios. For small α, correlations decay rapidly within the first coordination shell, while for large α oscillations in S₂ are visible over the whole range of r/d, which is here limited by r/d = 5, i.e., the width of the boundary region on top of the bulk region that restricts the maximum radius of the spherical shell used in eqn (6).

3.3.2 Bond orientational order parameters. The bond-orientational order metrics q_l and Q_l introduced by Steinhardt et al.⁴⁵ have most commonly been used to quantify translational order in disordered packings of spherical particles.^46–51 While Q_l is widely accepted as a well-defined parameter to measure global ordering in a packing, it has been suggested that the local order parameter q_l needs more caution to reliably identify local crystalline structures in these systems.^52,53 It was assumed that higher values of q₆ are associated with higher degrees of order⁴⁶ and averages 〈q₆〉 have been used to quantify the overall degree of order for disordered sphere packings.⁴⁸ However, it has been found that some local configurations of disordered sphere packings that are clearly non-crystalline have exhibited the same values of q₆ as hcp or fcc crystals.⁵² Therefore, in this study, we use recently introduced local order parameters defined by Eslami et al.⁵⁴ to improve the accuracy of determining local translational order in the dimer packings.

Steinhardt et al.⁴⁵ associated with every bond joining a particle and its neighbours a set of spherical harmonics:

where the Y_lm are spherical harmonics and θ_ij, ϕ_ij denote the polar and azimuthal angles which define the orientation of the vector (bond) pointing from the reference particle i to another particle j, see Fig. 9. NN(i) contains the set of neighbour indices for particle i, which are defined as those particles j that have at least one contact with i.

Fig. 9 Parametrization of the separation vector (bond vector) r_ij = r_j − r_i connecting the reference particle i (red) with j (yellow). The definitions of the polar and azimuthal angles, θ_ij and ϕ_ij, respectively, are indicated.

The local orientational order parameter q_l(i) of particle i is then defined as the following rotational invariant combination of q_lm:

Moreover, the global orientational order parameter Q_l is defined as

where

Recently, Eslami et al. introduced the local order parameters to improve the determination of liquid and different crystallized phases.⁵⁴ Starting from the q_lm of eqn (7), we first determine

where is the complex conjugate of q̂_lm(j) and q̂_lm(i) is defined as follows:

Then the order parameters are obtained by averaging over the first coordination shell of particle i:

The advantage of over q_l is that they can distinguish the liquid phase and different crystalline phases in a more accurate way.⁵⁴ They indicate in fact the correlation between the order in the first and the second coordination shell of a reference particle.⁵⁴ It has been observed that is large ≈1 for crystalline phases, while assumes values close to zero for disordered (liquid) phases, which thus allows to easily discriminate between such phases. On the other hand, the values of are sensitive to the crystal type, so is able to distinguish bcc, fcc, and hcp crystals.

We display the pairs for each dimer in the bulk region of the packing in Fig. 10 for various aspect ratios. By comparing these results to empirical data for liquid, bcc, hcp, and fcc phases of Lennard-Jones particles from ref. 54, we observe that the distributions at large aspect ratios (α > 1.4) are quite clearly in a liquid phase where −0.05 < < 0.3 and 0 < < 0.4. As the aspect ratio decreases, the region occupied by and expands and approaches the region occupied by the bcc/hcp crystal phases indicating the presence of a large proportion of dimers exhibiting some local translational order intermediate between a liquid and bcc/hcp crystalline order.

Fig. 10 The local order parameters and defined in eqn (13). Every data point corresponds to a dimer in the bulk region of the packing. The sketched regions for bcc, hcp, fcc, and liquid phases of Lennard-Jones particles are taken from ref. 54.

We also calculate the averages , and compare their values with the global order parameters Q₄, Q₆ for different aspect ratios, see Fig. 11. While Q₄ is close to zero for all aspect ratios, there is a slight increase in Q₆ for α < 1.4 implying some global ordering at small aspect ratios. In line with the observations in Fig. 10, we see that both and are non-zero and monotonically decreasing as α increases, whereby varies over a larger range than . For small aspect ratios, both averages are considerably larger than the corresponding averages of a fluid phase, which were determined as and in ref. 54. Overall, we observe that at large aspect ratios the packing is more translationally disordered than at small aspect ratios.

Fig. 11 The global bond orientational order parameters Q₄, Q₆ and the averages and vs. α. By comparison, and for the liquid phase of Lennard-Jones particles.⁵⁴

3.3.3 Radial distribution function. We calculate the radial distribution function g(r) to further examine the translational correlations between the dimers. The radial distribution function of the bulk dimers is determined aswhere n_i(r) denotes the set of particles in a spherical shell of width Δ(r) = 0.025d at a distance r from the centre of dimer i in the bulk, ρ is the particle number density, and V_shell(r) is the volume of the shell. As discussed for the orientational correlation function S₂(r), eqn (6), the restriction to a radial coordinate is only an approximation due to the fact the our packings are not rotationally invariant. As before the spherical shell can extend into the boundary region beyond the bulk. We plot g(r) as a function of r/d for various aspect ratios in Fig. 12. We see that for small aspect ratios g(r) exhibits the characteristic shape of sphere packings with a main peak at r/d = 1 and a split second-peak at r/d ≈ 1.7 and r/d ≈ 2.^{14,29,32,39,55} For larger aspect ratios, these sharp peaks broaden and reduce in height. These results are consistent with the variation of the bond orientational order parameters with the aspect ratio discussed above, where elongation in the dimers results in a reduction of translational correlations.

Fig. 12 The radial distribution function g(r) of the dimer packings, eqn (14), for different α. Inset: Enlargement of the regime r/d ∈ [1.125,3].

3.3.4 Bond angle distribution. We measure the probability for a dimer to have a contact at a particular direction relative to its long axis. For each dimer pair i, j, we determine the polar angle θ_ij and the azimuthal angle ϕ_ij of the bond vector r_ij = r_j − r_i in the reference frame of particle i, see Fig. 9. The probability density functions (PDFs) of θ_ij and ϕ_ij are shown for various aspect ratios in Fig. 13. It can be clearly seen from Fig. 13 that at small aspect ratios dimers have primarily contacts at θ_ij = 90°. As the aspect ratio increases, the narrow band around 90° widens and finally almost disappears at α = 2. For small aspect ratios, there are also symmetric secondary peaks visible at θ_ij = 30° and θ_ij = 150°, with all contacts occurring within the range θ_ij ∈ [30°,150°] up to α ≈ 1.4.

Fig. 13 PDFs of the polar and azimuthal angles θ_ij, ϕ_ij of the bond vectors r_ij for all neighbour pairs i, j and different aspect ratios.

To get a better insight into the origin of these structures, the PDFs of θ_ij, ϕ_ij are further refined according to the contact configuration type between neighbouring dimers, see Fig. 14–16. For aspect ratio α = 1.05 (Fig. 14), we see that for Types 2–5 only configurations with θ_ij ≈ 90° are possible due to the geometric constraint of these configuration types. The structure observed in the overall bond diagram at very small aspect ratios (Fig. 13a and b) is thus primarily due to Type 1 configurations and the peak at θ_ij ≈ 90°. For larger aspect ratios α = α_max = 1.4 and α = 2, the bands for Types 2–4 widen due to the increase in possible relative orientations that still satisfy the contact constraint (see Fig. 15 and 16). This excludes Type 5 configurations which are available only in a narrow range of possible polar angles by definition. As expected, Type 1 configurations with only a single contact point between neighbours, which thus least constrains the relative orientations, exhibit a wide band of possible polar angles at all aspect ratios, see Fig. 14(a), 15(a) and 16(a). Interestingly, this band still exhibits some structure, with a main peak at θ_ij = 90° and symmetric secondary peaks at θ_ij = 30° and θ_ij = 150° for both α = 1.05 and α = 1.4, which disappear for α = 2.

Fig. 14 PDFs of the polar and azimuthal angles θ_ij, ϕ_ij of the bond vectors r_ij for all neighbour pairs i, j with a specific contact type. Aspect ratio: α = 1.05.

Fig. 15 PDFs of the polar and azimuthal angles θ_ij, ϕ_ij of the bond vectors r_ij for all neighbour pairs i, j with a specific contact type. Aspect ratio: α = α_max = 1.4.

Fig. 16 PDFs of the polar and azimuthal angles θ_ij, ϕ_ij of the bond vectors r_ij for all neighbour pairs i, j with a specific contact type. Aspect ratio: α = 2.

4 Conclusions

One of the main results of our study is the identification of structural features that accompany the formation of the peak in the packing density of jammed dimer particles. In particular, we find that (i) the coordination number z_c; (ii) the fractions of Types 1–4 contact configurations; and (iii) the nematic order parameter S undergo rapid changes upon deforming spheres into dimers with aspect ratios up to α ≈ α_max, while further elongation of the dimers leaves these metrics largely unchanged. This highlights that the peak in the packing density of Fig. 3 arises due to microscopic re-arrangements up to α≈ α_max and subsequent excluded volume effects: the contact configurations remain statistically unchanged for α > α_max, but since the particles are longer the packing can sustain more empty space while being mechanically stable, in line with the phenomenological description of spherocylinder packings using the random contact equation, which predicts a decay ϕ_j ∼ 1/α.⁵⁶

Dimers are a convenient shape model, because their contact interactions can be easily implemented by overlapping spheres. As such they represent one of the simplest non-spherical and concave shapes. However, our analysis shows that such a particle model does not allow to resolve the contact configurations at very small aspect ratios when interactions are not truly hard. As such we are not able to probe in our simulations, e.g., the analytical predictions from effective medium theory on the contact number scaling for very small shape deformations.⁵⁷ The problematic double and cusp contacts should generally occur for shapes composed of overlapping (soft) spheres as used, e.g., in the optimization studies of ref. 23 and 24, which might prevent a detailed analysis of the contact properties of such simulated packings.

Our investigation highlights the competition between orientational and translational order as a result of elongation. While the translational order is larger for small aspect ratios, the elongation induces the dimers to have both more orientationally ordered local structures (with slight global oblate ordering) and less translational order akin to those of a liquid. Dimers at large aspect ratios thus exhibit structures that resemble a liquid crystal in terms of these metrics. Importantly, the structural features identified here might be specific to the gravitational packing protocol used and might not occur in dimer packings obtained with other packing methods such as energy minimization from a random initial configuration.²² Nevertheless, due to the simplicity of the protocol, which is also relevant in many real world scenarios, we expect our results to be significant to understand the packing density and structural properties of real granular matter composed of non-spherical particles.

Conflicts of interest

There are no conflicts to declare.

Appendix

A Calculation of the dimer volume

The overlap volume of the two constituent spheres of a dimer contains two equal spherical caps which lie above/below the plane through the cusp points at the dimer's centre, see Fig. 17. The volume of a spherical cap V_cap of height h is found as:where R is the sphere radius. The dimer volume V_α is then calculated by subtracting the overlap volume from the sum of its constituent sphere volumes as:

V_α = 2V_sphere − 2V_cap (16)

Fig. 17 The overlap volume of a dimer contains two equal spherical caps of height h (coloured in yellow).

B Algorithm for the identification of double and cusp contacts

Double and cusp contacts are identified by checking if there is any overlap between the circle enclosing the cusp on the dimer surface and a contacting sphere of its neighbouring dimer, see Fig. 18(a). This circle with centre c_c, radius r_c and unit normal w and a sphere with centre c_s, radius r_s are shown in Fig. 18(b). The next steps are followed for the identification:

Fig. 18 Detecting double and cusp contacts. (a) First, it is checked if there is any overlap between the black circle (dashed) enclosing cusp located on the yellow dimer's surface and the contacting sphere of the red dimer. If there is an overlap between the circle and the sphere, it is identified as a cusp contact. (b) 3D visualization of the circle and sphere interaction, it is determined if the plane of the circle cuts the sphere or not. (c) If the plane of the circle cuts the sphere, it forms a new circle (red) and then it is checked if there is overlap between the original circle and the new red circle.

1. The distance d_cs = |w·(c_c − c_s)| between the plane of the circle and the sphere's centre is calculated to check if the plane cuts the sphere or not. If d_cs > r_s then there is no intersection, so the plane passes above/below the sphere entirely.

2. If there is an intersection, i.e., d_cs < r_s, it will be between the original circle and a new one formed where this plane meets the sphere, with centre c_p = c_s + d_csw.

3. If d_cs = r_s then this is the sole point of intersection with the plane, otherwise a new circle with radius r_p occurs as displayed in Fig. 18(c), where . Then, the problem has been reduced to a circle–circle interaction.

4. If |c_p − c_c| < r_c + r_p, then there is overlap between the circle and the sphere, so the contact is identified as a cusp contact. If there is no overlap, then the contact is either a double contact or a Type 2 configuration.

5. To distinguish a double and a Type 2 configuration, two vectors v₁ and v₂ from the contacting sphere's centre to the centres of the constituting spheres of the reference dimer are determined as illustrated in Fig. 19. The projections of these two vectors onto the unit normal w of the circle enclosing cusp are determined and the directions of these projections are checked. If both of them have the same direction, the contact is identified as a double contact, otherwise it is regarded as a Type 2 configuration.

Fig. 19 Two vectors v₁ and v₂ from the contacting red sphere's centre to the centres of the constituting spheres of the yellow dimer are determined. The projections of these two vectors onto the unit normal w of the circle enclosing cusp are determined and the directions of these projections are checked. (a) If both of them have the same direction, it is identified as a double contact (b) otherwise it is regarded as Type 2 configuration.

C The order parameter χ

In ref. 10 the following order parameter has been introduced to measure the orientational order of prolate ellipsoidswhere β_i is the angle between the semi-major (long) axis of particle i and the ẑ-axis (gravity direction). Since the director identified with the Q-tensor in Section 3.3.1 is also aligned with the ẑ-axis, the expression for χ is the same as that for S, eqn (3), apart from the shift −π/2 in the argument of P₂. The parameter χ of eqn (17) thus takes values in the interval [−2,1]: when all particles are randomly oriented, χ = 0, while if all particles' long axes are oriented in the horizontal plane normal to the gravity direction χ = 1. When the long axes of particles are oriented along the gravity direction we have χ = −2. A plot of χ as a function of α for our dimer packing data is shown in Fig. 20.

Fig. 20 The orientational order parameter χ vs. α. Values of χ are shown averaged over 10 independent simulation runs for α ≥ 1.1 (dots), and for a single run for α < 1.1 (diamonds).

D Mapping between different contact configuration types

We introduce a heuristic method to re-assign configurations with double and cusp contacts to one of the Types 1, 2, and 4 configurations. The precise mapping depends on the number and the location of double and cusp contacts as summarized in Table 3. In general, double contacts are mapped to one contact point and cusp contacts to two. For Type 3 configurations, no double or cusp contacts have been found. For Type 5 configurations, two cusp contacts do occur, which leave the configuration as Type 5 after the mapping.

Table 3 Two-dimensional illustrations of configurations with double and cusp contacts. These configurations are re-assigned to Types 1, 2, and 4 as indicated in the table

Configuration type	Re-assigned configuration type
	Type 1	Type 2	Type 4
Type 2
	A double contact is counted as one contact point: two contact points are reduced to one.
Type 4
	Two overlapping double contacts are counted as one contact point: three contact points are reduced to one.	One double and one cusp contact (cusp 2 overlaps with the red sphere) are counted as two contact points: three contact points are reduced to two.
Type 5
	Two overlapping double contacts are counted as one contact point: four contact points are reduced to one.	Two distinct double contacts are counted as two contact points: four contact points are reduced to two.	One double contact (cusp 1 is not covered by one of the yellow spheres) and one cusp contact (cusp 1 overlaps with the other yellow sphere) are counted as three contact points: four contact points are reduced to three.

---

With this mapping, we count a smaller number of contact points and thus the average number of contacts z decreases. In fact, we obtain a rapid but smooth decrease of z as α → 1, whereby z approaches the corresponding value of spheres (Fig. 21). Resolving the contact counting by Types 1–5 configurations, we see that, as expected, the fraction of Type 1 configurations now increases for α < 1.05, while the fractions of Types 2, 4 and 5 configurations decreases in the same regime (Fig. 22). In fact, the adjusted counting of contact points leads to sharp peaks at α ≈ 1.05, i.e., at the aspect ratio at which double and cusp contacts start to occur, that appear unphysical.

Fig. 21 A double-logarithmic plot of z − z_svs. α − 1 for three different normal spring constants K_n. We define z_s as the contact number of the corresponding sphere packing, which approaches the isostatic value z_s = 6 as the particle hardness increases.

Fig. 22 The fractions of the contact configurations of Types 1–5 vs. the aspect ratio α. For α ≥ 1.05 the data shown is the same as in Fig. 7, but in the regime α < 1.05 (dashed lines) the contact counting has been adjusted by re-assigning configurations with double and cusp contacts to Types 1, 2, and 4 configurations as summarized in Table 3.

Acknowledgements

This research utilised Queen Mary's Apocrita HPC facility, supported by QMUL Research-IT http://doi.org/10.5281/zenodo.438045. We acknowledge the assistance of the ITS Research team at Queen Mary University of London.

Notes and references

1. S. Torquato and F. H. Stillinger, Rev. Mod. Phys., 2010, 82, 2633.
2. P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani and F. Zamponi, Annu. Rev. Condens. Matter Phys., 2017, 8, 265.
3. A. Baule and H. A. Makse, Soft Matter, 2014, 10, 4423.
4. A. Baule, F. Morone, H. J. Herrmann and H. A. Makse, Rev. Mod. Phys., 2018, 90, 015006.
5. A. Haji-Akbari, M. Engel, A. S. Keys, X. Zheng, R. G. Petschek, P. Palffy-Muhoray and S. C. Glotzer, Nature, 2009, 462, 773.
6. Y. Jiao and S. Torquato, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 041309.
7. P. Damasceno, M. Engel and S. Glotzer, Science, 2012, 337, 453.
8. R. Shepherd, J. Conrad, T. Sabuwala, G. Gioia and J. A. Lewis, Soft Matter, 2012, 8, 4795.
9. L. Liu, Z. Li, Y. Jiao and S. Li, Soft Matter, 2017, 13, 748.
10. B. J. Buchalter and R. M. Bradley, Europhys. Lett., 1994, 26, 159.
11. G. Delaney, J. Hilton and P. Cleary, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 051305.
12. A. Donev, I. Cisse, D. Sachs, E. Variano, F. Stillinger, R. Connelly, S. Torquato and P. Chaikin, Science, 2004, 303, 990.
13. W. Man, A. Donev, F. Stillinger, M. Sullivan, W. Russel, D. Heeger, S. Inati, S. Torquato and P. Chaikin, Phys. Rev. Lett., 2005, 94, 198001.
14. S. R. Williams and A. P. Philipse, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 051301.
15. J. Zhao, S. Li, R. Zou and A. Yu, Soft Matter, 2012, 8, 1003.
16. C. Abreu, F. Tavares and M. Castier, Powder Technol., 2003, 134, 167.
17. X. Jia, M. Gan, R. A. Williams and D. Rhodes, Powder Technol., 2007, 174, 10.
18. M. Bargiel, Computational Science-ICCS2008, 2008, vol. 5102, p. 126.
19. A. Wouterse, S. Luding and A. P. Philipse, Granular Matter, 2009, 11, 169.
20. A. V. Kyrylyuk, M. A. van de Haar, L. Rossi, A. Wouterse and A. P. Philipse, Soft Matter, 2011, 7, 1671.
21. S. Faure, A. Lefebvre-Lepot and B. Semin, Esaim: Proceedings, 2009, vol. 28, p. 13.
22. K. Shiraishi, H. Mizuno and A. Ikeda, J. Phys. Soc. Jpn., 2020, 89, 074603.
23. M. Z. Miskin and H. Jaeger, Soft Matter, 2014, 10, 3708.
24. L. Roth and H. Jaeger, Soft Matter, 2016, 12, 1107.
25. A. Philipse, Langmuir, 1996, 12, 1127.
26. S. Torquato, T. M. Truskett and P. G. Debenedetti, Phys. Rev. Lett., 2000, 84, 2064.
27. J. Gan and A. Yu, Powder Technol., 2020, 361, 424.
28. A. Donev, R. Connelly, F. Stillinger and S. Torquato, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 051304.
29. A. Wouterse, S. R. Williams and A. Philipse, J. Phys.: Condens. Matter, 2007, 19, 406215.
30. A. Wouterse, S. Luding and A. Philipse, Granular Matter, 2009, 11, 169.
31. C. Ferreiro-Córdova and J. S. V. Duijneveldt, J. Chem. Eng. Data, 2014, 59, 3055.
32. J. Zhao, S. Li, R. Zou and A. Yu, Soft Matter, 2012, 8, 1003.
33. L. Meng, Y. Jiao and S. Li, Powder Technol., 2016, 292, 176.
34. S. Sacanna, L. Rossi, A. Wouterse and A. Philipse, J. Phys.: Condens. Matter, 2007, 19, 376108.
35. Y. Han and M. Kim, Soft Matter, 2012, 8, 9015.
36. S. Plimpton, J. Comput. Phys., 1995, 117, 1.
37. LAMMPS, http://lammps.sandia.gov.
38. L. Silbert, D. Ertas, G. S. Grest, T. Halsey, D. Levine and S. Plimpton, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 64, 051302.
39. L. Silbert, D. Ertas, G. S. Grest, T. Halsey and D. Levine, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 031304.
40. A. Baule, R. Mari, L. Bo, L. Portal and H. Makse, Nat. Commun., 2013, 4, 2194.
41. K. Shiraishi, H. Mizuno and A. Ikeda, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2019, 100, 012606.
42. C. F. Schreck, N. Xu and C. S. O'Hern, Soft Matter, 2010, 6, 2960.
43. E. Azéma, F. Radja, B. Saint-Cyr, J.-Y. Delenne and P. Sornay, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 052205.
44. P. G. de Gennes and J. Prost, The Physics of Liquid Crystal, Clarendon Press, 2nd edn, 1995.
45. P. Steinhardt, D. Nelson and M. Ronchetti, Phys. Rev. B: Condens. Matter Mater. Phys., 1983, 28, 784.
46. A. Kansal, S. Torquato and F. Stillinger, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 041109.
47. T. Aste, M. Saadatfar and T. Senden, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 061302.
48. K. Lochmann, A. Anikeenko, A. Elsner, N. N. Medvedev and D. Stoyan, Eur. Phys. J. B, 2006, 53, 67.
49. A. Wouterse and A. Philipse, J. Chem. Phys., 2006, 125, 194709.
50. Y. Jin and H. Makse, Phys. A, 2010, 389, 5362.
51. W. Xu, Z.-Y. Sun and L. An, Eur. Phys. J. E: Soft Matter Biol. Phys., 2010, 31, 377.
52. S. Kapfer, W. Mickel, K. Mecke and G. Schröder-Turk, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 030301.
53. W. Mickel, S. Kapfer, G. Schröder-Turk and K. Mecke, J. Chem. Phys., 2013, 138, 044501.
54. H. Eslami, P. Sedaghat and F. Müller-Plathe, Phys. Chem. Chem. Phys., 2018, 20, 27059.
55. A. Donev, S. Torquato and F. Stillinger, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 011105.
56. A. Philipse, Langmuir, 1996, 12, 1127.
57. C. Brito, H. Ikeda, P. Urbani, M. Wyart and F. Zamponi, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 11736.

---