Disordered packings of binary mixtures of dimer particles

Disordered packings of non-spherical particles and their mixtures are abundant in nature, but have so far attracted only few systematic studies. Previous investigations of binary mixtures of specific convex shapes have established two generic properties: (i) the existence of a unique density maximum when shape or mixture composition of the two species are varied; (ii) the validity of an ideal mixing law indicating that the packing density is independent of the segregation state. These findings were so far only observed for mixtures of convex particles such as spherocylinders, ellipsoids, and spheres. Here, we investigate the packing properties of binary mixtures of frictionless dimer particles simulated by a gravitational pouring protocol in LAMMPS. Our results demonstrate the validity of (i, ii) also for such packings of non-convex particles. Moreover, we investigate the contact statistics of these packings to elucidate the microstructural features that underlie (i, ii). Our results show that the contact number per species also satisfies a simple mixing law and that similar microscopic rearrangements of contacts as in monodisperse dimer packings accompany the formation of the density peak in binary mixtures largely independent of the mixture composition.


Introduction
Jammed particle packings have been studied to understand the structures of amorphous materials such as powders, reinforcing fibres, granular matter, and glasses [1]. Most studies mainly focused on monodisperse packings of spherical and non-spherical particles for which a plethora of experimental and theoretical results are available [2][3][4][5]. Elongated non-spherical particles such as ellipsoids [6][7][8][9], spherocylinders [10][11][12][13][14][15][16], and dimers [17][18][19][20] exhibit a non-monotonic variation of the packing density upon deviation from the spherical shape, with a maximum at specific aspect ratios. Although unlike spheres these shapes include rotational degrees of freedom, they do not fully represent realistic materials yet since polydispersity in both shape and size is inevitable for particle aggregates in nature.
In fact, numerous studies of binary and polydisperse packings of spherical particles demonstrate that changing the size distribution of the particles improves the packing density [21][22][23][24][25][26][27][28][29][30][31]. On the other hand, there are only few systematic investigations of mixtures elucidating the effect of non-spherical particle shapes. Studies of jammed packings of spherocylinder-sphere mixtures with the same diameter report a density maximum occurring at the same aspect ratio of the spherocylinders as in monodisperse packings, regardless of the relative volume fraction of the two species [16,32,33]. The density maximum is also present under different conditions, such as when spherocylinders and spheres have equal volume or different diameters [14,34,35]. For binary mixtures of two species of spherocylinders with the same diameter but different aspect ratios α 1 , α 2 , the density maximum that arises upon varying α 2 while keeping α 1 fixed, always occurs at the same aspect ratio of α 2 [14,35,36]. This aspect ratio likewise agrees with value at which the packing density of the monodisperse spherocylinder packing is maximal.
Non-spherical particle mixtures also satisfy a remarkable empirical ideal mixing law, which states that the inverse packing density is a linear superposition of the inverse packing densities of pure (single species) phases weighted by their relative volume fraction [16,37,38]. For example, the total packing densityf of spherocylinder-sphere mixtures is then given by where X 1 denotes the relative volume fraction of spherocylinders with a monodisperse packing density of f 1 and likewise f 2 is the packing density of monodisperse spheres. For binary mixtures the ideal mixing law establishes a linear dependence of f -1 on the relative volume fraction of one species. Clearly, since the relative volume fraction is fixed by the setup, the packing density of the mixture is then fully determined by the packing densities of the pure phases and thus independent of the segregation state. This independence implies further that particle orientations are completely uncorrelated showing an interesting similarity with a plastic crystal [16].
All these results provide avenues to optimize the packing densities of granular materials by mixing nonspherical particles, which might be relevant in industrial applications. Therefore, it is a fundamental question whether the universal density maxima and the ideality in mixing exist also for binary mixtures of elongated nonconvex particles. In this study, we focus on dimer shaped particles, a simple non-convex particle, obtained by overlapping two identical spheres. We generate disordered packings of both dimer-sphere and dimer-dimer mixtures by molecular dynamics (MD) simulations using a gravitational pouring protocol in LAMMPS [39,40]. In addition to the packing density, we also investigate the contact statistics of the mixtures. In our recent work on the structural analysis of monodisperse packings of dimers, we could associate the emergence of the density maximum with microstructural re-arrangements that are manifest in the contact statistics of local configurations [20]. Whether such re-arrangements are also related to the behaviour of the packing density in binary dimer mixtures is thus an interesting open question.
The remainder of the paper is organized as follows. In section 2 we present the details of our simulation method with LAMMPS. In section 3 and section 4, we show results on our analysis of the packing density, contact and coordination numbers for dimer-sphere and dimer-dimer mixtures, respectively. Finally, we conclude in section 5 with a discussion of our results.

Simulation method
We generate disordered packings of frictionless dimer-sphere and dimer-dimer mixtures in three dimensions with MD simulations. A single dimer is formed by overlapping two identical spheres with diameter d, whereby the dimer's aspect ratio (length divided by width) α is varied in the range 1.05 α 2, see figure 1. We use a packing-generation protocol where N particles are poured under gravity into a three-dimensional box of side length ≈20d. The box is constrained in the ẑ-direction by a rough surface at the bottom and an open top with periodic boundary conditions in the xŷ -plane. In our packing protocol, the dimers are initially placed at random positions and with random orientations within a specified insertion region 30 − 40d above the bottom and then released to settle in the box under gravity. Similar gravitational pouring protocols have previously been employed to generate jammed disordered packings of monodisperse spheres [41], ellipsoids [6,7], and dimers [17]. Since we model dimers as overlapping spheres, the interaction between two dimers can be determined from the pairwise interaction of spheres, which we assume to follow a spring-dashpot model as in the studies of frictional sphere packings of [41,42] that used MD simulations. In this model, two contacting spheres i and j having positions r i and r j , respectively 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 on sphere i is = where F ij n,t are the normal and tangential contact forces, respectively, given as [42]: where 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 Δs t is the elastic tangential displacement. K n,t and γ n,t are the elastic and viscoelastic constants, respectively [42]. In a gravitational field = -g g ẑ, the total force F i tot and torque t i tot on sphere i is then given as: where the sum runs over all j spheres in contact with sphere i. The equations of motion are integrated with the MD solver LAMMPS [39,40]. It allows in particular to define rigid bodies like the dimers by fixing the distance between its two constituent spheres' centres. In every time step, each dimer particle's total force and torque are computed as the sum of the forces and torques on its constituent spheres. Then, the coordinates, velocities, and orientations of the constituent spheres are updated so that the dimer moves and rotates as a single entity.
Simulations are run until a static equilibrium is achieved when the kinetic energy per particle is less than 10 −8 mgd. The number of particles, material parameter values, and time step used in the simulations of dimersphere and dimer-dimer mixtures are given in table 1. Note that most of these values are chosen following the discussion in [41]. We run ten independent simulations for both dimer-sphere and dimer-dimer mixtures and average all data points in our plots over them.

Packing density
We investigate dimer-sphere mixtures for various relative dimer volume fractions X d defined as: where N i and V i are the number of particles and the particle volume of component i, respectively, for component i ≡ d(dimers) and i ≡ s(spheres). An example of a disordered solid of a dimer-sphere mixture is shown in figure 2.
In order to calculate the packing density of the disordered solid, we define a bulk region, see figure 3. We observe that our packing protocol can result in some crystallization at the bottom of the box, depending on many factors such as the box width, the time step and the pouring height. In order to have results unaffected by dimer-sphere 12,000-18,000 2 × 10 5 2/7 50 0.001 dimer-dimer 12,000 2 × 10 5 2/7 50 0.001 this crystallization, we exclude the particles within 5 − 8d from the box floor from our packing density calculation. To determine the density, we need to calculate the total volume occupied by the particles. The Voro ++ package in LAMMPS uses a conventional Voronoi tessellation that provides the Voronoi volume W l of each sphere in the packing [40]. The Voronoi volume W k of a dimer is obtained by summing the Voronoi volumes of its two constituent spheres. Note that the Voronoi volumes of the particles within 5d from the upper-most particles in the packing can not be determined accurately due to deficiencies in their neighbourhood. Therefore, we also exclude those particles from the bulk region. The total bulk volume occupied by N s spheres and N d dimers in the bulk is calculated as = å Here, the volume of a dimer, V d , is calculated by subtracting the overlap volume from the sum of its constituent sphere volumes [20]. To consider a dimer as a bulk particle, the centres of both constituent spheres should be within the bulk. We calculate all average quantities discussed in the following for the bulk particles only.
We plot the total packing density of the mixture as a function of the dimer aspect ratio α in figure 4(a) for different relative dimer volume fractions. There are two limiting cases: we obtain pure sphere and pure dimer packings when X d = 0 and X d = 1, respectively. It can be seen from figure 4(a) that f exhibits a non-monotonic behaviour with the dimer aspect ratio by yielding a maximum at α ≈ 1.4, as in the case of monodisperse dimer packings [17][18][19][20], independently of the relative volume fraction. We observe that the packing density monotonically grows upon the increase in the relative amount of dimers up to the absolute maximum (f = 0.707) that has been achieved in our previous study of monodisperse dimer packings using the exact same protocol [20]. The appearance of a unique maximum is in agreement with previous studies of packings of spherocylinder-sphere mixtures with the same diameter condition for various compositions [16,32,33] and have been explained by the competition between local caging (a short spherocylinder can be oriented to minimize the space left by its contacting neighbours) and excluded volume effects. Other studies simulated the packings of spherocylinder-sphere mixtures only for one fixed spherocylinder volume fraction and with different diameters, they still found a maximum in the packing density at one unique rod aspect ratio [14,34,35].
In order to investigate the validity of the ideal mixing law equation (1) in our dimer-sphere mixtures, we plot the inverse packing density as a function of the dimer volume fraction X d for several aspect ratios. As shown in figure 4(b), the different curves are indeed well described by a linear relationship whose validity has also been established statistically, see the discussion in appendix. This linearity suggests that the packing density of the mixture is independent of the segregation state, i.e., a completely mixed packing has the same volume as one consisting of two separate phases each composed of only dimers or spheres, respectively, as discussed in detail in [16].

Contact and coordination numbers
We measure the contact number z, the average number of contact points of a particle and the coordination number n, the average number of neighbours of a particle. Here, we define neighbours as particles sharing at least one contact point. Note that multiple contacts between neighbouring particles are available for concave shapes like dimers, so z n, whereas z = n for smooth convex shapes like spheres, ellipsoids, and spherocylinders. In our definition, two dimers A and B are in contact if the distance between the centre of the constituent sphere i of dimer A and the centre of the constituent sphere j of dimer B satisfies r ij d. We further define four types of contact numbers for the mixtures of dimers and spheres: dimer-to-dimer z dd , dimer-tosphere z ds , sphere-to-sphere z ss , and sphere-to-dimer z sd , which denote the average number of contact points of the former component with the latter one. The contact number of dimer particles z d and the contact number of spheres z s are found as, respectively: The overall contact number z is then defined as: Similarly, there are four types of coordination number: dimer-to-dimer n dd , dimer-to-sphere n ds , sphere-tosphere n ss , and sphere-to-dimer n sd , which denote the average number of neighbours of the former component with the latter one. The coordination number of dimer particles n d and the coordination number of spheres n s are found as, respectively: The overall coordination number n is then defined as: In figure 5(a), we show the overall contact number z as a function of the relative dimer volume fraction X d . As can be seen, for all aspect ratios, z monotonically increases from the contact number of monodisperse sphere packings, 6.14 to that of pure dimer ones, 10.28 as X d grows. These two limits for the contact number of monodisperse packings are slightly above the isostatic values of z = 6 and z = 10, generally observed for disordered sphere [5] and dimer packings [19,43,44], respectively. These differences are because of the soft particle interaction model used in our gravitational packing protocol. The overall coordination number n of the mixtures is also strongly dependent on the relative dimer volume fraction, as displayed in figure 5(b). The variation of n with the aspect ratio becomes noticeable and approaches that for monodisperse dimer packings as X d increases. The contact numbers of spheres and dimers, z s and z d are calculated separately and shown in figure 6 as a function of X d . Both z s and z d vary approximately linearly with X d exhibiting a positive slope for α 1.8 and a negative slope for α = 2. When α = 2, z s < 6 for large dimer volume fractions, indicating that the spheres have fewer contacts than required for mechanical stability in this regime. Similarly, while dimers of α = 2 have more contacts than in the monodisperse case (z d > 10.28), shorter dimers lose some of their contacts in the dimersphere mixtures. Overall, we see that the isostatic condition is generally violated for the individual components in the disordered packings of dimer-sphere mixtures. While one species is hyperstatic, the other one is always hypostatic. We also calculate the coordination numbers of spheres and dimers (n s , n d ) separately and display them as a function of X d in figure 7. The dependence on X d is also approximately linear with similar trends as for the contact numbers, but the change in slope from positive to negative occurs now already for α > 1.4.
The validity of a linear relationship of z s , z d , n s , and n d as a function of X d highlights that also the componentwise contact and neighbour numbers satisfy a simple superposition principle for mixing and are fully specified by the corresponding numbers of monodisperse packings.

Contact configurations
In order to better understand changes in the microstructure of the packing due to shape variations and mixing, we investigate five distinct contact configurations as defined in our previous work on monodisperse dimer packings according to the number of contact points shared by two neighbouring particles [20], see table 2. We count the number of contact types occuring in the mixture per component pair and calculate the fractions of these contacts as shown in figure 8. Here, the Type 1 sd,ds fraction refers to the fraction of Type 1 contacts among all contacts between dimers and spheres, i.e., the fractions of Type 1 sd,ds and Type 2 sd,ds add up to 1. Likewise, the fractions of the five different types of dimer-dimer contacts add up to 1. Figure 8 shows that these fractions change initially upon increasing the dimer aspect ratio α up to the region at which the density peak occurs (α ≈ 1.4), but remain approximately unchanged for α > 1.4. We also observe that the fraction of Type 1     Table 2. Two distinct contact configurations of sphere and dimer (sd, ds) and five distinct contact configurations of two neighbouring dimers (dd). We show illustrations for aspect ratio α = 2. The total number of contact points for each type is: one (Type 1), two (Type 2, 3), three (Type 4), four (Type 5).

Pair
Type 1  Type 2  Type 3  Type 4  Type 5 sd, ds dd contacts increases, while that of Type 2 contacts decreases as the mixture packs more dense for increasing α, which is somewhat counterintuitive, since Type 2 configurations are locally more compact, see table 2. These observations are analogous to the monodisperse case, which further supports the qualitative picture that the peak in the packing density arises due to the interplay of structural rearrangements for small α and subsequent excluded volume effects with unchanged structure. The latter means that if the structure is unchanged while α increases, there is simply more empty space in the packing and the packing density has to decrease. Surprisingly, all the different fractions in figure 8 show almost no variations with a change of X d . Indeed the fractions for the dimer-dimer contact configurations are almost identical to the monodisperse case. We believe that this is not a simple consequence of the normalization of these fractions. After all, it could have been expected that the fraction of Type 1 dimer-dimer configurations is different when there are a lot of dimers available as contacts than if there are few, but this is not the case.

Dimer-dimer mixtures 4.1. Packing density
We measure the bulk packing density f of binary dimer mixtures for various mixture compositions. The two species of dimers (dimer 1 and dimer 2) have different aspect ratios, α 1 and α 2 , respectively. The relative volume fraction of dimer 1, X 1 is defined as: where N i and V i (i = 1, 2) are the number of particles and the volume of the i-th component, respectively. We follow the same steps as in section 3 for the packing density calculation. The mixture packing density variation with α 1 and α 2 is shown in a heat map for three different X 2 values in figure 9(a)-(c). For all fixed aspect ratios of one component (e.g. α 1 ), f exhibits a non-monotonic relationship with the aspect ratio of the second species α 2 : it increases up to a peak and subsequently decays, whereby the peak always occurs at α 2 ≈ 1.4, i.e., at the aspect ratio at which the monodisperse dimer packing has its density maximum [20]. This behaviour is irrespective of the relative volume fractions of the two components. As in the case of dimer-sphere mixtures, the packing density of binary dimer mixtures never exceeds the maximum packing density observed for monodisperse dimer packings at f = 0.707 [20]. These results agree with the findings of disordered packings of binary spherocylinders [14,35,36]. We also plot the inverse packing density of the mixtures of dimer 1 with α 1 = 1.4 and dimer 2 with various aspect ratios α 2 as a function of X 2 and compare them with the results from the ideal mixing law equation (1) in figure 9(d). As can be seen, the curves exhibit a linear relationship with the volume fraction when α 2 < 2. However, we observe that the inverse density curve is slightly concave-upward for α 2 = 2, exhibiting deviations from the ideal mixing law. A statistical analysis confirming the validity of the ideal mixing law for α 2 < 2 can be found in appendix. For α 2 = 2 a perfect fit of the data can be obtained using a 4th order polynomial instead of a linear curve, see figure 14(b).

Contact and coordination numbers
For the contact and coordination number analysis, we focus on the mixtures of dimer 1 with α 1 = 1.4 and dimer 2 with various aspect ratios, α 2 . We use the same definitions for partial contact and coordination numbers of the two components as in section 3. We also measure the overall contact and coordination numbers of the mixtures. For all calculations, we consider only bulk particles.
We show the overall contact number z and the overall coordination number n as a function of the relative volume fraction of the second species in figure 10. We observe that the coordination number n is sensitive to both changes in α 2 and X 2 , while z is essentially constant with z = 10.3, the same value of monodisperse dimer packings, see figure 10(a) [20]. Therefore, disordered packings of binary dimer mixtures also satisfy the isostatic condition regardless of the relative volume fraction of the components. In this context, we should mention that previous studies of bidisperse and polydisperse sphere packings also found a constant mean contact number, irrespectively of the particle size distribution and the relative amount of different components [45][46][47][48][49].
Although the size effect is the only parameter in the sphere case whereas we take into account both shape and size effects in our packings when α 2 is varied, a constant mean contact number seems a generic result of mixing two components having the same isostatic value for their disordered monodisperse packings.
We calculate the contact and coordination numbers of the two components (z 1 , z 2 , n 1 and n 2 ) separately and display them as a function of X 2 in figure 11 and figure 12, respectively. All these numbers satisfy an approximate linear relationship as in the case of dimer-sphere mixtures. Figure 11 shows that the contact number of the longer dimer species is always hyperstatic with a higher coordination number than in the monodisperse case (X 2 = 0), while the shorter dimer species is always hypostatic. Both contact and coordination numbers exhibit a crossover from positive to negative slopes for α > 1.2. Comparing with the dimer-sphere case ( figure 6(b)), we see that dimers mixed with spheres need one to two fewer contacts at each aspect ratio than dimers mixed with another dimer species. On the other hand, the coordination numbers of the dimers are very similar (compare figure 7(b) and figure 12(b)) indicating that the extra contacts arise from multiple contacts between  neighbouring dimer pairs. The coordination number n 2 exhibits an intersection of the α 2 > 1.4 curves for X 2 = 0.9, whose origin is not clear.

Contact configurations
We also determine the fractions of the five contact configuration types of table 2 for each component (dimer 1 and 2) and display them as a function of the aspect ratio of the second component α 2 for three different volume fractions in figure 13. Comparing with the dimer-sphere case, the fractions of the dimer 2 contact configurations are very close to those of the dimers in dimer-sphere mixture and thus also to those of monodisperse dimer packings, but a slight dependence on X 2 can be observed, in particular for small α 2 . The Type 1 (Type 2) fractions of the dimer 1 contact configurations are considerably smaller (larger) than those of   spheres, but always larger (smaller) than those of dimer 2. The curves become slightly flatter when X 2 decreases, i.e., there are in particular more Type 1 and less Type 2 configurations for small α 2 , which is expected for the dimer 1 fractions, since they need to approach constants when X 2 → 0, but why those of dimer 2 change likewise is unclear. Overall, as in the dimer-sphere and monodisperse dimer case, we see that the fractions show the largest variation in the regime α 2 < 1.4 and remain approximately unchanged for α 2 1.4.

Conclusions
Our investigation shows that the packing densities of dimer-sphere and dimer-dimer mixtures exhibits a nonmonotonic variation in the packing density as the aspect ratio of one species is varied, confirming previous results of spherocylinder-sphere and ellipsoid-sphere mixtures. We also confirm the validity of the ideal mixing law equation (1) for both types of mixtures highlighting the independence of the packing density on the segregation state. Somewhat surprising is the observation that the packing density of dimer-sphere and dimerdimer mixtures is always below the maximum packing density of monodisperse dimers (with α = 1.4), while bidisperse spheres of different diameters, e.g., increase the packing density compared with monodisperse spheres. However, this behaviour trivially follows from the ideal mixing law, since the packing densities of the pure sphere or dimer phases are below the maximum and thus the total packing density as well, see equation (1). Another manifestation of this ideal mixing property is evident in the linear behaviour of the component-wise contact and coordination numbers, which is not observed for the corresponding total contact and coordination numbers. Our analysis of the contact configurations confirms the qualitative picture that the peak in the packing density arises due to the competition of locally optimal rearrangements and excluded volume effects, which is here manifest in the significant variation of the configuration statistics as the dimer is elongated until the maximal packing density is achieved at α = 1.4.
In terms of future work, it would be interesting to understand further what kind of observables in the mixture exhibit similar ideal mixing properties. One caveat when studying mixtures of non-spherical particles is the large parameter space, which complicates any systematic analysis. In our case, the fact that we keep the diameters of the spheres that constitute the dimers constant leads to both relative shape and size variation as the aspect ratio of one dimer component is varied. Disentangling the effect of shape and size variation as attempted e.g. in [36] could shed valuable further insight. The inverse packing density f -1 of the dimer-sphere mixtures as a function of the relative dimer volume fraction X d for several dimer aspect ratios (same data as shown in figure 4(b)). Solid lines are obtained from equation (1), and dashed lines are fitted to the data points with linear regression. Each shaded region shows the 99% confidence interval for the regression slope. (b) The inverse packing density f -1 of mixtures of dimer 1 with α 1 = 1.4 and dimer 2 with various aspect ratios α 2 as a function of X 2 (same data as shown in figure 9(d)). Solid lines are obtained from equation (1), and dashed lines are fitted to the data points with linear regression for α 2 < 2 (including the 99% confidence interval) and with a 4th order polynomial for α 2 = 2.