Emergent Geometry of Inhomogeneous Planar Crystals

Uniform triangular crystals are the ground state for particles that interact isotropically in two dimensions. However, when immersed in an external potential, for example, one arising from an electric field, a flow field, or gravity, the resulting phases are significantly distorted in a way reminiscent of conformal transformations of planar lattices. We study these “conformal crystals” using colloidal experiments and molecular dynamics simulations. By establishing a projection from these self-assembled inhomogeneous crystals to homogeneous crystals on curved surfaces, we are able to both predict the distribution of defects and establish that defects are an almost inevitable part of the ground state. We determine how the inherent geometry emerges from an interplay between the confining potential and the interparticle interactions. Using molecular dynamics simulations, we demonstrate the generic behavior of this emergent geometry and the resulting defect structures throughout a variety of physical systems.


I. INTRODUCTION
The work of M. C. Escher is often built on curiously distorted grids [1] [ Fig. 1(a)]. Similarly, repulsive particles confined by an external potential form distorted latticelike structures [ Fig. 1(c)] with large variations in size and orientation of identifiable local unit cells. This type of assembly occurs naturally in a diverse range of twodimensional systems, including electrons on the surface of liquid helium, quantum dots, one-component plasmas, and dusty plasmas [2][3][4][5][6][7][8][9][10], and has also found applicability in type-II superconductors [11][12][13]. Despite their ubiquity, understanding the principles that govern these structures has remained challenging because the large variations in scale and orientation are incompatible with typical treatments of condensed matter that rely on uniformity of the condensed phase.
Early theoretical studies of these structures leveraged the same mathematical transformations found in the works of Escher: conformal maps on the plane. These maps have the property of preserving angles so that perpendicular grid lines map to perpendicular grid lines. Applying these transformations to uniform triangular lattices produces "conformal crystals" [ Fig. 1(b)] that closely resemble the inhomogeneous self-assembled structures formed in a variety of systems including repulsive particles in a potential, foams, and phyllotactic patterns found in plants. This resemblance has stimulated the use of conformal maps to model such systems [14][15][16].
However, since conformal maps preserve angles, the resulting distorted crystals are everywhere locally sixfold coordinated. In contrast, inhomogeneous crystals found in self-assembled physical systems typically contain topological defects [3,4,[17][18][19][20][21]. Moreover, theoretical and computational studies of electrons confined in rotationally symmetric potentials show that topological defects are directly related to the generation of density variations in the resulting structures [3,4,18], raising the question of whether such systems can ever self-assemble into a perfect conformal crystal. Furthermore, while defects in ordinary crystals have proven crucial in understanding the material's order, mechanics, and phase transitions [22,23], the physics of topological defects in inhomogeneous systems remain far less understood.
Recent computational work on Coulomb interacting particles suggests that the total number of defects in inhomogeneous crystals can be related to that of crystals frustrated by Gaussian curvature [18]. In this work, we generalize this concept to show that relating nonuniformity to Gaussian curvature can be used to understand how defects are distributed throughout a variety of inhomogeneous systems. Using colloidal experiments and molecular dynamics simulations, we confine repulsive particle monolayers in spatially varying potentials. We find that the particles generically form defected, inhomogeneous crystals. Instead of starting from a perfect crystal in a plane, and distorting it using a conformal map, we gain insight by carrying out the inverse process: We conformally transform our nonuniform experimental structures into lattices of uniform density. We find that, unlike the transformations demonstrated in Figs. 1(a) and 1(b), which map one planar structure to another, making our experimental structures uniform requires that we map them onto a curved surface, where the physics of defects have been studied extensively [22,[24][25][26][27][28][29][30][31].
Taking advantage of the inherent connection between inhomogeneity, defects, and curvature has recently fostered unexpected insights into seemingly unrelated systems, such as twisted bundles of filaments and nonuniformly growing elastic sheets [32][33][34][35][36]. Similarly, we find that by viewing our flat, nonuniform structures as projections of uniform crystals on curved surfaces, we are able to predict the resulting defect distributions and gain insight into their geometric function and generic behavior.

II. EXPERIMENTAL SYSTEM
Our experimental model system consists of a monolayer of repulsive poly(methyl methacrylate) (PMMA) particles (∼1 μm diameter) that are bound to an interface between an oil (dodecane) phase and a water-glycerol mixture. We estimate that the total integrated curvature of the region the particles are confined in is on the order of 10 −3 and, thus, any interactions with the substrate curvature are negligible. An oscillating electric potential is applied to a sharp, needle-shaped electrode and placed above the interface [see Fig. 1(d)]. The resulting electric field gradient creates an approximately parabolic potential through dielectrophoresis (see Appendix A). By varying the applied potential, we can easily tune the confinement strength of the system. In this confinement, as shown in Fig. 1(c), the equilibrium particle distribution is inhomogeneous but appears to locally maintain a triangular crystal structure.

III. GEOMETRIC FRUSTRATION
In the absence of a spatially varying confining potential, the particles in our experiment will form a uniform triangular (Wigner) crystal everywhere (Fig. 8). If confined by a wall, the lattice spacing will be determined by the number of particles and the area of the box or, equivalently, the pressure supported by the edges of the box. Once the reference state is established, deviations from the Gallery. To the right of the artwork is the conformal grid, which bears striking similarity to the distorted grid used by Escher to create his piece (image from Ref. [1]). All M. C. Escher works © 2018 The M. C. Escher Company-the Netherlands. All rights reserved. Used by permission. www.mcescher.com All images may only be used unaltered, unchanged and not manipulated in any way. (b) Perfectly conformal lattice used in Ref. [14] to model repulsive particles in a gravitational field. (c) Inhomogeneous crystal of repulsive colloidal particles. (d) Experimental setup. Micron-sized PMMA colloids are confined at an oil-water interface. We estimate the total integrated curvature in the region where the particles are confined and, thus, any interactions with the substrate curvature are negligible. An external potential is created by applying an alternating voltage to a needle-shaped electrode above the sample. The voltage is grounded on a transparent indium-tin-oxide cover slip below the sample. The particles are fluorescently labeled and imaged from below. crystalline reference state, resulting in compression, shear, or topological defects, can be treated with elasticity theory [37,38]. A spatially varying confining force produced by an external potential will give rise to a spatially varying internal pressure and thus reference lattice spacing. The body may therefore, in principle, possess a nonuniform lattice spacing while maintaining vanishing elastic energy.
Analysis of our crystals, however, reveals many defect patterns, including "scars"-chains of alternating fivefold (positive topological charge) and sevenfold (negative topological charge) disclinations terminating in the bulk of the crystal [ Fig. 2(a)]. This is indicative of an incompatibility between our lattice spacing variation and the triangular structure. Since it is possible to have a defect-free structure in which locally triangular patches of lattice, each with their own lattice spacing, are stitched together seamlessly [see Fig. 1(b)], this incompatibility is not an inevitable consequence of the fact that there is a large variation in lattice spacing, but rather a consequence of how it varies globally. This is analogous to the frustration encountered when attempting to tile a curved surface with a uniform triangular crystal. In this case too, the frustration is not an inevitable consequence of curvature-a cylinder can be seamlessly tiled-but rather a consequence of how it is curved. If Gaussian curvature is present, frustration is inevitable. This is because placing a uniform crystal on a surface with Gaussian curvature will cause lattice lines, which would otherwise be parallel, to converge or diverge, distorting the structure from the crystalline reference state. The presence of topological defects allows these rows to terminate, enabling a uniform density crystal to conform to the surface [22,[24][25][26][27]29].

IV. CONFORMAL GEOMETRY, CURVATURE, AND DEFECTS
In inhomogeneous crystals, the lattice lines also converge and diverge. However, certain density variations can FIG. 2. Conformal geometry and defects. (a) After collecting microscope images of the equilibrated system, each particle is identified and the resulting structure is triangulated. Here, fivefold and sevenfold defects are colored green and red, respectively, while black particles are undefected. The shaded gray area represents the boundary region, which is not used for the analysis. In this region, the net defect charge drops to þ6, as demanded by topological constraints. (b) Conformal transformation between inhomogeneous and homogeneous structures. Above the flat, experimental structure is the corresponding transformed, curved crystal. (c) Average angle between bonds vs orientation of bisect. The angle between adjacent bonds is plotted against the orientation of the bisect of those two bonds with respect to the radial direction. We find that the angle of bonds is the value expected for a triangular lattice, irrespective of the orientation. In other words, adjacent bonds that are oriented in the azimuthal direction and those oriented in the radial direction both are separated by 60 degrees. (d) Average bond length vs orientation of bond. The length of bonds, normalized by the average local lattice spacing at that distance from the structure's center, is plotted against the orientation with respect to the radial direction. We find that bond lengths are uniform, irrespective of their orientation. cause the lattice lines to converge or diverge in ways that do not result in frustration and thus do not require defects. How then can the presence of frustration be visualized? We first note that despite the density variation, the local structure in our experiment is still conformal to a triangular crystal everywhere. This can be seen in Figs. 2(c) and 2(d), which show that the average angle between bonds and the local length of bonds are uniform, irrespective of their orientations. Therefore, there exists a conformal map that can transform our structure into a more crystalline one while preserving the local geometry. We use this map to create a clearer picture of the angular frustration intrinsic to the system.
Conformal transformations locally scale area elements by a factor Λ 2 ðr; ϕÞ. Two sufficiently close points which are initially spaced by ds are separated approximately by Λds after the transformation. Here, we want the spacing between particles to be uniform throughout the system. By measuring the experimental lattice spacing variation aðrÞ, we can thus deduce ΛðrÞ as Λ ¼ fa 0 /½aðrÞg, where a 0 is an arbitrary target uniform lattice spacing.
The metric of a sufficiently smooth surface can also be expressed as Λ 2 ds 2 , where ds 2 ¼ ðdr 2 þ r 2 dϕ 2 Þ corresponds to a flat metric [39]. When written in this form, the Gaussian curvature K of the surface can be computed in terms of Λ: Note that the spatial variations in lattice spacing do not automatically produce an underlying nonzero Gaussian curvature. Variations that satisfy a Laplacian Δln½ΛðrÞ¼0, such as those found in Figs. 1(a) and 1(b), are intrinsically flat.
By determining the appropriate scaling factor Λ [ Fig. 2(b)], one can establish the presence or absence of an underlying Gaussian curvature of the transformed, uniform crystal. We find that our experimental structures are always curved, and Fig. 2(c) shows the curved surface corresponding to a typical experimental structure. Note that there is an additional freedom in choosing the embedding of our curved surface; however, only the Gaussian curvature itself, and not the particular embedding, is used in our analysis.
Thus, such self-assembled inhomogeneous crystals are naturally modeled, not by maps from the plane to the plane, but from a curved surface to the plane, like stereographic projections used to make maps of the Earth. This makes curved-space crystallography the natural starting point for modeling the defect structures that are observed in selfassembled inhomogeneous crystals and provides a natural explanation for the presence of defect structures, such as scars and pleats, which enable the relief of angular strain induced by Gaussian curvature.
Equation (1) is identical to a quantity found in Ref. [40] known as the amount of "extra matter" in an elastic body. This quantity is derived from the nonmetricity tensor, the covariant derivative of the relaxed state's metric. Equation (1) is also identical to the relationship for the disclination charge in confined one-component plasmas derived in Ref. [4]. This relationship is derived from a counting argument for the Burgers vector density resulting from the lattice spacing variation. Using our mapping to a curved surface, we can understand the relationship between the disclination density and the frustration caused by lattice variation using known results from curved-space crystallography. On a curved surface, the sum of the internal angles in a geodesic triangle will differ from the flat space value by the amount of integrated Gaussian curvature in the interior of the triangle. This angular mismatch is accounted for by the presence of disclinations in the crystal. In particular, within a patch of smoothly curved surface, the integrated Gaussian curvature will be approximately proportional to the net defect charge [29].
To demonstrate this relationship for our inhomogeneous crystals, we integrate the effective Gaussian curvature within increasing radii from the center of the structure and compare the result to the net topological charge within those radii. Figures 2(f) and (9) show the relationship between the effective Gaussian curvature extracted from the experiment and the topological defect charge for a range of external potential strengths. We find the quantity shown in Eq. (1) effectively predicts the distribution of topological defects induced by density inhomogeneities in our experiment.
We note that the overall number of excess defects in the system will be fixed by topological constraints. In particular, our planar experimental structures must have an excess topological charge of þ6. We find that, while the integrated topological charge within the bulk of our system can exceed this amount by more than a factor of 2, the total charge meets the topological constraint because of disclinations located near the boundary.
The distribution of defects is predicted purely from geometric arguments. Thus, the analysis we have presented can easily be applied to a wide range of physical systems. To demonstrate the generality of these results, we perform molecular dynamics simulations [41-43] of particles interacting through a variety of interaction potentials within external confining potentials (see Appendix D for simulation details). In addition to systems with rotational symmetry, we examine potentials that contain translational symmetry, as shown in Fig. 3. Our results indicate that the effective integrated curvature robustly predicts the defect distribution across all systems we have investigated.
How is the effective Gaussian curvature determined by the physical parameters of the system? By tuning interaction and confinement potentials, we observe that the resulting structures vary in overall size by factors of 3 or more. Despite this, we find that the curvature distributions, when rescaling coordinates by the respective system sizes, collapse onto each other (Fig. 3).
To examine the generality of this scaling behavior, we relate the physical system parameters to the resulting curvature profiles. We first determine the lattice spacing variation by fixing the total chemical potential in the system: μðrÞ þ V ext ¼ E, where μ is the local chemical potential derived from the interaction potential, V ext is the external potential, and E is a constant. In the case of interparticle interactions of the form V int ∝ 1/r (monopoles), V int ∝ 1/r 3 (dipoles), and V int ∝ 1/r m (m > 3) confined in potentials of the form V ext ∝ r n (n > 0), we find that the resulting lattice spacing variations have the form aðrÞ ¼ a 0 /½1 − ðr/RÞ p 1/q , where the powers p and q depend only on m and n (see Appendix F). In this case, the integrated curvature is Likewise, for translationally symmetric confinements of the form V ext ∝ jyj n (n ≥ 2), for the above interactions, we obtain a similar form for the lattice spacing variation, aðyÞ ¼ a 0 /½1 − ðy/HÞ p 1/q , and the integrated curvature is FIG. 3. Structures of simulated confined crystals. (a) Rotationally symmetric potentials. An example low-energy configuration of Coulomb-repulsive particles confined inside the potential V ext ∝ r 2 is shown on the left next to the corresponding uniform curved crystal. On the right of the crystal is a plot of the integrated curvature for particles in rotationally symmetric parabolic potentials interacting via various inverse power-law potentials and the result after scaling by the overall system size. (b) Results for translationally symmetric potentials. We show an example low-energy configuration of 1/r 4 interacting particles inside of the potential V ext ∝ y 2 next to the corresponding uniform curved structures. The simulated system is carried out with periodic boundary conditions. Note that the entire crystal is not shown. Shown on the right is the integrated curvature profile for various inverse power-law interactions and the result after scaling by the overall system size. In this case, the integrated curvature scales linearly by the lateral system size, and thus, the integrated curvature must be scaled by the crystal aspect ratio. (c) Result for particles in a constant force field (V ext ∝ y) in the presence of a wall at the bottom o the crystal. Again, the simulated system is carried out with periodic boundary conditions, and the entire crystal is not shown.
We also examine the geometry found in Refs. [14,16]. Here, particles are placed in a constant force field (V ext ∝ y) towards a boundary wall [ Fig. 3(c)]. We find that while longer-ranged interacting particles (m ≥ 3) once again obey Eq. (3), Coulomb interacting particles do not (see Appendixes F and G).
For all systems obeying integrated curvatures of the forms in Eqs. (2) and (3), rescaling coordinates by the crystal size reveals that the curvature distributions collapse onto single curves that are determined solely by the powers of the internal and external potentials. In particular, all prefactors of the interaction and external potentials are absorbed into the system sizes R and H. This result suggests that the relative distribution of curvature in these systems is robust against physical parameters. In particular, changing the relative strengths of the external and interaction potentials modifies the crystal size but has little effect on the relative distribution of topological defects. We also explore the effect of modifying the interaction potential through the explicit addition of a length scale using a Yukawa potential. We observe that while this length scale may disrupt the above scaling, the relation between curvature and defect charge still holds (see Appendix H).
The above results reveal the extent to which the distribution of curvature and, thus, defects in confined systems may be manipulated. Examining Eq. (1) also reveals that the curvature can be tuned to change sign. A change in sign occurs when the derivative of the lattice spacing changes sign. Qualitatively, this will change the direction of polarization of dislocations and will cause the system to favor a net of sevenfold defects instead of fivefold defects (or vice versa). This can be seen, for example, by comparing Coulomb repulsive particles confined by hard walls (shown in Figs. 4 and 5) [4,18,20,21] with the structures shown thus far. In this case, the long-range interactions themselves create FIG. 4. Simulated hard-wall confined crystal with a corresponding negatively curved surface. (a) Circular wall boundary. An example low-energy configuration of Coulomb-repulsive particles confined inside a circular hard-wall boundary is shown next to the corresponding uniform curved crystal. Embedding the entire structure requires covering the surface multiple times (see Appendix E). In order to separate the different coverings, the surface is slightly distorted. (b) Integrated curvature, which collapses when rescaled by the system radius.
FIG. 5. Simulated hard-wall confined crystal with corresponding negatively curved surface. (a) Translationally symmetric hard-wall boundary. An example low-energy configuration of Coulomb-repulsive particles confined inside a hard-wall boundary is shown next to the corresponding uniform curved crystal. Embedding the entire structure requires covering the surface multiple times (see Appendix E). In order to separate the different coverings, the surface is slightly distorted. (b) Integrated curvature, which collapses when rescaled by the system radius. density inhomogeneities in the resulting Wigner crystal [4]. Again, we find that Eqs. (2) and (3) predict the effective integrated curvatures and, thus, the defect distributions. Figures 4 and 5 show the resulting structures and negatively curved surfaces.
In order to stabilize a defect-free structure, our results suggest that the integrated curvature must vanish everywhere in the domain. Equation (1) thus suggests that defects should be expected for all lattice spacing variations that do not have the form of a power law. In all the structures we have observed, the lattice spacings diverge at finite length scales and thus do not take on a power-law form. This suggests that, for potentials we have explored, topological defects should be present for all interaction and confinement strengths.

V. FINITE TEMPERATURE
The simulated results discussed above have informed our understanding of the minimum energy defect configurations of inhomogeneous crystals. In real assembly processes, however, the effects of temperature can also influence the excitation of defects [44]. This is the case, for example, in our experimental system.
In particular, we observe that as we reduce the confinement strength, the experimental structures are characterized by a growing disordered region at the edge, as shown in Fig. 6(f). To characterize this disorder, we use the magnitude of the orientational order parameter, jΨ 6 ðr m Þj ¼ ð1/N b Þj P N b n¼1 e 6iθ nm j, where N b is the number of bonds of a particle m located at r m , and θ nm is the angle of the bond joining particles m and n. Figure 6(a) shows that the disordered region grows in from the edge as the confinement strength is reduced. We observe that this effect is absent in zero-temperature simulated structures, irrespective of tuning.
What consequences will the observed inhomogeneous melting have on the distribution of curvature and, thus, defects? Surprisingly, we observe that the relationship between net defect charge and curvature persists, as shown in Fig. 2(f). In this figure, the partially melted crystals correspond to the brighter data points. This suggests that the notion of effective curvature may be applicable even in the melted or partially melted state.
While the balancing of net defect charge (fivefold minus sevenfold) and curvature is robust against the effects of temperature, partial melting does have consequences for the overall number of defects (fivefold plus sevenfold) in the system. In particular, in Fig. 6(b), we show the minimum number of defects required to balance the curvature and the measured number of defects. The minimum number of defects is calculated by relating the observed lattice variation to the dislocation density [4] (see Appendix C for more details). We find that the total number of defects grows away from the minimum required number as the confinement strength is reduced, a further characteristic of melting.

VI. CONCLUSION
We have used the inherent connection between topological defects and curvature to reveal generic behavior in confined two-dimensional crystalline systems. In particular, we show that inhomogeneity in density can frustrate a crystal in a manner similar to out-of-plane Gaussian curvature. For structures self-assembled in external potentials like the ones explored here, we find this frustration to be generic, suggesting they always have defects. Thus, our results suggest that it may not be possible for perfect conformal crystals to self-assemble. An interesting extension of these concepts may be the use of defects to engineer density variations, much like their use in engineering curvature [45][46][47]. Fellowship. L. R. G. acknowledges support from the National Research Council of Argentina, CONICET, ANPCyT, and Universidad Nacional del Sur.

APPENDIX A: EXPERIMENTAL DETAILS
Colloidal experiments are conducted in glass capillary channels (height ≈110 μm, width ≈3 mm, length ≈22 mm) as follows. First, water-glycerol droplets (50/50 wt %) are prepared in channels while in contact with air. All glassware is first base washed in order to remove any coatings and reduce the droplet contact angles. Next, fluorescent PMMA particles (diameter ≈1.5 μm), prepared using the methods of Refs. [48,49] and suspended in dodecane, are sent into the channel, leaving a monolayer of particles confined to the fluidfluid interfaces. We estimate that the total integrated curvature of the region the particles are confined in is on the order of 10 −3 and, thus, any interactions with the substrate curvature are negligible. The particles are viewed from below in fluorescence with an inverted microscope.
Because the particles interact isotropically, they naturally form a triangular lattice at high enough density in the absence of a spatially varying potential (Fig. 8).
Next, a needle-shaped electrode is placed on top of the capillary channel above a desired water-glycerol droplet. Alternating voltages between 1 and 5 kV are then applied at 100 Hz on the electrode using a high-voltage amplifier. A transparent indium-tin-oxide-coated cover slip is placed below the capillary and acts as the ground. After voltage is initially applied, the sample is left to equilibrate for at least a day before running experiments. After varying voltage during an experimental run, the sample is left to reequilibrate for 1-2 hours.
In order to approximately characterize the potential in which the particles are confined, the electric field profiles are numerically calculated. The calculation includes the glass boundary walls, electrode, ground, and each fluid phase with dielectric constants of each. The resulting external potential arises from a dielectrophoretic force F ∝ ∇jEðrÞj 2 . Thus, the external potential is given by V ext ∝ jEðrÞj 2 (Fig. 7).

APPENDIX B: EXPERIMENTAL ANALYSIS
After collecting microscope images of the equilibrated system, each particle is identified using standard routines [50] and the resulting structure is triangulated, allowing the determination of lattice spacings and particle coordinations. We note that all quantities are averaged over at least 20 images, taken over the course of several minutes after equilibration. In addition, lattice spacing measurements are azimuthally averaged.
Figure 2(f) shows the measured relationship between the integrated net disclination charge and integrated curvature. We note that we also directly compare the spatial variation of each quantity, as shown in Fig. 9.

APPENDIX C: MINIMUM DEFECT COUNT
To approximate the minimum number of defects, we calculate the dislocation density needed to support the observed lattice variation [51]. This is done by first computing the Burgers vector density [4]: Dividing the Burgers vector density by the distance between lattice lines (a ffiffi ffi 3 p /2), we obtain the density of dislocations. We then multiply this number by 2 and compare to the number of particles with a coordination number not equal to 6 to obtain Fig. 6(b).

APPENDIX D: SIMULATION DETAILS
Molecular dynamics simulations are performed using HOOMD-Blue [42,43,52] at constant volume. Particles were annealed from a disordered state by gradually lowering the system temperature until particles are no longer able to rearrange. Temperature ramping parameters vary between systems and need to be chosen independently between runs in order to ensure simulation stability and proper assembly. After this annealing procedure, the resulting configurations are relaxed using the Fast Inertial Relaxation Engine (FIRE) in order to obtain the final structures. Note that the result is thus a zerotemperature, approximately energy-minimized structure. Thus, other techniques, such as Monte Carlo methods, could also have been used.
In the case of rotationally symmetric external potentials, the simulation box size is chosen to be much larger than the resulting structures. In the case of translationally symmetric potentials, the lateral dimension of the simulation is treated with periodic boundary conditions, while the vertical direction is chosen to be much larger than the size of the resulting structures. Potentials that do not have their minimum in the interior of the simulation box need to include walls to prevent particles from escaping. In addition, walls are needed for cases in which no other confinement potential is present, as is the case for the negative effective curvature systems studied here. Walls are implemented with Lennard-Jones-type forces: where with parameters chosen such that the simulation remains stable and the wall does not directly interact with the bulk structure. In the case of the circular boundary, the circle was discretized with 24 straight walls.

APPENDIX E: SURFACE RECONSTRUCTION
We can write the metric associated with the conformal transformation as ds 2 ¼ Λ 2 ðrÞðdr 2 þ dθ 2 Þ. Interpreting this expression as the metric of a surface written in conformal (isothermal) coordinates, we can write the Gaussian curvature as KðrÞ ¼ f1/½Λ 2 ðrÞgΔ lnf1/½ΛðrÞg. FIG. 9. Experimental comparison of net disclination charge and integrated curvature distribution. The net disclination charge enclosed within varying radii is compared to the integrated curvature within that radii. The two quantities are found to be proportional at each radii. Here, each plot corresponds to a different-sized crystal, and each row of plots corresponds to a crystal with a different number of total particles, N.
Embeddability of the resulting surfaces in IR 3 is not automatically guaranteed. In the case of rotationally symmetric surfaces with negative Gaussian curvature everywhere, it can be shown that f½dρðr 0 Þ/ðdr 0 Þg 2 < 1 everywhere, and thus, the z coordinate is complex throughout the structure. In other cases, the resulting surface may only be partially embeddable. However, in these situations, the structure may be embedded by augmenting the procedure by transforming the angular variable θ → ϕ, where nθ ¼ ϕ, n > 1. The structure will thus cover the resulting surface n times. Figures 4 and 5 demonstrate this for the case of Coulomb interacting particles in a translationally symmetric hard-walled potential, in which the resulting surface has a catenoidlike geometry.
APPENDIX F: GROUND-STATE LATTICE SPACING

Power-law potentials
Suppose N repulsive particles are immersed in a rotationally symmetric external potential, V ext . We assume the system locally obeys the equilibrium condition: Here, nðrÞ is the local number density. The pressure P is computed from the internal energy U by P ¼ −dU/dV. Assuming local crystallinity, the lattice spacing is related to the number density as nðrÞ ¼ 1/V ¼ f2/½ ffiffi ffi 3 p a 2 ðrÞ. By relating the chemical potential μðrÞ to the pressure ndμ ¼ dP, we recast the condition in Eq. (F1) as where E is a constant throughout the system. Thus, given an interaction potential V int , one can compute the internal energy, and thus, the pressure and chemical potential μðrÞ in terms of the local density. Using Eq. (F2), the local density can then be determined as demonstrated below. Consider interaction potentials of the form V int ¼ ð1/r m Þ and external potentials of the form V ext ¼ γr n , where m > 0, γ > 0, and n > 1. If we assume that the neighborhood around each particle is crystalline, we can write down the pairwise interaction energy in terms of a multiple of the lattice spacing. The internal energy of a particle is then obtained from summing over the pairwise interaction between all the other particles in the system UðrÞ ¼ fk/½a m ðrÞg; where Note that this summation only converges for powers m > 2. Cases where m ≤ 2 must be treated separately (see Appendix F 3). Moreover, because of the varying lattice constant, the expression for the interaction constant shown in Eq. (F3), which assumes perfect crystal order, is only an approximation. The approximation is valid, however, when the lattice spacing change is small in the area in which the summation sufficiently approaches the convergent value.
In this case, Eq. (F2) takes the form Therefore, where R ¼ f½kðm þ 2Þ/ð2γa m 0 Þg ð1/nÞ is the resulting crystal radius. The initial lattice spacing a 0 is determined from the normalization condition Figure 10(b) shows the lattice spacing variation of a typical simulation compared to this calculation.

Translationally symmetric confinements
The above arguments can easily be applied to systems with translational symmetry. Consider, for example, potentials of the form V ext ¼ γjyj n and V int ¼ ð1/r m Þ, where m > 2, γ > 0, and n ≥ 1. The arguments in the previous section once again hold, resulting in the lattice spacing variation where H ¼f½kðmþ2Þ/ð2γa m 0 Þg ð1/nÞ . Figures 10(d) and 10(f) show example simulation lattice spacing variations compared to this calculation.

The Coulomb interaction
The approximation shown in Eq. (F3) is not valid in cases where m ≤ 2. In this case, the summation used to compute the interaction constant is not convergent. Therefore, in the case of Coulomb repulsion between particles, the local internal energy must take into account density variations over the entire system: In the case of rotationally symmetric parabolic external potentials, V ext ¼ γr 2 , the resulting lattice spacing variation can be shown as [2][3][4] aðrÞ where a 0 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ½ð2π 2 Þ/ð4 ffiffi ffi 3 p γRÞ q and R ¼ f½ð3πNÞ/ð8γÞg 1 3 are again given by the normalization condition Eq. (F6) [see Fig. 10(a)].
In the presence of only a hard-wall confinement, the distribution of particles also becomes inhomogeneous as a result of the long-range nature of the interactions. In this case, the resulting lattice spacing can be shown as [4] aðrÞ where a 0 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ½ð4πR 2 Þ/ð ffiffi ffi 3 p NÞ q is given by normalization Eq. (F6) [see Fig. 10(g)].
The distribution of particles can also be calculated in the case of 1D (translationally symmetric) confinements. For hard-wall confinements, the particle density has the form of the charge density on a conducting strip [53]: where y is the distance from the center of the strip in the finite dimension, H is the width of the strip, λ is the charge per unit length in the lateral dimension, and the lattice Fig. 10(h)]. In the case of 1D parabolic confinement, V ext ¼ γy 2 , the density variation has the form where ρ 0 ¼ ½ðγHÞ/π, H ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ½ð2λÞ/γ p , and λ is the charge per unit length in the lateral dimension [54,55]. See Fig. 10(c).
Finally, for linearly varying potentials in the presence of a hard wall, such as in the case of charged particles in a uniform electric field with a barrier or massive particles under the influence of gravity, the density takes the form of a charged conducting strip in a uniform external field: where H ¼ ½ð4NÞ/ðγΔxÞ is given by the condition that the density vanish at H [see Fig. 10(e)]. This result is derived from the solution for the potential near a charged conducting cylinder placed in a uniform electric field. Using the Joukowski transformation, the solution for the cylindrical case can be mapped to the case of a conducting strip of charge in a uniform field, giving Eq. (F13).  In cases where aðrÞ has the form shown in Eq. (F5), we have The crystal radius R in the above expression will, in general, vary between systems, depending on the particular system parameters. By rescaling by R, however, we see that the above form suggests that the integrated curvature profile depends only on the relevant powers of the interaction and external potentials, and not on the potential strengths γ and k. For Coulomb interacting systems, the lattice spacings are given by Eqs. (F9) and (F10). Using Eq. (G2) results in the same form of the integrated curvature as Eq. (G3), but with the parameters n and m redefined.

Translational symmetry
For translationally symmetric density variations, the integrated curvature takes the form FIG. 11. Yukawa interacting particles in a parabolic confining potential V ext ¼ γy 2 . (a) Integrated effective curvature predicts final defect distribution as in the case for pure power-law potentials. (b) Example particle configuration. (c, d) Scaling behavior of the integrated curvature profiles. With the additional length scale introduced by the screening length, our data suggest that the integrated curvature profiles no longer collapse when scaling the coordinates by the system size as shown in (d).
Here, Δx is the lateral size of the system. In the simulations presented here, periodic boundary conditions are used in the lateral dimension, and Δx is the distance between periodic boundaries. For the lattice spacing variation shown in Eq. (F7), the integrated curvature takes the form In this case, the integrated curvature scales with Δx. Thus, when the coordinates are rescaled by the system size, H, the integrated curvature depends only on the relevant powers of the interaction and external potentials, and not on the potential strengths γ and k.
Here, Coulomb interacting systems must again be treated separately. The lattice spacings are given by Eqs. (F11)-(F13). For all cases except Eq. (F13), using Eq. (G4) results in the same form of the integrated curvature as Eq. (G5), but with the parameters n and m redefined.
For the case of Coulomb particles in a uniform field, given by Eq. (F13), the integrated curvature takes the form In this case, the integrated curvature no longer collapses when rescaling coordinates by the system size. In particular, Eq. (G6) suggests that tuning each system parameter affects the resulting defect distribution. Another feature of this equation is a concentration of integrated curvature near the boundary, as is the case of hard-wall confined Coulomb systems.
FIG. 12. Yukawa interacting particles in a parabolic confining potential V ext ¼ γr 2 . (a) Example particle configuration. (b) Integrated effective curvature predicts final defect distribution as in the case for pure power-law potentials. (c)-(d) Scaling behavior of the integrated curvature profiles. With the additional length scale introduced by the screening length, our data suggest that the integrated curvature profiles no longer collapse when scaling the coordinates by the system size.