Crystalline Particle Packings on Constant Mean Curvature (Delaunay) Surfaces
Abstract
We investigate the structure of crystalline particle arrays on constant mean curvature (CMC) surfaces of revolution. Such curved crystals have been realized physically by creating chargestabilized colloidal arrays on liquid capillary bridges. CMC surfaces of revolution, classified by Charles Delaunay in 1841, include the 2sphere, the cylinder, the vanishing mean curvature catenoid (a minimal surface) and the richer and less investigated unduloid and nodoid. We determine numerically candidate ground state configurations for 1000 pointlike particles interacting with a pairwiserepulsive potential, with distance measured in 3dimensional Euclidean space . We mimic stretching of capillary bridges by determining the equilibrium configurations of particles arrayed on a sequence of Delaunay surfaces obtained by increasing or decreasing the height at constant volume starting from a given initial surface, either a fat cylinder or a square cylinder. In this case the stretching process takes one through a complicated sequence of Delaunay surfaces each with different geometrical parameters including the aspect ratio, mean curvature and maximal Gaussian curvature. Unduloids, catenoids and nodoids all appear in this process. Defect motifs in the ground state evolve from dislocations at the boundary to dislocations in the interior to pleats and scars in the interior and then isolated 7fold disclinations in the interior as the capillary bridge narrows at the waist (equator) and the maximal (negative) Gaussian curvature grows. We also check theoretical predictions that the isolated disclinations are present in the ground state when the surface contains a geodesic disc with integrated Gaussian curvature exceeding . Finally we explore minimal energy configurations on sets of slices of a given Delaunay surface and obtain configurations and defect motifs consistent with those seen in stretching.
pacs:
I Introduction
Much has been learned in the last decade about crystalline particle packings on a wide variety of curved twodimensional surfaces Bowick and Giomi (2009). The surfaces studied include the 2sphere (constant positive Gaussian curvature) Bowick et al. (2000, 2002); Bausch et al. (2003); Bowick et al. (2007), the 2torus (variable positive and negative Gaussian curvature with vanishing integrated Gaussian curvature) Giomi and Bowick (2008a, b); Pairam and FernandezNieves (2009); Yao and Bowick (2011), the paraboloid (variable Gaussian curvature and a boundary) Giomi and Bowick (2007); Bowick et al. (2008), the Gaussian bump (positive and negative Gaussian curvature) Sadoc and Charvolin (1989); Vitelli and Nelson (2004); Vitelli et al. (2006) and the catenoid minimal surface Irvine et al. (2010); Bowick and Yao (2011); Irvine and Vitelli (2012). Although the nature of condensed matter order on surfaces is applicable to many different physical settings the richest comparison between theoretical/numerical predictions and experiments has been made with micronscale colloidal systems. In these systems solid colloidal particles selfassemble at the interface of two distinct immiscible liquids, either in particlestabilized (Pickering) emulsions Binks and Horozov (2006); Pickering (1907); Dinsmore et al. (2002) or chargestabilized emulsions Leunissen et al. (2007a, b). The case of ordering on the 2sphere (the surface of a ball) is realized by particles selfassembling at the surface of a droplet held perfectly spherical by surface tension Bausch et al. (2003). The ordered configurations of particles may be imaged with confocal microscopy and the particles even manipulated with laser tweezers Irvine et al. (2012). In chargestabilized emulsions particles do not even wet the surface and so preserve the intrinsic shape of the droplet. The natural range in size of the droplets allows one to explore the effect of the Gaussian curvature, which varies as , on the ordering. The ground state has been found to exhibit qualitatively different types of defect arrays as the size changes Bowick et al. (2006); Bausch et al. (2003). Small droplets have 12 isolated 5fold coordinated disclinations (5s) whereas for larger droplets the isolated 5s sprout additional chains of dislocations (a dislocation is a tightly bound disclinationantidisclination or 57 pair) to become 12 grain boundary scars, each of which has a net disclination charge of (i.e. one excess 5) and freely terminates within the medium. Free termination of a grain boundary is possible in curved space because Bragg rows are curved along the converging geodesics of the 2sphere, eventually healing the mismatch in crystallographic axes found at the center of the grain boundary.
Recently a very rich and flexible experimental system has been developed Irvine et al. (2010) in which chargestabilized emulsions are created in the form of capillary bridges  these are structures in which a drop of liquid A, immersed in liquid B, spans the gap between two parallel flat glass surfaces. The surface separating the two liquids has the topology of a cylinder and necessarily has a constant mean curvature (CMC) determined by the Laplace pressure difference between the inside (A) and the outside (B) liquids. These CMC surfaces of revolution were classified in the 19th century by Charles Delaunay and come in 5 classes: the sphere, the cylinder, the catenoid, the unduloid and the nodoid Eells (1987); Delaunay (1841); Sturm (1841). Negative Gaussian curvature is expected to give rise to quite different structures from surfaces with positive Gaussian curvature Nelson (1983); Rubinstein and Nelson (1983); Bowick et al. (2000); Tarjus et al. (2012); Sausset et al. (2010); Irvine et al. (2010); Irvine and Vitelli (2012); Bowick and Yao (2011); Kusumaatmaja and Wales (2013). The simplest negative curvature surface is the constant negative curvature hyperbolic plane , the negative curvature analog of the 2sphere , but the isometric embedding of as a complete subset of Euclidean 3space is not differentiable Henderson and Taimina (2001) and may not be realizable physically. To illustrate one elementary feature of negative curvature surfaces we highlight in Fig. 1 a single geodesic triangle (in red) on a negative Gaussian curvature nodoid (blue)– the sum of the interior angles of this triangle is , which is a striking departure from the sum of a Euclidean triangle. Note also the diverging geodesics on the surface. The
In this paper we explore the structure of the crystalline ground state of particles strictly confined to a Delaunay CMC surface and interacting with a pairwiserepulsive shortrange power law potential. Such surfaces have varying negative Gaussian curvature and are thus technically more challenging to analyze and simulate numerically. We are particularly interested in the defect structure of the ground state and how distinctive defect motifs emerge as the integrated Gaussian curvature is varied within one class of CMC surface as well as the evolution of the ground state as one quasistatically varies the manifold by increasing the height of the capillary bridge: the sequence of Delaunay surfaces studied corresponds to the physical experiment of pulling the bounding plates slowly apart or pushing them slowly together and hence stretching or compressing the capillary bridge.
Ii Geometry of Delaunay surfaces
Delaunay surfaces are variationally determined by being constant mean curvature surfaces of revolution with minimum lateral area and fixed volume. The characterization that most directly gives rise to a simple and useful parametrization, however, is to consider them as the surfaces of revolution whose meridians are the roulettes of the foci of the conics, that is the curves traced by the foci of the conics as they roll on a conic tangent without slip (for a detailed study see Bendito et al. ()). Besides the extreme limiting cases of the sphere and the cylinder these surfaces are catenoids, unduloids and nodoids, generated by the roulettes of parabolas, ellipses and hyperbolas respectively.
The roulette of the focus of the parabola is the catenary, and the surface of revolution, the catenoid, is the best known nontrivial Delaunay surface. Its parameterization is (see Fig 2):
(1) 
where is the radius of the waist, and .
The parameterization of an unduloid is given by
(2) 
and are the semiaxes of the ellipse, and . If , the roulette generated by the closest focus to the tangent to the ellipse generates unduloids of the type shown in Fig. 3 (left) with negative Gaussian curvature. Taking, instead, , the same focus generates the unduloids with positive Gaussian curvature, such as the one shown in Fig. 3 (right). The same surfaces are obtained by using the other focus, but this time in the reverse order.
Similarly, the parameterization of a nodoid is given by
(3) 
where and are the semiaxes of the hyperbola, , and . The expressions in treat the roulette generated by the closest focus to the tangent to the hyperbola and yield nodoids with negative Gaussian curvature, as shown in Fig. 4 (left). Taking the roulette generated by the other focus yields the expressions and nodoids with positive Gaussian curvature, as illustrated in Fig. 4 (right).
(4) 
All the fundamental properties of the Delaunay surfaces can be obtained easily using these parameterizations. In particular, the mean curvature of catenoids, unduloids and nodoids is given by
(5) 
respectively, where is the semimajor axis of the related conic. The Gaussian curvatures of catenoids, unduloids and nodoids are given by
(6) 
Iii Numerical simulation: Finding a minimum energy configuration on a Delaunay surface
To determine candidate minimum energy configurations we use the Forces Method Bendito et al. (2007, 2009). The basic structure of the Forces Method is classical and explained in detail in Bendito et al. (2009). For completeness we give a brief description here. It can be viewed as a local relaxationgradientlike descent algorithm where each step consists of finding the update direction and the step size in a predefined way. It consists of four steps:

Choose a certain number of particles and an initial configuration for them

Update the positions of the particles in threedimensional space

Project the positions on the surface

Repeat steps 2 and 3 until a given threshold is reached
Briefly these are:
Initial configuration. The initial configuration is chosen to be as uniformly distributed on the given Delaunay surface as possible, with the restriction that the particles are not near the two boundaries. Many independent runs (order 10) are made starting from different initial configurations and the lowest energy configuration is selected.
Update and projection. For one particle
where , and are the position, update direction and step size at the th step, respectively. would be the new location of the particle if this location were on the surface. Since it is generically off the surface after the update the actual position is obtained by projecting onto the surface.
The update direction is in the direction of the net force on the particle following from the gradient of the potential. In this paper, as noted before, we use a power law potential. For a system of charge 1 particles, the potential energy is then given by , where is the potential created at by all the other particles. We denote by minus the gradient of in terms of the position of the th particle.
If the particles lie on a regular surface , then equilibrium is reached when the component of tangential to , , vanishes at . Then we choose as the step direction, where (we normalize by the maximum force). The magnitude of the step size is then , where the coefficient is a constant positive scalar. The existence of a minimum distance between particles makes it possible to adjust the step size to most efficiently access the various configurations that arise during the iterative process. The error at iteration is . The algorithm stops when reaches a certain prescribed threshold value .
The 1000 particle simulations are performed on an Intel core i7 processor at a speed of 36.5 iterations/sec. We simulate 23 distinct particle systems in all requiring a total of iterations for a total CPU time of 42 hours. ^{1}^{1}footnotetext: Coordinates and energies for the configurations analyzed in this paper may be obtained by contacting Mark Bowick (bowick@phy.syr.edu).
Iv Constant volume results
In this section we analyze the evolution of defect motifs on a capillary bridge as it is stretched while preserving the volume and the contact area with both bounding parallel flat plates (corresponding to the boundary disks of the Delaunay surface). At each stage of stretching the capillary bridge is a distinct Delaunay surface with different mean and Gaussian curvatures. Thus, in the experimental setting, the particles must continuously reequilibrate to a new ground state. Numerically we minimize the energy on each Delaunay surface separately. This analysis allows one to explore the structure of an entire sequence of Delaunay surfaces, one for each stage of the stretch. We study two types of constant volume stretching: one starts from a fat cylindrical capillary bridge and the other from a square cylindrical capillary bridge.
iv.1 Stretching from an initial fat cylindrical capillary bridge
The initial fat cylinder has contact radius and height , so both the volume and the aspect ratio (the ratio of the radius of the surface at an equatorial plane (waist) to the radius of the contact disk at the plates) are unity. In our numerical study starting from a fat cylindrical capillary bridge yields a rich family of surfaces that begins with the cylinder and subsequently deforms to a series of unduloids, the catenoid and then a series of nodoids. Thus all the classes of Delaunay surfaces, apart from the wellstudied 2sphere which can be generated by other simpler means, are explored in this single process! Minimum energy configurations for charged particles interacting via a Yukawa potential and confined to Delaunay surfaces were obtained in Kusumaatmaja and Wales (2013), but only for catenoids and unduloids. As is clear from Fig. 5, the predominant surfaces in stretching are nodoids. The proper treatment of nodoids is therefore essential to a comparison with the experiments of Ref.Irvine et al. (2010).
Planar sections that contain the common revolution axis of several surfaces from this family are shown in Fig. 5. The final unduloid has an aspect ratio .
We obtain minimum energy states for particles interacting via a pairwiserepulsive potential. Shortrange potentials (decaying faster than ) produce more uniform ground state configurations because longrange potentials (decaying slower than ) drive particles to the boundary circles to maximize the average separation. This can lead to defect structures that are influenced by the boundary – while of intrinsic interest as a boundary driven phenomenon this is not our primary focus in this paper. We found that longrange effects may also be reduced, but not eliminated, by adding a line of neutralizing charge along the central axis of the surface. Here we focus on the role of variable Gaussian curvature in the interior of Delaunay surfaces.
Various defect motifs emerge as the capillary bridge becomes higher and thinner (decreasing aspect ratio with our definition). Basically we identify the following sequence: at , individual dislocations (tightly bound 57 pairs) appear at the boundary, resulting from the repulsion between dislocations of identical orientation.
As decreases, these dislocations migrate to the interior of the surface and occasionally form multiple dislocation clusters. In Fig. 6(left, right) we show minimum energy configurations for and . The dual pentagons and heptagons shown here are obtained by first performing a Delaunay triangulation of the configuration and then connecting the barycenters of adjacent triangles. In Fig. 6(left) we also highlight the Bragg rows each side of two pairs of dislocations so that one can see the removal (insertion) of a halfline of particles characteristic of a dislocation (antidislocation). In Fig. 6(right) we highlight the Bragg rows surrounding two dislocations oriented almost tangentially to the boundary (and thus with Burgers vector almost parallel to the boundary).
At , pleats appear attached to the boundary, as shown in Fig. 7(left). Pleats are neutral grain boundaries with a 5fold disclination at one end and a 7fold disclination at the other end. A pleat can freely terminate at one or both ends and so they are easy to recognize.
Below pleats and dislocations proliferate both in the interior and at the boundaries (see Fig. 7(right)).
At one sees for the first time isolated 7fold disclinations (7s). In the figures these are seen clearly in the form of their dual heptagons. This is the true hallmark of the prevailing negative Gaussian curvature in the interior of the surface! 7fold disclinations may arise by the unbinding of a 7 from a pleat, which typically end with a 7 oriented towards the central interior of the surface. Indeed there are usually compensating scars with a net positive disinclination charge near an isolated 7 in this regime of aspect ratios. As the negative Gaussian curvature increases, however, isolated 7s are typically farther away from the positive scars that compensate them topologically.
In Fig. 8(left) we show a minimal energy configuration, with aspect ratio , containing an isolated 7fold disclination. We also show the geodesic circle, , with radius , that is the boundary of a disk whose integrated Gaussian curvature is . was determined by constructing a geodesic polygon with enough edges that the integrated curvature coincides with the total exterior angle deficit and applying the GaussBonnet theorem
(7) 
where is the Gaussian curvature and is the exterior angle deficit at vertex . Since disclinations naturally couple to Gaussian curvature in the effective Hamiltonian that controls their energetics Bowick et al. (2000); Bowick and Giomi (2009); Bowick and Yao (2011); Irvine et al. (2010) a natural criterion for the ground state to admit a defect with a net disclination charge is that there be a domain of the surface, centered on the defect, with integrated Gaussian curvature matching the disclination charge .
In Fig. 8(right) we show a configuration, with aspect ratio , containing a completely isolated 7 along with the geodesic circle centered on the 7 and enclosing an integrated Gaussian curvature of . We have also inscribed a 7sided geodesic polygon as determined by the triangulation. The exterior angles grow from at the 7fold disclination to at the inscribed geodesic heptagon. This nodoid has vanishing contact angle at the boundary, corresponding to complete wetting of the capillary bridge at the plates; it corresponds to the nodoid with parameter covering all . For smaller aspect ratio the contact angle grows to a maximum of .
As the capillary bridge is stretched to have a very narrow waist (the equatorial section of the Delaunay surface) with large negative maximal Gaussian curvature we observe a proliferation of 7fold disclination defects (heptagons). Two examples are displayed in Figs. 9(left, right). In both configurations there are a total of 11 isolated 7s plus scars, which have the same net disinclination charge. The 7s preferentially occupy the waist of the capillary bridge, whereas scars and pleats are found throughout the surface.
The sequence of defect motifs revealed by our simulations conforms remarkably closely to that found experimentally in Irvine et al. (2010) (see in particular their Fig. 4). The initial compressed bridge is free of defects. As the bridge is stretched, both experimentally and numerically, one observes the appearance of dislocations, neutral disinclination dipoles, polarized with the 7fold disclinations towards the maximally negatively curved neck (see Fig. 6 and Fig.4(i) of Irvine et al. (2010)). Further stretching leads to the appearance of pleats (see Fig. 7 and Fig.4(j) of Irvine et al. (2010)). Finally yet more stretching leads to the appearance of isolated 7fold disclinations and scars (see Fig. 8 and Fig.4(k) of Irvine et al. (2010)). These comparisons were also noted in Kusumaatmaja and Wales (2013) but we note that all the stretched manifolds displaying defects are nodoids (as clearly established in Fig. 5) whereas the simulations in Kusumaatmaja and Wales (2013) are only for catenoids and unduloids.
iv.2 Stretching from an initial square cylindrical capillary bridge
In this section we discuss the appearance of defects upon stretching through another family of Delaunay surfaces preserving volume and contact radius.
This family is generated by starting from a positive Gaussian curvature nodoid with aspect ratio and is designed to pass through the square cylindrical capillary bridge with and (aspect ratio ). It terminates in an unduloid with aspect ratio . Representative configurations for this sequence are shown in Fig. 10. The central cylinder exhibits very few defects. Increasing or decreasing the aspect ratio leads first to the appearance of pleats and then to longer pleats and scars. Geodesic discs with total curvature are shown surrounding scars and are completely contained in the interior of the surface. Note that the defect structure in Fig. 10(a) is a pleat with the top 7fold disinclination separated by one lattice spacing from the rest of the pleat.
V Minimum energy configurations on slices
For an isolated 7fold disclination (heptagon) to be present in minimum energy configurations of crystalline arrays on a negative curvature surface a reasonable criterion is that the integrated Gaussian curvature over the geodesic disc centered at the disclination be more negative than , the deficit angle corresponding to the topological charge of a 7fold disclination Bowick and Giomi (2009). To check that it is the integral of the Gaussian curvature, instead of the Gaussian curvature itself, that determines the character of defect motifs we examine minimum energy configurations on slices of a given Delaunay surface with steadily increasing values of , the maximal value of the meridian coordinate . Each slice has spanning the interval . The shape of a nodoid is completely determined by the parameters and . We merely restrict to a growing set of slices by varying . The first surface we investigate is the nodoid with and . This capillary bridge has maximal Gaussian curvature at the waist (equator) of . In Fig. 11(a) we take a slice with and the behavior is completely equivalent to a cylinder. For completeness we note that this surface encloses a volume , has a lateral area and an equilibrium potential energy .
In Fig. 11(b) we show a slice with . For this slice, dislocations are found in the interior just as in the case of stretching within the fat cylinder family. The evolution of defect motifs continues with the appearance of pleats and, in Fig. 11(c), where , disclination charge scars sitting inside discs of integrated Gaussian curvature . The first appearance of an isolated 7fold disclination (heptagon) is for a slice with . In Fig. 11(d) we show the geodesic disc of geodesic radius , the same value as in the previous slice, because the surfaces are the same except for their height. This last surface encloses a volume of and has a lateral area .
As we increase the parameter , the Gaussian curvature in the waist is strongly reduced and yet the sequence of defect motifs is essentially the same. To be concrete we show the family of slices with and . Now the Gaussian curvature in the waist is .
In Fig. 12(a) we display a slice with and the pattern is again equivalent to a cylinder. Its volume is and the lateral area is . In Fig. 12(b) we again observe the formation of dislocations in the interior. In Fig. 12(c) we have enclosed a dislocation with a polygon to show the associated Burgers’ vector. In Fig. 12(d) we see an isolated 7fold disclination with a nearby dislocation.
Finally in Fig. 13 we include the evolution of defects for the nodoid with parameters and , when increases and the slices contain most of the nodoid. In both cases there are 11 total defects counting isolated 7s and scars. In Fig. 13 (right) the isolated 7s are more concentrated in the waist.
This we find the full range of phenomena explored experimentally in Irvine et al. (2010) by simulating the full range of Delaunay surfaces generated in stretching. We emphasize that catenoids and unduloids alone are insufficent in making a proper comparison with experiment Kusumaatmaja and Wales (2013) since the majority of surfaces encountered in stretching are nodoids with negative mean curvature as opposed to the positive mean curvature of unduloids.
The presence of isolated 7fold disclinations (or equivalently their dual heptagons) or scars in the ground state offers new possibilities for supramolecular chemistry via functionalization of crystalline arrays on these surfaces, provided one can chemically detect the 7fold disclination (or scars) and attach ligands there just as smectic disclinations are functionalized by placeexchange reactions for nanoparticles coated with an equal mixture of short and long chain alkanethiols in the work of DeVries et al. DeVries et al. (2007). Isolated 7fold discinations, scars and pleats may also be sources of material weakness or perhaps improved performance through capturing material dislocations and thus limiting plastic deformation.
Acknowledgements
This research work was partly supported by the Spanish Research Council under project MTM201019660, by the National Science Foundation through award DMR0808812 and by the Soft Matter Program at Syracuse University. MJB would like to thank David Nelson and William Irvine for many stimulating discussions on the nature of condensed matter order on capillary bridges over a period of several years.
References
 Bowick and Giomi (2009) M. J. Bowick and L. Giomi, Adv. Phys. 58, 449 (2009), eprint http://arXiv.org/abs/0812.3064.
 Bowick et al. (2000) M. J. Bowick, D. R. Nelson, and A. Travesset, Phys. Rev. B 62, 8738 (2000), eprint http://arXiv.org/abs/condmat/9911379.
 Bowick et al. (2002) M. Bowick, A. Cacciuto, D. R. Nelson, and A. Travesset, Phys. Rev. Lett. 89, 185502 (2002), eprint http://arXiv.org/abs/condmat/0206144.
 Bausch et al. (2003) A. Bausch, M. Bowick, A. Cacciuto, A. Dinsmore, M. Hsu, D. Nelson, M. Nikolaides, A. Travesset, and D. Weitz, Science 299, 1716 (2003), eprint http://arXiv.org/abs/condmat/0303289.
 Bowick et al. (2007) M. Bowick, H. Shin, and A. Travesset, Phys. Rev. E 75, 021404 (2007), eprint http://arXiv.org/abs/condmat/0610819.
 Giomi and Bowick (2008a) L. Giomi and M. J. Bowick, Phys. Rev. E 78, 010601 (2008a), eprint http://arXiv.org/abs/0801.3484.
 Giomi and Bowick (2008b) L. Giomi and M. J. Bowick, Eur. Phys. J. E 27, 275 (2008b), eprint http://arXiv.org/abs/0807.4538.
 Pairam and FernandezNieves (2009) E. Pairam and A. FernandezNieves, Phys. Rev. Lett. 102, 234501 (2009).
 Yao and Bowick (2011) Z. Yao and M. J. Bowick, Eur. Phys. J. E 34, 1 (2011), eprint http://arXiv.org/abs/1011.3437.
 Giomi and Bowick (2007) L. Giomi and M. Bowick, Phys. Rev. B 76, 054106 (2007), eprint http://arXiv.org/abs/condmat/0702471.
 Bowick et al. (2008) M. J. Bowick, L. Giomi, H. Shin, and C. K. Thomas, Phys. Rev. E 77, 021602 (2008), eprint http://arXiv.org/abs/0709.2731.
 Sadoc and Charvolin (1989) J.F. Sadoc and J. Charvolin, Acta Crystallogr. Sec. A 45, 10 (1989).
 Vitelli and Nelson (2004) V. Vitelli and D. R. Nelson, Phys. Rev. E 70, 051105 (2004), eprint http://arXiv.org/abs/condmat/0406328.
 Vitelli et al. (2006) V. Vitelli, J. B. Lucks, and D. R. Nelson, Proc. Natl. Acad. Sci. 103, 12323 (2006), eprint http://arXiv.org/abs/condmat/0604203.
 Irvine et al. (2010) W. T. M. Irvine, V. Vitelli, and P. M. Chaikin, Nature 468, 947 (2010).
 Bowick and Yao (2011) M. J. Bowick and Z. Yao, Europhys. Lett. 93, 36001 (2011), eprint http://arXiv.org/abs/1012.4378.
 Irvine and Vitelli (2012) W. T. M. Irvine and V. Vitelli, Soft Matter 8, 10123 (2012), eprint http://arXiv.org/abs/1206.4159.
 Binks and Horozov (2006) B. P. Binks and T. S. Horozov, Colloidal particles at liquid interfaces (Cambridge University Press, 2006).
 Pickering (1907) S. Pickering, J. Chem. Soc. 91, 2001 (1907).
 Dinsmore et al. (2002) A. Dinsmore, M. F. Hsu, M. Nikolaides, M. Marquez, A. Bausch, and D. Weitz, Science 298, 1006 (2002).
 Leunissen et al. (2007a) M. E. Leunissen, A. van Blaaderen, A. D. Hollingsworth, M. T. Sullivan, and P. M. Chaikin, Proc. Natl. Acad. Sci. USA 104, 2585 (2007a).
 Leunissen et al. (2007b) M. E. Leunissen, J. Zwanikken, R. van Roij, P. M. Chaikin, and A. van Blaaderen, Phys. Chem. Chem. Phys. 9, 6405 (2007b).
 Irvine et al. (2012) W. T. M. Irvine, M. J. Bowick, and P. M. Chaikin, Nature Mater. 11, 948 (2012).
 Bowick et al. (2006) M. J. Bowick, A. Cacciuto, D. R. Nelson, and A. Travesset, Phys. Rev. B 73, 024115 (2006), eprint http://arXiv.org/abs/condmat/0509777.
 Eells (1987) J. Eells, Math. Intell. 9(1), 53 (1987).
 Delaunay (1841) C. Delaunay, J. Math. Pures et Appl. Sér. 1, 6, 309 (1841).
 Sturm (1841) M. Sturm, J. Math. Pures et Appl. Sér. 1, 6, 315 (1841).
 Nelson (1983) D. R. Nelson, Phys. Rev. B 28, 5515 (1983).
 Rubinstein and Nelson (1983) M. Rubinstein and D. R. Nelson, Phys. Rev. B 28, 6377 (1983).
 Tarjus et al. (2012) G. Tarjus, F. Sausset, and P. Viot, Adv. Chem. Phys. 148, 251 (2012), eprint http://arXiv.org/abs/1005.2684.
 Sausset et al. (2010) F. Sausset, G. Tarjus, and D. R. Nelson, Phys. Rev. E 81, 031504 (2010), eprint http://arXiv.org/abs/0911.0387.
 Kusumaatmaja and Wales (2013) H. Kusumaatmaja and D. J. Wales, Phys. Rev. Lett. 110, 165502 (2013), eprint http://arXiv.org/abs/1303.6786.
 Henderson and Taimina (2001) D. W. Henderson and D. Taimina, Math. Intell. 23(2), 17 (2001).
 (34) E. Bendito, M. Bowick, and A. Medina, “Delaunay surfaces” (in preparation).
 Bendito et al. (2007) E. Bendito, A. Carmona, A. M. Encinas, and J. M. Gesto, J. Comp. Phys. 225, 2354 (2007).
 Bendito et al. (2009) E. Bendito, A. Carmona, A. Encinas, J. Gesto, A. Gómez, C. Mouriño, and M. Sánchez, J. Comp. Phys. 228, 3288 (2009).
 DeVries et al. (2007) G. A. DeVries, M. Brunnbauer, Y. Hu, A. M. Jackson, B. Long, B. T. Neltner, O. Uzun, B. H. Wunsch, and F. Stellacci, Science 315, 358 (2007).