Using particle shape to control defects in colloidal crystals on spherical interfaces

Abstract

Spherical particles confined to a sphere surface cannot pack densely into a hexagonal lattice without defects. In this study, we use hard particle Monte Carlo simulations to determine the effects of continuously deformable shape anisotropy and underlying crystal lattice preference on inevitable defect structures and their distribution within colloidal assemblies of hard rounded polyhedra confined to a closed sphere surface. We demonstrate that cube particles form a simple square assembly, overcoming lattice/topology incompatibility, and maximize entropy by distributing eight three-fold defects evenly on the sphere. By varying particle shape smoothly from cubes to spheres we reveal how the distribution of defects changes from square antiprismatic to icosahedral symmetry. Congruent studies of rounded tetrahedra reveal additional varieties of characteristic defect patterns within three, four, and six-fold symmetric lattices. This work has promising implications for programmable defect generation to facilitate different vesicle buckling modes using colloidal particle emulsions.

---

Introduction

Colloidal emulsions are found in everyday products from cosmetics to paint.^1,2 Pickering-Ramsden emulsions are a subset of this class of materials that are stabilized by the interfacial adsorption of colloidal particles rather than a molecular surfactant, such as block-copolymer stabilized polymersomes and micellar droplets.^3–6 These colloidal emulsions have novel applications including encapsulating payloads for controlled drug delivery,^7,8 facilitating interfacial bio-catalysis,⁹ and controlling swelling dynamics to drive particle assembly.¹⁰ They also serve as an exciting option for non-toxic and renewable emulsions with high stability, as seen in studies of cellulose particle-stabilized emulsions and foams that are shelf stable on the order of months.^11,12 Their novel qualities arise from highly tunable properties (e.g. mechanical stability, droplet size, and porosity) that are dependent on parameters such as solvent choice,¹³ particle characteristics (size,¹³ material, shape,^14–16 charge,¹⁷ and surface coverage¹¹). Many of these properties are relevant for meeting cargo release,^18,19 retention,²⁰ and efficiency criteria²¹ in drug delivery applications. Examples of precise pore size control, as seen in polystyrene particle stabilized droplets, are achieved by sintering, adsorption of polymers, or particle aggregation, leading to a broad range of tunable elastic moduli and breaking forces.^22,23

Particle ordering in Pickering-Ramsden emulsion droplets is mirrored in many natural systems such as spherical monolayers of epithelial cells,²⁴ virus capsids,^25–29 and pollen grains.³⁰ Models of these systems often use spherical particles constrained to an emulsion interface. The shape and pair-wise interactions of these spherical particles drive lattice structure, local symmetry, and defect morphology on the two-dimensional surface, mirroring the impact of shape and interactions in bulk crystal structure in three-dimensions.^31,32 On flat interfaces, spheres pack most efficiently with hexagonal ordering, where each particle i has the expected coordination number of neighbors c_i = c_hex = 6. Defects in hexagonal packings take on individual topological charge q_i = c_hex − c_i, and form as isolated point disclinations or as = ±1 dislocation pairs. On non-flat spaces, such as the sphere and topologically equivalent convex polyhedra, the Euler characteristic³³ is given as χ = V − E + F where V, E, and F are the number of vertices, edges, and faces, respectively, of the confining surface. Given the Euler characteristic χ = 2 of a spherical droplet, topology requires the preferred hexagonal packing to integrate topological defects with a total surface charge of . The required total topological charge manifests through the formation of dislocations, disclinations, or long defect “scars” of connected disclinations.^29,34 This is readily seen in the familiar patterning of a soccer ball, which has twelve five-fold + 1 charge defects distributed with icosahedral symmetry over the surface, or in the assembly of protein subunits leading to the distinctive icosahedral buckling behavior of many virus capsids.^27,35–37 Lattice symmetry and defect structure are fundamentally related on spherical interfaces, as illustrated in molecular models of Hertzian spheres that assemble into hexagonal and simple-square lattices with competing defect motifs³⁸ and simulations of non-overlapping square tetratic molecules.³⁹

Introducing shape anisotropy has proven to be a strong design parameter to affect assembly and packing behavior.^40,41 Shape anisotropy affects ordering phenomena in flat space, as seen in an experimental study of superballs that showed a continuous deformation from a simple-square (cubic) to hexagonal (fcc) lattice^42,43 and in numerous simulation studies of hard polyhedra,faceted spheres, and a myriad of other shapes that used shape-induced valence to drive the system towards various crystal structures.^44–55 Studies of Pickering-Ramsden emulsions stabilized by rods,⁵⁶ cubes,¹⁴ peanuts,⁵⁷ and Janus spheres⁵⁸ have all reported novel ordering and demonstrated the feasibility of combining shape anisotropy and interfacial curvature. Further simulation work has shown that rods of varying asphericity adapt to spherical interfacial curvature readily, packing with nematic order with defects distributed on the vertices of a tetrahedron or about the arc of a great circle when they are orientationally locked^59–61 and show complex two-step melting behavior when free to rotate over a range of number densities.⁶²

In this paper, we significantly expand on these previous studies by investigating two shape families of rounded polyhedra – cubic and tetrahedral – constrained to a spherical surface. We focus our simulations on hard particles that form ordered, entropy-maximizing structures.⁴¹ We demonstrate how topology and local order conspire to dictate defect morphology and distribution about the interface. We quantify order for a variety of lattice types over a range of particle densities, sphere interface radii, and particle roundedness. Furthermore, we show the different ways in which these systems resolve the geometric frustration imposed by the incompatible topology.

Methods

We ran hard particle Monte Carlo (HPMC) simulations of n ∈ {1000, 1500, 2000, 2500} particles, with their centroid bound to the surface of a sphere with radius R_S. The constraining sphere surface radius R_S is coupled with the particle number-area density ρ_n = n/(4πR_s²), as described pictorially in Fig. 1a. Particle shape varies continuously between a sphere and a polyhedron (cube or tetrahedron) via a rounding parameter s = R_R/(R_R + R_PI) ∈ [0,1], given the polyhedron insphere radius R_PI and rounding radius R_R; s = 0 indicates no rounding of the underlying shape (sharp) and s = 1 indicates an isotropic sphere,⁶³ as seen in Fig. 1b. While superball schemas are commonly used to describe and interpolate shape between the cube and sphere, the rounding radius method used here gives similar self-assembly behavior and is computationally less expensive.⁴⁸ We clarify values of s between polyhedral families by subscripts, s_cube and s_tetrahedron. All particle volumes are normalized to the unit sphere (a sphere with radius σ = 1), .

Fig. 1 Shape definition and constraining radius: (a) Particle centroids are constrained to move on a spherical interface with radius R_S. (b) Rounded polyhedra are defined by the union of an underlying polyhedron (cube or tetrahedron) and a rounding sphere. The rounding parameter s uses the polyhedron insphere radius R_PI and a sphere with rounding radius R_R to define particle shape, seen here for a head-on view of a rounded cube. (c) Shape descriptive vectors are defined for both the cube and tetrahedron. These include face normals (f⃑), vertices (v⃑), edge normals , and edge tangents . These descriptive vectors are used in the calculation of various orientation dependent order parameters for systems seen within this work.

To obtain particle assemblies on a sphere of radius R_S, particles are first placed on a sphere with radius 5R_S. The sphere is slowly reduced, over ∼10⁶ simulation MC sweeps, to a target radius R_S and corresponding number density ρ_n ∈ [0.215, 0.287]. We simulated systems with particles in the cube family for 10⁷ MC trial moves, while tetrahedral particle systems were run for up to 4 × 10⁷ trial moves, with an acceptance ratio of 33%. We found that these simulation lengths were sufficient to equilibrate the particle positions and orientations, by observing autocorrelation functions decaying to zero for the relevant system order parameters for the total number of “ordered” particles over the simulation length. The breadth of system parameters studied was chosen to demonstrate the robustness of our analysis methods and the general effect on the resulting assemblies of imposed topological constraints. Each combination of these parameters was simulated with 10 replicates each for n = 1000 and 5 replicates for n > 1000.

All simulations were performed with an add-on to the HOOMD-blue simulation toolkit (version 3.10.0)⁶⁴ that replaces translational moves in 3D Euclidean space with translational moves in S². This modified HPMC algorithm is modeled after a schema described for simulations in hyperspherical geometry^65,66 and provides a robust and efficient particle displacement method. Particle translation and rotation are uncoupled, allowing for full three-dimensional rotational degrees of freedom of the particles. All particles interact purely via excluded volume interactions.

We used the freud⁶⁷ python library to analyze our simulations and the signac⁶⁸ software package for data management. Simulation snapshots are visualized using Ovito visualization software.⁶⁹ There are numerous sophisticated order parameters available to quantify local order in the assemblies.^63,70–72 We found that the order parameters described below, while relatively simple, are fast, easily interpretable and sufficient for our purposes. For comparison to we applied dimensionality reduction methods (principal component analysis) to our systems and found it yields no additional information; (see supplementary S5).

All order parameters are normalized such that a value of 1 denotes perfect agreement with the desired order and 0 denotes low agreement. Although each order parameter is normalized to the same range, they use parameter specific cutoffs to denote particles within their respective lattice motif. These cutoffs are informed by probability distributions. Detailed explanation of order parameter calculations and justifications for cutoffs are included in the SI Section S2. The variety of per-particle order parameters described are averaged over the first nearest-neighbor shell of the particle using a distance cutoff to give a subscript nn. This cutoff was calculated to be between the first two peaks of the geodesic radial distribution function g(θ). The geodesic distance in radians is the length of the shortest curve that connects two points on the sphere surface and is size-agnostic. We use g(θ) for simple comparison of translational order across different sphere radii, and symmetry measurements for the distribution of defects on the surface. Using local symmetry and shape vector information, Fig. 1c, we construct four different order parameters to describe the four crystalline motifs found in this study:

1. The hexagonal motif: the hexagonal order parameter, |ψ₆|_nn, quantifies local six-fold orientational symmetry via a per-particle value of . We assert that hexagonal packing describes particles with |ψ₆|_nn ≥ 0.80 and consider particles with |ψ₆|_nn < 0.80 to be defects. This cutoff is in agreement with previous studies of hard-sphere-like colloidal particles of silica and polymethylmethacrylate.⁷³

2. The face-aligned motif: we quantify the face-to-face alignment of rounded cubes by the order parameter f_nn. A high value of f_nn implies the following. First, one face of the cube f⃑ is relatively parallel to the normal vector out of the spherical interface, making the cube effectively “face-up” relative to the surface. Second, for each nearest neighbor within a cutoff distance face normals that are roughly tangential to the spherical interface are all in line with one another subject to the symmetry of the cube θ_ij = 〈f_i,f_j〉%(π/2). The order parameter f is then This asserts that face-up cubes are also aligned face to face with their neighbors, rotated to be in the same tangential reference plane. We use a cutoff of f_nn ≥ 0.75 to describe highly face-aligned rounded cubes, and consider particles with f_nn < 0.75 to be defects.

3. The honeycomb motif: The honeycomb order parameter, hc_nn, captures the high local three-fold orientational symmetry of neighboring tetrahedra. This motif has a restriction that particles must be vertex-up or face-up, corresponding to the basis positions of the honeycomb point set. We calculate the approximate orientation of each particle as either “face-up”, “edge-up”, or “vertex-up” using f⃑, , v⃑ relative to the surface normal. For each particle either “vertex-up” or “face-up” we calculate given its three nearest neighbors are of the opposite orientation. We use a cutoff of hc_nn ≥ 0.75 to describe the honeycomb motif, and consider particles with hc_nn < 0.75 to be defects. This motif is similar to the honeycomb phase of tetrahedron superlattices on substrates.⁷⁴

4. The woven motif: the woven order parameter, w_nn, for rounded tetrahedra is described using a two-particle arrangement of “edge up” tetrahedra rotated by π/2 relative to each other. This is calculated between “edge-up” tetrahedra with θ_ij = 〈e_‖i,e_‖j〉%(π/2) and the woven parameter is then . This order parameter captures the motif's distinct “interwoven” appearance and tetratic symmetry. We use a cutoff of w_nn ≥ 0.65 to describe the woven motif, and consider particles with w_nn < 0.65 to be defects.

Description of defect morphology. We expand the concept of defect scars to examine the shape of the defects for spheres with s = 1, cubes with scube = 0, and tetrahedra with s_tetrahedron = 0. We binarize the systems using |ψ₆|_nn < 0.80, f_t,nn < 0.75, and w_nn < 0.65 for the sphere, cube, and tetrahedron systems, respectively. We do not consider defect morphologies of honeycomb dominated systems due to a lack of isolated defect regions. We construct a connected graph with these defect particles as nodes and connect these nodes based on a distance cutoff corresponding to the first neighbor peak within the geodesic radial distribution function. For each component of the graph (isolated defect) the defect length d[σ] is defined as the longest shortest path. The shortest longest path is formulated as such: for a graph compute the shortest connected path between all pairs of nodes i and j, the maximum of these values is then the final value of shortest longest path. We study defect region length as the system size, n ∈ 1000, 1500, 2000, 2500, increases over a range of ρ_n.

Results and discussion

We constructed phase diagrams for assemblies of n = 1000 particles for our two shape families of rounded polyhedra as shown in Fig. 2. The various phases reported in the figure demonstrate our ability to change the self-assembled crystal structure by varying particle shape and roundedness s.

Fig. 2 Phase diagrams for systems with n = 1000 particles over a range of particle roundedness s and number densities ρ_n. For cubic rounded polyhedra the face-aligned motif dominates at low rounding values s_cube, gradually shearing as the rounding of the shape increases, and exhibits high hexagonal ordering as s_cube approaches 1. Striped regions at s_cube = 0.3 and ρ_N ≥ 0.272 indicate regions where both hexagonal and face-aligned motifs are observed. When s = 1 (middle), idealized spheres assemble the hexagonal motif. At high values of s_tetrahedron the hexagonal motif continues to dominate. With intermediate rounding, tetrahedra assemble the honeycomb structure before giving way to the woven structure as s → 0. The four distinct motifs (right) are color coded to match the assembly phase diagram, with motif insets corresponding to the appropriate black-filled marker. These specific motifs correspond to the following values of (shape, s, ρ_n, marker): (sphere, 1.0, 0.258, yellow circle), (cube, 0.0, 0.282, dark blue diamond), (tetrahedron, 0.2, 0.272, light blue downward triangle), (tetrahedron, 0.0, 0.258, red upturned triangle) from left to right, top to bottom. The white filled marker at s_cube = 0.6 and ρ_n = 0.287 indicates an ordered jammed state rather than self-assembled structure. Jammed states are where there was no initial disordered condition. Particles cannot reconfigure to assemble and are ordered only as an effect of “jamming” during initialization system compression.

Results for spheres: the phase diagram in Fig. 2 reveals that the dominant lattice order for assembles of n = 1000 spheres is hexagonal across a range of number densities ρ_n ≥ 0.239, with a simulation image and unwrapped Mercator projection in Fig. 3a and b. The ordering of the hexagonal motif is supported by distinct peaks in the radial distribution function g(θ) seen in Fig. 3c. The hexagonal lattice is disrupted by defect scars that each have the expected +1 overall topological charge, and for sufficient ρ_n we identify 12 distinct scars with low |ψ₆|_nn. These 12 distinct defect regions are distributed with icosahedral symmetry, as seen in Fig. 3d. Peaks in the geodesic radial distribution function calculated between defects, g_DD(θ), appear at geodesic distances θ^I₁, θ^I₂, and θ^I₃. These distances are the geodesic distances between 12 points distributed on a sphere with icosahedral symmetry. As we increase system size by increasing n at constant ρ_n the defect scar length d[σ] increases as seen in Fig. 3e and f. This is in agreement with previous studies of spheres on the surface of a sphere³¹ and serves as a baseline for our studies of cubic and tetrahedral rounded polyhedra.

Fig. 3 Assemblies of spheres: (a) A representative simulation image of n = 1000 spheres and (b) its Mercator projection at a density of ρ_n = 0.258 and s = 1, colored by |ψ₆|_nn. (c) The translational order of the hexagonal motif is showcased in the geodesic radial distribution function g(θ) at these same parameters. (d) The distribution of defect regions with low |ψ₆|_nn is icosahedral across densities, seen by peaks in g_DD(θ) at geodesic distances corresponding to the vertices of an icosahedron. (e) The average length of defect regions d[σ] over a range of densities ρ_n and particle number n. The length of the defect reliably increases as a function of n for sufficiently high rho_n, elongating into commonly seen “scars”. (f) This increase in d[σ] is shown in representative defect images at n = 1000 (left) and 2500 (right).

Results for rounded cubes: we investigate the family of rounded cubes by decreasing s_cube, with constant n = 1000. For highly rounded cubes with s_cube ≥ 0.4 a hexatic rotator phase dominates, similar to experimental and simulated systems of rounded cubes⁵⁵ on flat substrates and superballs⁴³ in bulk. This phase has high hexatic order for ρ_n ≥ 0.244 and little orientational correlation between particles, as indicated by the high and low global averages of 〈|ψ₆|_nn〉 and 〈f_nn〉, respectively. The crystallization transition gradually increases to a maximum value ρ_n = 0.272 at s_cube = 0.3. At this rounding value the particles form both hexagonal and face-aligned motifs, indicating a coexistence of spatially disparate domains within one system. The boundary between the two phases can have a region of disordered particles, a region of particles that have high orientational and rotational order, or can be a sharp transition between hexagonal and face-aligned, see SI Fig. S21. As the particle shape approaches a sharp cube with s_cube = 0, face-aligned order becomes dominant and gives rise to a simple square lattice with tetratic symmetry and high orientational order, shown in Fig. 2 and 4a–c. Furthermore, the onset number density for crystallization rapidly decreases to ρ_n = 0.239 at s_cube = 0.

Fig. 4 Assemblies of rounded cubes: (a) a representative simulation image of n = 1000 cubes and (b) its Mercator projection at a density of ρ_n = 0.282 and s = 0, colored by f_nn. (c) The simple-square translational order of the face-aligned motif is showcased in the geodesic radial distribution function g(θ) at these same parameters. (d) The distribution of eight defect regions with low f_nn is square antiprismatic, as seen in peaks in g_DD(θ) at , θ^AP₂, and θ^AP₄. There are no peaks at and θ^C₅ that would indicate a cubic distribution. (e) In systems of cubes with s_cube = 0, defect lengths d[σ] can be seen to generally decrease as ρ_n increases, across all values of n. (f) The length d[σ] is consistent across n, shown in representative defect images at n = 1000 (left) and 2500 (right).

From the simulation image in Fig. 4a and its Mercator projection in Fig. 4b, we see isolated defect regions colored by low values of f_nn. The preferred local symmetry around particles changes from hexagonal to tetragonal, leading to an expected coordination number of c_i = c_square = 4. From the topological expectations that arise from the Euler characteristic χ = 2 and the dominant tetragonal symmetry, we rationalize the symmetric distribution of eight three-fold symmetric q = +1 topological charges, Fig. 4b, similar to k = −1/2 defects seen in nematic liquid crystals.⁷⁵ Spatial correlations between face-aligned defect particles as quantified by g_DD(θ) in Fig. 4d show an organized distribution of defects about the surface. We hypothesize that there are two candidates for the distributed symmetry of these defects, one aligning to the vertices of a cube and the other to a square antiprism. Geodesic distances for these symmetries are denoted by θ_i, i ∈ [1, 2, 3, 4, 5], and the superscript denotes whether this distance is attributed to cubic (C) symmetry, square antiprismatic (AP) symmetry, or both (*). Peaks in g_DD(θ) show high correlations corresponding to ,θ^AP₂, and θ^AP₄, demonstrating that the distribution aligns more closely with the symmetry of the square antiprism, reminiscent of a study of hard squares on the sphere by Li et al.³⁹ This distribution of defects has more space between defects and is the solution to the Thompson problem for eight repulsive points distributed about a sphere surface.

The three-fold nature of the defect region is showcased in Fig. 4e, appearing through a cooperative loss of orientational order and attributed to topological resolutions to lattice incompatibility. The defect has either one or three particles centered about the geometric center of the defect; leading to the two peaks in d[σ] shown in the SI, Fig. S13. These central particles help stabilize the defect and rotate more freely compared to particles on the periphery or outside of the defect, rather than forming as a result of jamming between nearby misoriented grains. As particle number n increases, the absolute length of the defect d[σ] remains constant, but decreases with increasing ρ_n as seen in Fig. 4f.

We observed no smooth transition from five-fold defect scars to localized three-fold symmetric defects with corresponding icosahedral and square antiprismatic defect distribution symmetry as s_cube decreases. Instead, we observed a dissolution of one defect morphology and a subsequent emergence of the counter morphology. The dominant translational order transition between hexagonal and simple-square is different to those of superball assemblies in bulk⁴² where lattice transformation occurs via a smooth variable angle shear square transition as the roundedness of the particles decreases. Confined to the surface of a sphere, bounded system size and required topological defects hinder the system's ability to smoothly interpolate between these two ends of the spectrum. We also note the monotonic increase in allowable ρ_n as s_cube decreases, which follows from the high maximum number density ρ_n = 0.385 given in the flat R = ∞ limit with defect-free simple-square packing.

Results for rounded tetrahedra: we repeated the simulations with the family of rounded tetrahedra by starting from sphere particle systems of size n = 1000 and decreasing s_tetrahedron. Rounded tetrahedra show very distinct phase regions as a function of s_tetrahedron; we observed no gradual decrease of one order parameter corresponding to an increase in another. This leads to a phase diagram with large regions that are geometrically inaccessible for rounded tetrahedra but accessible to rounded cubes.

As we begin to decrease s_tetrahedron the hexatic rotator phase quickly gives way to a disordered phase. This rotator phase is disallowed at high ρ_n and moderate s_tetrahedron due to the effective radius of the particle as it rotates. The circumsphere radius of rounded tetrahedra grows much more quickly for tetrahedra as compared to cubes, and correspondingly the maximally allowed ρ_n for perfect hexagonal packing sharply decreases with decreasing s_tetrahedron. Decreasing s_tetrahedron further, we see the emergence of a honeycomb phase, Fig. 5a–c. The first neighbor shell includes three particles that align face to face with the center particle, and are at distances larger than the minimal steric distance. For example, given ρ_n = 0.263 and s_tetrahedron = 0.3, the average Euclidean distance to the first nearest neighbors is 1.754, compared to a minimal Euclidean distance of 1.434. This phase does not span the surface at any simulated ρ_n, as evidenced by the lack of long range correlations in g(θ), Fig. 5d, and have only ∼40% surface coverage in the best case. Consequentially, the assembled structures lack clearly isolated defect regions and have no spatial correlations in their defect distribution, Fig. 5d, even though they share three-fold symmetry with the hexagonal motif. Well-defined defects are unstable and fluid in position. The system in Fig. 5a and b has exceptional honeycomb coverage compared to the majority of replicate simulations.

Fig. 5 Assemblies of Rounded Tetrahedra: (a) A representative simulation image of n = 1000 tetrahedra and (b) its Mercator projection at a density of ρ_n = 0.267 and s_tetrahedraon = 0.2, colored by hc_nn. (c) The translational order of the honeycomb motif is showcased in the geodesic radial distribution function g(θ) at these same parameters. (d) In honeycomb systems with s = 0.2 there is no regular distribution of defects about the surface, confirmed by a flat g_DD(θ). (e) A representative simulation image of n = 1000 tetrahedra and (f) its Mercator projection at a density of ρ_n = 0.263 and s_tetrahedraon = 0.0, colored by w_nn. (g) The simple-square translational order of the woven motif is showcased in the geodesic radial distribution function g(θ) at these same parameters. (h) In woven systems with s_tetrahedraon = 0, eight isolated defects appear only at a high density ρ_n = 0.263. The spatial distribution of defects show little to no correlation across the surface for all densities ρ_n < 0.263, with the first trough in g_DD(θ) indicating only that defect regions are spatially isolated from each other. (i) Tetrahedra with s_tetrahedraon = 0 form a large variety of defects, with broadly distributed lengths d[σ], seen for ρ_n = 0.263. (j) A range of exemplary defects are shown over the range n = 1000, 1500, 2000, and 2500 (left to right).

As s_tetrahedron approaches 0 the woven motif, Fig. 5e and f, appears for a small range of number densities 0.248 ≤ ρ_n ≤ 0.263. This motif has tetratic symmetry and forms a simple square lattice as evidenced by the peaks in Fig. 5g, bearing a marked similarity to the simple square lattice of face-aligned cubes in Fig. 4c. When rounding is finely varied within the region 0 < s_tetrahedron < 0.1 in 0.01 increments (see SI Fig. S10c and b), we see that the woven motif is highly sensitive to particle rounding and persists for a small window of s_tetrahedron < 0.06 before the honeycomb motif returns. For ρ_n > 0.263 and s_tetrahedron = 0 the surface area constraint prevents further assembly. The woven motif forms a tightly bound network, maintaining favorable face-to-face alignment with four neighbors. The minimum distance between particles in this motif with s_tetrahedron = 0 would be 1.644σ. Yet, at s_tetrahedron and ρ_n = 0.263 the average distance to the first neighbor shell is 1.914σ, larger than the minimum distance due to warping of the woven motif due to the surface curvature.

For sharp tetrahedra with s_tetrahedron = 0, a single, surface-spanning woven cluster forms only at a density of ρ_n = 0.263. Although the woven motif has tetratic symmetry, defects rarely exhibit the topologically expected three-fold symmetric morphology and even more rarely form eight ideally isolated defects. The instance of clearly isolated defects given in Fig. 5e and f is exceptional. In spite of this, peaks in g_DD(θ) show that defects in systems of woven tetrahedra display little to no long range correlations, Fig. 5h. The defect-ridden assembly of the woven phase is attributed to the lack of reconfigurability of the lattice. In comparison to the face-aligned simple square lattice seen in cubes, the woven motif is translationally and orientationally rigid in its construction; this prevents the resolution of misalignment between two grains of woven tetrahedra by sliding. From an elastic perspective, the woven tetrahedra lattice cannot locally splay nor bend to a large extent without generation overlaps or by lowering the number density so much that the entire motif dissolves. As the number of particles n increases the defects in ordering continue to lack a clear morphology, and the distribution of defect lengths d[σ] remains wide, Fig. 5i. Some of these defects are highly localized, but a variety of morphologies are present, including elongated scars and three-fold chiral spirals, shown in Fig. 5j We believe this is due to the difference in crystallization pathway, woven tetrahedra grains grow from nucleated “seeds” that spread across the surface, and once a particle has crystallized its position and orientation relative to its neighbors remains largely unchanged. As the grains continue to grow they come into contact with one another and create larger, irregular grain boundaries and irregular enclosed defects regions. This is contrasted by the coordinated formation of defects in systems of cubes.

Conclusions

We have presented work that seeks to provide insight into the effect of spherical topology on the self-assembled ordering of rounded polyhedra and the distribution of topological defects. Hard spheres exhibit hexagonal order, and distribute defects with icosahedral symmetry. By continuously changing the particle shape to cubes, the entropic preference for face-to-face alignment drives the system towards simple square order, with a notable lack of a smooth shear-square transition. Nonetheless, the topological incompatibility of the local tetragonal symmetry seen in cubic rounded polyhedra and the surface drives the formation of eight three-fold symmetric isolated defects, distributed with approximately square antiprismatic symmetry about the surface. By using shape as a continuous tunable parameter, we have effectively bridged the gap between what is known for disk and square packing on the sphere. For our other shape family, as we morph spherical particles into tetrahedra we see the appearance of the honeycomb and woven motifs. These self-assembled motifs are highly dependent on utilizing the full three-dimensional rotational degrees of freedom of the tetrahedron, unlike cubes. For n = 1000 particles, both motifs give rise to incoherent and non-distinct defect regions. As we increase n for sharp tetrahedra we continue to see a wide distribution of defect morphologies, only rarely observing distinctly three-fold symmetric defects. Overall, the reconfigurability of the assembled lattice impacts how evenly topological strain is localized and distributed. Overall, tetrahedra assemblies are significantly more sensitive to rounding, a great insight for experimentalists who want to synthesize specific structures.

In systems of anisotropic particles at interfaces, entropic interactions are only one of many important contributions. Even particle orientation and depth of the particle centroid within the interface become functions of particle wettability and fluid–fluid surface tension.^76–78 Theoretical studies of cubic particle assemblies at flat fluid–fluid interfaces have predicted orientations relative to the interface normal (face, edge, or vertex up) with varying frequency.⁷⁹ Interparticle interactions dictated by orientationally-dependent induced capillary forces between particles can lead to a variety of non-trivial assemblies, as observed in experimental systems of cuboidal hematite particles.⁸⁰ Incorporating each of these phenomena, penetration depth, orientation, and capillary interactions, would require a complex feedback loop so complex and highly system dependent that it warrants an independent study. A similar orientational constraint also occurs for crystals of anisotropic particles confined within a spherical droplet.^81–83 It has been shown that the curvature of the droplet surface imparts a fixed, preferred orientation to each particle at the boundary, which in turn affects the subsequent inner layers of particles. As a result, the positional ordering of the outermost boundary layer is dominated by the orientational constraints, leading to a different set of surface ordering compared to the freely rotating particles in our simulations. Again, while this is most evident in the tetrahedral case, it is not a given that cubes would preferentially align in a face-up configuration.

As the separation between defects drives buckling transitions in icosahedral virus capsids, our observation of eight consistent three-fold defects opens avenues for new buckling phenomena. Vesicles with a faceted cuboid morphology have been observed in phospholipid liposomes by enhancing intra-membrane attraction forces to overcome surface tension, yet there is variation in aspect ratio of these vesicles.⁸⁴ Furthermore, complicated buckling modes have been shown to appear in simulated and natural systems, such as the infolding of desiccated pollen grains, preventing vesicle collapse and allowing for later re-hydration.^85–89

We have shown that shape anisotropy is a strong design parameter to influence defect morphology and distribution. Our findings within this minimal entropic model show that defect distribution varies in non-trivial ways and the morphologies that have been shown in the paper are unlike those in systems of spherical particles. The total number of isolated defects and their characteristic topological charge vary as a function of particle number, intrinsic shape, and particle rounding, showing that these tunable parameters dictate particle order and defect distribution and separation. By showing the self-assembly behavior of crystal motifs and the distribution of defects in hard particle systems to be designable and robust we hope to inspire new directions in colloidosome design.

Author contributions

GNJ and PWAS were both heavily involved in the conceptualization of the study and software development. GNJ performed the simulations and analyzed the data. All authors contributed to the investigation and the writing of the manuscript. SCG directed the research.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data and software for this article, including HOOMD-blue source code and python code for data and analysis generation are available at Deep Blue Data at https://doi.org/10.7302/aagt-dh84.

Additional analysis and simulation data supporting this article has been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5sm01271f.

Acknowledgements

We would like to acknowledge interesting discussions on this work with Prof. Jonathan Selinger, Prof. Xiaoming Mao, and Dr Nan Cheng.

Research supported in part by the Center for Complex Particle Systems (COMPASS) a National Science Foundation Science and Technology Center, Division of Materials Research Award #DMR 2243104. GNJ acknowledges support from the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Department of Energy Computational Science Graduate Fellowship under Award Number DE-SC0022158. This work used Anvil at Purdue RCAC through allocation DMR 140129 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. Computational resources and services were also provided by Advanced Research Computing (ARC), a division of Information and Technology Services (ITS) at the University of Michigan, Ann Arbor.

Disclaimer: “This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof”.

References

1. E. Guzmán, F. Ortega and R. G. Rubio, Cosmetics, 2022, 9, 68.
2. F. B. De Carvalho-Guimarães, K. L. Correa, T. P. De Souza, J. R. Rodríguez Amado, R. M. Ribeiro-Costa and J. O. C. Silva-Júnior, Pharmaceuticals, 2022, 15, 1413.
3. W. Ramsden, Proc. R. Soc. London, 1904, 72, 156–164.
4. S. U. Pickering, J. Chem. Soc., Trans., 1907, 91, 2001–2021.
5. D. E. Discher and F. Ahmed, Annu. Rev. Biomed. Eng., 2006, 8, 323–341.
6. S. J. Hunter and S. P. Armes, Langmuir, 2020, 36, 15463–15484.
7. S. Simovic, N. Ghouchi-Eskandar and C. Prestidge, J. Drug Delivery Sci. Technol., 2011, 21, 123–133.
8. Z. Xu, T. Hueckel, W. T. M. Irvine and S. Sacanna, Chem. Mater., 2023, 35(16), 6357–6363.
9. Z. Sun, U. Glebe, H. Charan, A. Böker and C. Wu, Angew. Chem., Int. Ed., 2018, 57, 13810–13814.
10. L. Tran and K. J. M. Bishop, ACS Nano, 2020, 14, 5459–5467.
11. I. Kalashnikova, H. Bizot, B. Cathala and I. Capron, Langmuir, 2011, 27, 7471–7479.
12. H. A. Wege, S. Kim, V. N. Paunov, Q. Zhong and O. D. Velev, Langmuir, 2008, 24, 9245–9253.
13. F. Qi, J. Wu, G. Sun, F. Nan, T. Ngai and G. Ma, J. Mater. Chem. B, 2014, 2, 7605–7611.
14. M. Pang, A. J. Cairns, Y. Liu, Y. Belmabkhout, H. C. Zeng and M. Eddaoudi, J. Am. Chem. Soc., 2013, 135, 10234–10237.
15. X. Cai, X. Deng, Z. Xie, S. Bao, Y. Shi, J. Lin, M. Pang and M. Eddaoudi, Chem. Commun., 2016, 52, 9901–9904.
16. X. Yue, X. Liu, N. Yan and W. Jiang, Chem. Commun., 2020, 56, 12737–12740.
17. T. Nallamilli, E. Mani and M. G. Basavaraj, Langmuir, 2014, 30, 9336–9345.
18. J. S. Sander and A. R. Studart, Langmuir, 2013, 29, 15168–15173.
19. S. Li, B. A. Moosa, J. G. Croissant and N. M. Khashab, Angew. Chem., Int. Ed., 2015, 54, 6804–6808.
20. A. San Miguel and S. H. Behrens, Soft Matter, 2011, 7, 1948–1956.
21. X. Huang, W. Liu, Z. Zheng, Y. Zhang, X. Zhao, Z. Li, W. Gao, P. Xie, Y. Pan and Q. Sun, Colloids Surf., A, 2024, 702, 135107.
22. M. F. Hsu, M. G. Nikolaides, A. D. Dinsmore, A. R. Bausch, V. D. Gordon, X. Chen, J. W. Hutchinson, D. A. Weitz and M. Marquez, Langmuir, 2005, 21, 2963–2970.
23. A. D. Dinsmore, M. F. Hsu, M. G. Nikolaides, M. Marquez, A. R. Bausch and D. A. Weitz, Science, 2002, 298, 1006–1009.
24. D. S. Roshal, K. Azzag, E. Le Goff, S. B. Rochal and S. Baghdiguian, Sci. Rep., 2020, 10, 7652.
25. A. Klug and D. Caspar, Advances in Virus Research, Elsevier, 1961, vol. 7, pp. 225–325.
26. M. Martín-Bravo, J. M. G. Llorente, J. Hernández-Rojas and D. J. Wales, ACS Nano, 2021, 15, 14873–14884.
27. R. Zandi, B. Dragnea, A. Travesset and R. Podgornik, Phys. Rep., 2020, 847, 1–102.
28. R. Twarock and A. Luque, Nat. Commun., 2019, 10, 4414.
29. T. Einert, P. Lipowsky, J. Schilling, M. J. Bowick and A. R. Bausch, Langmuir, 2005, 21, 12076–12079.
30. E. Katifori, S. Alben, E. Cerda, D. R. Nelson and J. Dumais, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 7635–7639.
31. R. E. Guerra, C. P. Kelleher, A. D. Hollingsworth and P. M. Chaikin, Nature, 2018, 554, 346–350.
32. J. E. Bond, A. J. Yeh, J. R. Edison and M. A. Bevan, Soft Matter, 2024, 20, 6371–6383.
33. P. H. J. Pedersen, SPHM, 1982, vol. 8, pp. 1–17.
34. A. R. Bausch, Science, 2003, 299, 1716–1718.
35. M. J. Bowick and L. Giomi, Adv. Phys., 2009, 58, 449–563.
36. J. Wang, C. F. Mbah, T. Przybilla, B. Apeleo Zubiri, E. Spiecker, M. Engel and N. Vogel, Nat. Commun., 2018, 9, 5259.
37. B. W. Clare and D. L. Kepert, J. Math. Chem., 1991, 6, 325–349.
38. H. Xie, W. Liu, Z. Lu, J. Z. Y. Chen and Y. Li, Nano Lett., 2025, 25, 1193–1198.
39. Y. Li, H. Miao, H. Ma and J. Z. Y. Chen, Soft Matter, 2013, 9, 11461.
40. S. C. Glotzer and M. J. Solomon, Nat. Mater., 2007, 6, 557–562.
41. P. F. Damasceno, M. Engel and S. C. Glotzer, Science, 2012, 337, 453–457.
42. Y. Zhang, F. Lu, D. van der Lelie and O. Gang, Phys. Rev. Lett., 2011, 107, 135701.
43. R. Ni, A. P. Gantapara, J. De Graaf, R. Van Roij and M. Dijkstra, Soft Matter, 2012, 8, 8826.
44. A. Haji-Akbari, M. Engel, A. S. Keys, X. Zheng, R. G. Petschek, P. Palffy-Muhoray and S. C. Glotzer, Nature, 2009, 462, 773–777.
45. E. R. Chen, M. Engel and S. C. Glotzer, Discret. Comput. Geom., 2010, 44, 253–280.
46. E. R. Chen, D. Klotsa, M. Engel, P. F. Damasceno and S. C. Glotzer, Phys. Rev. X, 2014, 4, 011024.
47. F. Smallenburg, L. Filion, M. Marechal and M. Dijkstra, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 17886–17890.
48. C. Avendaño and F. A. Escobedo, Soft Matter, 2012, 8, 4675.
49. M. Dijkstra, Proc. Int. Sch. Phys. “Enrico Fermi”, 2013, 184, 229–244.
50. E. G. Teich, G. Van Anders, D. Klotsa, J. Dshemuchadse and S. C. Glotzer, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, E669–E678.
51. J. A. Anderson, J. Antonaglia, J. A. Millan, M. Engel and S. C. Glotzer, Phys. Rev. X, 2017, 7, 021001.
52. W. Shen, J. Antonaglia, J. A. Anderson, M. Engel, G. van Anders and S. C. Glotzer, Soft Matter, 2019, 15, 2571–2579.
53. D. Wang, M. Hermes, S. Najmr, N. Tasios, A. Grau-Carbonell, Y. Liu, S. Bals, M. Dijkstra, C. B. Murray and A. Van Blaaderen, Nat. Commun., 2022, 13, 6001.
54. A. K. Sharma and F. A. Escobedo, J. Chem. Phys., 2023, 158, 044502.
55. Y. Zhong, T. C. Moore, T. Dwyer, A. Butrum-Griffith, V. R. Allen, J. Chen, Y. Wang, F. Cheng, S. C. Glotzer and X. Ye, Nat. Chem. Eng., 2024, 1, 532–541.
56. D. Rajendra, J. Mandal, Y. Hatwalne and P. K. Maiti, Soft Matter, 2023, 19, 137–146.
57. T. G. Anjali and M. G. Basavaraj, Phys. Chem. Chem. Phys., 2017, 19, 30790–30797.
58. X.-C. Luu and A. Striolo, J. Phys. Chem. B, 2014, 118, 13737–13743.
59. F. Smallenburg and H. Löwen, J. Chem. Phys., 2016, 144, 164903.
60. M. A. Bates, J. Chem. Phys., 2008, 128, 104707.
61. H. Shin, M. J. Bowick and X. Xing, Phys. Rev. Lett., 2008, 101, 037802.
62. J. Mandal, C. Dasgupta and P. K. Maiti, Soft Matter, 2025, 21, 6458–6468.
63. R. Van Damme, G. M. Coli, R. Van Roij and M. Dijkstra, ACS Nano, 2020, 14, 15144–15153.
64. J. A. Anderson, J. Glaser and S. C. Glotzer, Comput. Mater. Sci., 2020, 173, 109363.
65. D. W. Sinkovits and E. Luijten, Nano Lett., 2012, 12, 1743–1748.
66. P. W. Schönhöfer, K. Sun, X. Mao and S. C. Glotzer, Phys. Rev. Lett., 2023, 131, 258201.
67. V. Ramasubramani, B. D. Dice, E. S. Harper, M. P. Spellings, J. A. Anderson and S. C. Glotzer, Comput. Phys. Commun., 2020, 254, 107275.
68. C. S. Adorf, P. M. Dodd, V. Ramasubramani and S. C. Glotzer, Comput. Mater. Sci., 2018, 146, 220–229.
69. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.
70. A. Haji-Akbari and S. C. Glotzer, J. Phys. A: Math. Theor., 2015, 48, 485201.
71. D. Fijan, M. R. W. Rashidi and S. C. Glotzer, Quantifying Local Point-Group-Symmetry Order in Complex Particle Systems, arXiv, preprint, arXiv:2509.12665, 2025, https://arxiv.org/abs/2509.12665, [physics].
72. F. Camerin, S. Marín-Aguilar and M. Dijkstra, Nanoscale, 2024, 16, 4724–4736.
73. Y. L. Wu, D. Derks, A. Van Blaaderen and A. Imhof, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 10564–10569.
74. S. Zhou, J. Li, J. Lu, H. Liu, J.-Y. Kim, A. Kim, L. Yao, C. Liu, C. Qian, Z. D. Hood, X. Lin, W. Chen, T. E. Gage, I. Arslan, A. Travesset, K. Sun, N. A. Kotov and Q. Chen, Nature, 2022, 612, 259–265.
75. X. Tang and J. V. Selinger, Soft Matter, 2017, 13, 5481–5490.
76. L. Botto, E. P. Lewandowski, M. Cavallaro and K. J. Stebe, Soft Matter, 2012, 8, 9957.
77. T. G. Anjali and M. G. Basavaraj, J. Colloid Interface Sci., 2016, 478, 63–71.
78. H. Kumar, V. R. Dugyala and M. G. Basavaraj, Langmuir, 2021, 37, 7295–7304.
79. C. Anzivino, G. Soligno, R. van Roij and M. Dijkstra, Soft Matter, 2021, 17, 965–975.
80. T. G. Anjali and M. G. Basavaraj, Langmuir, 2017, 33, 791–801.
81. D. Wang, M. Hermes, R. Kotni, Y. Wu, N. Tasios, Y. Liu, B. De Nijs, E. B. Van Der Wee, C. B. Murray and M. Dijkstra, et al., Nat. Commun., 2018, 9, 2228.
82. R. S. Skye, E. G. Teich and J. Dshemuchadse, Soft Matter, 2022, 18, 6782–6790.
83. S. Schyck, J.-M. Meijer, L. Baldauf, P. Schall, A. V. Petukhov and L. Rossi, JCIS Open, 2022, 5, 100037.
84. P. Bakardzhiev, A. Forys, B. Trzebicka, T. Andreeva and S. Rangelov, Nanoscale, 2021, 13, 15210–15214.
85. H. Yuan and M. Olvera De La Cruz, Phys. Rev. E, 2019, 100, 012610.
86. G. Vernizzi, R. Sknepnek and M. Olvera De La Cruz, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 4292–4296.
87. A. Božic and A. Šiber, Biophys. J., 2022, 121, 782–792.
88. D. A. Matoz-Fernandez, S. Li, M. O. D. L. Cruz and R. Sknepnek, PyMembrane: A flexible framework for efficient simulations of elastic and liquid membranes, preprint, arXiv:2308.12754, 2023, https://arxiv.org/abs/2308.12754, [cond-mat].
89. C. Waltmann, A. Shrestha and M. Olvera De La Cruz, Phys. Rev. E, 2024, 109, 054409.

---