Andraž
Gnidovec
*a,
Anže
Božič
b and
Simon
Čopar
a
aUniversity of Ljubljana, Faculty of Mathematics and Physics, SI-1000 Ljubljana, Slovenia. E-mail: andraz.gnidovec@fmf.uni-lj.si
bDepartment of Theoretical Physics, Jožef Stefan Institute, SI-1000 Ljubljana, Slovenia
First published on 29th September 2022
Packing problems are abundant in nature and have been researched thoroughly both experimentally and in numerical models. In particular, packings of anisotropic, elliptical particles often emerge in models of liquid crystals, colloids, and granular and jammed matter. While most theoretical studies on anisotropic particles have thus far dealt with packings in Euclidean geometry, there are many experimental systems where anisotropically-shaped particles are confined to a curved surface, such as Pickering emulsions stabilized by ellipsoidal particles or protein adsorbates on lipid vesicles. Here, we study random close packing configurations in a two-dimensional model of spherical geodesic ellipses. We focus on the interplay between finite-size effects and curvature that is most prominent at smaller system sizes. We demonstrate that on a spherical surface, monodisperse ellipse packings are inherently disordered, with a non-monotonic dependence of both their packing fraction and the mean contact number on the ellipse aspect ratio, as has also been observed in packings of ellipsoids in both 2D and 3D flat space. We also point out some fundamental differences with previous Euclidean studies and discuss the effects of curvature on our results. Importantly, we show that the underlying spherical surface introduces frustration and results in disordered packing configurations even in systems of monodispersed particles, in contrast to the 2D Euclidean case of ellipse packing. This demonstrates that closed curved surfaces can be effective at introducing disorder in a system and could facilitate the study of monodispersed random packings.
In particular, packing of anisotropic hard objects in Euclidean geometries has been extensively studied during the last two decades, with development of theoretical models9–11 and research on different particle shapes,12,13 effects of confinement,14–17 and new methods to describe structural properties of jammed states18,19 continuing to this day. In comparison to packings of spherical particles, anisotropy introduces additional degrees of freedom into the system and in turn allows for denser packing configurations for all aspherical shapes in three dimensions,10,20 with asphericity in general influencing the packing density in a non-monotonic way.21 One of the prototypical cases for studying the effects of anisotropy on particle packing are systems of ellipsoids, the most simple aspherical shapes that are routinely used to model particles in a variety of soft matter, granular, and molecular systems. They are known to improve the density of jammed disordered states22–25 and can form unusually dense crystal packings,26,27 with numerous investigations highlighting the mechanical properties of jammed configurations of hard ellipsoids and ellipses.19,24,28–32 Additionally, ellipses were found to maximize the packing fraction of random sequential adsorption in two dimensions among the set of shapes of smoothed n-mers and spherocylinders.33
Besides particle shape, confinement also significantly affects packing properties—most notably its density—both for spherical particles34,35 and anisotropic particles such as rods and ellipses.16,17,36,37 In two dimensions, one can consider systems on a sphere as a peculiar case of confinement that arises due to the compactness of the spherical surface, replete with obligatory topological defects.38,39 The question of optimal packing of circles on the sphere was first posed by Tammes already in 193040 and has been well-researched since,41–44 with studies of jamming having been extended to different surfaces with positive curvature.45 However, investigations into the packing of anisotropic objects on spherical surfaces have to date been mostly limited to systems of hard rods.46,47 There, it was observed that in contrast to their optimal flat-plane packing, short rods adapt to the spherical geometry more efficiently than both spheres and longer rods, demonstrating the significant interplay between curvature and particle shape.
Only very recently have some works started to explore also the packings of ellipsoids on curved surfaces such as sphere.48 Packing of ellipses on the sphere is fundamentally different compared to their 2D packings in Euclidean space. Monodispersed systems of ellipses always crystallize in the plane when simulated by conventional packing algorithms and their packing fraction ϕC will be the same as for the hexagonal packing of circles (ϕC ≈ 0.9069), with ellipse and circle packings directly connected by an affine transformation.49 Analysis of disordered ellipse packings in flat plane therefore requires introducing bidispersity into the system in order to suppress crystallization. However, the closed surface of the sphere and its topological requirement for lattice defects preclude any periodic packing solutions, so in this sense, every configuration of ellipses on the sphere can be considered globally disordered.
In this work, we study a natural extension of the Tammes problem from circles to purely 2D spherical ellipses of arbitrary aspect ratios. The ellipses are defined as having a constant sum of geodesic distances to two foci, which also corresponds to having elliptical orthogonal projections onto a plane.50 We first introduce the model and explain the packing simulation procedure used to generate jammed configurations of spherical ellipses. We then present packing results, demonstrating that the packing fraction depends non-monotonously on the ellipse aspect ratio and can change with system size. We further analyze the generated packings by calculating the structure factor and tensor harmonic expansion to better understand the ordering of spherical ellipses. Finally, we compare our results to previous investigations of similar systems.
A = 2πR2(1 − Λ0(ψ,k)), | (1) |
![]() | (2) |
It is important to consider that the shape of spherical ellipses depends on their size α and therefore on the number of particles in a random close packing (RCP). The definition of a spherical ellipse can be satisfied by taking the intersection of a sphere and an elliptical cylinder with semiaxes x1 = Rsin
α and x2 = R
sin
α/ε. The details of this description are given in ref. 50. As α → π/2, the semi-major axis of this ellipsoid approaches R, resulting in a sharper and sharper vertex of the spherical ellipse (see Fig. 1(b)). At α = π/2, the spherical ellipse becomes the shape of a lemon wedge. This shape-changing effect is intrinsic to the purely 2D spherical ellipse model and is more pronounced at larger α (lower N in RCPs) and higher ε. As an alternative to keeping ε of the spherical ellipse constant for different particle sizes α, one could conserve the aspect ratio x1/x2 of the elliptical cylinder that defines the ellipse; however, this does not address the shape-changing effect. Additionally, in this case small changes in cylinder semiaxes close to x1 = R result in large changes of α, which is unfavorable for the iterative packing simulation described below. However, in physically relevant systems, N will be large enough for these shape effects to be negligible.
In contrast to the Euclidean space where the bulk of packing research is focused on understanding system properties in the thermodynamic limit N → ∞, systems on a sphere are intrinsically of a finite size. The characteristic length scale, given by the semi-major axis a, and the sphere radius R can be varied independently, thus changing the effect of curvature on system behavior. We therefore focus on studying the interplay between finite-size effects and curvature that is most prominent at smaller system sizes (N ≤ 1000). In the thermodynamic limit, one should recover the 2D Euclidean results, as discussed in Section 4.
As packings in nature are not necessarily ordered in the most dense configuration, we do not focus solely on finding the optimal packing solution but rather on analyzing the general properties of jammed configurations of ellipses. That is, we adopt an ensemble-based approach to analyzing packings where the packing simulation is run repeatedly from random initial conditions to gather a statistical sample of jammed configurations with different packing densities, weighted by the probability of occurrence. The generated packings can be related to the concept of RCP, where “randomness” stems from initialization and is not necessarily quantified for end configurations.52–54 This is fundamentally different from the approach based on individual configurations, used to study the mathematically rigorously defined maximally random jammed (MRJ) packings, determined as packings that minimize a chosen order metric (among other requirements).55 As we discuss in Section 4, finding an appropriate order metric to quantify randomness in packings of spherical ellipses proves to be elusive because it is not clear what the ordered phase should be in the first place.
We generate the packings by a compression–decompression algorithm similar to the method used previously to study RCP of spheres and ellipses.29,34,52,56 The procedure requires a measure of overlap to be defined as the system energy; however, conventional overlap parameters used for measuring distances between ellipsoids in Euclidean geometries57–59 do not work for purely 2D ellipses on a sphere. We therefore use a recently developed algorithm by Gnidovec et al.50 to determine the overlap function λ between two spherical ellipses. This algorithm is based on calculating the eigenvalues of a linear interpolation between quadratic forms that define spherical ellipses, and returns λ < 1 if two ellipses overlap. The energy for a pair of ellipses is then defined as E(λ) = 1/λ + λ − 2 for overlapping ellipses and E(λ) = 0 otherwise. Note that the first nonzero term in the expansion of energy around λ = 1 on the overlapping side is the harmonic term (1 − λ)2.
The packing simulation is initialized with random positions and orientations of ellipses on the sphere at a packing fraction ϕ ≈ 0.5. The size of the particles α is then increased step-wise and the configuration is relaxed using energy minimization after each iteration to eliminate overlaps. As we approach a jammed configuration, step size is decreased and, additionally, a decrease of ellipse sizes is permitted if overlaps cannot be removed by translations and rotations of particles. We stop the simulation once the energy falls in the interval E0 < E < 2E0, with E0 = 10−8. It follows that the packing fraction of generated configurations is precise within 10−4, which is sufficient for our analysis. Note that the initial particle growth rate can be set arbitrarily in the algorithm and can impact the density of generated packings. We choose a value of Δϕ ∼ 10−3, which is small enough to prevent significant “unphysical” rearrangements that follow large overlaps of particles while also providing satisfactory time performance of the simulation.
Our simulation procedure technically considers spherical ellipses as soft particles, and some caution is required before classifying the resulting states as equivalent to jammed configurations of hard particles. It can be shown, however, that configurations with locally minimal energy always correspond to jammed states for repulsive interparticle potentials.60,61 What is more, in the limit of zero external forces and pressures on the system (as is the case for our method), the resulting configurations can be closely related to collectively jammed packings of hard ellipses.24,52,62 Here, we adopt the classification of jamming from ref. 63. Note that the more demanding jamming categories, namely, the strictly jammed and metric jammed definitions (the latter proposed by Burke45), cannot be applied to spherical systems, as the underlying surface cannot be deformed. We also expect that alternative algorithms for generating packings of hard particles, such as the Lubachevsky–Stillinger algorithm64 or the Zinchenko algorithm,65 would produce qualitatively similar packings.62
![]() | ||
Fig. 2 (a) Average packing fraction 〈ϕ〉 and (b) number of contacts 〈Z〉 as functions of the ellipse aspect ratio ε. Dotted curves do not represent a model fit and serve only as guides to the eye, with the shadowed areas around the curves indicating the scale of errorbars. The detailed distribution of packing densities around the plotted averages is shown in Fig. 4. Inset in panel (b) shows the dependence of 〈Z〉 on contact tolerance δλ. A small number of particles, called rattlers, can be underconstrained and able to move within the cage of the surrounding particles. These are also included in the calculation of 〈Z〉. |
Jammed configurations of particles must also satisfy force and torque balance that is already implied in our packing generation algorithm. As a requirement to achieve mechanical stability, the isocounting conjecture was proposed, which states that in static packings of frictionless particles, the number of contacts between particles needs to be equal to the number of degrees of freedom (DOF) in the system. It follows that the mean number of contacts per particle is Ziso = 2df − 2/N, where df is the number of DOF per particle. The second term is a consequence of confinement and finite size (configurations are invariant to two global rotations), and becomes small for large systems. On a spherical surface and in 2D in general this results in Ziso = 4 for circles and Ziso = 6 for ellipses. While random packings of spheres in a flat space are indeed isostatic, packings of anisotropic particles such as ellipses and ellipsoids are stable also when hypoconstrained with Z < Ziso, as has been demonstrated both in simulations and experiments.22,24 The missing contacts were shown to belong to quartic modes in the spectrum of the dynamical (Hessian) matrix, describing collective rotational motions of particles.24,29
In Fig. 2(b), we show the dependence of the mean number of contacts per particle 〈Z〉 on the aspect ratio ε, averaged again over all generated configurations at each point in the parameter space. In the simulated packings, contacts are determined between pairs of ellipses if |λ − 1| < δλ, as λ represents the area scaling factor required to reach exact tangency. The value of tolerance is chosen to be δλ = 10−4, where the curves of 〈Z〉 as a function of (δλ) are almost flat over the entire range of simulated ε, as shown in the inset of Fig. 2(b) and similar to what has been observed for ellipsoid packings.68 Our results show that for most circle packings on a sphere we have 〈Z〉 < 4, but there are exceptions with higher contact numbers, such as N = 30, where the average configuration is hyperconstrained (〈Z〉 > 4). As particle shape starts to deviate from a perfect circle, the mean contact number steeply increases as the introduced additional degrees of freedom in the system need to be balanced to reach stability of the configuration. It reaches a maximum at ε ≈ 2, which notably differs from the value where the peak in the packing fraction is attained (ε ≈ 1.5). In 3D packings of ellipsoids, these two extrema are closer together.22 Another difference with 3D results is the subsequent decrease in the mean contact number, as it stays mostly constant for ellipsoid packings after reaching the maximum value. This decrease at high ε could be connected to local crystallization in packings on the sphere (ordered local domains) where many particles only have four contacts. We discuss this more thoroughly in the ESI† (see Fig. S1). Finally, the system with N = 30 particles exhibits consistently lower mean number of contacts than other presented cases at all values of ε > 1. For such small system sizes, the effect of confinement to a compact surface that lowers the number of degrees of freedom in the system (2/N term) becomes significant. Additionally, at high ε, the interplay between ellipse shape and curvature that drives the increase in packing fraction also promotes parallel ordering of neighboring ellipses and consequently domain formation, which amplifies the effect observed at large system sizes.
We further study the dependence of packing fraction on system size in Fig. 4, where we use boxplots to show how the results of distinct simulation runs are distributed. Orange lines represent the average values at the respective points. For circular particles at ε = 1, 〈ϕ〉 is non-monotonic, with large variations especially at small system sizes where curvature and compactness of the spherical surface can lead to highly frustrated local arrangements of circles. These variations are not surprising—the best packing configurations closely follow the Tammes solutions (densest circle packings, marked by red crosses in Fig. 4), which exhibit the same property.42 Note that all generated packings at N = 20 and most of those at N = 30 return the respective Tammes configurations. A slight increasing trend in 〈ϕ〉 can be observed as we move towards large system sizes. The packing fraction for circles on the sphere should converge to the planar case ϕC ≈ 0.9069; however, it has been demonstrated that this convergence is very slow.69 In general, the distribution of packing fraction within the ensemble is centered below the maxima dictated by the Tammes solutions and becomes approximately Gaussian for N ≳ 80, which is found to be true also for distributions at other values of ε (with a few exceptions, all distributions pass the normality test). As the aspect ratio ε is increased above 1, any notable differences in packing densities between configurations with different N begin to disappear. This is a direct consequence of adding additional degrees of freedom to the system, resulting in ellipse packings that are less frustrated than the packings of circles. Further details on the transition from circles to ellipses are available in ESI.† For aspect ratios around the peak of the packing fraction (ε ≈ 1.5), the packing efficiency remains approximately constant with respect to N, as indicated already in Fig. 2(a). Increasing the aspect ratio further, the effects of particle curvature at small system sizes that increase the packing fraction (as already mentioned in the discussion of Fig. 2) become evident. Whereas the decrease of 〈ϕ〉 with N is relatively gradual at ε = 5, it is significantly more pronounced at ε = 10 where the case for N = 20 satisfies the condition perfect for packing (the average generated packing has a packing fraction of ϕ = 0.994 (Fig. 3(a))). Average packing fraction then steeply decreases with N until the slope flattens at around N = 50 and continues to decrease more gradually.
![]() | (3) |
![]() | ||
Fig. 5 Heatmap of the decimal logarithm of the spherical structure factor [eqn (3)] for N = 100 at different aspect ratios ε and wave numbers ![]() ![]() ![]() ![]() |
In Euclidean space, jamming is conjectured to be related to hyperuniformity, i.e., the suppression of density fluctuations at large distances described by a vanishing Euclidean structure factor S(k) as |k| → 0 or—for some systems of aspherical particles—vanishing spectral density.70–72 The notion of hyperuniformity was recently expanded to spherical surfaces, where assumes the role of k, and by connecting the spherical structure factor to the spherical cap density, one can show that the hyperuniformity criterion on the sphere corresponds to the limit
.39Fig. 5 demonstrates that this condition is satisfied up to a small residual on the order of 10−3, at least for the range of aspect ratios shown in the figure. Some larger deviations were observed at high ε, particularly at ε = 10 as seen also in the inset of Fig. 5, where this residual is between one and two orders of magnitude higher (more significantly so at small N). We also calculated the asimptotic expansion coefficients AN for the cap number variance as another way to establish hyperuniformity on the sphere.39,73 The coefficients AN are on the order of 10−3 or smaller for all ensemble samples in the inset of Fig. 5, including ε = 10, which demonstrates that number density fluctuations vanish for larger cap sizes and the systems can be considered hyperuniform. The different extents to which the limit
is satisfied for different N and ε is similar to deviations from perfect hyperuniformity that were observed also in Euclidean jammed states and are believed to be related to both imperfections in packings (rattler particles)74 as well as finite size effects, which cannot be avoided on the closed surface of the sphere.
Besides the harmonic expansion of the ellipse positions, further insight into orientational ordering of ellipses in close packings can be obtained by expanding the tensor field , where unit vectors
k describe the ellipse orientations, over two orthogonal tangential tensor spherical harmonics on the sphere. Details of this expansion are given in ESI.† In contrast to the scalar structure factor obtained from the positional order, the two tensor structure factors, defined in an analogous fashion from tensor expansion coefficients, do not vanish at small
for any value of ε (Fig. S4 in ESI†). Instead, they have a peak at
= 2 (the smallest
where the expansion is defined), which shows that even though the ellipses are not positionally ordered, there exists some trend in the average orientational order over long length scales. The peak is higher at large ε, showing that its magnitude indicates the preference for parallel ellipses during the densification process of the packing simulation, as parallel local structures enable higher packing densities. See ESI† for further tensor expansion results and clarifications.
On the other hand, the dependence of packing density on the aspect ratio ε that we obtain for monodispersed spherical ellipses strongly resembles the curve for the bidispersed planar case, not only qualitatively but also quantitatively. The diamond-marked yellow curve in Fig. 6 shows the results obtained by Delaney et al.,66 and other studies of 2D ellipse packings have also observed similar densities.24,31 Just like in the packings of spherical ellipses, the best packing is achieved for ε ≈ 1.5, with the peak being slightly higher in the Euclidean plane. What is more, the decrease in the packing efficiency when ε increases is comparable between the two cases. The strong quantitative resemblance could indicate that there exists a fundamental connection between random packings of ellipses on flat and curved surfaces, even though the mechanism that prevents crystallization is different. As expected, both results also show higher packing densities compared to the elusive monodispersed MRJ packing of hard disks with ϕmono = 0.826, generated by an algorithm designed specifically to avoid crystallization.76
Finally, exploiting curvature as a means of introducing disorder into a system raises the question of the behavior of spherical packings at large system sizes. In the limit of N → ∞ and, consequently, vanishing curvature, one would expect the system to form large domains with crystalline order and with defects concentrated in the grain boundaries, effectively transitioning to the (ordered) planar case where the packing fraction is independent of the ellipse aspect ratio. This is indeed the case for spherical packings of circles, where this transition is continuous as the density of grain boundaries very slowly decreases with increasing system size.69 Similar behavior would be expected for systems of spherical ellipses; however, given the limited range of N studied in our simulations, we are unable to exclude the possibility that the transition from random to crystalline order appears faster around some specific system size. This problem is left for future investigations as it would require a significant computational effort.
In real systems, building blocks of anisotropic shapes are often preferred over spherical ones: for instance, diatom frustules of various adnate diatom species are elliptical with mean aspect ratio ε ≈ 1.7,77 elongated fat particles have been shown to improve the stability (crystallization) of ice cream,78 and using elongated carrier particles in dry powder inhalers can considerably increase the amounts of drug delivered to the lower airway region,79 to name just a few. On spherical surfaces, elliptical particles such as polystyrene ellipsoids can also be used to stabilize Pickering emulsions more efficiently than spherical particles.80,81 They form a dense packing at the spherical liquid–liquid interface, providing a possible realization of packings studied in this work. Similarly, elliptical particles could be used to manipulate properties of liquid marbles.82 In biological systems, our work could help with understanding the packing of adsorbed elliptically-shaped BAR proteins on lipid vesicles.83 These proteins are intrinsically curved and resemble the shape of geodesic spherical ellipses which makes our model preferable compared to the 3D ellipsoid alternatives. However, in all these examples, the interactions between particles and building blocks are more complex than the hard-core repulsion studied in this work. Our model is furthermore purely two-dimensional, which is more suited for thin systems that adhere to the surface, such as structures within biological membranes, but can lead to some discrepancies when describing real systems with three-dimensional ellipsoidal building blocks. We also do not take into account the possible deformations of the underlying spherical surface, or any reshaping of the particles themselves. Each system of interest must therefore be evaluated with regard to whether the geodesic ellipse is a sufficient approximation to the actual shape for the purposes of the investigation.
Some additional challenges still remain. As an ordered reference configuration for the arrangement of ellipses on the sphere does not exist, we are unable to quantify disorder in the system. This essentially precludes us from a detailed analysis of individual configurations, also in terms of finding the MRJ packing density. Additionally, the system sizes studied in this work are relatively limited due to the computational costs involved. In future investigations, more efficient packing algorithms could enable the study of large domain formation that is expected to occur with increasing system size. The model could also be expanded either to non-hard interactions between particles where penetration is possible or to include long-range interactions such as multipolar electrostatic interactions. These more complex models could be used to better describe the behavior of various colloidal or biological systems where building blocks are preferably oblong, or even to enable the design of packings with a desired local structure.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2sm00624c |
This journal is © The Royal Society of Chemistry 2022 |