key: cord-0297566-xviaq6qr
authors: Bowick, Mark J.; Giomi, Luca
title: Two-Dimensional Matter: Order, Curvature and Defects
date: 2008-12-16
journal: nan
DOI: 10.1080/00018730903043166
sha: 8e72c04f5585dc7a4fbfb66c215a116087087b9d
doc_id: 297566
cord_uid: xviaq6qr

Many systems in nature and the synthetic world involve ordered arrangements of units on two-dimensional surfaces. We review here the fundamental role payed by both the topology of the underlying surface and its detailed curvature. Topology dictates certain broad features of the defect structure of the ground state but curvature-driven energetics controls the detailed structured of ordered phases. Among the surprises are the appearance in the ground state of structures that would normally be thermal excitations and thus prohibited at zero temperature. Examples include excess dislocations in the form of grain boundary scars for spherical crystals above a minimal system size, dislocation unbinding for toroidal hexatics, interstitial fractionalization in spherical crystals and the appearance of well-separated disclinations for toroidal crystals. Much of the analysis leads to universal predictions that do not depend on the details of the microscopic interactions that lead to order in the first place. These predictions are subject to test by the many experimental soft and hard matter systems that lead to curved ordered structures such as colloidal particles self-assembling on droplets of one liquid in a second liquid. The defects themselves may be functionalized to create ligands with directional bonding. Thus nano to meso scale superatoms may be designed with specific valency for use in building supermolecules and novel bulk materials. Parameters such as particle number, geometrical aspect ratios and anisotropy of elastic moduli permit the tuning of the precise architecture of the superatoms and associated supermolecules. Thus the field has tremendous potential from both a fundamental and materials science/supramolecular chemistry viewpoint.

More than 200 years ago, in his treatise on the resistance of fluids, d'Alembert wrote: "Geometry, which should always follow physics when used to describe nature, sometimes commands it" [1] . Since then the physics community has explored the power of geometry not only to describe, but also to explain structures and their properties. In the past 20 years soft condensed matter physics has provided many examples of how the geometry of matter is not a quiescent background for some microscopic degrees of freedom, but instead plays a major role in determining structural and mechanical properties and designing the phase diagram of materials such as colloids, liquid crystals, membranes, glasses and carbon nanostructure.

Geometric models of condensed matter systems have been developed for a wide class of materials since the pioneering work of Bernal and Finney [2, 3, 4, 5, 6] . In a series of classic papers they suggested that several properties of liquids have their geometric counterpart in randomly packed arrays of spheres. The difference in density between the solid and the liquid phase of a simple monoatomic substance, for instance, is approximately the same between periodically and randomly packed hard spheres (roughly 15% − 16%). Also the radial distribution function of randomly packed spheres corresponds well with that determined by X-ray and neutron diffraction for rare-gas liquids. After Bernal, a significant amount of work has been done on random close packing and, even though the legitimacy of the notion of random close packing itself has been questioned frequently in recent years [7] , it is now established that many features of the liquid state have in fact a purely geometrical nature.

After the discovery of icosahedral order in metallic glasses [8, 9, 10] the idea of geometrical frustration (the geometric impossibility of establishing a preferred local order everywhere in space, see Sec. 2), became a fundamental concept for the characterization of amorphous solids. Farges and coworkers [11, 12, 13] were the first to show, by electron diffraction experiments and computer simulations, that the first atoms of small aggregates of rare gases condensed in ultra-high vacum form regular tetrahedra, which later organize in the form of small icosahedral clusters.

Since icosahedra do not fill three-dimensional Euclidean space R 3 , the structure resulting from the aggregation of these icosahedral building blocks does not exhibit long range translational order. The lack of crystallization in covalent glasses is also rooted in the geometrical frustration associated with the constant coordination number of their constituents. Tetravalent monoatomic materials, for instance, cannot form a constant angle between bonds incident at the same atom and organize in a regular network at the same time. In multiatomic glass-forming materials the situation is more involved and the route to the formation of amorphous structures is related to the fact that crystallization would require complex activated phenomena and too large a decrease in entropy.

A breakthrough in the geometrical description of amorphous solids came in 1979 when Kléman and Sadoc first observed that a number of continuous random lattices can be classified as specific mappings of ordered lattices in spaces of constant curvature onto R 3 [14] . This idea was inspired by the observation that whereas regular tetrahedra do not fill three-dimensional Euclidean space, they do regularly tile the three-dimensional sphere S 3 (the manifold described by the equation 4 i=1 x 4 i = R 2 in R 4 ) on which they build a regular polytope (Schläfli symbol {3, 3, 5}) with 120 vertices of coordination number 12. This idea, later developed by several others (see Sadoc and Mosseri [15] , Steinhardt et al. [16] , Nelson [17] , Kléman [18] ), paved the way for a new approach to spatial disorder based on the interplay between order and geometry on three-dimensional manifolds of constant Gaussian curvature.

In-plane order on two-dimensional manifolds has been the subject of much research since the discovery of the ordered phases L β and P β of phospholipid membranes and is now a rich and mature chapter of condensed matter physics [19] . After the seminal work of Nelson and Peliti [20] on shape fluctuations in membranes with crystalline and hexatic order, much work has been done elucidating the intimate relation between in-plane order and the geometry of the underlying substrate with many striking results and even more open questions. The fundamental role of topological defects in two-dimensional systems, first elucidated in a series of pioneering papers by Berezinskii [21] , Kosterlitz and Thouless [22, 23, 24] , is enhanced in the presence of Gaussian curvature in the underlying medium, and leads to structures in the ground state that would normally be highly suppressed in flat systems. The goal of this article is to review the most recent developments in the study of the ground state properties of two-dimensional order on curved media; that is the structure and the mechanics of ordered phases on two-dimensional substrates equipped with non-zero Gaussian curvature, in a regime where thermal fluctuations are negligible in comparison with other energy scales in the system. While focusing on ground states we note that finite temperature physics on curved spaces is a rich source of open problems.

This review is organized as follows. In Sec. 2 we discuss the concept of geometrical frustration and review some fundamentals of the elasticity of topological defects in flat and curved spaces. The study of crystalline and orientational order on the sphere is the basis of most of our current knowledge of order in curved space and will be reviewed in Sec. 3. Sec. 4 is dedicated to crystalline order on surfaces with variable Gaussian curvature and boundary. The existence of defective ground states in toroidal monolayers, with intrinsic crystalline or hexatic order, is a recent development in the study of order on curved surfaces and is reviewed in Sec. 5. Finally, in Sec. 6 , we discuss some current and potential applications of defective structures to materials science and nano-engineering.

Before discussing more technical aspects of the physics of ordered structures on curved surfaces we want to recall some salient features of physical systems with in-pane order and spatial curvature.

Amphiphilic membranes are thin sheets (50 − 100Å) of amphiphilic molecules immersed in a fluid and organized in the form of a bilayer (see Fig. 2 ). The most common constituents of biological membranes are phospholipids consisting of a polar head group and a hydrophobic tail made up of two fatty acyl chains (see Fig.  1 ). Tails have typical length 14 − 20 carbon atoms and regulate the thickness and the stability of the bilayer. The polar head group contains one or more phosphate groups -POOH-O-R. Most phospholipid head groups belong to the phosphoglycerides, which contain glycerol joining the head and the tail. Examples of phosphoglycerides include phosphatidylcholine (PC), phosphatidylethanolamine (PE) and phosphatidylserine (PS). These are distinguished by the residue R carried by the phosphate group (choline: R=-CH 2 -CH 2 -N + -(CH 3 ) 3 and ethanolamine: R=-CH 2 -CH 2 -NH 2 ). The fatty acyl chain in biomembranes usually contains an even number of carbon atoms. They may be saturated (neighboring C atoms are all connected by single bonds) or unsaturated (some neighboring C atoms are connected by double bonds). At low temperature pure phospholipids crystallize and form a bilayer with all the tails in the trans configuration and the heads parallel to the bilayer surface and firmly linked into a lattice (mostly by hydrogen bonds). At higher temperature the order is disrupted and the bilayer becomes fluid. In solution lipid bilayers can be found in a number of phases with differing degrees of order among the hydrocarbon chains. For phospholipids in the PC family these phases are usually called L α , L β and P β . In the L α phase, the tails are liquid and disordered. L β is a gel-like phase in which hydrocarbon tails are ordered and diffusion is severely reduced. The L β → L α transition is usually referred to as main transition. The main transition temperature T m increases with the length of the molecules due to a stronger Van der Waals attraction between adjacent molecules. The degree of unsaturation of the hydrocarbon tails also affects T m . An unsaturated double bond can produce a kink in the alkane chain, disrupting the regular periodic structure and thus lowering the transition temperature. For bilayers in the PC family, T m ranges form −60 to 80 • C , depending on the tail length and the number of double bonds. The order of the hydrocarbon chains also implies a larger thickness of the bilayer. The two lamellar phases can be separated by an intermediate "rippled" phase P β in which the bilayer exhibits an undulated structure and almost solidlike diffusion properties. Hydrocarbon chains can also appear tilted with respect to the bilayer plane. Tilted phases are generally denoted as L β and P β . There are in fact several L β phases characterized by differing amounts of tilt and in-plane orientational order [25] .

Upon changing their concentration, amphiphiles in solution aggregate in a wide variety of structures in addition to bilayers. Above the critical micelle concentration (which is of the order of 10 −3 mol/ ) spherical micelles appear (see Fig. 2 ). Their formation occurs more readily for single-chain amphiphiles (e.g. monoglycerids) and is favored by the presence of large head groups. At higher amphiphile concentrations spherical micelles are replaced by non-spherical ones and eventually by cylindrical rods. Spherical and cylindrical micelles can themselves organize in higher order structures such as cubic lattices or hexagonally packed rod piles. More exotic phases can be obtained by adding oil to the solution. Once the oil is dispersed in water, amphiphiles can form a monolayer across the water-oil interface and self-assemble in complex tubular structures known under the common name of plumber's nightmare [26] .

The name colloidosome was coined by Dinsmore et al. to indicate microcapsules consisting of a shell of coagulated or partially fused colloidal particles surrounding a liquid core [27] . Because of their controllable size, elasticity and permeability, colloidosomes have been recognized to form a promising class of "soft devices" for the encapsulation and delivery of active ingredients with a variety of potential applications for the development of novel drug and vaccine delivery vehicles and for the slow release of cosmetic and food supplements. Their major advantage relies on the fact that the permeability of the shell depends mainly on the size of the gaps between neighboring colloidal particles which can be tuned by controlling their size, interactions and degree of fusion. Velev et al. [28, 29, 30] were the first to report a method for the preparation of colloidosomes by templating octanolin-water emulsions stabilized by latex particles and subsequently removing the octanol core by dissolution in ethanol. Structures of similar architecture have been obtained by templating water-in-oil emulsions [31, 32] . Multilayer shells consisting of alternating positive and negative polyelectrolytes and/or nanoparticles have also been prepared by using layer-by-layer assembly techniques, with the final hollow shells being obtained by removal of the central, sacrificial colloidal particles [33, 34] . Loxley and Vincent [35] developed a new way of preparing polymeric capsules with liquid cores based on a phase separation of the polymer within the templated emulsion. The colloidosomes produced by Dinsmore et al. [27] were obtained by the assembly of latex particles into shells around water-in-oil emulsion drops, followed by thermal fusion of the particles in the shell and centrifugal transfer into water Figure 3 . Scanning electron microscope image of a colloidosome from Dinsmore et al. [27] . The colloidosome is composed of 0.9 µm diameter polystyrene spheres sintered at 105 • C. The close-up on the left shows the effect of sintering at the contact points of neighboring spheres. through a planar oil-water interface (see Fig. 3 ).

The coverage of emulsified droplets by colloidal particles takes place by selfassembly. The particles dispersed in the fluid spontaneously adsorb on the interface provided the surface energy between the fluid on the inside and the outside of the droplet (σ i,o ) is larger than the difference of those between the particles and the internal fluid (σ p,i ) and the particles and the external fluid (σ p,o ). Thus σ i,o > |σ p,i − σ p,o |. A similar mechanism is used in Pickering emulsions, which are stabilized by surface adsorption of colloidal particles. Once adsorbed at the interface, interacting particles distribute evenly, assuring a full and uniform coverage of the droplet. The interactions depend on the type of colloidal particles used as well as the liquids. Coated polymethylmethacrylate (PMMA) or polystyrene spheres, for instance, acquire a permanent electric dipole moment at the interface between the two fluids, probably because of the dissociation of charges on the hydrated surface similar to what happens at a water-air interface [36] . The resultant dipolar interaction stabilizes the particles and allows full coverage of the droplet.

The colloidal particles comprising the shell are then locked together to achieve the desired permeability and robustness of the colloidosome. Several techniques are available to achieve this. By sintering polystyrene particles at a temperature slightly above the glass transition (T g ≈ 100 • C ), it is possible to achieve a partial fusion of neighboring particles at the contact points. This process also allows one to control precisely the size of the gaps between particles and therefore the permeability of the colloidosome. Another method consists of binding the particles with a polyelectrolyte of opposite charge which can bridge neighboring particles and immobilize them at the interface. Particle locking clearly enhances the toughness of the colloidal shell and increases its rupture stress. For sintered polystyrene particles the latter can be tuned within the range 1 − 100 MPa. Colloidosomes locked with polyelectrolytes are even more deformable and can withstand strains of order 50 % before rupturing.

The crystalline arrangement of charged colloids on a hemispherical droplet has been recently studied by Irvine and Chaikin who fabricated colloidal suspensions of PMMA spheres at the interface between water and cyclohexyl bromide (CHB) [37] . Because of the large difference in dielectric constants, ions from the oil strongly partition into the water phase [38] . When nearly 100% PMMA particles are dispersed into the oil phase they form a Wigner crystal and a monolayer near the interface separated by a zone depleted of particles. These last two features are due the ion partitioning that provides the water phase with mobile charges. The net charge of water can be controlled through the pH, shifting the hydrolysis and ion partitioning equilibria. An electrically neutral water droplet acts as a conductor and attracts PMMA particles at the interface through an image charge mechanism. As a result the particles are barely wet by the water phase and organize in a perfect monolayer at (and not across) the interface.

Viral capsids are protein shells that enclose the genetic material of a virus and protect it from enzymatic digestion. Capsid proteins are expressed from the DNA or RNA genome of the virus and in physiological conditions self-assemble in very efficient structures which can withstand high forces (their Young modulus is ∼ 2 MPa) and at the same time effectively disassemble to allow the viral genome to be released in the host cell. Most viral capsids have spherical or rod-like shape, but less standard shapes, such as conical or toroidal, also occur.

The crystallographic structures of spherical viruses have been extensively investigated and, thanks to the modern techniques of X-ray spectroscopy and cryotransmission electron microscopy, form part of the core knowledge of modern virology. In most cases the capsid proteins are grouped in subunits called capsomers, oligomers made of either five (pentamer) or six (hexamer) proteins. Spherical viruses typically posses icosahedral symmetry with twelve pentamers located at the vertices of a regular icosahedron. The number of hexamers that complete the capsids is given by 10(T − 1), where T , the triangulation number, takes values from a sequence of "magic numbers" (i.e. T = 1, 3, 4, 7 . . .) associated with the lattice structure of the capsid, as brilliantly explained by Caspar and Klug (CK) in a seminal 1962 [39] (the CK construction of icosahedral lattices will be reviewed in Sec. 3). The diameters of spherical viruses span the range from 10 to 100 nm. While small capsids are almost perfectly spherical, large viruses, such as the bacteriophage HK97 or the phycodnavirus, typical exhibit a faceted geometry with nearly flat portions separated by ridges and sharp corners corresponding to the twelve pentamers. This morphological difference was explained by Lidmar et al [40] as a buckling transition resulting from a balance between the stretching energy associated with the pentamers in capsomer lattices and the bending elasticity of the viral capsid.

Non-icosahedral capsids of spherocylindrical shape are common among bacteriophages such as some T -even phages as well as the φCBK and the φ29. In this case the capsid has the form of a cylindrical tube composed of a ring of hexamers closed at the ends by two half-icosahedral caps. This structure is also found in a vari- and the cowpea mosaic virus (CPMV) (on the right). Both have 60 capsomers, but the latter has a marked faceted geometry. From VIPERdb [46] . ant of the T = 7 papovavirus and can be induced in other icosahedral viruses by point mutation in the capsid proteins [41, 42] . Of special interest are polymorphic viruses, which can appear in either spherical or spherocylindrical conformations. Polymorphism has been observed in the polyoma/SV40 animal virus [43] and the cowpea chlorotic mottle virus (CCMV) [44] and, for the latter case, appears to be related to the pH and salt concentration of the environment. The human immunodeficiency virus (HIV) also shows broad polymorphism in its capsid shape, including cone-like structures in addition to tubular and spherical ones and has been well studied in recent years [45] .

The science of carbon nano-materials has experienced a period of phenomenal growth since carbon nanotubes (CNTs) were found by Iijima in 1991 [47] and, since then, a large number of similar structures, including helix-shaped graphitic nanotubes [48] , nanotori [49] , graphitic nanocones [50] and nanoflowers [51] have been reported in the literature. The excitement in this field stems from certain exceptional properties that make CNTs potentially useful in many applications in nanotechnology, electronics, optics and other fields of materials science. They exhibit extraordinary strength, unique electrical properties and are efficient conductors of heat. Several methods for the preparation of graphitic nanomaterials have been developed, including laser pyrolysis, arc discharge, and electron irradiation. Recently, metal-catalyzed methods have also been used to synthesize carbon nanomaterials.

Most single wall nanotubes (SWNT) have diameter of about 1 nm, while the length is often of the order of microns. The lattice structure of a SWNT can be obtained from that of a graphene plane by assigning a pair of indices (n, m) which specify how the graphene lattice is rolled up into a seamless cylinder (see Sec. 5). n = 0 nanotubes are referred to as "zigzag", while n = m tubules are called "armchair". The generic name "chiral" is used otherwise [52] . In terms of tensile deformations, SWNT are the stiffest materials known with a Young modulus in the range 1-5 TPa and a tensile strength of 13-53 GPa. This strength results from the covalent sp 2 bonds formed between the individual carbon atoms.

The lattice structure of a nanotube strongly affects its electrical properties because of the interplay between the unique electronic structure of graphene and the tubular geometry. For a given (n, m) structure, a nanotube can be a conductor (if n = m), a small-gap semiconductor (if n − m is a multiple of 3) or a standard semiconductor (otherwise). Thus all armchair (n = m) nanotubes are metallic, and nanotubes with (n, m) equal (5,0), (6, 4) , (9,1), etc. are semiconducting.

Topological defects may occur on the side-wall of carbon nanotubes and related materials in the form of atomic vacancies, 5−7 dislocations and Stone-Wales defects [53] (i.e. quadrupoles consisting of two pairs of 5-membered and 7-membered rings) and are believed to significantly change the mechanical and transport properties of carbon nanotubes. The latter, in particular, were suggested to serve as possible nucleation centers for the formation of dislocations in the original ideal graphene network and constitute the onset of possible plastic deformations [54] . Electronic transport is affected by the presence of defects by lowering the conductivity through the defective region of the tube. Defect formation in armchair-type tubes (which can conduct electricity) can cause the region surrounding that defect to become semiconducting. Furthermore, single monoatomic vacancies induce magnetic properties. Phonon scattering by defective regions heavily affects the thermal properties of carbon nanotubes and leads to an overall reduction of the thermal conductivity. Phonon transport simulations indicate that substitutional defects such as nitrogen or boron will primarily lead to scattering of high frequency optical phonons. Larger scale defects, however, such as Stone Wales defects, cause phonon scattering over a wide range of frequencies, leading to a greater reduction in thermal conductivity [55] . On the other hand, defective regions appear to be natural places for chemical functionalization. Numerical simulations have shown how the presence of Stone-Wales defects considerably enhances the adsorption of carboxyl groups (COOH) which can then bind to molecules with amide and ester bonds [56] .

The notion of geometrical frustration was introduced to describe situations where certain types of local order, favoured by physical interactions, cannot propagate throughout a system [57] . The expression was used for the first time by Toulouse in 1977 [58] to describe certain particular magnetic systems with nearest-neighbours interactions which cannot all be satisfied simultaneously. A textbook example of frustration in magnetic models is represented by a system of Ising spins on a triangular lattice with antiferromagnetic bonds: while a perfect antiferromagnetic alignment would minimize all terms in the Ising Hamiltonian, such an alignment is not allowed by the topology of the underlying lattice so that for any triangular plaquette there is always at least one unsatisfied bond (see Fig. 6 ). This concept can be extended naturally to any system where interactions impose a local order, but the most favoured local configuration is geometrically incompatible with the structure of the embedding space.

Two-dimensional manifolds equipped with some microscopic field for which a notion of local order can be defined unambiguously provide a paradigm for systems exhibiting geometrical frustration. Consider for instance an assembly of identical particles interacting with a spherically symmetric pair potential V ij = V (|x i −x j |), with x i the position vector of the ith particle in a suitable coordinate system. In flat two-dimensional space, particles almost always pack in triangular lattices, unless the interaction potential is carefully tuned to select some some other lattice topology. Endowing the medium with a non-planar topology introduces frustration in the sense that the energetically favoured 6−fold orientational order can no longer be established everywhere in the system. Such geometrical frustration arises at the microscopic level from the celebrated Euler theorem of topology which relates the number of vertices V , edges E and faces F of any tessellation of a 2−manifold M :

where χ is the Euler characteristic of M . If M is an orientable closed surface, one 

(Center) Vector field of a sphere. As a consequence of the Poincaré-Hopf theorem, a vector field must vanish at least at two points, corresponding to the north and south pole of the sphere in this example. (Right) Triangulation of the sphere. As required by the Euler theorem, any triangulation of the sphere must contain some vertices of coordination number different from six (i.e. twelve 5−fold vertices in this case).

can show that χ is determined uniquely by an integer g ≥ 0, called the genus of M , which represents the number of "handles" of M ; namely χ = 2(1 − g). Two orientable closed surfaces with the same genus (thus the same Euler characteristic) are homeomorphic: they can be mapped into one another without changing their topological properties. In a surface with boundary, the Euler characteristic is given by

where h is the number of boundaries or "holes" of M . Thus a sphere, which has no handles nor boundary (g = h = 0) has χ = 2, while the embedded torus (g = 1 and h = 0) has χ = 0. A disk, on the other hand, has χ = 1 (g = 0 and h = 1) and is topologically equivalent to a sphere with one hole. In Secs. 3-5 we discuss ordered structures on three important topologies: the sphere, the disk and the torus.

In the case of two-dimensional crystals with 6−fold local order Eq. (1) can be rephrased in a form that is particularly useful in describing the presence of defects in the lowest energy state by defining a topological charge as the departure from the ideal coordination number of a planar triangular lattice: q i = 6 − c i , with c i the coordination number of the ith vertex. Now, consider a tessellation in which each face is an n−sided polygon and let k faces meet at each vertex. Since each edge is shared between two faces and links two vertices it follows that:

where V k is the number of vertices of degree k. For a triangulation n = 3. From Eq. (1) it follows then:

In the case of a sphere, with χ = 2, Eq. (2) implies any triangulation contains defective sites such that the total topological charge of the lattice is Q = 12. This can be achieved, for example, by incorporating twelve 5−fold disclinations (with q = 1) in a network of 6−fold coordinated sites like in a common soccer ball. These twelve disclinations, which would be suppressed in the lowest energy state of a planar crystal, result from the geometrical frustration associated with the topology of the sphere. Eq. (2) is a special case of geometrical frustration on 2-manifolds where local orientations are defined modulo π/3. More generally one can consider a p−atic director field for which local orientations are defined modulo 2π/p. The topological charge of a disclination in this case is q = ∆θ/ 2π p , where ∆θ is the angle the director rotates in one counter-clockwise circuit of any closed contour enclosing the defect. Eq. (2) becomes then:

where N is the total number of defects. The ratio k = q/p is commonly referred to as the winding number of the disclination. In the case of a simple vector field, for instance, p = 1 and Eq. (3) corresponds to the well known Poincaré-Hopf theorem according to which the sum of the indices of all the isolated zeros of a vector field on a oriented differentiable manifold M is equal to the Euler characteristic χ of M . Thus for a sphere a vector field must have at least one sink and one source, each having topological charge one, while on a torus (χ = 0) a vector field can be defect free. A nematic director n, on the other hand, has p = 2 (i.e. physical configurations are invariant under inversions n → −n). Thus disclinations with ±1 topological charge correspond to configurations where the director rotates ±π in one circuit enclosing the defect. These elementary disclinations have semi-integer winding number k = ±1/2. As a consequence of Eq. (3) the total topological charge of a nematic texture on the sphere is Q = 4, corresponding for instance to the typical baseball texture consisting of four q = 1 (k = 1/2) disclinations located at the vertices of a regular tetrahedron [59, 60] .

The equilibrium structure of two-dimensional locally ordered systems on curved substrates depends crucially on the existence and arrangement of the defects. The seminal work of Kosterlitz, Thouless, Halperin, Nelson and Young (KTHNY) on defect-mediated melting in two dimensions makes it clear that it is useful to employ a theoretical framework where the fundamental objects in the system are the defects themselves, while treating the microscopic constituents within a continuum elastic theory. This approach has the advantage of far fewer degrees of freedom than a direct treatment of the microscopic interactions and allows one to explore the origin of the emergent symmetry observed in non-Euclidean ordered structures as a result of the interplay between defects and geometry. This interplay is one of the fundamental hallmarks of order on two-dimensional manifolds and leads to universal features observed in systems as different as viral capsids and carbon macromolecules. In this section we will briefly review some concepts of differential geometry, mostly to establish notation. In the next section we will review some fundamentals of the elasticity of defects in two dimensions. The coupling mechanism between curvature and defects will be introduced in Sec. 2.4.

Points on a two-dimensional surface S embedded in R 3 are specified by a three-dimensional vector R(x) as a function of a two-dimensional coordinate x = (x 1 , x 2 ). For each point of S we define three vectors:

and

where ∂ i = ∂/∂x i . The vectors g i belong to T R S, the tangent space of S at R whereas n is a normal vector. Note that while n is a unit vector, g i are generally not of unit length. The metric, or first fundamental form, of S is defined as:

where g ij is the metric tensor:

The dual tensor is denoted as g ij and is such that:

with δ j i the Kronecker symbol. This allows us to introduce contravariant tangentplane vectors g i = g ij g j , satisfying g i · g j = δ j i . Any vector v on the tangent plane can be expressed as a linear combination of basis vectors g i or

The extrinsic curvature of the surface S is encoded in the second fundamental form b ij (also known as the extrinsic curvature tensor):

The eigendirections of b ij at a given point correspond to the principal curvature directions of S at that point and the associated eigenvalues κ 1 and κ 2 are the extremal (or principal) curvatures. The mean curvature H and the Gaussian curvature K are defined as the sum and the product of the principal curvatures:

where ij is the dual of the Levi-Civita tensor, whose components are given by:

where g = det g ij . The contravariant form is ij = ij /g and satisfies ik jk = δ j i . Since ij v i v j = 0, where v i is any contravariant vector, it follows that the vector ij v i is perpendicular to v j . Thus inner multiplication by ij rotates a vector by π/2. The covariant derivative of a vector field v in the ith coordinate direction is defined as usual by:

where Γ k ij is the Christoffel symbol:

Both the metric and the Levi-Civita tensor are invariant under parallel transport. This translates into:

Much of the elastic theory of defects, either in flat or curved systems, relies on the calculation of the Green function of the Laplace operator. On a generic 2−manifold the latter obeys:

where ∆ is the Laplace-Beltrami operator:

and δ the delta-function:

The Stokes theorem is frequently invoked when calculating elastic energies of defects. In covariant form it can be stated as follows: given a vector field v on a 2−manifold M with boundary ∂M , the following identity holds:

For consistency we will adopt covariant notation throughout this review. When discussing planar systems, in particular, we have: b ij = H = K = 0 and the elements of the metric tensor are given in Cartesian coordinates by g xx = g yy = 1, g xy = g yx = 0 and in polar coordinates by g rr = 1, g φφ = r 2 , g rφ = g φr = 0. The definition of orientational order on a surface clearly requires a nonambiguous notion of angular distance between vectors on the same tangent plane. This is traditionally achieved by introducing a pair of orthonormal vectors e α (α = 1, 2) called vielbein (note that the canonical coordinate vectors g i are generally neither orthonormal nor orthogonal) so that:

where δ αβ is the usual Kronecker symbol in the indices α and β. A vector field v = v i g i can be expressed alternatively in the basis e α :

Since the coordinates v α are locally Cartesian, δ αβ = δ αβ and there is no distinction between upper and lower Greek indices: v α = v α . Vielbein are constructed to be invariant under parallel transport, and thus ∇ i (e α ) j = 0 and:

where Ω iαβ is the so called spin connection, and is given by:

Taking the derivative of the left equation in (17) one immediately sees:

from which it follows that Ω iαβ = −Ω iβα . In two dimensions this permits the parametrization of the spin-connection Ω iαβ by a single covariant vector Ω (i.e. the antisymmetry under exchange of the Greek indices leaves only d 2 (d − 1)/2 = 2 independent components in d = 2 dimensions).

Thus

with αβ = δ 1 α δ 2 β − δ 1 β δ 2 α the antisymmetric symbol. It is well known how the curvature of a manifold manifests itself when a vector is parallel transported around a closed path. Taking an infinitesimal square loop of sides dx and dy, the parallel transported vector v differs from the original vector v by an amount:

where R k lij is the Riemann tensor given in two-dimensions by:

In the orthonormal basis e α , Eq. (19) reads:

where R ijαβ is the curvature tensor associated with the spin-connection Ω iαβ :

In two dimensions, using Eq. (20) , one can write:

which, combined with Eqs. (21) and (18), yields which also implies:

The XY −model is the simplest setting where particle-like objects emerge from a purely continuum theory [22, 23, 24] . The order parameter is the angular field θ ∈ [0, 2π], which may represent the orientation of two-dimensional vectors S = S(cos θ, sin θ) or the phase of a complex field ψ = |ψ|e iθ . The interaction that tends to minimize the spatial variations of the order parameter results from the continuum free energy:

Despite its simplicity this model successfully describes several aspects of the physics of vortices in superfluid 4 He or thin superconducting films, where the angle θ is identified with the phase of the collective wave function. Eq. (24) also describes nematic liquid crystals in the limit of equal splay and bending moduli. This approximation, however, favors q = ±2 (k = ±1) disclinations rather than the more natural q = ±1 (k = ±1/2) disclinations of nematics and is not particularly suitable to describe the ground state. At T > 0, thermal fluctuations drive the two elastic constants to the same value at long wavelengths, so that there is a unique Kosterlitz-Thouless transition temperature [61] . As discussed by Deem [62] , the essential effect of unequal elastic constants is to create a distinct long-range contribution to the core energy of each defect. More generally Eq. (24) can be considered as the simplest phenomenological free energy describing p−atic order (assuming θ defined modulo 2π/p). Calling v = ∇θ, the minimization of the free energy (24) leads to the equilibrium condition:

In the presence of a number of disclinations of topological charge q α , θ changes by (2π/p) α q α in one circuit around any contour enclosing a total topological charge α q α :

Using Stokes' theorem, Eq. (26) can be translated into the requirement:

where:

is the topological charge density. Using standard manipulations (see for example [63] ) a vector field v satisfying Eq. (25) with the constraint (26) can be found in the form:

where G L (x, y) is the Laplacian Green function. Using Eq. (28) in Eq. (24) leads to the well known expression:

Like charged particles, disclinations interact via a Coulomb potential, which in two dimensions is proportional to the logarithm of the distance in the plane. As noted, this framework is valid for both nematic liquid crystals in the one elastic constant approximation and superfluids. An important difference between these two systems lies in the choice of the boundary condition for the field θ. Nematogens are typically forced to be normal to boundary of the substrate and this implies a constraint for θ, while there is no such constraint for 4 He films since the wave function is defined in a different space from that to which the superfluid is confined. This difference is crucial on a curved substrate (see Ref. [64] for a detailed review of this topic). Eq. (29) represents the elastic energy associated with the distortion introduced by defects in the far field, where the elastic variables change slowly in space. This expression breaks down in the neighborhood (or core) of a defect, where the order parameter is destroyed and the actual energetic contribution depends on microscopic details. In order to describe a defective system at any length scale, Eq. (29) must be corrected by adding a core energy F c representing the energetic contribution within the core of a defect, where standard elasticity breaks down. A detailed calculation of the core energy requires some microscopic model and is usually quite complicated. Nonetheless its order of magnitude can be estimated by writing:

where a is the core radius, corresponding to the short distance cut-off of the elastic theory, and f c is some unknown energy density independent of a. f c can then be estimated by minimizing the total energy of the system with respect to a. For a planar system with a single disclination of winding numer k, the total energy can be easily calculated from Eqs. (29) and (30):

Minimizing Eq. (31) with respect to a, one finds f c = K A k 2 /(2a 2 ) from which:

The quantity c is the energy needed to increase the number of defects by one, regardless of its location. The elasticity of dislocations and disclinations in solids resembles in many respects that of vortex lines in the XY −model. In two dimensional elasticity a pure in-plane deformation is encoded in a displacement field u i , i = 1, 2, which maps any point in the system x to:

where g i is a suitable basis in the coordinates of the undeformed system. If there are no defects, the displacement field is a single-valued mapping of the plane into itself. Topological defects introduce an incompatibility in the displacement field, in the sense that u i is not a single-valued mapping anymore. The elastic stress in the region surrounding a defect appears at the macroscopic level from the Hooke's law of elasticity:

where Y and ν are the two-dimensional Young modulus and Poisson ratio respectively and σ ij is the stress tensor. Eq. (34) can be obtained, for example, by minimizing the elastic energy:

where λ and µ are the Lamé coefficients in two dimensions:

In the absence of body forces, the force balance equation requires σ ij to be divergence free:

The strain tensor u ij represents the variation in the first fundamental form of the surface due to the deformation field (33), namely:

In the presence of a dislocation line L the function u becomes multivalued so that, while traversing any closed counterclockwise loop C containing L:

where b is the Burgers vector representing the amount by which the image of a closed loop under the mapping (33) fails to close in presence of a dislocation line [65, 66] . If L is an isolated straight dislocation with origin at the point x 0 , Eq. (38) implies:

Generally speaking, maps such as that of Eq. (33), induce variations in both the length and the orientation of an infinitesimal distance vector dx in the deformed medium. This statement can be clarified by writing:

is the infinitesimal rotation tensor induced by the deformation Eq. (33) . In two dimensions the latter can be conveniently written in the form:

with:

Since the application of ij to a vector rotates the vector by π/2 clockwise, the coefficient Θ in Eq. (42) can be interpreted as the average rotation of an infinitesimal line element under a distortion of the form Eq. (33) . In presence of an isolated disclination the field Θ becomes multi-valued. If s is the angular deficit associated with the disclination, integrating along a circuit enclosing the disclination core yields:

which can be rephrased as a statement about commutativity of partial derivatives:

where x 0 is the position of the disclination core. In a lattice with n−fold rotational symmetry, s = (2π/n)q, where q is the topological charge introduced above. A comparison between Eq. (27) and Eq. (44) clarifies the common identification of Θ with the bond angle field of a crystal [67] . Eqs. (39) and (44) provide a set of relations between the fundamental features of topological defects in crystals (i.e. the Burgers vector b and the topological charge q) and the partial derivatives of the displacement field u. These relations can be used to derive an equation for the elastic stress arising in a system as a consequence of the defects. A simple and elegant way to achieve this task in the context of planar elasticity relies on the parametrization of the stress tensor via a single scalar field χ known as the Airy stress function (see for example [68] ). Taking advantage of the commutativity of partial derivatives in Euclidean space, one can write

so that the force balance condition (36) is automatically satisfied. Applying the operator ik jl ∇ k ∇ l to both sides of Eq. (34) and using Eq. (45) one finds the following fourth order Poisson-like problem:

where ∆ 2 is the biharmonic operator and η is the defect charge density:

where s α denotes the topological charge of a disclination at x α and b β the Burgers vector of a dislocation at x β . Let's now turn our attention to the case of a planar crystal with local 6−fold orientational symmetry populated with N disclinations of density:

Assuming a free boundary (i.e. all the components of the stress tensor are zero along the boundary and thus χ = ∇χ = 0), the stress function χ can be expressed in the Green form:

where G 2L (x, y) is the biharmonic Green function. The elastic energy of the system is given by Eq. (35) and (48) in the form:

Eqns. (24) and (49) are the fundamental equations in the elastic theory of defects in two-dimensional planar p−atics and crystals. In the next section we will see how a non-zero Gaussian curvature in the underlying medium affects these energies by effectively screening the topological charge of the defects.

The elasticity of topological defects on curved substrates equipped with local orientational order was first considered by Nelson and Peliti in the context of hexatic membranes [20] . It is now standard knowledge in the statistical mechanics of membranes that long range forces appearing as a consequence of local orientational order crucially affect the behavior of membranes at finite temperature. The stiffness associated with orientational correlations leads to an enhancement of the bending rigidity that counteracts the thermal softening that occurs in fluid membranes. Furthermore, for planar membranes, the stabilizing effect induced by the orientational stiffness opposes the entropy-driven tendency to crumple, causing a transition between a flat and crumpled phase at T > 0. In this review we focus on the ground state properties of ordered structures on curved surfaces and we refer the reader to the literature for a discussion of finite temperature physics [19] . In this section we will show how an underlying non-zero Gaussian curvature couples with the defects by screening their topological charge. As a result, topological defects, whose existence would be energetically suppressed if the same system was lying on a flat substrate, instead proliferate. The universality of such a curvature screening mechanism is remarkable in the sense that it occurs in a conceptually identical fashion in both p−atics and solids, although the kernel of the elastic interactions between defects differs. As a consequence of curvature screening topological defects organize themselves on a rigid surface so as to match the Gaussian curvature of the substrate. On the other hand, if the geometry of the substrate is allowed to change (for instance by lowering the bending rigidity) the system eventually enters a regime where the substrate itself changes shape in order to accommodate some preferential in-plane order. The latter is the fundamental mechanism behind the buckling of crystalline membranes [69] and is believed to be the origin of the polyhedral geometry of large spherical viruses such as bacteriophages [70] .

In the case of manifolds equipped with a pure rotational degree of freedom like p−atics, the curvature affects the elastic energy through the connection which determine how vectors change when parallel transported. This statement can be clarified by rewriting the elastic energy (24) in the form:

where m is a p−atic director field which can be conveniently expressed in a local orthonormal frame e α (α = 1 2): m = cos θ e 1 + sin θ e 2 , with θ defined modulo 2π/p. Since ∂ i m α = − αβ m β ∂ i θ, using the properties of vielbein outlined in Sec. 2.2, we have:

The elastic energy (50) thus becomes

Using Eq. (23), the vector field Ω can be expressed as:

where G L (x, y) is again the Green function of the Laplace-Beltrami operator (14) . Substituting Eqns. (53) and (28) in the expression for the elastic energy (52) we obtain:

As anticipated, the Gaussian curvature of the underlying substrate couples with the defects by screening their topological charge. Eq (54) is identical to the Coulomb energy of a multi-component plasma of charge density η in a background of charge density −K. Like particles in a plasma, we can expect defects to screen each other's charge, thus forming clusters of zero net charge, and matching the charge distribution of the surrounding background. This implies that disclinations will be attracted to regions of like-sign Gaussian curvature. Crystalline surfaces (i.e. non-Euclidean crystals) differ from manifolds equipped with a p−atic director field because bonds, whose local orientation is encoded in the field θ in p−atics, can now be compressed and sheared to relieve part of the elastic stress due to the geometrical frustration provided by the embedding manifold. A first attempt at describing the elasticity of defects on a two-dimensional curved crystal was made by Dodgson [71] , who studied the ground state of the Abrikosov fulx lattice in a model thin-film superconductor on a sphere (subject to a field radiating from a magnetic monopole at the center) and found evidence for twelve 5−fold disclinations at the vertices of an icosahedron in a otherwise 6−fold coordinated environment (see Sec. 3 for a detailed review of spherical crystals). Later, Dodgson and Moore [72] proposed adding dislocations to the ground state of a sufficiently large spherical vortex crystal to screen out the strain introduced by the twelve, topologically required, disclinations. A general framework to describe the elasticity of defects in non-Euclidean crystals was proposed by Bowick, Nelson and Travesset (BNT) in 2000 [73] . The BNT model, obtained from the covariantization of the elastic energy (49) , relies on the following expression for the elastic energy of a collection of disclinations in a triangular lattice of underlying Gaussian curvature K:

where η(x) is the topological charge density (47) and G 2L (x, y) the Green function of the covariant biharmonic operator on the manifold. The origin of the coupling between topological charge and curvature, in this case, is rooted in a profound result of discrete geometry originally due to Descartes which can be considered the oldest ancestor of the Gauss-Bonnet theorem of differential geometry. Let P be a convex polyhedron and define the angular deficit of a vertex v of P as:

where α i (v) with i ∈ [1, c(v)] are the angles formed by all the faces meeting at v. Descartes' theorem states that the sum of the angular deficits of a convex polyhedron is equal to 4π. More generally:

which is exactly the Gauss-Bonnet theorem for the case in which the Gaussian curvature is concentrated in a finite number of points (vertices) rather than smoothly Figure 9 . (Color online) Parallel transport of a vector around a corner of a cube. As a consequence of the discrete Gaussian curvature, the vector is rotated by π/2 after parallel transport.

distributed across the whole surface. The deficit angle k(v) is thus the discrete analog of the Gaussian curvature. This analogy is not limited exclusively to Eq. (56). Consider, for example, the corner of a cube and imagine parallel transporting a vector v around a closed loop encircling the corner (see Fig. 9 ). After parallel transport, the vector has rotated by k(v) = π/2 with respect to its original orientation. Analogously, a vector parallel transported around a closed loop on a surface rotates by an angle

where the integral is over the portion of the surface enclosed by the loop. Now, on a triangulated surface α(v) ≈ π/3 and:

Thus the source term η − K figuring in Eq. (55) corresponds to the difference between the pre-existing Gaussian curvature of the manifold with local 6−fold orientational order and the additional Gaussian curvature induced by a disclination of topological charge q. A more formal discussion can be found in Ref. [74] . As in the case of a purely rotational degree of freedom, disclinations in non-Euclidean crystals organize so as to match the Gaussian curvature of the underlying medium while maximizing their reciprocal distance as a consequence of the repulsive interaction between like-sign defects. In the next three sections we will see how the formalism outlined here applies to three physically relevant examples of manifolds: the sphere, the topological disk and the torus.

3. Crystalline and nematic order on the sphere

Experimental realizations of spherical crystals are found in emulsions of two immiscible fluids, such as oil and water, stabilized against droplet coalescence by coating droplets of one phase (say water) with small colloidal particles [75, 76] . The colloidal "armor plating" of these Pickering emulsions also plays a role in colloidosomes, colloid-coated lipid bilayer vesicles used for encapsulation and delivery of flavors, fragrances and drugs [27] . Identical micron-sized particles tend to crystallize under typical experimental conditions. The defect structure of crystalline ground states on a sphere is important in this context since defects influence the strength of the colloidal armor.

Defects play an essential role in describing crystalline particle packings on the sphere. At least twelve particles with 5-fold coordination (i.e., 12 disclinations) are required for topological reasons, and like the 5-fold rings in carbon fullerenes, one might expect that the energy would be minimized if the disclination positions approximated the vertices of a regular icosahedron. This expectation, which also plays a role in geodesic domes and in the protein capsomere configurations of spherical virus shells [39, 77] , is nevertheless violated when the shells are sufficiently large and disclination buckling [40] out of the spherical environment is suppressed by surface tension.

To understand the structure of spherical crystals, we must find the ground state of M particles distributed on a surface Σ and interacting with, say, pairwise repulsive potentials. For simplicity consider power law potentials V (r) = e 2 /r γ , where e is a generalized electric charge. The energy of M particles at positions r( ), with = (l 1 , l 2 ) ∈ Z⊗Z, becomes 2E 0 = = e 2 /|r( )−r( )| γ . This covers most of the interesting physical systems. The case γ = 1 corresponds to the pure 3D Coulomb potential, appropriate for charged colloids. This is known as the Thomson problem in the physics and mathematics literature and has been widely studied for over 100 years [78, 79, 80] . The case γ = 3 corresponds to a dipole interaction, appropriate for neutral colloids at the interface between two liquids. The different dielectric constants of the two liquids leads to an asymmetric distribution of charge on the colloids and a net dipole moment. The case γ = 12 is the repulsive part of the Lennard-Jones potential and is the important piece of the interaction for driving crystallization. Finally the limit γ → ∞ is an extreme short range interaction that selects the pair of particles with the minimum distance. The equivalent packing of spherical caps was introduced in 1930 by the biologist Tammes [81] and connects the present problem to the rich physics and mathematics of packing [82, 83] . The packing of particles on surfaces of nontrivial topology is relevant to a broad range of physical, chemical and biological systems in addition to its intrinsic mathematical interest. An almost literal realization of the Thomson problem is provided by multi-electron bubbles [84, 85] . Electrons trapped on the surface of liquid helium have long been used to investigate two dimensional melting [86, 87] . Multi-electron bubbles result when a large number of electrons (10 5 −10 7 ) at the helium interface subduct in response to an increase in the anode potential and coat the inside wall of a helium vapor sphere of radius 10−100µm. Typical electron spacings, both at the interface and on the sphere, are of order 2000Å, so the physics is entirely classical, in contrast to the quantum problem of electron shells which originally motivated J.J. Thomson [78] . Information about electron configurations on these bubbles can, in principle, be inferred from studying capillary wave excitations [88, 89] . Similar electron configurations should arise on the surface of liquid metal drops confined in Paul traps [90] .

A Thomson-like problem also arises in determining the arrangements of the protein subunits which comprise the shells of spherical viruses [39, 40, 45, 77, 91] . Here, the "particles" are clusters of protein subunits arranged on a shell. This problem of protein arrangements was solved in a beautiful paper by Caspar and Klug [39] for intermediate values of R/a, where R is the sphere radius and a is the mean particle spacing. Caspar and Klug constructed icosadeltahedral particle packings characterized by integers m and n, which provide regular tessellations of magic numbers M = 10(n 2 + nm + m 2 ) + 2 particles on the sphere. The Caspar-Klug tessellations provide an excellent starting point for finding low energy particle configurations on the sphere for intermediate values of

Particle numbers M that are not magic numbers can be accommodated by introducing vacancies or interstitials into these regular packings.

Other realizations of Thomson-like problems include regular arrangements of colloidal particles in colloidosome cages [92, 93, 94, 95] proposed for protection of cells or drug-containing vesicles [27] and fullerene patterns of spherical carbon atoms [96] . An example with long range (logarithmic) interactions is provided by the Abrikosov lattice of vortices which would form at low temperatures in a superconducting metal shell with a large monopole (solenoid tip) at the center [72] .

The problem of the best possible packing on spheres also has applications in the micro-patterning of spherical particles [97] , relevant for photonic crystals, and in understanding the structure of Clathrin cages, responsible for the vesicular transport of cargo in cells [98] . Crystalline domains covering a fraction of the sphere are also of experimental interest. In the context of lipid rafts [99] , confocal fluorescence microscopy studies have revealed the coexistence of fluid and solid domains on giant unilamellar vesicles made of lipid mixtures. The shapes of these solid domains include stripes of various widths and orientations [100, 101, 102, 103] .

Consider a collection of classical point particles constrained to a frozen (nondynamical) two-dimensional surface Σ embedded in three-dimensional Euclidean space. The particles interact through a general potential defined in the three dimensional embedding space or solely within the 2D curved surface itself. We focus primarily on the potential

Here, e is an "electric charge" such that if r is some quantity with dimensions of length then e 2 /r γ has dimension of energy. The case γ = 1 corresponds to the Coulomb potential in three dimensions. Although we do not treat this problem explicitly here, the replacement

allows us to treat the two dimensional Coulomb potential by taking the limit γ → 0,

The electrostatic energy of a system of M particles at positions r( ), interacting via Eq.(58), with = (l 1 , l 2 ), l 1 , l 2 ∈ Z, becomes

Note that with this definition the power law interaction acts across a chord of the sphere, as would be the case for electron bubbles in helium. Although we focus here on curved crystals, there are some quantities which are insensitive to the curvature of the surface, and the simpler geometry of the plane can be used to compute them. The following two subsections hence focus on planar crystalline arrays of particles interacting via the potential Eq.(58).

The electrostatic energy Eq. (61) and the corresponding elastic tensor, from which follows the elastic constants of the system, may be explicitly computed for crystalline orderings of particles in a triangular lattice. For any non-compact surface, like the plane, the energy Eq. (61) is divergent for all γ ≤ 2. If γ > 0, the divergence comes exclusively from the zero mode G = 0 associated with the thermodynamic limit of infinite system size. This term (which would be subtracted off if a uniform background charge were present) can be isolated by setting |G| ≡ ε 1 for this mode. The final result for the ground state energy reads

A C is the area of the unit cell of the triangular Bravais lattice (A C = √ 3 2 a 2 ) and Γ is the Euler Gamma function. The coefficient σ is a sum over Misra functions. The coefficient θ(γ) parameterizes the nonsingular part of the energy; its dependence on the exponent γ is shown in Table 1 . This negative quantity parameterizes the binding energy of the triangular lattice after the positive "zero mode" contribution is subtracted off. For γ = 1, we have a two dimensional "jellium" model. In the case of a sphere no neutralizing background is necessary and the energy is rendered finite by the compact nature of the sphere itself. The maximum distance between points on the sphere provides a natural infrared cut-off.

For small displacements of the particles from their equilibrium positions, one has

where u( ) is a small displacement of the particle in the plane of the surface from its equilibrium position r( ), and therefore a tangent vector to the surface. The elastic tensor Π ij ( , ) is defined as the leading term in an expansion of Eq. (63),

In deriving Eq.(63), we assume a constraint of fixed area per particle, enforced by a uniform background charge density or boundary conditions. This eliminates the 

The coefficients η(γ) and ρ(γ) depend only on the potential. In Table 1 , some values of the coefficients for a range of potentials with 0 < γ < 2 are listed.

When the deviations from the ground state are small, the long wavelength lattice deformations may be described by a continuous Landau elastic energy of the form (35) . The elastic tensor Eq. (64), within Landau elastic theory, is then

A comparison with Eq. (65) immediately yields an explicit expression for the elastic constants of the crystal

where Y is Young's modulus. The result λ = ∞ is equivalent to a divergent compressional sound velocity as p → 0 and for γ = 1 is just a statement of the incompressibility of a two-dimensional Wigner crystal. Alternatively, we can allow for wavevector-dependent elastic constants µ(p) and λ(p) in Eq. (66) . In this case . λ(p) diverges as p → 0:

while lim p→0 µ(p) is given by Eq. (67).

Spherical crystals have many properties not shared by planar ones, one of the most remarkable being that there is an excess of twelve positive (5−fold) disclinations. These disclinations repel, and the simplest spherical crystals will be those having the minimum number of defects (i.e. twelve) located at the vertices of an icosahedron. Triangular lattices on the sphere with an icosahedral defect pattern are classified by a pair of integers (n, m), as illustrated in Fig. 10 . The path from one disclination to a neighboring disclination for an (n, m) icosadeltahedral lattice consists of n straight steps, a subsequent 60 • turn, and m final straight steps. The geodesic distance between nearest-neighbor disclinations on a sphere of radius R is d = R cos −1 (1/ √ 5). The total number of particles M on the sphere described by this (n, m)-lattice is M = N nm , where

Such (n, m) configurations were once believed to be ground states for relatively small numbers (M ≤ 300, say) of particles interacting through a Coulomb potential [104, 105, 106, 107, 108] . The energy of discrete particle arrays described by Eq. (61) can be evaluated by starting with some configuration close to an (n, m) one and relaxing it to find a minimum. It is found that the (n, m) configurations are always local minima. Whether these icosahedral configurations are global minima as well will be analyzed later. Results for the energy E(M ) are shown in the inset to Fig. 11 . From Fig. 11 it is clear that energies grow very fast for increasing volume. More interestingly, the (n, 0) and (n, n) configurations show a growing difference in energy for increasing volume size, implying that the energy of icosahedral configurations does not tend to a universal value for large numbers of particles but rather remains sensitive to the (n, m) configuration, a result also noted by other authors [105] . Further insight comes from investigating the distribution of energy. Plots of the local electrostatic energy, the electrostatic energy at point x on the sphere, as defined in Eq. (61) are shown in Fig. 12 [109] . From Fig. 12 it should be noted that the triangles obtained by the Voronoi-Delaunay construction, after minimization of the potential Eq. (61), are very close to equilateral. The distribution of the local energies for the (n, 0) and (n, n) configurations are very different. The (n, 0) configuration shows maximum energies along the paths joining the defects. The (n, n) configuration, on the other hand, has its maximum energies along the directions defined by the triangles formed by three nearest neighbor disclinations. The size of these regions of differing electrostatic energy turns out to scale with system size, making it very plausible that there might be small differences in the energy per particle for (n, 0) and (n, n) configurations in the limit n → ∞.

Spherical substrates provide the simplest example of the problem of crystals on curved surfaces. The study of spherical crystals is simplified by two important properties: there is a unique scale with dimensions of length, the radius R, and there is a fixed excess disclinicity of twelve following from the Euler theorem (1)

The free energy Eq. (55), applied to the sphere, is tractable analytically because the inverse biharmonic operator on a sphere of radius R can be computed explicitly. It is shown in [73] that the Green's function for the biharmonic, in spherical coordinates (θ, φ), has the following simple form on a unit sphere:

where β is the geodesic distance between two disclinations located at (θ a , φ a ) and

The total energy of a spherical crystal with an arbitrary number of disclinations follows from Eq. (55) and Eq. (72) and has the simple form [73] 

where {θ i , φ i } i=1,··· ,N are the coordinates of N defects and we restrict ourselves to 5−fold (q i = +1)and 7−fold (q i = −1) defects. The quantity E 0 is the zero point energy and is defined in Eq. (62) . Although 5 and 7-fold disclinations will in general have different core energies [110] , we assume equal core energies here for simplicity. What matters for our calculations in any case is the dislocation core energy E d , which we take to be

The value of the Young's modulus and the flat space ground state energy E 0 have been computed in Sec. 3.2.1. When the sphere radius R is large compared to the particle spacing a, we can use flat space values of Y and the flat space energy E 0 (M ) associated with a finite number of particles M . To obtain the leading terms in the expansion of the ground state energy for large but finite M , the precise compactification of the plane employed is irrelevant -it may be achieved by periodic boundary conditions, for example. For a sufficiently large plane the finite size effects will be negligible. The density σ of particles is then M divided by the total surface of the compact plane, taken to be the surface area S of a sphere of radius R,

From Eq. (67) the expression for the Young's modulus suitable for M particles on a spherical crystal of radius R with 0 < γ < 2 is then

One remaining detail is the divergent contribution to the energy E 0 in Eq. (62) .

Since the divergent part comes solely from the zero mode, the spatial variations in the density of the actual distribution are irrelevant. It may therefore be computed for a uniform density of charges. The divergent part is identical to the energy of a constant continuum of charges as described by the density Eq. (75) . We now evaluate this divergent part of the energy on a sphere, instead of a plane.

The divergent part has thus been regularized, and the energy is finite and welldefined for all M < ∞. Note that for the case γ < 2,

Hence E D is not simply a function of the particle density M/S, as one would expect for a short range interaction.

Upon substituting the elastic constant of Eq. (76) into Eq. (74) , one arrives at

where E 0 is defined in Eq. (62) and the function C(i 1 · · · i N ) depends on the position i 1 = (θ 1 , φ 1 ) etc. of the N disclination charges and is universal with respect to the potential. The total energy of a spherical crystal, including the contributions to E 0 is then

Note that the leading correction to the zero mode energy proportional to M 2 varies as M 1+γ/2 , and depends both on the flat space function θ(γ) and on the C-coefficient

associated with a particular configuration of disclinations. Note that the core energies contribute to the second sub-leading coefficient. For short-range potentials, such as γ > 2, the ground energy is extensive, and the leading term varies as M 1+γ/2 . The extensive nature of the M 1+γ/2 term becomes clear upon noting that

where a is the particle spacing. Comparison with Eq. (74) shows that the dimension of Young's modulus Y arises solely from the lattice constant a and the electric charge e, consistent with elastic constants arising from physics on the scale of the lattice constant in an essentially flat geometry. This observation is now generalized to the rest of the couplings discussed in the previous section. For a fluid membrane not on a frozen topography, Helfrich terms arising from the mean curvature H as well as the Gaussian curvature can be important. Their scaling behavior is given by

Both terms would therefore contribute to the same order in the M expansion, although the last term is purely topological. For crystals embedded in a frozen topography we expect an expansion along the lines of Eq. (79),

The nonextensive term a 0 M 2 arises from the long range interactions. The next extensive contribution comes from the interaction between Gaussian curvature and defects as well as the extensive energy per particle in flat space. Bending rigidity contributions are higher order in 1/M and can be absorbed into a redefinition of the disclination core energy. Core energies also depend on non-universal details of the short-distance physics.

The results presented so far are strictly for systems at zero temperature. In systems with short range interactions, the elastic constants can be strongly temperature-dependent. An extreme example is hard disks of radius a 0 , which may be viewed as a limiting case of a power law potential of the form V (r) 0 (a 0 /r) γ , with γ → ∞. In this case, the elastic constants are strictly proportional to temperature. It is straightforward, however, to adapt the techniques presented here to the simpler problem of short range pair potentials.

The configuration on a sphere with the minimum number of charge ±1 defects is twelve +1 (5−valent) disclinations, which minimize their energy by sitting at the vertices of an icosahedron Y. The energies of such configurations will be computed for the discrete spherical tessellations described in Sec. 3.2.3 and compared with the predictions of continuum elastic theory, as illustrated in Fig. 13 . It is well established that for sufficiently large values of M configurations with more than 12 disclinations (i.e., those with "grain boundary scars") have lower energies [73, 111] . 

It is of interest, however, to study simple icosahedral configurations for large M , as metastable states with a well defined energy.

Within the continuum elastic theory it can be shown that twelve disclinations at the vertices of an icosahedron minimize the energy [73] when no further defects are allowed. The C-coefficient of Eq. (79) for this configuration of defects has been computed in [73] 

Y here stands for a particle configuration with 12 defects at the vertices of an icosahedron. Using the energy of Eq. (79), the coefficient a 1 (γ, Y) appearing in the expansion of Eq. (83) may be computed, with the results shown in Fig. 14 and Table 2 . From the results described in Sec. 3.2.3, the a 1 coefficient may be extrapolated to very large numbers of particles using the expansion derived from Eq. (83) . Indeed, as shown in Fig. 15 , plots of

vs 1/M are linear, with a slope that determines a 1 (γ) and an intercept related to the higher order core energy-like contribution. The results of these extrapolations are shown in Table 2 . The agreement between the continuum elastic theory and the explicit computation for the (n, n) configuration is remarkable, holding to almost five significant figures. For the (n, 0) lattice there is agreement to four significant figures. This agreement is even more striking when it is recalled that the a 1 -coefficient is obtained after subtraction of the term a 0 (γ)M 2 , as illustrated in Fig. 11 . Furthermore, in the range from γ = 0.125 to γ = 1.875, all the significant digits vary and yet the accuracy of the calculation is virtually independent of γ.

The a 1 -coefficient computed within the continuum elastic approach above does not depend on the icosadeltahedral class (n, m). Results from the direct minimization of particles do, however, show a weak dependence (in the 4th significant digit) on the particular (n, m) configuration, as is apparent from Fig. 11 and Table 2 . It should be noted that the discrepancy from the continuum result has a well defined sign, and is therefore reasonably attributed to a term not present in the energy expansion. 

When the number of particles is extremely large, the minimum energy configurations can be approximated by a closed analytical form, upon assuming a continuous distribution of defects. The formal elimination of the geometric frustration introduced by the Gaussian curvature may be formulated as a concrete set of equations in the case of the sphere.

We shall use the identity

which follows from the topological constraint Eq. (71) . Provided a disclination configuration exists such that

for each (l ≥ 1, m), the disclination density completely screens the Gaussian curvature. A configuration of defects satisfying Eq. (87) is an absolute minimum of the elastic energy, a result easily understood by writing the energy in the form

where the zero point energy E 0 in Eq. (62) is kept. A configuration satisfying Eq. (87)) will be denoted by G. For this hypothetical configuration, the C-coefficient in Eq. 79 vanishes, although there is now a large contribution (linear in R) from the dislocation core energies represented by the last term of Eq. (88). The G configuration may be characterized more explicitly. It consists of a density of dislocations that converges to the local Gaussian curvature. It can be shown that upon approximating the dislocations (each regarded as a disclination dipole with spacing a) as a continuum distribution, this dislocation density for a sphere

The summation here runs over the six coordinates of the northern hemisphere of an icosahedron, (0, 0) and (θ Y , 2πk/5), where θ Y = arccos(1/ √ 5) and α k is the angle θ relative to a coordinate system with the north pole located at (θ Y , 2πk/5) for k = 1, · · · 5. This angle is given implicitly by

The implicit form of Eq. (89) can be further simplified

Close to one of the 12 disclinations with charge s = 2π/6 Eq. (89) predicts a singularity in the dislocation density [74] b ≈ s 2πRa

.

For small angles, close to each disclination, there is a short-distance singularity

in agreement with known results in flat space. Eq. (89) represents a continuous distribution of dislocations, and neglects both dislocation discreteness and their mutual interactions. It represents six families of dislocations with azimuthal Burgers vectors associated with antipodal pairs of the 12 original disclinations in the icosahedron. When discreteness and interactions are taken in account, we expect these dislocations to condense into grain boundary arms, containing with quantized Burgers vectors and variable spacing in the radial direction [73, 112, 113] . The total number of discrete arms remains, therefore, the variable that needs to be determined for a discrete solution of the Thomson problem.

Within the continuum elastic approach, the dominant configurations for a small number of particles are 12 defects with an icosahedral symmetry [73] . We have just seen, however, that adding a continuous distribution of dislocations, as might be appropriate when the particle number is large, can more efficiently screen the Gaussian curvature on a sphere. The natural problem then becomes to determine the precise structure of the defect arrays for intermediate numbers of particles when the discreteness of interacting dislocations is taken into account.

We note first that the particular arrangement of defects dominating in this regime will not be fully universal. The particular array structure favored can vary from system to system with fixed particle number, depending, e.g., on details such as the dislocation core energy. This result may be traced back to the M -expansion of Eq. (83), in which the sub-leading terms which depend on non-universal properties influence the dominant terms for finite values of M . Some typical defect configurations obtained by direct minimization of particles on the sphere are shown in Fig.  16 and show incipient scars, already at 500 particles (in [73] the minimum number of particles where scars are systematically found is predicted to be around 400). By using the geometrical model described here, where the energy is parameterized just by a Young's modulus and a dislocation core energy [73, 92] one can simulate larger particle numbers and one obtains results as in Fig. 17 . Note the occurrence of low energy configurations with scars (m = 2) in one instance and pentagonal buttons (m = 5) in another. The dislocation spacing decreases the further a dislocation is from the central disclination.

An overview of results involving grain boundary scars is presented in Fig. 18 . If a disclination is placed on a perfect crystal, no additional defects will appear if the disclination is located on the tip of a cone with total Gaussian curvature equal to the disclination charge. If a disclination is forced into a flat monolayer, then m low angle grain boundaries, with constant spacing between dislocations as shown in Fig. 18 and grains going all the way to the boundary, will be favored (see [112] for a detailed discussion). In the intermediate situation where a finite Gaussian curvature is spread over a finite area, as in the case of a spherical cap, a disclination arises at the center of the cap and finite length grain boundaries stretched out over an area of (π/3)R 2 with variable spacing dominate, again as illustrated in Fig. 18 .

Additional results may be obtained for the number of arms within the grain boundary, the actual variable spacing between dislocations within the grain and the length of the grains as a function of the number of particles.

When grain boundary scars appear, one can estimate the number of excess dislocations which decorate each of the 12 curvature-induced disclinations on the sphere using ideas from Ref. [73] . This estimate is in reasonable agreement with experiments probing equilibrated assemblies of polystyrene beads on water droplets [92] . Consider the region surrounding one of the 12 excess disclinations, with charge s = 2π/6, centered on the north pole. As discussed in Ref. [73] , one expects the stresses and strains at a fixed geodesic distance r from the pole on a sphere of radius R to be controlled by an effective disclination charge

Here the Gaussian curvature is K = 1/R 2 and the metric tensor associated with spherical polar coordinates (r, φ), with distance element ds 2 = d 2 r + R 2 sin 2 (r/R)d 2 φ, gives √ g = R sin(r/R). Suppose m grain boundaries radiate from the disclination at the north pole. Then, in an approximation which neglects interactions between the individual arms, the spacing between the dislocations in these grains is [73] 

which implies an effective dislocation density

This density vanishes when r → r c , where

which is the distance at which the m grain boundaries terminate. The total number of dislocations residing within this radius is thus

As discussed in Ref. [114] , it is also of interest to consider 2π disclination defects (appropriate to crystals of tilted molecules [115] ) on the sphere. The icosahedral configuration of 12 s = 2π/6 disclinations is now replaced by just two s = 2π disclinations at the north and south poles. Using the approximation discussed above, it is straightforward to show that the density of dislocations in each of m (noninteracting) grain boundary arms now reads

This density vanishes at r c = (π/2)R, corresponding to a hemisphere of area on the sphere for each cluster of arms. It is of considerable interest to repeat the above calculation for a square lattice, as found for example in the protein surface layers (s-layers) of some bacteria [116, 117] . In this case the basic disclination has s = 2π/4. The effective dislocation density becomes

This density vanishes when r → r c , where

which is the distance at which the m grain boundaries terminate. The longer angular length of square lattice scars reflects the larger initial disclination charge (90 • ) that must be screened. The total number of dislocations residing within this radius is thus

Thus the angular length of scars and the total number of excess dislocations is a measure of the underlying topology of the lattice tiling the sphere.

Interstitial fractionalization on the sphere

Consistent with theoretical expectations, experiments on particle-coated water droplets in oil [92] reveal that the twelve excess disclinations sprout grain boundary "scars" for sufficiently large R/a. When triangulations of microscopic particle packings are used to reveal the local coordination number, these grain boundaries appear as additional dislocations, i.e., 5-7 pairs, arrayed around an unpaired 5, in a pattern such as 5-7-5-7-5-7-5-7-5. Although the critical value of R/a above which grain boundaries appear depends on microscopic details, both theoretical estimates and experiments indicate that this instability arises as soon as R/a 5 − 6, i.e., when the total number of particles exceeds several hundred. Thus, unlike crystals in flat space, dislocations arrayed in grain boundaries are an intrinsic part of the ground state. These grain boundaries can, moreover, stop and start freely on the sphere, unlike their flat space counterparts. Such terminations occur naturally (and with low energy cost) because crystalline grains rotate under parallel transport due to the nonzero Gaussian curvature of the sphere.

If disclinations and dislocations are crucial for understanding spherical crystallography, what can we say about vacancies and interstitials, which are well known to play a key role in conventional crystals [118, 119] . It is natural to introduce vacancies and interstitials in an attempt to understand spherical particle packings that deviate from the sequence of "magic numbers" N nm = 10(n 2 + nm + m 2 ) + 2. It is tempting to ignore the instability to grain boundary scars for large N nm and regard the commensurate (n, m) lattice as an interesting metastable state. It would then be natural to introduce vacancies and interstitials to describe candidate ground state packings for N nm ± t particles, where t = 1, 2, ... and much less than the distance to the next magic number. There is, however, another surprise in store: in contrast to flat space, where vacancy and interstitial defects are stable and well defined, one finds that interstitials and vacancies are typically ripped apart into dislocation fragments by the strain fields of nearby 5-fold disclinations. The dislocations then combine with some of the excess 5's to form defect clusters such 5-7-5. Thus, vacancies and interstitials lose their integrity via fragmentation in spherical crystals, and mediate formation of small grain boundary scars. Figure 19 . (Color online) Initial icosahedral configuration of an (8,3) spherical crystal lattice with 12 disclination defects. This lattice is chiral, and is arranged such that 8 steps along a Bragg row, followed by 3 steps to the right along a Bragg row at a 120 • angle, connect neighboring 5-fold disclinations. Distinct initial locations for interstitials are shown as triangular plaquettes labeled by characters and numbers.

We are interested in studying the effect of inserting an interstitial or a vacancy into a regular (n, m) icosahedral lattice by adding or removing a single particle, giving rise to particle numbers N nm ±1 not falling in the classification of icosadeltahedral lattices and in studying the relation of interstitials to the extra dislocation defects (scars) found in spherical crystals above a critical system size [73, 92] . Refs. [120, 121, 122, 123, 124, 125, 126, 127] study configurations and energies of interstitial and vacancy defects and their energetics in triangular lattices in flat space.

The presence of excess dislocation defects in the ground state of spherical crystals is dramatically illustrated by the following numerical experiment: we start with a regular icosadeltahedral tessellation of the sphere -say an (8, 3) , corresponding to N 83 = 972 (Fig. 19 ). This may be done with the applet located at [109] using the Construct (m, n) algorithm. Although the true ground state for 972 particles on the sphere with most pair potentials has additional dislocation defects (i.e. tightly bound pairs of 5-and 7-coordinated particles) arrayed in grain boundary scars [73] , the regular icosadeltahedral lattice is a local minimum from which it is difficult to escape without the addition of thermal noise. In fact it is a major challenge to find fast and reliable algorithms to locate the true ground state (global minimum) in this problem with its complex energy landscape. Now add a single particle to the lattice at the center of mass of a spherical triangle whose vertices are 3 nearest-neighbor 5-fold disclinations (using commands shift + click). The self-interstitial so formed is then relaxed by a standard relaxation algorithm, with sufficient thermal noise to allow dislocation glide over the Peierls potential [93] . One immediately finds that an interstitial is structurally unstable. In a few time steps it morphs into a complex of dislocations with zero net Burgers vector. The most common structure observed is a set of three dislocations, with Burgers vectors perpendicular to a line joining each 5-7 pair, inclined at 120 • angles to each other. Eventually the interstitial complex is ripped apart entirely, as illustrated schematically in Fig. 20 (see also Fig. 21 ). Intermediate configurations and final states as the dislocations glide apart will be classified later. Most often three separate dislocations are formed which each glide toward a 5-fold disclination. The end result is the formation of a "mini-scar" (a 5 − 7 − 5 grain boundary) at each of the vertex 5s. Subsequent removal of a particle to restore the particle number to the original 972 and relaxing still leaves scars with total energy lower than the starting configuration with 12 isolated 5's. This observation confirms that scars are definitely lower energy states and not simply artifacts of the relaxation algorithm.

The above phenomenon of low-temperature (T 0) unbinding of dislocations by spatial curvature is a curved space analog of melting at finite temperature. The extended nature of fractionated interstitials (each separating dislocation component involves an extra row of particles) means that they cannot be treated as small perturbations from the initial spherical crystal with particle number N nm .

Let's return to the specific case of the N = 973 particle configuration generated by an interstitial inserted into one triangular plaquette of a regular (8, 3) lattice of N nm = 972 particles with the requisite 12 disclination defects (5-fold coordinated particles) at the vertices of a regular icosahedron, as shown in Fig. 19 . The spherical crystal is distorted by the additional particle -the local configuration adopted by the interstitial changes as the crystal relaxes toward a lower energy state. As in the case of planar lattices, one also finds here that various interstitial defect configurations appear, such as the twofold symmetric interstitial I 2 , the threefold symmetric interstitial I 3 , and the fourfold symmetric interstitial I 4 (see Figs. 21, 22, and 23). The most common intermediate complex formed by the interstitial is threefold symmetric in the rough form of a triangular loop composed of three dislocations with radially oriented Burgers vectors. All of the configurations adopted by an interstitial prior to unbinding can be described as a set of dislocations with zero net Burgers vector.

In marked contrast to interstitial defects in a planar crystal, interstitial defect configurations in a spherical crystal are metastable states with characteristic decay processes. As we shall see, the instability is caused by interactions with the inevitable disclinations associated with the nonzero Gaussian curvature of the sphere. A representative evolution of an interstitial inserted at the center of three neighboring disclinations in an (8, 3) icosadeltahedron is shown in Fig. 21 . After some local relaxation the interstitial configuration denoted I 3 is formed, as shown in Fig. 21 (a) . We also show there the construction of a triangular contour surrounding the original interstitial defect. The presence of an interstitial follows because the contour encloses 7 particles instead of 6, the number appropriate to a perfect triangular lattice. During the annealed relaxation illustrated in Figs. 21(b) through (d), the dislocations which are bound together in the initial interstitial unbind into individual dislocations and subsequently glide towards nearby disclinations, eventually forming minimal 5-7-5 grain boundary scars. If one thinks of the initial dislocations as internal degrees of freedom within the interstitial, one could say that one-third of an interstitial is present in each mini-scar in the final state and thus that the interstitial demonstrates 1/3 fractionalization. For other initial conditions the fractionation is into two dislocations, each representing 1/2 of the original interstitial. The instability of interstitial defects in curved crystals may be studied via continuum elasticity theory by calculating the interaction energies between defects at each stage of the relaxation process of Fig. 21 [119] . We note that the triangular plaquette around the initial defect has been deformed such that it conforms as closely as possible to the regular triangular lattice during the relaxation process. Its deformation reveals the passage of the escaping dislocations.

In this section, we discuss the energy of an interstitial defect in a spherical crystal. Consider N point particles constrained to lie on the two-dimensional surface of a unit sphere. The energy of N particles interacting through a generalized Coulomb potential within the curved surface is given by

where x is the position of the particle in three dimensions and γ is an integer. For a flat triangular lattice with periodic boundary conditions, the interstitial defect energy at constant density was defined in [122, 123] as

where E relaxed is the relaxed energy of the rearranged lattice of N particles with the interstitial defect in the area A, and E per is the energy of the perfect crystal at the same areal density N/A. In curved space, however, the definition of a "perfect crystal" is more subtle, since disclination defects resulting from the Gaussian curvature and the topology are inevitable. We will take as a reference crystal the (n, m) icosadeltahedral configurations corresponding to triangular tessellations of a magic number of particles N nm = 10(n 2 + m 2 + nm) + 2. Once an interstitial or vacancy is added to such a (n, m) configuration, we are no longer at a magic number of particles since these are quite sparsely distributed. We thus need to define the energy of the perfect crystal. Here we define the energy of the interstitial (vacancy) defect at constant density in the spherical crystal as

where E local measures the energy of the relaxed interstitial while the constituent dislocations are still bound and E * annealed is the minimum energy of all possible final states attained after annealed relaxation leading to interstitial fractionalization. This definition will be more explicitly discussed in the following section (see Table 5 ). We note that both E local and E * annealed are measured at the same areal density (N nm ±1)/A, where ±1 correspond to an interstitial (vacancy) respectively. The lowest relaxed energy E * annealed plays the role of the energy of the perfect lattice in the planar case at the density of (N nm ± 1)/A.

Numerical measurements of E local and E * annealed for the power-law potentials with γ = 1, 3, 6 and 12 have been performed by adding one interstitial at the center of a spherical triangle formed by three nearest-neighbor disclinations in the (8, 3) lattice (the location represented by C in Fig. 19 ) [119] . E local is measured by quenching the system at the moment just prior to the fractionation of the interstitial into individual dislocations [(a) in Fig. 21 ]. The results are reported in Table 4 . Table 4 . The lowest local and annealed relaxed energy with the central interstitials created by putting a particle at C in Fig. 19 of the (8, 3) lattice, for γ = 1, 3, 6, and 12. The differences between two relaxed energies are calculated. Because the particles are embedded in a sphere of unit radius with our conventions, near-neighbor particle spacings are of order a ∼ N −1/2 1, leading to a strong γ dependence in the total energy given by Eq. (103) 

By adding a particle in different plaquettes within the large spherical triangle of Fig. 19 , one can explore the dependence of the final state on the initial interstitial location. In contrast to the case for planar crystals, both the location and orientation of the interstitial defect relative to nearby disclinations influences the resultant configuration and its corresponding evolution, leading to distinct relaxed configurations. The insertion of a particle at the center of the large spherical triangle leads to an I 3 -type initial configuration, whereas adding the interstitial to the plaquette along an edge results in an I 2 -type initial configuration.

During the relaxation process, one also finds that the dislocation complex representing an interstitial can rotate so as to reorient the Burgers vectors so that constituent dislocations can glide towards a nearby disclination and bind to it. This phenomenon is especially noticeable if one places one extra particle slightly off the absolute center C, such as at the locations C or 0 in Fig. 19 .

In Figs. 22 and Fig. 23 this phenomenon is illustrated explicitly. In Fig. 22(a) , an interstitial initially placed at the position C in Fig. 19 morphs quickly to an I 3 configuration. In (a), however, the orientations of the dislocations within I 3 are not appropriate for dislocation glide to the surrounding disclinations since the 5 end of the dislocations point towards rather than away from these 5-fold disclinations, Table 5 . The three classes of final annealed state depending on the initial interstitial location. The relaxed energies are measured for the power law potential γ = 6. Here E annealed with 3 scars corresponds to E * annealed .

Initial Location E annealed Annealed State C, C , 7, 8, 10 9.60570 ×10 8 3 scars (5-7-5) 0, 2, 3, 4, 6, 9 9 .61011 ×10 8 2 scars (5-7-5) 1, 5 9.61062 ×10 8 2 scars (5-7-5-7-5) causing them to be repelled. Glide-induced fractionalization is therefore prohibited in this orientation. Remarkably, though, the entire complex of dislocations can change its orientation by a transition through an intermediate I 2 configuration (shown in Fig. 22(b) ) and subsequently to a second I 3 configuration [ Fig. 22(c) ]. The final I 3 configuration is rotated by 60 • with respect to the first I 3 and can now fractionate analogously to an interstitial with initial position C in Fig. 19 .

One also find rotational reorientation of an interstitial defect with an I 4 -type intermediate state, as shown in Fig. 23 . This particular relaxation process reveals an interesting feature of dislocation dynamics on a curved surface. The I 3 configuration generated after the intermediate [see Fig. 23 (c)] now has one dislocation with its glide plane such that it can glide head-on into a vertex disclination. This disclination absorbs the dislocation but hops over one lattice spacing to accommodate the curved space Burgers vector. The other two dislocations end up bound in the form of mini-scars. In this case then the interstitial has fractionated into 2 rather than 3 parts and one say that there has been 1/2 fractionization of the interstitial. Absorption and emission of dislocations by 5-fold disclinations are somewhat analogous to absorption and emission of vacancies and interstitials by dislocations (allowing dislocations to climb), a phenomenon well-known in flat space.

What is the dependence of the final state on the initial location of the interstitial? The distinct initial conditions shown in Fig. 19 lead to three different final annealed states, as summarized in Table 5 . The final state with three mini-scars of the form 5-7-5 has the lowest energy of all final states and provides a measure of E * annealed . It is also informative to track the position dependence of the interstitial energy after it relaxes. Numerical measurements of local relaxed energy for the interstitial as a function of initial location are presented in Table 6 . In most cases the initial I 3 complexion undergoes transitions to more stable interstitial configurations, except the configuration that starts from the very center location C in Fig. 19 . For the initial conditions, C, C , 0, and 3, I 3 is the most stable. Note that the interstitial created at the exact center C has the lowest defect energy, while the one nearest to the disclination from the location 10 requires the largest energy for interstitial defect formation. Fig. 19 ) mediated by the transition: I 3 → I 4 → I 3 . Table 6 . The energy of interstitial defects created at different initial positions within the spherical crystal. The energies shown are for a power law potential with γ = 6.

n

Transition 

The case of nematic order on the sphere is also of great interest from both a fundamental and materials science viewpoint [60, 128] . The free energy for nematogens on the sphere is given by

where β ij is the geodesic angle between defects i and j with winding numbers n i and n j , and E(R) is the defect self-energy. Nematic spheres might be made by coating a droplet with gemini lipids, ABA triblock copolymers or nanorods. Microfluidic techniques for creating a thin spherical shell of liquid crystal in a double emulsion have been explored in [129, 130] . Although the tetrahedral array of four type-1/2 nematic disclinations might be unstable to Boojum formation (type-1 defects) for spheres immersed in a nematic bulk fluid [131] , Huber and Stark [132] have shown that nematic wetting for spheres immersed in an isotropic fluid is a promising alternative route to a spherical nematic. An important motivation for exploring nematic order on spherical topologies is the search for superatoms of low valency. We have seen that spherical crystals have 12 distinguished regions, each region being marked by the presence of either an isolated 5-disclination or an extended defect array like a scar. By functionalizing these regions one could hope to engineer 12-valent supermolecules. In particular one might be able to take advantage of the unusual presence of 5-7 pairs (dislocations) in scars to convert scars to sticky zones of the spherical superatom. Twelve is very likely too high a coordination number, however, to produce viable molecules or bulk materials. How can we change the coordination number of spherical superatoms? The solution lies in changing the local order on the sphere. An ordered state on the sphere with p-fold symmetry will have 2p net topological defects by Eq. (3). The crystalline case discussed so far corresponds to p = 6. Vector-like order on the sphere corresponds to p = 1 and the defects correspond to index-one vorticesone can also visualize the source and sink of fluid flow on a sphere [see Fig. 24(a) ] or imagine the base and crown bald spots of hair combed on a sphere. Thus a 2-sphere decorated with magnetic (or any other vector) field lines is a candidate for a valence-2 superatom. Consider now nematic order: p = 2. In this case the vortex configurations above each break up into 2 type-1/2 nematic disclinations, leading to a total of 4 defects located at the corners of a tetrahedron (in the one Frank constant approximation) [59] (see Fig. 24(b) ). This raises the possibility of making tetravalent superatoms with sp 3 type bonding [128] .

The precise geometry of the defect structure is important for determining what kind of chemical linkers could be attached to the defects to facilitate the creation of large scale structures like a diamond-type lattice of linked tetravalent colloids. Diamond lattices are desirable for photonic crystals (optical analogs of semiconductors) since they are known to possess a large photonic band gap [133] .

Although the elementary type-1/2 disclinations favor a tetrahedral geometry if splay and bend deformations are degenerate in energy this changes when elastic anisotropy is present. In the limit of pure splay or pure bend the energy of a +1 defect becomes degenerate with two 1/2 defects. Thus two +1 defects sitting at the north and south pole can be separated be cutting the sphere on a great circle joining the two defects and rotating by an arbitrary angle α. This gives rise to a set of four 1/2 disclinations lying on the same great circle and therefore coplanar in the embedding space [134] . Thus the tetrahedral configuration crosses over to a great-circle configuration as one increases the anisotropy of the elastic constants. This explicitly demonstrates that defect positions can be controlled by varying the elastic anisotropy. Thus the structure of directional bonding for the associated supermolecules is controllable. Under specific experimental conditions amphiphilic molecules in solution, such as lipids or amphiphilic block copolymers, self-assemble in a spectacular variety of shapes including spherical and cylindrical micelles, vesicles and lamellae, together with more complex geometries such as uni-and multilamellar vesicles, onion vesicles, toroidal and cage-shaped micelles. The exact shape and size of these structures has been observed to depend on both molecular details, such as molecular size, the hydrophilic/hydrophobic ratio and molecular stiffness, and collective parameters, such as the concentration or the diffusivity of the molecules. As noted previously, amphiphilic membranes can exist in different thermodynamic states depending on the degree of orientational and positional order and their precise molecular constituents. The L β and P β phases, which occur in lipid membranes featuring the phosphatidylcholine (PC) group, have been found, in particular, to exhibit in-plane orientational correlations extending over 200Å, one order of magnitude larger than the typical spacing between PC groups [25, 135] . Because of thermal fluctuations, such membranes are typically corrugated, possessing spatially inhomogeneous Gaussian and mean curvature. Furthermore it is likely for lipid membranes to have pores providing a passage from the exterior to the interior of a vesicle. It is therefore natural to ask how variable Gaussian curvature and the presence of one or more boundaries affects the phenomenology reviewed thus far for the sphere.

In this section we discuss crystalline order for two important surfaces with variable Gaussian curvature and (possibly) boundary, namely the bumpy surface that is obtained by revolving the graph of a Gaussian function around its symmetry axis (i.e. a Gaussian bump) and the paraboloid of revolution. The former can be thought as a gentle deformation of a plane and thus can serve as a playground to analyze the onset of behavior not found in planar systems; the latter is possibly the simplest two-dimensional Riemannian surface having variable Gaussian curvature and boundary and provides a setting that is simple enough to carry out a full analytical treatment and analyze also large curvature regimes. Furthermore, since paraboloidal shapes naturally occur across the air/liquid interface of a fluid placed in a rotating cylindrical vessel, a direct physical realization of paraboloidal crystals can be constructed by assembling monodisperse objects on the surface of a rotating liquid. In Sec. 4.5 we review a simple experiment by Bowick et al [136] in which such a macroscopic model for a paraboloidal crystal is constructed by assembling a two-dimensional soap bubble "raft" on the air/liquid interface of a water-soap solution, thus extending the classic work of Bragg and Nye on planar bubble rafts.

Purely orientational order on the Gaussian bump has been recently reviewed by Turner et al [64] and won't be discussed here.

Before analyzing the ground state structure of a crystalline paraboloid and Gaussian bump, it is useful to review the geometry of surfaces of revolution. The notion of conformal mapping of Riemannian surfaces will also be used in the following and will be briefly reviewed in this section with special attention to its application to the calculation of Green functions on simply connected 2−manifolds.

A surface of revolution M is a surface obtained by revolving a two-dimensional curve around an axis. The resulting surface therefore always has azimuthal symmetry. The standard parametrization of a surface of revolution is:

where r ∈ [0, R] (with R possibly infinite) and φ = [0, 2π). The metric of the surface (107) is given by:

where the prime indicate a partial derivative with respect to r. The Gaussian curvature is a function of r only and is given by:

In the presence of a boundary ∂M the condition (2) for the total topological charge of any triangulation on M reads:

where V ∂M and V M are the number of vertices on the boundary and the interior of the manifold respectively (with V = V ∂M + V M the total number of vertices). q i,∂M = 4−c i is the topological charge of a vertex of coordination number c i located on the boundary, where the coordination number of a perfect triangular lattice is four rather than six. A surface of revolution with boundary ∂M = {r = R}×[0, 2π] is homeomorphic to a disk (g = 0 and h = 1) and has therefore χ = 1 and total topological charge Q = 6. This topological constraint can be satisfied, for instance, by placing six isolated 3−fold disclinations along the boundary and keeping the interior of the surface defect-free or by placing a 5−fold disclination in the interior and the remaining five 3−fold disclinations along the boundary. Let us consider now a generic Riemannian surface M and two curves on S intersecting at some point x 0 . The angle between the two intersecting curves is, by definition, the angle between the tangents to these curves at x 0 . A mapping of a portion S of a surface onto a portion S * is called conformal (or angle-preserving) if the angle of intersection of every arbitrary pair of intersecting arcs on S * is the same as that of the corresponding inverse images on S at the corresponding point (see for example [137] ). It is not difficult to prove that a mapping from a portion S of a surface onto a portion S * is conformal if and only if, when on S and S * the same coordinate systems have been introduced, the coefficients g * ij and g ij of the metric tensor of S * and S are related by:

with w a positive function of the coordinates x = (x 1 , x 2 ). Eq. (111) implies indeed that the angle between any pair of intersecting curves is the same in S * and S.

Isometries are a special case of conformal mappings where w = 1 and the mapping is both distance and angle-preserving. A special type of conformal mapping is that of a portion S of a surface into a plane. This may be accomplished by introducing a set of coordinates u = (u 1 , u 2 ) such that:

Coordinates (u 1 , u 2 ) satisfying Eq. (112) are called isothermal (or conformal). In general, any simply connected Riemannian manifold with a C ∞ −smooth metric ds 2 can be equipped with a set of local isothermal coordinates. This important result can be stated by saying that any simply connected Riemannian manifold is locally conformally equivalent to a planar domain in two-dimensions. In conformal coordinates, the Gaussian curvature reads:

Conformal mapping is the fundamental tool behind the celebrated uniformization theorem according to which, every simply connected Riemannian surface is conformally equivalent to the unit disk, the complex plane or the Riemann sphere. This theorem, first proved by Koebe and Poincaré independently in 1907, extends the Riemann mapping theorem for simply connected domains in the complex plane to all simply connected Riemannian surfaces and provides an insightful classification scheme. Its formidable power lies in the fact that the mapping that allows one to transform a generic surface into a simpler "irreducible" one is not an arbitrary homeomorphism but is conformal, and thus preserves part of the geometrical structure of the original manifold. In the following we will see how this feature has important consequences in the elastic theory of defects on curved surfaces. The uniformization of Riemannian surfaces is historically the first example of geometrization. In the case of manifolds of higher Hausdorff dimension, and 3-manifolds in particular, the latter program has become, following Thurston, one of the most challenging and fascinating chapters of modern geometry. A bounded Gaussian bump and a paraboloid of revolution are both conformally equivalent to the unit disk D of the complex plane. Calling z = e iφ , the new metric will be:

The conformal factor w can be found by equating the metrics (108) and (114) . This yields:

with and r related by the differential equation

whose solution is given by:

The sign of the exponent and the integration constant in Eq. (117) can be tuned to obtain the desired scale and direction of the conformal map. The calculation of the Green function of the Laplace and biharmonic operator is considerably simplified once a surface of revolution has been endowed with a local system of isothermal coordinates. In this case it is easy to show the Laplace-Beltrami operator ∆ g takes the form:

where ∆ is now the Laplacian in the Euclidean metric tensor:

with determinant γ. Since the determinant of the metric tensor is rescaled ( √ g → w √ γ) under conformal mapping, the Laplace and biharmonic equation for the Green function become:

where δ(z, ζ) is the standard delta function at the point z = ζ of the unit disk. Eq. (120a) is now the standard Laplace-Green equation. Its associated Dirichlet problem has the familiar solution:

Eq. (120b) is known as the weighted biharmonic Green equation. The uniqueness of its solution requires imposing both Dirichlet and Neumann boundary conditions:

where ∂ n(z) denotes the derivative with respect to the variable z along the normal direction at ∂D. Its solution can be expressed in integral form as

where H(σ, ζ) is a harmonic kernel that enforces the Neumann condition. Such a function depends on the form of the conformal weight w. For radial weights w = w 0 (|z| 2 ), such as those obtained by conformally mapping a surface of revolution, H(σ, ζ) has been calculated explicitly by Shimorin [138] :

where:

and the coefficients c n are given by:

The most natural way of introducing a non-vanishing Gaussian curvature on an initially flat medium is to gently deform the medium at one point in such a way that the curvature introduced by this deformation dies off at infinity. If the initial planar domain is the entire Euclidean plane R 2 , the Euler characteristic is zero and, independently of the value of the Gaussian curvature, an ordered phase embedded on it is not topologically required to contain defects. Such a construction is clearly ideal to detect the onset of structural behavior not occurring in the ground state of a planar system such as the appearance of dislocations and disclinations. Vitelli et al [139, 140] analyzed crystalline and p−atic order on the bumpy surface obtained by revolving the graph of a Gaussian function about its symmetry axis. The resulting Gaussian bump has parametrization:

with r ∈ [0, R] and φ ∈ [0, 2π). In this parametrization the metric tensor g ij (with determinant g) and the Gaussian curvature are given by:

where α = h/r 0 is the aspect ratio of the bump and (r) is given by

It is instructive to verify that the Euler characteristic χ vanishes when the boundary radius R is set to infinity. Employing the Gauss-Bonnet theorem one has

where κ g , the geodesic curvature of the circular boundary C R of radius R, is given by

If the bump is unbounded the second term in Eq. (130) disappears and, in the limit R → ∞, one has 130) gives:

so that the terms proportional to 1/ (R) cancel each other, yielding χ = 1. Because the infinite bump has χ = 0 disclinations must appear in pairs and dislocations must have total Burgers vector zero:

Vitelli et al showed that topological defects appear in the ground state when the aspect ratio α exceeds a critical value α c . Because of the topological constraint, defects appear initially in the form of a pair of unbound dislocations, roughly located in the region where K = 0. Upon increasing the aspect ratio, more dislocations appear. This mechanism clearly resembles the defect proliferation that occurs in two-dimensional melting and suggests an interpretation of the curvature as a local effective temperature. At the onset of defect proliferation, inter-defect interactions are negligible compared to the interactions with the "smeared out" topological charge associated with the curvature of the underlying medium. The dislocation unbinding mechanism is then described by a non-local function of the Gaussian curvature representing the defect-curvature interaction part of the elastic energy of Eq. (54):

where F int is the curvature-defect interaction part of the elastic energy and ϕ(x) can be interpreted as a geometric potential associated with the Gaussian curvature of the embedding manifold: Thus

where U (y) is a harmonic function which enforces the boundary conditions. The Laplacian Green function has the form (121) with the conformal distance ρ = |z| given here by

For a single dislocation of Burgers vector b, F int has the form [139] :

where λ = r/r 0 and Λ = R/r 0 . The first term in Eq. (136) corresponds to the R → ∞ geometric potential, while the second is a finite size correction arising from a circular boundary of radius R. A plot of the function (136) is shown in Fig. 26 . The profile of the function F int can be undestood by regarding a dislocation as bound pair of disclinations of opposite topological charge. Each disclination interacts with a potential of the form (133) . For small r, positive (negative) disclinations are attracted (repelled) by the center of the bump. As a consequence disclinations experience a force which increases linearly with r. Thus if the positive disclination in the dipole is closer to the top, it will experience a force that is opposite and slightly less than that acting on the negative disclination that is further from the top. As a result an effective "tidal" force will push the dislocation downhill. For large r, however, the geometric potential saturates and the attractive force exerted on the positive disclination takes over and drags the dislocation toward the center of the bump. The minimum of the elastic energy corresponding to the equilibrium between these two competing forces is obtained for λ ≈ 1.1. The origin of these forces is the Peach-Koehler force f k = kj b i σ ij acting on a dislocation of Burgers vector b i in an external stress field σ ij . In the case of a two-dimensional crystal on a substrate with preexisting Gaussian curvature, dislocations couple with the internal stress due to the curvature of the medium, with similar effects. Unbound dislocations occur in the ground state of a crystalline Gaussian bump when the aspect ratio α exceeds a critical value. For nearly flat landscapes the energetic cost of a pair of unbound dislocations is larger than that arising from the distortion of the medium and the system is favored to be defect free. Upon increasing the aspect ratio, however, the resulting elastic strain can be partially relieved by introducing a pair of dislocations with equal and opposite Burgers vector. The transition occurs when the energy gain from placing each dislocation in the minimum of the potential energy F int outweights the total work needed to tear them apart plus the core energies 2E d . The resulting critical aspect ratio is given by:

. The number of unbound dislocations increases with increasing aspect ratio and elementary 5−7 dislocations cluster in more complicated structures with zero net Burgers vector. Fig. 27 shows a plot of the critical aspect ratio α c as the function of the dimensionless parameter r 0 /b from Ref. [139] , as well as density plots of the strain energy corresponding to a defect-free bump (c) and a representative configuration featuring two unbound dislocations (b). The interaction between defects and curvature also has some remarkable stabilizing effects on defect dynamics. Dislocation dynamics consists of two distinct processes: glide and climb. Glide is motion along the direction of the Burgers vector. It requires only local rearrangement of atoms and thus has a very low activation energy, making it dominant at low temperatures. Climb consists of motion in the direction perpendicular to the Burgers vector. It requires diffusion of vacancies and interstitials and is usually suppressed therefore relative to glide. As a consequence of the underlying curvature, however, a gliding dislocation experiences a "recalling" force f recall ≈ k d |y|, where y is the transverse displacement and k d is a position-dependent effective spring constant. At leading order in y and α, the latter is given by:

This effective recalling force is not due to any external field nor to the interaction of dislocations with other defects, but exclusively to the coupling of the gliding dislocation with the curvature of the substrate. As expected, the effective spring constant (138) vanishes for planar crystals (i.e. α → 0). Since Y b 2 can be hundreds of k B T , at finite temperature the harmonic potential 1 2 k d y 2 associated with the recalling force raises the activation energy of thermally induced dislocation glide. The harmonic recalling potential and the effective spring constant (138) are plotted in Fig. 28 , together with numerical data obtained from a fixed connectivity harmonic model [139] .

Spatial curvature also provides an effective potential in the thermal diffusion of interstitials and vacancies. These can be constructed by grouping three dislocation dipoles. The elastic energy associated with a single interstitial/vacancy can be derived from Eq. (132) in the form

where V = ∆ϕ and Ω is the area excess or deficit associated with the defects. As a result, interstitials tend to climb to the top of the bump while vacancies are pushed into the flat regions. Like a disclination of positive topological charge, an interstitial is attracted to regions of positive Gaussian curvature while a vacancy is attracted to regions of negative Gaussian curvature. The crystalline structure arising on a Gaussian bump has been investigated numerically by Hexemer et al [141] . By means of smart Monte Carlo (SMC) simulations, the authors analyzed the ground state configuration of a system of N point-like particles (with N up to 780) interacting on a Gaussian bump of variable aspect ratio with a Yukawa potential of the form U (r) = exp(−κr)/r, with r the Euclidean distance between two particles in R 3 . This potential approximately describes the interaction between charged particles in solution with counter-ions, with κ −1 being the Deybe-Hückel screening length. Starting from an initially defect-free lattice and relaxing it with order 10 5 SMC iterations, Hexemer et al observed the appearance of defects for increasing aspect ratio. Fig. 29 shows the arrangement of 780 particles for three different aspect ratios. As predicted by the elastic theory, the proliferation of defects starts with the appearance of a pair of isolated dislocations ( Figs. 29a and b) . They are sitting at λ = 1.46 from the center and are rotated by 180 • , giving rise to a configuration with total Burgers vector close to zero. Upon increasing the aspect ratio, a third dislocation is observed (Fig. 29c  and d) . The three dislocations are rotated by 120 • with respect to each other. The green and blue circles in Fig. 29c and d mark the region of zero Gaussian curvature and the minimum of the geometric potential at 1.1 r 0 . As shown in the figure, for α = 0.95 two of the three dislocation are not sitting at 1.1 r 0 as one might expect. One is deep inside the area of positive Gaussian curvature while the other is in the region of negative Gaussian curvature. Hexemer et al attributed this discrepancy to finite-size effects. Fig. 29e and f show the lowest energy state on a bump of aspect ratio α = 1.58. In this situation dislocations appear arranged in the form of grain boundaries. This phenomenon has a simple interpretation. The total elastic energy of a collection of defects on a curved surface consists of three parts: the curvature-defect interaction discussed above, which tends to localize the defects where the geometric potential is minimal, a defect-defect interaction which is repulsive for like-sign defects and attractive for defects of opposite sign and an additional contribution due to the curvature alone. For a small number of dislocations, inter-defect interactions are negligible compared with curvature-defect interactions and dislocations tend to gather where F int is minimal. As shown from the data in Fig. 30 , however, the total number of dislocations in the bump is linearly proportional to the aspect ratio. When the total number of dislocations on the bump exceeds some threshold value the repulsive interaction takes over and dislocations start spreading. Dislocations with parallel (antiparallel) Burgers vector are organized in such a way as to maximize (minimize) their reciprocal distance while still taking advantage of the Gaussian curvature screening. The competition leads to the appearance of the Y −shaped grain boundaries shown in Fig. 29 . The branching allows the dislocations to lower their potential energy by keeping a large number of dislocations close to the 1.1 r 0 circle.

A paraboloid of revolution is a simple surface with both variable Gaussian curvature and boundary. Its standard parametrization is given by

with r ∈ [0, R] and φ = [0, 2π]. Here h is the height of the paraboloid and R the maximum radius. In the following we will call κ = 2h/R 2 the normal curvature of the paraboloid at the origin. The metric tensor g ij and the Gaussian curvature are given respectively by:

The problem of finding the optimal arrangement of disclinations in the ground state of paraboloidal crystals has been considered by the authors of this review article [142, 143] . As in the case of a Gaussian bump, when the maximal Gaussian curvature exceeds some critical value (depending on κ and R) defects proliferate in an initially defect-free configuration. Unlike the Gaussian bump, however, the Gaussian curvature on a paraboloid is strictly positive and only vanishes at infinity.

Positive isolated disclinations are thus energetically preferred to dislocation dipoles and the defect proliferation mechanism consists of a "migration" of one of the six topologically required +1-disclinations from the boundary to the origin of the paraboloid. Dislocations, on the other hand, appear in the system clustered in the form of grain boundary scars in the high density regime as in spherical crystals.

Since the Gaussian curvature is not constant throughout the manifold, however, this transition is preceded by a regime in which isolated disclinations and scars coexist in the crystal. Let ρ(x) be the effective topological charge density of a system of N disclinations on a background of Gaussian curvature K(x):

If free boundary conditions are choosen, the elastic energy (55) can be written as

where Γ(x) = ∆χ(x), and χ(x) satisfies the inhomogeneous biharmonic equation

with boundary conditions

where ν i is the ith component of the tangent vector ν perpendicular to the boundary. In the parametrization (140), the normal vector ν is simply given by g r /|g r |, with g r the basis vector associated with the coordinate r. The stress function Γ(x) can thus be expressed as

G L (x, y) is given in (121) and U (x) is a harmonic function on the paraboloid that enforces the Neumann boundary conditions in Eq. (145) . The conformal distance on the paraboloid is given by

Using Eqns. (147) and (121) in Eq. (143), the elastic energy F el of a collection of N disclinations in a paraboloidal crystal can be expressed as:

where the first term represents the bare contribution of the defects to the energy density and the second corresponds to the screening effect of the Gaussian curva- ture. Explicitly:

where r = |x| and

is a normalization constant depending on boundary radius R and the ratio κ. Fig.  31a shows a plot of the screening function Γ s (r) for different values of κ ∈ [1, 2]. As expected, the contribution due to Gaussian curvature is maximal at the origin of the paraboloid and drops to zero at the boundary. The calculation of the harmonic function U (x) requires a little more effort. If the crystal was defect-free (or populated by a perfectly isotropic distribution of defects) the function U (x) would be azimuthally symmetric and constant on the boundary. By the maximum principle of harmonic functions, U (x) would then be constant on the whole manifold and depend only on κ and the radius R: U (x) = U κ,R . This constant can be determined by integrating ∆χ(x) = Γ(x) and imposing the second boundary condition in Eq. (145) . This gives:

where A is the area of the paraboloid:

As shown in Figure 31b , the value of U κ,R quickly approaches the linear regime as the size of the radius increases:

Then, for a defect-free configuration, the contribution of the boundary to the energy density is a constant offset that persists even for large radii. In the presence of disclinations, on the other hand, the function χ(x) is no longer expected to be azimuthally symmetric and the harmonic function U (x) will not be constant throughout the paraboloid. In this case U (x) can be expressed in the integral form

where H(x, y) is the harmonic kernel (124) . In the regime of large core energies F c F el , the creation of defects is strongly penalized and the lattice necessarily has the minimum number of disclinations allowed by the topology of the paraboloidal substrate. From symmetry considerations, we might expect the optimal distribution of defects to consist of b +1−disclinations arranged along the boundary at the base vertices of a b-gonal pyramid and a b-fold apex (of topological charge q 0 = 6 − b) at the origin. The homogeneous boundary conditions adopted require the first term in Eq. (148) to vanish at the boundary. In the minimal energy configuration then, the system has the freedom to tune the total number of defects along the boundary to minimize the elastic energy Eq. (143) for any given value of the ratio κ. This behavior is exclusive to manifolds with boundary and doesn't have any counterpart in crystals on compact surfaces like the sphere and the torus.

We will label a pyramidal configuration by Y b , where b denotes the number of base +1−disclinations. The coordinates (r, φ) of the vertices are given by

Using the Euler theorem one can show that it is possible to construct infinite families of polyhedra with the symmetry group C bv from the pyramidal backbone Y b . The number of vertices is given by:

where n is a positive integer which represents the number of edges (not necessarily of the same length) of the polyhedron which separates two neighboring disclinations. In the following we will refer to these polyhedra with the symbol Y b,n . Fig.  32 illustrates two Y b,n lattices for the cases b = 4 and n = 7 (with V = 113), and b = 5 and n = 10 (V = 276). By a numerical minimization of the energy Eq. (143) one can establish that the Y b are indeed equilibrium configurations for b ∈ [3, 5] , for some range of the parameters κ and R. The cases of b = 5, 6 are particularly significant because they are characterized by an equal number of defects (N = 6) of the same topological charge (q = 1). The two configurations will be associated therefore with the same core energy F c and this introduces the possibility of a structural transition between Y 5 and Y 6 governed by the curvature ratio κ and the boundary radius R. For fixed R and small values of κ, the 6−fold symmetric configuration Y 6 is the global minimum of the free energy Eq. (143). For κ larger than some critical value κ c (R), however, the Y 6 crystal becomes unstable with respect to the 5−fold symmetric configuration Y 5 . A numerical calculation of the intersection point between the elastic energies of Y 5 and Y 6 for different values of κ and R allow us to construct the phase diagram shown in Fig. 33 . The word "phase" in this context refers to the symmetry of the ground state configuration as a function of the geometrical system parameters κ and R. In principle, if we keep increasing the curvature we might expect the crystal to undergo a further transition to the Y 4 phase. In this case, however, the core energy will also increase by a factor 4/3 and so this is not generally possible in the regime in which F c F el . For intermediate regimes (i.e. F c ∼ F el ), Y 5 → Y 4 and Y 4 → Y 3 transitions are also possible. The critical value of the parameters κ and R, however, is not universal and will depend on the precise values of the core energy and the Young modulus.

When the core energy F c is small, the elastic energy Eq. (143) can be lowered by creating additional defects. Let us assume that a fivefold disclination is sitting at the point x 0 = (r 0 , φ 0 ). We can introduce a notion of distance on the paraboloid by setting up a system of geodesic polar coordinates (s, ϕ) with origin at x 0 . We expect that the stress introduced by the defect is controlled by an effective disclination charge inside a circular domain C L of geodesic radius L:

where q = π/3 is the charge of the isolated defect and the integral measures the screening due to the total Gaussian curvature within the domain. The metric tensor and the Gaussian curvature of a generic Riemannian manifold can be expressed in geodesic polar coordinates in the form (see for example Do Carmo [144] ):

where G = g ϕ · g ϕ . Furthermore, an expansion of the metric around the origin (0, ϕ) yields:

For small distance from the origin, Eq. (157) becomes:

The right hand side of Eq. (160) is a very general expression for the effective disclination charge at small distance and doesn't depend on the embedding manifold. If a grain boundary is radiating from the original disclination, we expect the spacing between consecutive dislocations to scale like a/q ef f , with a the lattice spacing [73] . When q ef f → 0 + the dislocation spacing diverges and the grain boundary terminates. Since the Gaussian curvature is not constant, the choice of the origin (i.e. the position of the central disclination along the grain boundary) affects the evaluation of q ef f . One can identify upper and lower bounds by observing that:

Unlike the case of surfaces of constant Gaussian curvature, the phase diagram for paraboloidal crystals consists of three regions separated by the curves:

When L − L(K min ) → 0 + , the effective disclination charge goes to zero and the distance between two consecutive dislocations diverges at any point. On the other hand if L − L(K max ) → 0 − , the disclination charge will prefer to be delocalized in the form of grain boundary scars. For L(K max ) < L < L(K min ) the paraboloid will be equipped with regions where the Gaussian curvature is high enough to support the existence of isolated disclinations as well as regions where the screening due to the curvature is no longer sufficient and the proliferation of grain boundary scars is energetically favored. This leads to a three region phase diagram in which the regime of isolated disclinations is separated from the delocalized regime of scars by a novel phase in which both isolated disclinations and scars coexist in different parts of the paraboloid according to the magnitude of the Gaussian curvature. It is useful to measure the distance L c in terms of the lattice spacing a and rephrase Eq. (162) as a condition on a (or equivalently on the number of vertices V ). To do this we note that in order for the domain C L to completely screen the topological charge of the shortest scar possible (i.e. 5 − 7 − 5), the geodesic radius L has to be large enough to enclose the entire length of the scar. Calling the geodesic distance associated with a single lattice spacing a, we will then approximate L ∼ 3 . This leads to the following expression for the lattice spacing a at the onset of scar formation:

The lattice spacing a can be approximately expressed as a function of the number of vertices of the crystal by dividing the area A of the paraboloid by the area of a hexagonal Voronoi cell of radius a/2 with a 2 ≈ A/ Fig. 34 for the case R = 1. The two phase boundaries that separate the isolated defects (ID in the plot) regime from the coexistence regime and this one from the grain boundaries phase correspond to the solutions for K 0 = K min and K 0 = K max respectively. Fig. 35 shows the Voronoi lattice and the Delaunay triangulations obtained from a numerical minimization of a system of V classical particles (up to V = 200) constrained to lie on the surface of a paraboloid and interacting with a Coulomb potential of the form U (r) = 1/r, with r the Euclidean distance between two particles in R 3 . Simulations are carried out using a parallel implementation of Differential Evolution (DE) [142, 145] . In all the systems observed disclinations always appear clustered in either grain boundary scars or dislocations with the exception of isolated +1−disclinations which appear in the bulk as expected from the curvature screening argument discussed above. The complex aggregation of defects along the boundary together with the presence of negatively charged clusters indicates that the effect of the boundary, in the case of relatively small systems like the ones simulated, is more drastic than predicted by the homogeneous boundary conditions Eq. (145) . Even in the computationally expensive case of V = 100, the distance between the origin and the boundary of the paraboloid is only four lattice spacings. In this situation we expect the distribution of particles along the boundary to play a major role in driving the order in the bulk.

For larger systems, such as V = 200 (top of Fig. 35 ), the behavior of the particles in the bulk is less affected by the boundary and the crystalline order reflects more closely the free-boundary problem previously presented. A comparison of the lattices in the first two rows of Fig. 35 , in particular, reveals substantial agreement with the scenario summarized in Fig. 34 . For κ = 0.8 and V = 200, the defects are all localized along the boundary with the exception of one length-3 scar in the bulk at distance r ≈ 0.63 from the center. For κ = 1.6, the pattern of defects in the bulk is characterized by the coexistence of an isolated +1−disclination at the origin and a length-5 W −shaped scar displaced along a parallel one lattice spacing away from the central disclination. Apart from the evident difficulty in comparing the struc-tures of small systems with those predicted from continuum elasticity theory, this behavior is consistent with the simple picture sketched in the phase-diagram of Fig.  34 . The local 5−fold symmetry at the origin of the κ = 1.6 configuration, compared with 6−fold symmetry for κ = 0.8, suggests, as in the case of spherical crystals [94] , that the complicated structure of defect clusters appearing in large systems is the result of the instability of the simpler Y b,n configurations from which they partially inherit their overall symmetry. A more accurate numerical verification of our theory remains a challenge for the future.

The symmetry of the configurations presented in Fig. 35 deserves special attention. As for any surface of revolution, the circular paraboloid possesses the symmetry group O(2) of all rotations about a fixed point and reflections in any axis through that fixed point. Any given triangulation of the paraboloid may destroy the full rotational symmetry completely or just partially, leaving the system in one of the following two subgroups: the pyramidal group C nv or the reflection symmetry group C s . In general we found the latter symmetry group for system sizes up to V = 200 particles. The symmetry for larger system sizes is under investigation.

Some sixty years ago Bragg and Nye used bubble rafts to model metallic crystalline structures [146] . A carefully made assemblage of bubbles, floating on the surface of a soap solution and held together by capillary forces, forms an excellent two-dimensional replica of a crystalline solid, in which the regular triangular arrangement of bubbles is analogous to the close packed structure of atoms in a metal. Feynman considered this technique to be important enough that the famous Feynman lectures in physics include a reproduction of the original Bragg-Nye paper in its entirety [147] . Bubble rafts can be made easily and inexpensively, equilibrate quickly, exhibit topological defects such as disclinations, dislocations and grain boundaries, and provide vivid images of the structure of defects. Bubble raft models have been used to study two-dimensional polycrystalline and amorphous arrays [148] , nanoindentation of an initially defect-free crystal [149] , and the dynamic behavior of crystals under shear [150] . In this section we review the experimental realization of a bubble-raft model for a paraboloidal crystal done by Bowick et al [136] by assembling a single layer of millimeter-sized soap bubbles on the surface of a rotating liquid, thus extending the classic work of Bragg and Nye on planar soap bubble rafts.

As we mentioned in the introduction, paraboloidal shapes naturally occur across the air/liquid interface of a fluid placed in a rotating cylindrical vessel. It is simple to verify that the height z of such an interface above the xy−plane of a system of Cartesian coordinates is given by

where ω is the angular velocity and g the gravitational acceleration and C = D − ω 2 R 2 /4g is an integration constant following from the requirement that the volume of the liquid during rotation be equal to the volume of the liquid at rest (D is its height) [151] . Thus the normal curvature at the origin is given by κ = ω 2 /g. Fig. 36 shows a computer reconstruction from Ref. [136] of two bubble rafts with κ = 0.15 cm −1 and κ = 0.32 cm 1 respectively. Soap bubbles have diameter a = 0.84(1) mm with monodispersity ∆a/a ≈ 0.003. The images shown in Fig.  36 are taken with a CCD digital camera rotating along with a cylindrical vessel of radius R = 5 cm and Delaunay triangulated to determine the lattice topology. The degree of order of the crystalline raft can be characterized by the translational and orientational correlation functions g(r) and g 6 (r) [110] . The former gives the probability of finding a particle at distance r from a fixed particle at the origin. The function is normalized with the density of an equivalent homogeneous system in order to ensure g(r) = 1 for a system with no structure. Interactions between particles build up correlations in their position and g(r) exhibits decaying oscillations, asymptotically approaching one. For a two-dimensional solid with a triangular lattice structure the radial correlation function is expected to exhibit sharp peaks at r/a = √ n 2 + nm + m 2 = 1, √ 3, 2, 2 √ 3 . . ., while the amplitude of the peaks decays algebraically as r −η with η = 1/3 (dashed line in Fig. 37 ). Within the precision of the data, the positional order of the paraboloidal crystals assembled with the bubble raft model reflects this behavior, although with more accurate measurements one might see dependence of the exponent η on the curvature (because of the proliferation of defects with curvature).

The orientational correlation function g 6 (r) is calculated as the average of the product ψ(0)ψ * (r) of the hexatic order parameter over the whole sample. For each bubble (labeled j) that has two or more neighbors, ψ j (r) = (1/Z j ) Zj k=1 exp(6iθ jk ), where Z j is the number of neighbors of i and θ jk is the angle between the j − k bonds and a reference axis. One expects g 6 (r) to decay exponentially in a disordered phase, algebraically in a hexatic phase and to approach a finite asymptote in the case of a crystalline solid. In the systems studied g 6 (r) approaches 0.8 within 5-6 lattice spacings.

Of particular interest is the structure of the grain boundaries appearing in the paraboloidal lattice for different values of the curvature parameter κ. Grain boundaries form in the bubble array during the growing process as a consequence of geometrical frustration. As noted, any triangular lattice confined in a simply connected region with the topology of the disk is required to have a net disclination charge Q = 6. In the absence of curvature, however, the elastic stress due to an isolated disclination is extremely high and defects are energetically favored to cluster in the form of a grain boundary consisting of one-dimensional arrays of tightly shows that they may deform slightly to better fill space, whereas the computer reconstruction of the lattice (right) uses perfect spheres of uniform size. From Ref. [136] .

bound (5, 7)−fold disclination pairs. In a planar confined system, grain boundaries typically span the entire length of the crystal. In the curved case, however, they can appear in the form of scars carrying a net +1 topological charge and terminating in the bulk of the crystal.

Prominent examples of grain boundaries are visible in the two lattices shown in Fig. 36 . For a gently curved paraboloid (with κ ≈ 0.15 cm −1 ), grain boundaries form long (possibly branched) chains running from one side to the other and passing through the center. As the curvature of the paraboloid is increased, however, this long grain boundary is observed to terminate in the center (see Fig. 36d ; a closeup version of this image is seen in Fig. 38 ). For R = 5 cm, the elastic theory of defects predicts a structural transition at κ c = 0.27 cm −1 in the limit of large core energies. In this limit the creation of defects is strongly penalized and the lattice has the minimum number of disclinations required by the topology of the embedding surface. In a low curvature paraboloid (κ < κ c ) these disclinations are preferentially located along the boundary to reduce the elastic energy of the system. When the aspect ratio of the paraboloid exceeds a critical value κ c (R), however, the curvature at the origin is enough to support the existence of a 5−fold disclination and the system undergoes a structural transition. In the limit of large core energies, when only six disclinations are present, such a transition implies a change from the C 6v to the C 5v rotational symmetry group.

Together with the theoretical argument reviewed in the previous section, these experimental observations point to the following mechanism for scar nucleation in a paraboloidal crystal. In the regime in which the creation of defects is energetically inexpensive, geometrical frustration due to the confinement of the lattice in a simply connected region is responsible for the formation of a long side-to-side grain boundary. When the curvature of the paraboloid exceeds a critical value, dependent on the radius of the circular boundary, the existence of a +1 disclination near the center is energetically favored. Such a disclination serves as a nucleation site for 5-7 dislocations and the side-to-side grain boundary is replaced by a terminating center-to-side scar.

Circa twenty years afer their first observation in partially polymerized diacetylenic phospholipid membranes by Mutz and Bensimon [152] , self-assembled toroidal aggregates are now considered the progenitors of a magnificent cornucopia of complex structures, also featuring branched network and micellar surfaces of high genus. The existence of such complex structures has become an experimental fact thanks to the enormous work in the past decade on the study of self-assembly of amphiphilic compounds such as lipids, surfactants and amphiphilic block copolymers. After the work of Eisenberg and coworkers [153, 154, 155] it became clear that block copolymers in particular afford access to a variety of complex structures spanning an unexpectedly vast range of topologies and geometries. More recently, several experimental studies have been performed on di-and triblock copolymers with the intent of unraveling the origin of morphological complexity in copolymer surfactants and some possible pathways for micelle and vesicle formation have been proposed.

Jain and Bates, for instance, reported the formation of several non-simplyconnected micellar structures from the self-assembly of diblock copolymer poly(1,2butadiene-b-ethylene oxide) (PB-PEO) [156] (see Fig. 39 ). These polymers selfassemble in Y −shaped junctions and these form the building blocks for more complex micellar structures via re-assembly. Later, Pochan et al. [157] found that almost all of the microstructures assembled from poly(acrylic acid-b-methyl acrylateb-styrene) (PAA99-PMA73-PS66) triblock copolymers are ringlike or toroidal micelles. The route to toroidal micelles in block copolymers, however, is believed to be different from that to more complex network structures suggested by Jain and Bates, since residual Y −shaped aggregates are very rarely found in the sample. According to the authors of Ref. [157] , the formation of toroidal micelles has to be attributed primarly to the collapse of cylindrical micelles. On the other hand, mere end-to-end connection of cylindrical micelles doesn't appear to be the exclusive ring-forming mechanism in triblock systems since the average circumference of the self-assembled tori appears smaller than the contour length of the average cylindrical micelles. It seems to be a more complicated process, possibly involving interactions and exchange of matter between neighboring cylinders. Toroidal micelles have also been observed in recent experiments by Kim et al [158] from the self-assembly of amphiphilic dumbbell molecules based on a aromatic rod segment that is grafted by hydrophilic polyether dendrons at one end and hydrophobic branches at the other end. Molecular dumbbells dissolved in a selective solvent self-assemble in an aggregate structure due to their amphiphilic character. This process has been observed to yield coexisting spherical and openended cylindrical micelles. These structures, however, change slowly over the course of a week to toroidal micelles which thus appear more stable. The formation of ringlike structures was explained in this case as result of the coalescence of spherical micelles occurring to reduce the contact between hydrophobic segments and water molecules. An alternative pathway for the self-assembly of toroidal structures in copolymer solutions has been proposed and numerically tested by He and Schimd [159] . In this pathway, the micelles do not coalesce, but simply grow by attracting copolymers from the solution. Once a critical micelle size is exceeded, copolymers start to flipflop in such a way that the micelle core becomes itself solvent-philic (semi-vesicle state). Finally the solvent diffuses inside the core and the semi-vesicle swells into a vesicle.

The existence of toroidal aggregates in amphiphilic compounds was suggested theoretically for smectic-C (SmC) membranes, based on the argument that orientational order is not frustrated on a surface with zero Euler characteristic and thus a toroidal topology may be energetically favoured in p−atic membranes with large order parameter coupling compared with the bending rigidity [59, 160] . This hypothesis was tested by Evans in 1995 with the result that a toroidal topology is indeed prefered over the spherical one for wide range of geometrical and mechanical parameters in defect-free p−atic tori. The role of disclinations in p−atic tori was investigated systematically by Bowick, Nelson and Travesset nine years after the original work of Evans with the conclusion that, even if not required by topological constraints, unbound disclinations can be energetically convenient even in very dense systems [161] . The precise number of unbound disclinations is controlled primarly by the aspect ratio of the torus. The existence of defects in the ground state of an ordered phase on the torus is crucial for toroidal crystals where it leads to some unique structural features as well as a spectacular example of a curvature-driven transition to a disordered, liquid-like, state in the limit the aspect ratio approaches to one [162, 163] . In the remainder of this section we will review these three examples of order on embedded tori. 

The standard two-dimensional axisymmetric torus embedded in R 3 is obtained by revolving a circle of radius R 2 about a coplanar axis located at distance R 1 ≥ R 2 from its center. Choosing the symmetry axis as the z direction of a Cartesian frame, a convenient parametrization can be found in terms of the polar angle ψ along the revolving circle and the azimuthal angle φ on the xy−plane:

where ψ , φ ∈ [−π, π). These coordinates satisfy the Cartesian equation

In this parametrization the metric tensor g ij and Gaussian curvature K are given by

The Gaussian curvature is therefore positive on the outside of the torus, negative on the inside and zero along the two circles of radius R 1 at ψ = ±π/2 (see Fig.  40 ). Moreover K is maximally positive along the external equator at ψ = 0 and maximally negative on the internal equator at ψ = ±π. The Gauss-Bonnet theorem requires the total topological charge and the integrated Gaussian curvature to be zero on the torus:

A global measure of the curvature of the embedded torus is provided by the aspect ratio r = R 1 /R 2 . Our discussion will be limited to the case r ≥ 1. A "fat" torus with r = 1 is obtained by taking R 1 = R 2 and is characterized by a singularity in the Gaussian curvature along the internal equator at ψ = ±π. The case r < 0 corresponds to a self-intersecting torus whose symmetry axis lies in the interior of the revolving circle. The "skinny" torus limit, r → ∞, can be obtained either by taking R 2 → 0 with finite R 1 or by letting R 2 stay finite and taking R 1 → ∞. The Gaussian curvature diverges in the first case and goes to zero in the second. Neither case, however, reflects a real physical situation since the area A = (2π) 2 R 1 R 2 of the torus becomes zero or infinite respectively. Non-axisymmetric tori can be generated in various ways by deforming the surface of revolution (165) . A special and physically relevant transformation consists of an inversion wrt the unit sphere, a translation about a vector β and a second inversion. Thus a point R on the surface is mapped onto:

This map is conformal and transforms an axisymmetric shape to a non-axisymmetric one if the vector β is not parallel to the symmetry axis [164] . As β is increased in magnitude, the asymmetry increases until, in the limit, the torus becomes a perfect sphere with an infinitesimal handle. The physical importance of mapping (170) relies on the fact that it leaves the bending energy κ d 2 x H 2 invariant with significant consequences for the equilibrium shape and stability of fluid toroidal membranes.

As in any closed manifold the traditional Green-Laplace equation doesn't have a solution on the torus. An alternative Green function, can be obtained from the modified equation:

As usual the calculation of the Green function can be simplified considerably by conformally mapping the torus to a domain of the Euclidean plane via a suitable system of isothermal coordinates. Intuitively the torus is conformally equivalent to a rectangular domain described by a system of Cartesian coordinates. To make this explicit, one can equate the metric of the torus in the coordinates (ψ, φ) to a conformally Euclidean metric in the coordinates (ξ, η):

where w is a positive conformal factor. Taking η = φ and w = (R 1 + R 2 cos ψ) 2 , the coordinate ξ is determined by the differential equation:

where r = R 1 /R 2 , the aspect ratio of the torus, may be taken greater than or equal to one without loss of generality. Choosing the plus sign and integrating both sides of Eq. (172) we find:

Taking ψ ∈ [−π, π), the integral (173) yields:

In the transformed coordinate system (ξ, η) the modified Green-Laplace equation reads:

where ∆ and δ are now the Euclidean Laplacian and delta function. The function G L (x, y) can be expressed in the form:

where G 0 (x, y) is the Laplacian Green function on a periodic rectangle and the angular brackets stand for the normalized integral of the function G 0 (x, y) with respect to the dotted variable:

Analogously the function G 0 (· , · ) is given by

and ensures the neutrality property:

The modified Laplacian Green function on a periodic rectangle of edges p 1 and p 2 can be conveniently calculated in the form:

where u λ is the eigenfunction of the Laplace operator with periodic boundary conditions:

such that:

.

In Cartesian coordinates the eigenfunctions are simple plane waves of the form:

where λ n and µ m are given by:

and the eigenvalue λ is given by:

Calling for simplicity x = (x, y) and y = (ξ, η), the function G 0 is given by:

Summing this series eventually leads to [163] :

where ϑ 1 (u|τ ) is the Jacobi theta function [165, 166] defined as:

with q = exp(iπτ ), z = x + iy and ζ = ξ + iη.

There has been much effort in recent years to shed light on possible pathways in the self-assembly of toroidal micelles and vesicles from block copolymer solutions and other amphiphilic compounds. In particular, the spontaneous formation of structures of genus g ≥ 1 from preexisting spherical objects, a process appearing in several different scenarios proposed in the literature, has attracted the most attention and debate because of its exotic character. The simplest question one can ask in this context is whether a vesicle can be energetically favored to change its topology from spherical to toroidal once the total surface area, and therefore the number of constituent molecules in the vesicle, is specified. Evans addressed this problem in 1995 [160] and showed how such a transition between genera is indeed possible in fluid membranes and is even enhanced if the membrane is endowed with in-plane orientational order. Although oversimplified if compared to the great complexity of the real self-assembly mechanism, Evans' calculation provides insight into the delicate problem of stability of toroidal vesicles as well as a good starting point for our discussion of order on the torus. Let θ be the local orientation of a p−atic director field on a surface as defined in Sec. 2. A standard orientational order parameter is given by scalar field:

where · deontes a thermal average. The total elastic energy of the vesicle consists of a pure bending term F b , describing the elasticity of the membrane in the liquid state, and a term F p associated with the internal p−atic order. Thus F = F p + F b , with:

where κ and κ g are the bending and Gaussian rigidity respectively 1 and the operator ∇ − ipΩ is obtained from the covariant derivative of a p−atic tensor order parameter expressed in a local orthonormal frame as described in Sec. 2.4, with Ω the covariant vector of Eq. (53) that parametrizes the spin-connection. The free energy (186a) is invariant under rotations of the local reference frame by an arbitrary angle χ:

which is a typical gauge transformation. The gradient term in Eq. (186a) is the same as in Eq. 52, upon identifying K A = p 2 C. As observed by Park et al. [168] , Eqns. (186a) and (186b) closely resemble the Landau-Ginzburg Hamiltonian of a superconductor in an external magnetic field:

When exposed to an external magnetic field, a super-conducting material can undergo a second order mean-field transition from a metal to a super-conductor characterized by an Abrikosov lattice of vortices whose density is determined by the temperature and the applied magnetic field H. The magnetic field, in particular, is conjugate to the vortex number N v since:

where L is the length of the sample in the direction of H and φ 0 = hc/2e is the flux quantum. In this context the vector potential associated with the applied magnetic field is replaced by the covariant vector Ω whose curl is the Gaussian curvature. Morover, on a closed surface:

Thus a p−atic phase on a closed 2−manifold is analogous to an Abrikosov phase with a fixed total vorticity rather than a fixed applied magnetic field.

The transition to an ordered phase in controlled by the parameter τ in Eq. (186a). For τ above a critical value τ c (τ c = 0 on a flat surface), ψ = 0 and the system is in the isotropic phase. In this regime the only contribution to the energy is determined by the shape of the vesicle as expressed by the bending term (186b). For τ < τ c , on the other hand, the gauge symmetry is spontaneously broken and ψ = 0, with the exception of a number of isolated defective points. Within a meanfield approach, thermal fluctuations can be neglected and the preferred configuration of the system corresponds to the minimum of the free energy. A ground state configuration for the p−atic order parameter can be found by rewriting the gradient energy term in Eq. (186a) as

with D k = ∂ k −ipΩ k , and expressing ψ in the basis of eigenfunctions of the operator in parentheses: ψ = n a n ϕ n (190) with ϕ n satisfying the Hermitean equation

with the normalization condition

The complex coefficients a k are then determined by minimization of the free energy F . The eigenfunctions are sometimes referred to as "Landau levels" and are typically degenerate. If ψ is expanded in this complete set of eigenfunctions then, in meanfield, the partition function is dominated by configurations involving only Landau levels with the lowest eigenvalue λ 0 . Taking lowest levels only, the p−atic free energy (186a) becomes

To find the lowest energy configuration of the field ψ one now has to solve the eigenvalue problem (191) to determine the lowest eigenvalue λ 0 and minimize the free energy with respect to the parameters of the linear combination (190) . Notice that the mean field transition has been lowered from τ c = 0 to τ c = −4πλ 0 C/A. Although λ 0 = 0 on a flat surface, it is finite and positive-definite on a surface with finite Gaussian curvature [160] . For a sphere of unit radius parametrized in standard spherical coordinates (θ, φ), the eigenvalue problem (191) is solved by λ 0 = p and ϕ n of the form [169] ϕ n = 2p + 1 4π(p + n)!(p − n)! sin p+n θ 2 cos p−n θ 2 e inφ (194) for integer values of n between −p and p. These functions have 2p zeros, corresponding to topological defects of unit charge and winding number 1/p. Because of the lowest Landau level approximation used to derive Eq. (193) , eigenfunctions (194) are degenerate to quadratic order in ψ. This implies the freedom to place defects anywhere on the sphere with no additional energy cost of order ψ 2 . The ψ 4 term, on the other hand, lifts this degeneracy and make the defects repel. For axisymmetric tori the eigenvalue problem (191) was solved numerically by Evans with the result shown in Fig. 41 . Eigenvalues are plotted as a function of the aspect ratio r of the torus for various p−atic textures labelled p n with 0 ≤ n ≤ p. Eigenfunctions ϕ n have 4n zeros corresponding to 2n disclinations of topological charge q = 1 (winding number 1/n) and 2n disclinations of topological charge q = −1 (winding number −1/n). Except for vector order (p = 1), defective configurations are energetically favored for small aspect ratios and the number of disclination pairs increases with r as a consequence of the larger (in magnitude) Gaussian curvature.

Carrying out the integrals in Eq. (186a) and (186b) in the lowest Landau level approximation, the total elastic free energy of spherical and toroidal vesicles can be written as

where J is the minimal value of the integral of |ψ MF | 4 and ψ MF = p n=−p a n ϕ n sphere a + ϕ p + a − ϕ −n torus is the linear combination of eigenfuctions that are degenerate to quadratic order in ψ. The second term in Eqs. (195a) and (195b) disappears in the isotropic phase when τ > −4πλ 0 C/A due to the intrinsic orientational order on the manifold. Even in this case, Eqs. (195a) and (195b) reveal the existence of a transition line between spherical and toroidal shape. In the absence of in-plane orientational order the bending energy term in Eq. (195b) is minimized by the so called Clifford torus with r = √ 2. In this case F torus = 4π 2 κ, which is smaller than F sphere = 4π(2κ+κ g ) when κ g /κ > π − 2. Toroidal vesicles with r = √ 2 are therefore energetically favored in the fluid phase for large values of κ g /κ. This argument cannot be invoked to explain the complex self-assembly of toroidal structures from homogeneous solutions mentioned in the introduction, but does provide a simple (but non-trivial) example of a toroidal shape being energetically preferred to a spherical one.

The critical value of κ g /κ of the genera transition is lowered for vesicles possessing in-plane p−atic order. This is displayed in the phase diagrams of Fig. 42 for the case of vector (p = 1), nematic (p = 2), tetradic (p = 4) and hexatic (p = 6) order. The lines, whose position depends on κ g /κ, separate spheres (above the line) from tori (below the line). Even for simple vector order, stable toroidal vesicles exist as well for κ g /κ < π − 2. For κ g /κ > π − 2, furthermore, transition lines are closed, with spheres on the inside and tori on the outside. In the shaded regions nonaxisymmetric tori are favored over symmetric ones. The hexatic phase diagram (bottom right corner of Fig. 42 ) deserves special attention. The shaded region is split in two parts. That on the right contains spheres if κ g /κ < π − 2 and nonaxisymmetric tori otherwise, while that on the left contains spheres above the transition line and non-axisymmetric tori below. The white region of the phase diagram is also divided in two. The top left region of the diagram behaves as described above, while the other part is characterized by closed lines separating spheres (outside) and small tori (inside). These last class of tori exhibit ten pairs of q = ±1 disclinations with positive disclinations distributed on the external equator of the torus and negative disclinations along the internal equator. This last feature is an important property sheared by toroidal objects with in-plane order and will be clarified in the following sections. The results reviewed here are valid within the mean-field approach and the lowest Landau level approximation, which both hold deep in the ordered phase that we are most interested in. At higher temperatures Evans proved the approximations are still valid away from the transition line and for C/A κ/k B T 1 [160] . The latter condition is generally fulfilled by several systems (i.e. κ/k B T = 1 − 10 for lipid bilayers).

Evans' analysis, summarized in the previous section, indicates that disclinations, even if not required by topological constraints, can nonetheless appear in the ground state of p−atic tori as a consequence of the coupling between in-plane orientational order and spatial curvature. The occurrence of defects in the ground state of toroidal hexatic vesicles was systematically investigated by Bowick, Nelson and Travesset based on the formalism outlined in Sec. 2 [161] . Before embarking on a detailed analysis of the elasticity of defects in hexatic tori, it is instructive to obtain a rough estimate of how many pairs of ±1-disclinations would be required to achieve perfect screening of the background topological charge associated with the Gaussian curvature of the torus. Consider a wedge of angular width ∆φ on the outside wall of positive Gaussian curvature. The net curvature charge associated with this region is

Upon equating s ef f to 2π/6, the charge of a single disclination, one finds that ∆φ = 2π/12, independently of R 1 and R 2 . Thus, 2π/∆φ = 12 positive disclinations would be required to completely compensate the negative curvature of the inner wall. This simple argument neglects core energies and interactions between disclinations, effects which will cause the preferred number of defect pairs to be less than twelve. The total elastic energy of a toroidal hexatic vesicle containing N disclinations of topological charge q i (i = 1 . . . N ) takes the form 197) where the first two terms represent the contributions due to the pair interaction between defects and the interaction of defects with the topological charge of the substrate associated with the Gaussian curvature. The third term in Eq. (197) is the spin-wave part of the frustrated hexatic energy while the last two terms represents the bending energy and the defect core energy respectively. The defect- Figure 43 . (Color online) The total energy (in units of K A /2 of hexatic tori of aspect ratio r = √ 2 (left) and r = 2.6926) for varying number of defects. The bending energy at fixed r is subtracted off. The disclination core energy is set to 0.1K A , which is 0.2 in the above units.

defect interaction potential has been calculated in Ref. [161] :

where α is related to the polar angle ψ on a cross-section of the torus by: 199) and η is the Dedekind eta function:

The one-body interaction between defects and curvature, on the other hand, is expressed via the potential energy

For opposite sign disclinations, the defect-defect interaction, determined by the function Q(x i , x j ), is attractive for all separations at constant φ. If only this term were present, the attraction would bring both charges as close as possible, binding all disclinations into dipoles which have a higher energy than a defect-free configuration. Thus if no other terms were present, the ground state configuration would be defect free. The defect-curvature interaction term L(x i ), however, favors the appearance of additional defects. This term acts like an electric field pulling the positive (negative) disclinations into regions of positive (negative) Gaussian curvature.

The elastic energy (197) was analyzed numerically in Ref. [161] for various aspect ratios and defect core energies. A plot of the energy of configurations obtained by placing a ring of (N/2) (+1)-disclinations equally spaced along the same parallel on the outside of the torus and a second ring of (N/2) (−1)-disclinations on the inside, is shown in Fig. 43 for various N and variable angular separation ∆ψ between the rings. The disclination core energy is taken to be c = cK A , with c a numerical constant taken equal to 0.1 in the plots. The energy at the maximum separation (∆ψ = π) first decreases and then increases with N . The optimal number of defect pairs N/2 is 6 and 7 (the latter curve is not shown in the plot) for r = √ 2 and r = 2.6926 respectively.

The existence of defects in the ground state of a toroidal vesicle with hexatic order is also affected by the total number of molecules N forming the hexatic phase. In Ref. [161] it was found that defects disappear in the limit N → ∞, but are present for numbers of molecules as large as N = 10 10 . This number is several orders of magnitude larger than the typical number of molecules of biological vesicles such as red-blood cells (i.e. N ∼ 10 8 ). To make this estimate, one considers a pair of opposite sign disclinations that have been pulled apart from the circle of zero Gaussian curvature to the equators in the regions of like sign Gaussian curvature; thus ψ + = 0 and ψ − = π. Upon approximating the defect-pair energy by its flat space value, the total energy reads:

where a 0 is the lattice spacing. Eq. (202) changes sign when

Taking

leads to the conclusion that defects are favored for:

where the last identity has been obtained by taking c /K A = 0.1. This result establishes that defects are present in the ground state of a hexatic tours for any fixed number of molecules provided the torus is sufficiently fat. For the energetically favored Clifford torus with r = √ 2, Eq. (204) predicts a critical number of molecules to be order N ∼ 10 10 .

In absence of defects the ratio between the hexatic stiffness K A and the bending rigidity κ dictates the optimal shape of the toroidal vesicle. Taking q i = 0 in Eq. (197) one obtains in this case:

The first term represents the energetic cost associated with the distortion of the hexatic director field due to the Gaussian curvature alone. In the limit of large and small hexatic stiffness Eq. (205) is minimized by:

which provides a compelling example of the interplay between order and geometry: if the stiffness associated with the intrinsic hexatic order is much smaller than the bending rigidity, the Clifford torus is the optimal geometry; if on the other hand, the hexatic stiffness dominates, then a thin torus, similar to a bicycle tire, is optimal.

Crystalline assemblages of identical sub-units packed together and elastically bent in the form of a torus have been found in the past ten years in a variety of systems of surprisingly different nature, such as viral capsids, self-assembled monolayers and carbon nanomaterials. In the introduction we mentioned the self-assembly of toroidal micelles and vesicles from homogeneous solutions of amphiphilic molecules such as oligomers of aromatic compounds or block copolymers. Toroidal geometries also occur in microbiology in the viral capsid of the coronavirus torovirus [170] . The torovirus is an RNA viral package of maximal diameter between 120 and 140 nm and is surrounded, as other coronaviridae, by a double wreath/ring of cladding proteins. Carbon nanotori form another fascinating and technologically promising class of toroidal crystals [49] with remarkable magnetic and electronic properties. The interplay between the ballistic motion of the π electrons and the geometry of the embedding torus leads to a rich variety of quantum mechanical properties including Pauli paramagnetism [171] and Aharonov-Bohm oscillations in the magnetization [172] . Ring closure of carbon nanotubes by chemical methods [173] suggest that nanotubes may be more flexible than at first thought and provides another technique of constructing carbon tori. In this section we review some recent developments in the study of the geometry and the elasticity of toroidal crystals. Additional details can be found in Refs. [162, 163] .

Toroidal crystals with local 6−fold order can be described as toroidal polyhedra with all triangular faces. After the discovery of carbon nanotubes in 1991 and the subsequent theoretical construction (later followed by the experimental observation) of graphitic tori, this class of polyhedra has drawn considerable attention and many possible tessellations of the circular torus have been proposed by the community [174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184] . Here we review the construction of a defect-free triangulated torus and we show how the most symmetric defective triangulations can be generally grouped into two fundamental classes corresponding to symmetry groups D nh and D nd respectively.

The structure of a triangulated cylinder can be specified by a pair of triangular lattice vectors c and t, called the chiral and translation vector respectively, which together define how the planar lattice is rolled up. In the canonical basis a 1 = (1, 0) and a 2 = ( 1 2 , √ 3

2 ), the vector c has the form The translation vector t, on the other hand, can be expressed as an integer multiple

of the shortest lattice vector e t perpendicular to c. The vector e t is readily found to be of the form e t = (n + 2m)a 1 − (2n + m)a 2 (n + 2m : 2n + m) ,

where (a : b) denotes the greatest common divisor of a and b and enforces the minimal length. The three-dimensional structure of the torus is obtained by connecting the edge OT of the rectangle in Fig. 44 (left) to O C and OC to O T . The edge OT is then mapped to the external equator of the torus while the edge OC to the φ = 0 meridian. The resultant toroidal lattice has characteristic chirality related to the initial choice of the vector c. In the nanotubes literature armchair refers to the lattice obtained by choosing n = m, zigzag to that obtained for m = 0 and chiral to all other lattices. An example of an (n, m, l) chiral torus is shown in Fig. 44 (right) for the case n = 6, m = 3 and l = 6. The chirality is extremely important in graphitic carbon nanotube or nanotori, where it determines whether the electronic behavior of the system is metallic or semiconducting. By Euler's theorem one can prove that the number of triangular faces F and the number of vertices V of a triangular toroidal lattice is given by

Denoting A R the area of the rectangle with edges c and t and A T the area of a fundamental equilateral triangle, the number of vertices of a defect-free toroidal triangulation is then:

The planar construction reviewed above allows only lattices with an even number of vertices. Defect-free toroidal deltahedra with an odd number of vertices are also possible and their construction is generally achieved by assembling congruent octahedral building blocks. We refer the reader to Ref. [185] for an comprehensive review of the topic.

The embedding of an equal number of pentagonal and heptagonal disclinations in the hexagonal network was first proposed by Dunlap in 1992 as a possible way to incorporate positive and negative Gaussian curvature into the cylindrical geometry of carbon tubules [174] . According to the Dunlap construction the necessary curvature is incorporated by the insertion of "knees" (straight cylindrical sections of the same diameter joined with a kink) in correspondence with each pentagonheptagon pair arising from the junction of tubular segments of different chirality (see Fig. 46 ). In particular, a junction between a (n, 0) and a (m, m) tube can be obtained by placing a 7−fold disclination along the internal equator of the torus and a 5−fold disclination along the external equator. Since the radii of the two segments of a junction are different by construction, the values of n and m are commonly chosen to minimize the ratio |c (n,0) |/|c (m,m) | = n/ √ 3m. By repeating the 5 − 7 construction periodically it is possible to construct an infinite number of toroidal lattices with an even number of disclinations pairs and dihedral symmetry group D nh (where 2n is the total number of 5−7 pairs, Fig. 45 top and Fig. 47 ). The structure of the lattice is described by the alternation of two motifs with crystalline axes mutually rotated by 30 • as a consequence of the connecting disclination. One of the fundamental aspects of Dunlap's construction is that all the disclinations are aligned along the two equators of the torus where the like-sign Gaussian curvature is maximal. As we will see below, this feature makes these arrangements optimal in releasing the elastic stress due to curvature. Other defective triangulations with D nh symmetry group have been proposed in Ref. [163] and won't be discussed here.

Another class of crystalline tori with dihedral antiprismatic symmetry D nd was initially proposed by Itoh et al [177, 178, 179] shortly after Dunlap. Aimed at reproducing a structure similar to the C 60 fullerene, Itoh's original construction implied ten disclination pairs and the point group D 5d . In contrast to Dunlap tori, disclinations are never aligned along the equators in antiprismatic tori, instead being staggered at some angular distance δψ from the equatorial plane. Hereafter we will use the symbols TPn and TAn to refer to toroidal lattices with 2n disclination pairs and symmetry group D nh and D nd respectively. A systematic construction of defected triangulations of the torus can be achieved in the context of planar graphs [188, 189, 190] . A topological embedding of a graph in a two-dimensional manifold corresponds to a triangulation of the manifold if each region of the graph is bounded by exactly three vertices and three edges, and any two regions have either one common vertex or one common edge or no common elements of the graph. The simplest example of toroidal polyhedra with D nd symmetry group, featuring only 5−fold and 7−fold vertices, can be constructed by repeating n times the unit cell of Fig. 48a . These toroidal antiprisms have V = 4n vertices and can be obtained equivalently from the edge skeleton of a n−fold antiprism by attaching at each of the base edges a pentagonal pyramid and by closing the upper part of the polyhedron with n additional triangles. By counting the faces one finds F = 5n + 2n + n = 8n from which V = 4n. The simplest polyhedron of this family has V = 12 and D 3d symmetry group (see top left of Fig. 49 ) and corresponds to the "drilled icosahedron" obtained by removing two parallel faces of an icosahedron and connecting the corresponding edges with the six lateral faces of an antiprism with triangular base (i.e. a prolate octahedron). Starting from this family of toroidal antiprisms a number of associated triangulations having the same defect structure can be obtained by geometrical transformations such as the Goldberg inclusion [39, 191, 192] . Such transformations, popularized by Caspar and Klug for the construction of the icosadeltahedral structure of spherical viruses [39] , consist in partitioning each triangular face of the original graph into smaller triangular faces in such a way that old vertices preserve their valence and new vertices have valence six. The partition is obtained by specifying two integer numbers (L, M ) which define how the original vertices of each triangle are connected by the new edges so that the total number of vertices is increased by a factor T = L 2 + LM + M 2 .

A general classification scheme for D nd symmetric tori was provided by Berger and Avron [189, 190] in 1995. Their scheme is based on the construction of unit graphs comprising triangular tiles of different generations. In each generation, tiles are scaled in length by a factor 1/ √ 2 with respect to the previous generation. This rescaling approximates the non-uniformity of the metric of a circular torus.

In the past few years, alternative constructions of triangulated tori have been proposed as well as novel geometrical and graph-theoretical methods to express the coordinates of their three-dimensional structures (see for example Kirby [180] , László at al [181, 182] , Diudea et al [183, 184] ). Here we choose to focus on the defect structure associated with the two most important class TPn and TAn with groups D nh and D nd .

The total free energy of a toroidal crystal with N disclinations can be expressed as usual as

where Y is the two dimensional Young modulus and

The defect part of the stress function Γ(x) can be found by integrating the Green function (183) and takes the form

where Li 2 is the usual Eulerian dilogarithm and

The function Γ s (x) representing the stress field due to the Gaussian curvature of the torus, on the other hand, is given by

A derivation of the functions Γ d (x, x k ) and Γ s (x) is reported in Ref. [163] and won't be repeated here.

To analyze the elastic free energy (209) we start by considering the energies of two opposite sign disclinations constrained to lie on the same meridian. The elastic free energy of this system is shown in Fig. 50 as a function of the angular separation between the two disclinations. The energy is minimized for the positive (5−fold) disclination on the external equator (maximally positive Gaussian curvature) and the negative (7−fold) disclination on the internal equator (maximally negative Gaussian curvature). The picture emerging from this simple test case suggests that a good ansatz for an optimal defect pattern is a certain number p of equally spaced +1 disclinations on the external equator matched by the same number of equally spaced −1 disclinations on the internal equator. We name this configuration with the symbol T p , where p stands for the total number of disclination pairs.

where the two pairs of numbers specify the (ψ, φ) coordinates of the positive and negative disclinations respectively. A comparison of the energy of different T p configurations, as a function of aspect ratio and disclination core energy, is summarized in the phase diagram of Fig. 51 . We stress here that only T p configurations with p even have an embedding on the torus corresponding to lattices of the TP p 2 class. Nevertheless a comparison with p−odd configurations can provide additional information on the stability of p−even lattices. For small core energies, moreover, thermally excited configurations with a large number of defects and similar p−polar distributions of topological charge are expected to exhibit an elastic energy comparable in magnitude with that of these minimal constructions. The defect core energy has been expressed here in the form

The core energy c of a single disclination depends on the details of the crystalforming material and the corresponding microscopic interactions. A simple phenomenological argument (see for example Ref. [193] ) gives

where a the lattice spacing. Taking a 2 = A/ √ 3 2 V , with A the area of the torus, yields

For a system of order V = 10 3 subunits, then, the dimensionless core energy on the left hand side of Eq. (216) is of order 10 −5 . This estimate motivates our choice of the scale for c /(AY ) in Fig. 51 . For dimensionless core energies below 4·10 −5 and aspect ratios r between 3.68 and 10.12 the ground state structure is the TP5 lattice corresponding to a double ring of +1 and −1 disclinations distributed on the external and internal equators of the torus as the vertices of a regular decagon (the T 10 configuration). The TP5 lattice has dihedral symmetry group D 5h . That this structure might represent a stable configuration for polygonal carbon toroids has been conjectured by the authors of Ref. [186] , based on the argument that the 36 • angle arising from the insertion of ten pentagonal-heptagonal pairs into the lattice would optimize the geometry of a nanotorus consistently with the structure of the sp 2 bonds of the carbon network (unlike the 30 • angle of the 6−fold symmetric configuration originally proposed by Dunlap) . In later molecular dynamics simulations, Han [194, 195] found that a 5−fold symmetric lattice, such as the one obtained from a (9,0)/(5,5) junction (see Fig. 47 ), is in fact stable for toroids with aspect ratio less then r ∼ 10. The stability, in this case, results from the strain energy per atom being smaller than the binding energy of carbon atoms. Irrespective of the direct experimental observation of such disclinated toroidal crystals, which is still open, we show here how continuum elasticity predicts that a 5−fold symmetric lattice indeed constitutes a minimum of the elastic energy for a broad range of aspect ratios and defect core energies.

For small aspect ratios the 5−fold symmetric configuration becomes unstable and is replaced by the 9−fold symmetric phase T 9 . As we mentioned, however, this configuration doesn't correspond to a possible triangulation of the torus. It is likely that the ground sate in this regime consists of ten skew disclination pairs as in the antiprismatic TAn lattice. The latter can be described by introducing a further degree of freedom δψ representing the angular displacement of defects from the equatorial plane:

A comparison of the TP5 configuration and the TA5 configuration is shown in Fig. 52 for different values of δψ. The intersection points of the boundary curves with the δψ−axis has been calculated by extrapolating the (r, δψ) data points in the range δ ∈ [0.07, 0.8] with ∆(δψ) = 2.5π · 10 −3 . For small δψ and r ∈ [3.3, 7.5] the prismatic TP5 configuration is energetically favored. For r < 3.3, however, the lattice undergoes a structural transition to the TA5 phase. For r > 7.5 the prismatic symmetry of the TP5 configuration breaks down again. In this regime, however, the elastic energy of both configurations rapidly rises because of the lower curvature and defects disappear.

In the regime of large particle numbers, the amount of curvature required to screen the stress field of an isolated disclination in units of lattice spacing becomes too large and disclinations are unstable to grain boundary "scars" consisting of a linear array of tightly bound 5 − 7 pairs radiating from an unpaired disclination [73, 142] . In a manifold with variable Gaussian curvature this effects leads to a regime of coexistence of isolated disclinations (in regions of large curvature) and scars. In the case of the torus the Gaussian curvature inside (|ψ| > π/2) is always larger in magnitude than that outside (|ψ| < π/2) for any aspect ratio and so we may expect a regime in which the negative internal curvature is still large enough to support the existence of isolated 7−fold disclinations, while on the exterior of the torus disclinations are delocalized in the form of positively charged grain boundary scars.

This hypothesis can be checked by comparing the energy of the TP5 lattice previously described with that of "scarred" configurations obtained by decorating the original toroid in such a way that each +1 disclination on the external equator is replaced by a 5−7−5 mini-scar. The result of this comparison is summarized in the phase diagram of Fig. 53 in terms of r and the number of vertices of the triangular lattice V (the corresponding hexagonal lattice has twice the number of vertices, i.e. V hex = 2V ). V can be derived from the angular separation of neighboring disclinations in the same scar by approximating V ≈ A/A V , with A V = √ 3 2 a 2 the area of a hexagonal Voronoi cell and a the lattice spacing. When the aspect ratio is increased from 1 to 6.8 the range of the curvature screening becomes shorter and the number of subunits required to destroy the stability of the TP5 lattice decreases. For r > 6.8, however, the geodesic distance between the two equators of the torus becomes too small and the repulsion between like-sign defects takes over. Thus the trend is inverted.

Some example of toroidal lattices obtained from numerical simulations are shown in Fig. 54 . The configurations are found by optimizing a system of V point-like particles constrained to lie on the surface of a torus and interacting via a pair potential of the form U ij = 1/|x i − x j | 3 where | · | denotes the Euclidean distance in R 3 . The choice of the cubic potential is motivated here by the so called "poppy seed bagel theorem" [80] , according to which the configuration of points that minimizes the Riesz energy E = i<j 1/|x i − x j | s on a rectifiable manifold of Hausdorff dimension d is uniformly distributed on the manifold for s ≥ d. In the case of a torus of revolution this implies that for small s the points are mostly distributed on the exterior of the torus (the interior becomes completely empty in the limit V → ∞). As s is increased, however, the points cover a progressively larger portion of the surface. The distribution becomes uniform for s ≥ 2. On the other hand, since the number of local minima of the Riesz energy increases with s, it is practical to choose a value not much larger than two. The choice s = 3 has the further advantage of modelling a real physical system of neutral colloidal particles assembled at an interface [36] and is therefore suitable for direct comparison with experiments on colloidal suspensions. The lowest energies found, as well as the number of defects in the corresponding configuration, are reported in Table 7 . The complete set of data produced in these simulations together with a collection of interactive 3D graphics for each low energy configuration studied can be found on-line [196] .

We have seen that disclination defects, forbidden in the lowest energy state of a planar crystal, may be energetically favored on a substrate of non-vanishing Gaussian curvature. It is therefore natural to ask whether large curvature can completely destroy crystalline order by driving the proliferation of a sufficiently high density of defects. The resulting state would be amorphous. The problem of generating amorphous structures by tiling a two-dimensional curved space with identical rigid subunits has drawn attention over the years, particularly through the connection to the structure of such disordered materials as supercooled liquids and metallic glasses. Since the work of Frank [197] the notion of geometrical frustration arises frequently in investigations of supercooled liquids and the glass transition. A paradigmatic example is represented by the icosahedral order in metallic liquids and glasses which, although locally favored, cannot propagate throughout all of threedimensional Euclidean space. A two-dimensional analog, consisting of a liquid of monodisperse hard disks in a 2-manifold of constant negative Gaussian curvature (the hyperbolic plane) was first proposed by Nelson and coworkers in 1983 [17, 198] . In such a system the impossibility of covering the entire manifold with a 6-fold coordinated array of disks mimics many aspects of the geometrical frustration of icosahedral order in three dimensions. In all these models of geometrical frustration, however, the origin of the disorder is primarily due to the short-range nature of the potential between the subunits. In a more realistic setting, part of the frustration is relieved by the fact that hexagonal unit cells can compress in order to match the underlying geometry.

The embedding of a triangular lattice on an axisymmetric torus, provides a particularly suitable playground to study curvature-driven disorder. When r → 1 the Gaussian curvature on the inside of the torus grows like 1/(r − 1) and diverges on the internal equator at ψ = π. We thus expect a high density of defects in the vicinity of the curvature singularity and a resultant loss of the local 6−fold bond Table 7 . Low energy configuration for a selected number of toroidal lattices with aspect ratios r = 3, 4, 6 and 20. For each aspect ratio the table displays the number of particles V , the lowest energy found and the number of k−fold vertices V k with k = 4-8. orientational order. In this regime the system will have crystalline regions on the outside of the torus and amorphous regions near the curvature singularity.

In this section we this claim can be supported using the elastic theory of continuous distributions of edge dislocations on a "fat" torus. Our argument is based on the following construction. As a consequence of the curvature singularity the surface area of an arbitrary wedge of angular width ∆φ becomes smaller and smaller as the sectional angle ψ increases and vanishes at ψ = π. If a defect-free lattice is embedded on such a wedge, Bragg rows will become closer and closer as the singularity is approached with a consequent rise in the elastic energy (see Fig. 55 ). An intuitive way to reduce the distortion of the lattice is to recursively remove Bragg rows as one approaches the point ψ = π (see Fig. 56 ). This is equivalent to introducing a growing density of edge dislocations. This dislocation "cloud" will ultimately disorder the system by destroying the local 6−fold bond orientational order. One might therefore view the curvature as playing the role of a local effective temperature which can drive "melting" by liberating disclinations and dislocations. In two-dimensional non-Euclidean crystals at T = 0, however, the mechanism for dislocation proliferation is fundamentally different from the usual thermal melting. While the latter is governed by an entropy gain due to unbinding of dislocation pairs, the amorphization at T = 0 is due to the adjustment of the lattice to the geometry of the embedding manifold via the proliferation of defects and the consequent release of elastic stress. A similar phenomenon occurs in the disorder-driven amorphization of vortex lattices in type-II and high-T c superconductors [199, 200, 201] . Since the shrinking area per plaquette on the inside of the torus necessitates a high density of dislocations we may approximate the dislocation cloud in this region by a continuous distribution of Burgers vector density b. Minimizing the elastic energy with respect to b yields a variational equation from which the optimal dislocation density can be calculated as a function of the ratio d /(Y R 2 ) between the dislocation core energy d and elastic energy scale Y R 2 with R = R 1 = R 2 .

As a starting point, we calculate the Green function G L (x, y) in the fat torus limit r → 1 (i.e. κ → ∞ and ω → 0, see Eq. (174)). The conformal angle ξ in this limit is lim r→1 ξ = tan ψ 2 ,

and to leading order of κ we have

To handle the limit of the Jacobi theta function we can take u = ∆z/κ, q = e iπτ = e − 2π κ and calculate the limit q → 1. This can be done by using the modular transformation properties of Jacobi functions [165, 166] : .

Taking the logarithm and neglecting irrelevant constant terms, we obtain

which finally leads to

with z = tan(ψ/2) + iφ. With the Green function in hand, we can calculate the effect of the curvature singularity at ψ = π on the distribution of defects. Let b be the Burgers vector density of the dislocation cloud. Hereafter we work in a local frame, so that

with g i = ∂ i R a basis vector in the tangent plane of the torus. R is a threedimensional surface vector parameterizing the torus and its Cartesian components are given in Eq. (165) . The quantity b has to be such that

with b D the total Burger's vector in a generic domain D. Because on a closed manifold dislocation lines cannot terminate on the boundary, extending the integration to the whole torus we have:

Since the basis vectors g i in Eq. (220) have the dimension of length, contravariant coordinates b i have dimensions of an inverse area. Assuming all defects to be paired in the form of dislocations (i.e. q i = 0 everywhere), the total energy of the crystal reads

where d is the dislocation core energy. The function Γ(x) encoding the elastic stress due to the curvature and the screening contribution of the dislocation cloud obeys

where ∇ i is the usual covariant derivative along the coordinate-direction i and i k is the Levi-Civita antisymmetric tensor on the torus:

The stress function Γ(x) can be expressed in the form Γ(x) = Γ d (ψ, φ) − Γ s (ψ) with Γ s (ψ) Y = log 1 2(1 + cos ψ)

The integral term in Eq. (224b) can now be reduced to a more friendly functional of b, suitable for a variational approach. Given the azimuthal symmetry we assume that all dislocations are aligned along b = b φ g φ .

Γ d (ψ) Y = 1 2π π −π dψ √ g b φ (ψ ) ψ + sin ψ + π sgn(ψ − ψ ) .

Substituting Eq. (225) and (224a) in Eq. (222) and minimizing with respect to b φ , the problem can be converted, after some algebraic manipulations, to the following Fredholm equation of the second kind:

where λ = d /(Y R 2 ), B(ψ) = R 2 (1 + cos ψ) 2 b φ (ψ) and the kernel K(ψ, ψ ) is given by:

K(ψ, ψ ) = 1 4(1 + cos ψ ) π −1 (ψ + sin ψ )(ψ + sin ψ) + |ψ − ψ | + 2 cos ψ + ψ 2 sin |ψ − ψ | 2 . 

where Cl 2 is the Clausen function (see Ref. [202] , pp. 1005-1006) defined as

As previously noted the dislocation core energy is d is much smaller than the elastic energy scale Y R 2 . Eq (226) is then suitable to be solved in powers of the dimensionless number λ: B(ψ) = B 0 (ψ) + λB 1 (ψ) + λ 2 B 2 (ψ) + · · ·

The corrections to the zero-order term B 0 (ψ) can be calculated recursively by solving a set of Fredholm equations of the first kind:

The function B 0 (ψ) associated with the Burgers vector density of the dislocation cloud in the limit λ → 0, on the other hand, can be calculated directly from Eq. (223) by setting the effective topological charge density on the right hand side to zero:

from which one obtains B 0 (ψ) = sin ψ and b φ = sin ψ R 2 (1 + cos ψ) 2 .

(230)

The Burgers vector density b φ obtained by a numerical solution of Eq. (226) is shown in Fig. 57 for different values of λ. The Burgers vector density is measured in units of R −2 . The function b φ has cubic singularities at ψ = ±π and is approximately zero on the outside of the torus. The solid blue curve in Fig. 57 represents the zeroth order solution of Eq.(230). Now, in the theory of dislocation mediated melting a system at the solid liquid phase boundary is described as a crystalline solid saturated with dislocations. In three-dimensions, in particular, there is a strong experimental evidence of the existence of a critical dislocation density at the melting point ρ(T m ) ≈ 0.6b −2 where b is the length of the length of the smallest perfect-dislocation Burgers vector [203] . Several theoretical works have motivated this evidence both for three-dimensional solids and vortex lattices in super conductors [200] . On the other hand, given the existence of such a critical density, its value can be empirically used to determine whether a system is in a solid or liquid-like phase in the same spirit as the Lindemann criterion. With this goal in mind we can calculate the dislocation density by requiring |b| = ρa with ρ the density of single lattice spacing dislocations. This yields:

Solving ρ(ψ, V )a 2 = 0.6 as a function of ψ and V we obtain the diagram of Fig. 58 .

As expected the inside of the torus contains an amorphous region whose angular size decreases with the number of vertices V as a consequence of the reduction of the lattice spacing.

In this article we have reviewed recent progress toward understanding the ground state properties of two-dimensional ordered phases on substrate of non-zero Gaussian curvature. Curved phases of matter are found in a broad class of systems comprising materials as different as viral capsids and carbon nanotori. The geometry of the substrate gives rise to novel defective structures that would be energetically prohibitive in the ground state of conventional flat systems. Defects may appear showing the two polar defects that must exist to allow the alternation of concentric rings. (Right) TEM images of a chain obtained when 11-mercaptoundecanoic acid functionalized nanoparticles are reacted with a water phase containing divalent 1,6-diaminohexane in a two-phase reaction. Scale bar 50 nm. From [204] .

in curved space either because they are required by the topology of the underlying substrate (as in the case of spherical crystals and nematics) or because they are favored by the curvature itself. The purely curvature driven energetic nature of the latter phenomenon is of fundamental interest. On the material science side, non-Euclidean systems provides a promising route to constructing arrays of nanoparticles via the chemical functionalization of topological defects created on the surface of the particles by coating them with an ordered monolayer. Candidates for such coatings include triblock coplymers, gemini lipids, metallic or semiconducting nanorods and conventional liquid crystal compounds. If the induced order is vector-like or striped, then there must necessarily be two defective regions corresponding to the bald spots exhibited by a combed sphere or the source and sink of fluid flow confined to a sphere. Not only are the defects physically and mathematically distinguished, since the condensed matter order vanishes there, but it also turns out that chemical processes can detect the defects and insert linker molecules at the precise defect locations. This creates the possibility of efficiently making functionalized nanoparticles with a precise valence and corresponding directional bonding. The nanoparticles themselves can then be linked into well-defined two and three dimensional arrays. Control over the valence and the geometry of the directional bonding can be achieved by varying the nature of the ordered monolayer.

This scheme has been demonstrated recently in the work of DeVries et al. [204] in which gold nanoparticles are coated with two species of naturally phase separating ligands (see Fig. 59 ). Phase separation on the spherical gold surface translates to a striped arrangement of alternating ligands. Such a spherical smectic has topological defects at the north and south pole of the nanopartcle. These may be functionalized by attaching thiol-based linker molecules at the defects via place-exchange reactions. This creates a divalent gold nanoparticle with directional bonds 180 degrees apart. These in turn can be linked to create polymers and free-standing films.

Alembert Essai d'une nouvelle théorie de la résistance des fluides

Statistical Mechanics of Membranes and Surfaces, 2nd World Scientific

Physics of amphiphilic layers

Proc. Natl. Acad. Sci

VIPERdb2: an enhanced and web API enabled relational database for structural virology

Dresselhaus Physical properties of carbon nanotubes

Sadoc and R. Mosseri Geometrical Frustration

Principles of Condensed Matter Physics

Nabarro Theory of Crystal Dislocations

Lothe Theory of Dislocations

On the origin of number and arrangement of the places of exit on pollengrains

The pursuit of perfect packing

Sloane Sphere packings, lattices and groups

Solid State Comm

Proc. Natl. Acad. Sci. USA 96

Nelson Defects and Geometry in Condensed Matter Physics

Mermin Solid state physics

Differential Geometry

Proc. Natl. Acad. Sci. USA

Differential geometry of curves and surfaces

Spiegel Applied Differential Equations

String Theory

The Coronaviridae Plenum Press

Lavrentovich Soft Matter Physics, an Introduction

Toroidal Single Wall Carbon Nanotubes in Fullerene Crop Circles

Toroidal lattices database

Stegun Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables

We thank V. Vitelli, A. Hexemer, W. Irvine, P. Chaikin, A. Lucas and A. Fonseca for sharing with us both images and data presented in this review. We thank D. R. Nelson, H. Shin and A. Travesset for long-standing collaborations and M. Forstner for discussions of amphiphilic membranes.