Defective nematogenesis: Gauss curvature in programmable shape-responsive sheets with topological defects

Flat sheets encoded with patterns of contraction/elongation morph into curved surfaces. If the surfaces bear Gauss curvature, the resulting actuation can be strong and powerful. We deploy the Gauss-Bonnet theorem to deduce the Gauss curvature encoded in a pattern of uniform-magnitude contraction/elongation with spatially varying direction, as is commonly implemented in patterned liquid crystal elastomers. This approach reveals two fundamentally distinct contributions: a structural curvature which depends on the precise form of the pattern, and a topological curvature generated by defects in the contractile direction. These curvatures grow as different functions the contraction/elongation magnitude, explaining the apparent contradiction between previous calculations for simple +1 defects, and smooth defect-free patterns. We verify these structural and topological contributions by conducting numerical shell calculations on sheets encoded with simple higher-order contractile defects to reveal their activated morphology. Finally we calculate the Gauss curvature generated by patterns with spatially varying magnitude and direction, which leads to additional magnitude gradient contributions to the structural term. We anticipate this form will be useful whenever magnitude and direction are natural variables, including in describing the contraction of a muscle along its patterned fiber direction, or a tissue growing by elongating its cells.


Introduction
Differential growth and differential muscular contraction underpin many impressive and important shape transformations found in biology, 1 ranging from embryonic gastrulation and limb formation to gut peristalsis and the beat of a heart. The key idea is that, with the right spatial patterning, a limited set of local movesisotropic growth, directional growth, or uniaxial muscular contraction -can be choreographed into a complex and reliable global shape transformation. Correspondingly, there is considerable interest in developing spatially-programmable shape-changing materials, 2,3 both to elucidate the basic geometric and mechanical principles of differential shape change, [4][5][6] and for use in soft machines 7,8 and deployable structures. 9 Several different programmable shape-changing materials have been developed, each responding to different stimuli and offering a different palette of local shape changes. Recent work has highlighted patterned (isotropic) swelling in hydrogels, [10][11][12][13] patterned contraction/elongation in liquid crystaline elastomers/ glasses (LCE/Gs) subject to heat or light, [14][15][16][17] and patterned contraction in ''baromorphs'' subject to inflation. 18,19 A natural distinction thus arises between a patterned magnitude of a locally isotropic shape change, and a patterned direction of a fixed magnitude shape change (Fig. 1). Returning to biology, although evolution is surely able to pattern magnitudes, 1,20 most biological tissues are anisotropic and also exhibit patterned directions, be they patterns of elongational growth 21,22 or patterns of muscular contraction.
In each of the above responsive materials, shape change is programmed into an initially flat sheet, which, on activation, morphs into a curved surface. Strikingly, such transformations are often impossible with passive sheets, generating timeless struggles for tailors, mapmakers, architects, gift-wrappers and manufacturers. As immortalized in Gauss's Theorema Egregium, 23,24 the key difficulty is geometric: the Gauss curvature of a surface (calculated as the product of the two principal curvatures, K = 1/(R 1 R 2 )) cannot be modified without changing the in-surface distance between points, encoded by the surface's metric. Thus a flat sheet of paper can be bent into a cylinder, but cannot be coerced into a sphere. Active shape-changing sheets, like their biological counterparts, can side-step this geometric constraint as their programmed shape-change does indeed change the in-surface distance between points. 10,25,26 Directional growth offers a unique additional possibility: patterns with topological defects. Such directional defects are familiar from liquid nematics 27 and have inspired many patterns encoded in LCE/Gs. 26,28,29 Indeed, LCE sheets encoded with concentric circles of contraction, forming a +1 nematic defect, have become a ubiquitous test-case for new patterning techniques. 15,28,[30][31][32][33] Such sheets morph dramatically into cones ( Fig. 1) demonstrating that topological defects offer a way to produce sharp features with divergent Gauss curvature. 26 Conversely, Gauss's geometric coupling of bend and stretch make curved surfaces strong and mechanically favourable. For example, if a flat sheet is bent in one direction, it cannot then bend in the transverse direction without stretch. 34 This curvature-induced rigidity, which is familiar from corrugated cardboard and pizza consumption, 35 leads to the general result that, whilst flat sheets buckle in a pure bending mode at a characteristic compressive force F p t 3 , curved shells require both stretch and bend to buckle, leading to a much higher buckling force F p t 2 dictated by the harmonic mean of the bending and stretching moduli. [36][37][38] Consequently, when an active sheet morphs into a Gauss-curved surface, the transition is mechanically strong: the LCE defect cones will lift 2500 times their own weight 32 as they activate. Patterns of Gauss curvature are also an essential component of inverse design: in order to program a sheet to morph into a desired shape it is necessary (though not sufficient) that the programmed shape change produces the desired Gauss curvature. 11,16,28,[39][40][41] In this paper, we use the Gauss-Bonnet theorem to calculate the distribution of Gauss curvature in an initially flat sheet programmed with a directional pattern of shape change. Our calculation is greatly facilitated by the use of the natural coordinate system of a nematic field, recently introduced by Niv and Efrati. 42 The Gauss-Bonnet approach allows us to construct a single result that unifies previous work on distributed Gauss curvature 39,40,43,44 and sharp points with divergent Gauss curvature. 26,45 Our unified result reveals a delicate interplay between topology and spatial patterning in the resultant Gauss curvature, and allows us to calculate the Gauss curvature encoded in topological defects with charges other than +1. We extend and verify these results by conducting numerical shell calculations on sheets encoded with higher-order defects to clarify the full 3D form of the surface that emerges. Finally, we present an analogous result for sheets encoded with directional growth in which both the magnitude and direction vary spatially in plane. Although such patterning is yet to be demonstrated in liquid crystalline solids, we anticipate that it will soon (either via local crosslink density or via control of the imprinted order parameter), and has already been achieved by evolution in many different biological contexts.
2 Curvature induced in flat sheets with a programmed direction of contraction/elongation We first consider a planar sheet, in which each point is programmed with a planar direction n = (cos c,sin c) along which the material will contract/elongate by a factor of l 8 upon activation, accompanied by, a sympathetic change of l > in the perpendicular direction. In an LCE/G, n would correspond to the nematic director (alignment direction), and, as indicated in Fig. 1, the response on heating/illumination is a large contraction, 0.25 o l 8 o 1, while the lateral expansion, l > = l 8 Àn , is determined by the opto-thermal Poisson ratio, which is strictly n = 1/2 in incompressible elastomers, but can be as high as 2 in photo-glasses. In baromorphs, n is the programmed pneumatic channel direction, and on inflation the channel does not change in length, l 8 = 1, but ideally contracts laterally by l > = 2/p = 0.63. . ..
As shown in Fig. 2, we may also define an orthogonal dual director on the sheet with cc + p/2, such that n* = (Àsin c,cos c). Then, in the undistorted sheet, the infinitesimal length element dl = (dx,dy) has length dl 2 = dl T Idl, which, on activation, becomes dl A 2 = dl T (l 8 2 nn + l > 2 n*n*)dl dl T % adl.
The activated sheet (indicated via subscript A) must deform into a surface following the programmed metric, % a. Suitable  An initially flat patch of nematic elastomer encoded with a spatially varying director n (left) will, upon activation, contract everywhere along n and expand along n*, causing it to morph into a curved surface (right). The nematic pattern defines a natural orthogonal coordinate system for the patch, (u,v), with u-lines along n and v-lines along n*.
This journal is © The Royal Society of Chemistry 2020 Soft Matter spatial variation in n will generate a spatially varying metric bearing Gauss curvature, guaranteeing that the resultant surface will also be curved. Our task is to compute the intrinsic geometry of the resultant surface, as characterized by its geodesic and Gauss curvatures. The metric is invariant under n -Àn, which is an inevitable consequence of the quadrupolar nature of nematic order. However, even if the underlying order really were vectorial, the metric, as a quadratic form, would nevertheless have this quadrupolar symmetry. Thus nematic rather than vectorial patterns are the natural language of directional shape change, and nematic defects, including half integer defects, could be found in any directional shape-changing system.
Previous authors have applied the Theorema Egregium to directly calculate Gauss curvature from gradients of the metric, yielding 40,43,44 This form has been used successfully to design surfaces with uniform finite Gauss curvature (such as spherical caps) and arbitary surfaces of revolution, 46,47 and underpins recent work on inverse design. 40 However, if we attempt to apply the above derivative formula at a topological defect in the pattern we encounter a problem: the derivatives are divergent. To clarify this, let us first consider two well-understood +1 defect nematic director patterns: those with the director in concentric circles around a single point, and those with the director emanating radially from a single point, each shown in Fig. 1. If a disk is encoded with concentric circles and then activated, it forms a new surface where a circle with original radius r has new (in material) radius rl > 4 r, and new circumference 2prl 8 o 2pr. This geometric contradiction is resolved by the disk morphing into a cone, with half-angle sin À1 (l 8 /l > ). The sides of the conical surface are Gauss flat everywhere, K A = 0, in agreement with a direct application of eqn (1). However, at the tip of the cone the Gauss curvature is clearly infinite. Furthermore, this tip of infinite Gauss curvature makes a finite contribution to the integrated (aka total) curvature, Ð K A dA A ¼ 2p 1 À l k l ? À Á 4 0, as may be computed by regularizing the tip with a spherical-cap: the cone's Gauss curvature is a delta function at the tip, with an infinite value but a finite integral. In contrast, the radial defect pattern creates a surface where the circumferences are too long for the radii, which buckles out of plane into an ''anticone'' shape resembling a saddle. This surface is also Gauss flat, except for an infinite negative curvature at the center which, since this is equivalent to exchanging l 8 and l > , gives a finite integrated curvature 26,48 With suitable care, differentiation can yield delta functions; for example in 3D electrostatics we routinely write r 2 r À1 = À4pd 3 (r). Thus one might reasonably hope to derive these cone/anticone tip curvatures via a careful application of eqn (1). However, any such attempt is doomed to fail. Eqn (1) has the general property that the Gauss curvature encoded by an orthogonal dual pattern (cc + p/2) is the negative of that of the original pattern, K A -ÀK A . In contrast, although the cone and anticone are orthogonal duals, and do have opposite signs of integrated Gauss curvature, the magnitudes are not equal.
A new approach is provided by the Gauss-Bonnet formula, 24 which relates the integrated Gauss curvature over a surface S, the geodesic curvature (k g ) of the surface's boundary, and the topological classification of the region via its Euler Characteristic w(S): ð S KdA þ ð @S k g ds ¼ 2pwðSÞ: Geodesic curvature, like Gauss curvature, is an intrinsic property of a path on a surface, and can be computed directly from the metric. We may use Gauss-Bonnet to calculate Ð KdA for any patch of surface from the geodesic curvature of the patch boundary. Gauss-Bonnet thus provides a rigorous basis for the assignation of finite integrated curvature to the tip of an (anti)cone, since it can be inferred from the finite geodesic curvature of an encircling path. Furthermore, we may slit a cone along a straight generator and unroll it isometrically into the plane to form a sector of a disk. Applying Gauss-Bonnet in the unrolled state, we then see that the tip's finite integrated curvature is exactly equal to the angular deficit of the sector, which becomes an angular surplus in the anticone case.
In general, any such point with finite integrated curvature is associated with a sharp (aka singular) point in the surface with a discontinuous surface normal. Unlike many other sharp points, such as those created in paper-folding origami, these sharp points cannot be removed by isometric bending deformations, as integrated curvature is a property of the metric alone. We thus describe them as intrinsically sharp points of the surface, and interpret the finite integrated curvature as an angular deficit/surplus associated with the point itself. Away from such points, where K is finite, we may apply Gauss-Bonnet to an infinitesimal patch (over which K is effectively constant) to infer the value of K itself. We will term the finite K in such regions distributed Gauss curvature.
Returning to nematogenesis, our general approach is thus to consider a region of patterned nematic sheet, bounded by a curve along which the nematic director is smooth, compute the geodesic curvature along this path, and then use Gauss-Bonnet to deduce the integrated Gauss curvature within. This approach will work even if the region contains topological defects, and whether or not the activated region contains intrinsically sharp points. Since we may choose any patch of surface, we will thus be able to identify the value of K A in regions with distributed finite curvature, and the finite integrated curvature Ð K A dA A at intrinsically sharp points where K A itself is infinite.

Geodesic curvature
The first step is to compute the geodesic curvature in the activated surface, of a path r(l) = (x(l),y(l)), defined in the flat sheet. In general this is a very involved calculation, but given an orthogonal coordinate system (u,v), in which the metric will take the simplified form diag (E,G), the geodesic curvature for

Soft Matter
This journal is © The Royal Society of Chemistry 2020 the unit-speed curve (u(l),v(l)) takes the simple form ( 24 p. 296) sometimes known as Liouville's formula: Above, f is the angle between the curve tangent and the u-line tangent in the tangent plane of the surface (as shown in Fig. 2). Following Niv and Efrati, 42 we thus move into the nematic pattern's natural coordinate system (u,v) in which u-lines (i.e., v = const.) are tangent to n and v lines (u = const.) are tangent to n* (Fig. 2). Since this parameterization is locally orthogonal, the metric of the flat (non-activated) sheet will be diagonal, where the scalar fields Z(u,v) and b(u,v) are computed from the geometry of the pattern in question. 42 Furthermore, given the principal stretches on activation are also along n and n*, the u and v lines will also be orthogonal on the deformed surface, and the (diagonal) activated metric is simply 40 These activated u and v lines thus also describe a nematic field (and its orthogonal dual) on the deformed surface. Given the activating deformation is continuous, this nematic field has the same topology as the reference state field, including any topological defects. A direct application of eqn (2) reveals that the integrated geodesic curvature along a path (u(l),v(l)) in the activated surface is Examining Fig. 2, the angle f between the path and a u-line in the unactivated flat (aka reference) state is given by This allows us to compute the f A term in the geodesic curvature as a reference-state integral The second term of eqn (4) can also be evaluated in the reference state by recognizing that dl = Z du n + b dv n*: Along a u line or a v line we have f A 0 = 0, as the angle between the curve tangent and the coordinate line is fixed. Setting l 8 = l > = 1 to interrogate the non-activated flat sheet, we see that the curvature of a u-line is ÀZ v /(Zb), and the curvature of the v-line is b u /(Zb). However, these curvatures are simply the 2D bend and splay (i.e., 2D curl and divergence) of the flat nematic pattern, 42 which allows us to eliminate the unknown fields Z and b from eqn (5), to get a coordinate-free result for the activated geodesic curvature along any path: This integrated geodesic curvature is appropriate for Gauss-Bonnet. To evaluate the local value, k gA one may equate integrands of eqn (8), whilst accounting for the infinitesimal

Gauss curvature away from topological defects
We now apply Gauss-Bonnet to a closed loop built from alternating u and v line segments. Along each segment df A = df = 0, so the only contribution to k gA comes from the final term in eqn (8), whilst at each corner there is a concentrated contribution Ð df A ¼ AEp=2, since the u and v lines meet perpendicularly even in the activated surface. In the simplest case, we combine two u lines and two v lines to form a quadrangular region, as shown in Fig. 2, giving four +p/2 contributions in the Gauss-Bonnet formula, which, given the region is topologically a disk with w = 1, exactly cancel the 2pw. ð The value of this integral is unchanged by an equal rotation of each vector. Taking a clockwise p/2 rotation, the line element rotates to point along the loop's outward normal, dlvdl, so the integral is re-cast as a boundary flux. This rotation also transforms n -Àn* and n*n, yielding: We further simplify this expression by introducing the bend and splay vectors, b = bn* and s = sn, respectively the curvature vector of u lines and minus the curvature vector of v lines. Unlike the director fields themselves, these are true nematic objects (invariant under n -Àn) and lead to a particularly natural expression: In a region without intrinsically sharp points, we may now apply the divergence theorem to identify the local value of the Gauss curvature Although this was derived for a quadrangular region bounded by u and v lines, it applies to any shape of region, as can be seen by tiling with infinitesimal quadrangles. Expanding the divergence, and recognizing that dA A = l 8 l > dA, we find that the local value of Gauss curvature is K A = l > 2 (Àb 2 + (n*Ár)b) + l 8 À2 (Às 2 À (nÁr)s), in exact agreement with 42 and. 40 Given the dependence on and l 8 and l > , this form doesn't immediately resemble eqn (1), although direct substitution of (smooth) b, s, n and n* reveals the two do indeed agree.

Gauss curvature with topological defects
However, the above analysis fails if there is a topological defect within the region because, as seen in Fig. 3, the simplest region is no longer quadrangular: it will require different numbers of sides depending on the topological charge, m, enclosed. Working in the reference domain, the tangent vector must wind by 2p around the closed loop. Of these, it will wind by 2pm in the u and v line segments, and hence the contribution from the corners must be 2p(1 À m). This introduces an additional 2pm into eqn (11) where the m i are the topological charges enclosed. This form can be simplified by noting that rc = bn + sn*, and hence 2pm ¼ so we can integrate the b À s component of our previous result, to get This formula is our main result for the Gauss curvature. We will call the first term the structural curvature, since it depends on the local structure of the pattern, while the second term is naturally called the topological curvature. The topological curvature contributes to any patch containing the defect, so it is a finite contribution to the integrated curvature at the point of the defect. Importantly, this does not imply that all defects with the same charge encode intrinsically sharp points with the same integrated curvature, as the structural term may also make a contribution potentially canceling or even reversing the topological one. If we are away from any defect, we can again apply the divergence theorem and equate integrands to compute the local distributed Gauss curvature. Recognizing that dA A = l 8 l > dA, this gives which manifestly has the same dependence on l 8 and l > as eqn (1), and matches exactly upon substitution for b and s. Alternatively, one may use the gradient result to bring all the topological terms into the integral, which leads to a particularly compact form for the Gauss curvature: In both cases, the new forms are equivalent to the original result for distributed Gauss curvature in the absence of defects, but also give the correct result for a region containing defects provided the director is smooth on the boundary.

Constant-rotation m = 1 defects
We first test our new form on the familiar constant-rotation +1 defect patterns 26 which, in plane polar coordinates (r,y) have the form c = y + a and hence rc = e y /r, n = cos(a)e r + sin(a)e y , and n* = Àsin(a)e r + cos(a)e y . The director curves thus form logspirals with constant angle a between the director and the radial direction, such that a = 0 is a radial defect and a = p/2 produces circles. In this case the bend and splay are given by b ¼ sin a r n Ã ; s ¼ cos a r n: In the particular case of a circular pattern, the splay vanishes and b = Àe r /r, while a radial pattern has vanishing bend and s = e r /r. Applying eqn (14) to a circular patch centered on the origin, we can calculate the integrated Gauss curvature as Since this is independent of the radius of the patch (and the patterns have rotational symmetry), we conclude that there is no distributed Gauss curvature, but there is an intrinsically sharp point at the origin with finite integrated curvature. Equivalently, we can reach this conclusion directly by applying the divergence theorem to eqn (14) and recalling that rÁ(e r /r) = 2pd 2 (r) and rÁ(e y /r) = 0, so rÁb = À2pd 2 (r)sin 2 a and rÁs = 2pd 2 (r)cos 2 a.
We can compute the same result from eqn (13) by noting Either way, setting a = 0, p/2 we recapitulate the familiar results for anticones/cones formed from radial/circumferential patterns respectively. However, for an intermediate angle, cos(2a c ) = (l 8 À l > )/(l 8 + l > ), the integrated curvature at the origin is zero, so there is no intrinsically sharp point in the activated surface; in fact, given the distributed curvature is also zero everywhere in the flanks, the activated nematic sheet will be completely flat. 26 The second approach clarifies that the integrated curvature at the origin contains a topological contribution, which is identical for all m = 1 defects, and a pattern-dependent structural contribution p cos(2a) which is not. The structural contribution reverses sign under taking the orthogonal dual, while the topological contribution does not, leading circular and radial patterns to have different but not exactly opposite curvature. Perhaps surprisingly, the structural term contributes a finite integrated curvature at the origin, despite the distributed curvature being zero throughout the flanks.

Integrated curvature at a general nematic defect
We next consider the neighborhood of a general defect with topological charge m. Constructing a polar coordinate system centered on the defect, the director will have the form c = y + a(y,r), where a is now a varying function that will wind m À 1 times on a path around the origin. Computing the bend and splay, we now get b ¼ sin a r 1 þ @a @y þ @a @r cos a s ¼ cos a r 1 þ @a @y À @a @r sin a; so, defining a 0 (y) = a(y,0), near the defect we have Considering a small circular patch at radius r, we can again use eqn (13), to find the finite integrated curvature at the defect point, where the angle brackets indicate an angular average around the defect. This result is strongly reminiscent of eqn (16), and again is the sum of a structural part, which depends on the particular pattern of the defect, and a topological part which depends only on its charge. Given that the factor 1 À l ? l k 1 À l k l ?
is strictly negative, we conclude that positive charge defects are inclined to produce negative Gauss curvature, and vice versa. Furthermore, given À1 o hcos(2a 0 )i o 1, the structural contribution provides an identical range of curvatures in all defect charges, while the topological term produces an offset to this range proportional to the defect charge. In particular, a higher-order defect will have zero integrated curvature, and hence will not generate an intrinsically sharp point in the surface, if which, for a fixed value of l 8 /l > , will be possible to satisfy over a range of m centered on zero.

Constant-speed higher-order defects
The simplest manifestations of higher-order defects are those with constant rotational speed, which, in plane polar coordinates (r,y) have the form c = my + g, and hence a(y,r) = a 0 (y) = (m À 1)y + g. These patterns are also of particular interest because they have been experimentally realized. 29 In this case we have simply Away from the origin we have distributed Gauss curvature which arises, via the structural curvature, from the divergence of this term: This distributed curvature vanishes for m = 1, as is familiar for log-spiral patterns. In contrast, all other constant-speed defects produce surfaces with distributed Gauss curvature K A p cos(2a) = cos(2(m À 1)y + 2g), in agreement with the calculations in ref. 43 and 45.

Paper Soft Matter
Open Access Article. Published on 01 October 2020. Downloaded on 11/9/2020 6:13:20 PM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

View Article Online
This journal is © The Royal Society of Chemistry 2020

Soft Matter
The integrated curvature at the defect point itself is characterized by eqn (20). Along an infinitesimal circle around the defect we have hcos(2a 0 )i = 0 for all m a 1, so the structural term vanishes leaving only the topological term. For m = 1 the angle a is constant, so hcos(2a 0 )i = cos(2a) and both terms contribute. Within the circle we thus have where the Kronecker d 1m equals one when m = 1 and zero otherwise. Constant speed defects thus always encode intrinsically sharp points, with the sole exception of m = 1 defects at the critical angle a c discussed previously. This result highlights a profound difference between m = 1 and m a 1 constant-speed defects. For m = 1, the structural curvature gives zero distributed curvature, but produces a finite contribution to the integrated curvature at the point of the defect. For m a 1 the structural curvature encodes a pattern of distributed Gauss curvature that diverges in strength towards the defect point, but makes no contribution to the integrated curvature at the point itself because it cancels under integration due to its angular variation. A (heuristic) differential perspective on this difference is that taking the divergence in eqn (22) also produces a singular term at the origin, (l > À2 À l 8

À2
)pm cos(2a)d 2 (r). Then, integrating the combined (singular + distributed) divergence over an infinitesimal disk, the cos(2a) variation causes the resulting structural curvature to vanish for all m a 1; but the singular term leads to a finite structural contribution for m = 1. This reflects a key symmetry difference in the patterns: adding an equal increment to a everywhere fundamentally changes the structure of an m = 1 pattern, but simply rotates an m a 1 pattern about the origin. Correspondingly, such an increment may change the integrated curvature at the defect point for m = 1, but not for m a 1.

Numerical calculation of activated surfaces for constant speed higher-order defects
Constant-speed higher-order defects thus offer a clear test of our main result, since the structural term accounts for the distributed curvature while the topological term accounts entirely for the finite integrated curvature at the defect point, and hence the intrinsically sharp point in the activated surface. We thus conduct numerical calculations using a bespoke shell-elasticity code (MorphoShell) to independently quantify their Gauss curvature and reveal their full activated shapes. Our approach closely follows the spirit of ''non-Euclidian shells'' with programmed metrics, 5 including both stretching and bending energies, with some specializations for incompressible rubber sheets.
In more detail, our elastic shell calculations describe the current shape of the sheet via the (3D) position x(r) of the material originally at the (2D) position r in the flat unactivated state. Using Latin indices for 3D and Greek for 2D, one can then compute the local deformation gradient F ia = q a x i , and hence the local metric of the deformed sheet as dx T dx = dr T F T ÁFdr dr T adr.
If this current metric differs from the activated target metric % a = (l 8 2 nn + l > 2 n*n*) then the sheet must pay a stretching energy which is simply the familiar stretching energy for a thin membrane of incompressible Neo-Hookean elastomer, with shear modulus m and (unactivated) thickness t. This energy is minimized (and vanishes) if the sheet achieves an isometry, a = % a. However, minimizing this stretching energy alone quickly reveals infinite numbers of non-smooth and unphysical isometries. To obtain sensible and physical surfaces, one must also include a sub-dominant bending energy, that penalizes extrinsic curvatures, encoded via the surface's second fundamental form, B ab = (q a q b x)ÁN, where N is the current surface normal. Following, 5 we use the standard bending energy for an incompressible sheet, which is minimized when B = 0, (i.e., flat), as expected when the encoded metric does not vary through the thickness. MorphoShell represents the deforming sheet via an unstructured triangulated mesh. The current metric for each triangle, a, requires first derivatives and is estimated from the unique linear deformation that describes the current positions of the triangle's three nodes: a standard constant-strain finite-element approach. The current second fundamental form, B, requires second derivatives and hence is estimated from the unique quadratic deformation consistent with the positions of six 'patch' nodes close to the triangle's centroid. Full details of these estimates are given in the ESI. † The activated form of the surface is then given by minimizing the total energy, E S + E B , over current node positions via either damped Newtonian dynamics or gradient descent.
The resultant surfaces for a range of 10 higher-order constantspeed defects are shown in Fig. 4. Since the defects form complicated 3D surfaces, additional viewpoints are provided in Fig. S1 and supplementary videos M1-M10 (ESI †). Our calculations are for free floating films, and have a large (rubber-like) actuation strains, leading to dramatic high-amplitude shapes. The results clearly exhibit the azimuthally oscillating distributed Gauss curvature given by eqn (22), leading to increasingly flowery surfaces at higher-orders, and good qualitative agreement with the experiments in ref. 29. As familiar from arrays of cones, 15 many of the defect shells can pop into different stable configurations. The configurations shown in Fig. 4 were selected as likely global energy minima, and some examples of alternatives are given in Fig. S2 (ESI †).
To verify our predicted structural and topological contributions to the Gauss curvature, we first compute the Gauss curvature of the activated surfaces, K A (r,y), via the angular deficit at each node in the mesh (see ESI † for details). In Fig. 5 we plot the local activated Gauss curvature around a ring of fixed large reference radius, K A (194t,y), for the m = +5/2 defect. Since this radius is far from the intrinsically sharp point at the origin, the local Gauss curvature stems only from the distributed contribution, and shows excellent agreement with the magnitude and form predicted in eqn (22).

Soft Matter Paper
Open Access Article. Published on 01 October 2020. Downloaded on 11/9/2020 6:13:20 PM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Soft Matter
This journal is © The Royal Society of Chemistry 2020 However, since the topological curvature is concentrated in a delta function at the origin, it cannot be interrogated by plotting K A directly. Rather, in the spirit of Gauss-Bonnet, in Fig. 6 we plot the integrated Gauss curvature within a reference-state disk centered on the origin, Ð r 0 Ð 2p y¼0 K A dA A , as a function of integration disk radius. In a perfect isometry the only contribution to the integrand would be a delta function of concentrated curvature at the center, leading the integral to have a constant value, mp 1 À l ? l k 1 À l k l ?
, stemming entirely from the topological curvature. In practice, the bending energy prohibits a singular curvature at the origin, so this contribution is 'smeared out' over a finite non-isometric region around the origin. However, as one moves away from the origin, all the programmed curvature is indeed accounted for, and each integral asymptotes to the expected topological constant. We emphasize both the quantized nature of the ladder of asymptotes shown in Fig. 6, which reflects their topological origin, and also the different scalings with l 8 and l > of the curvatures in Fig. 5 and 6, which highlight the distinction between the topological and structural contributions to the curvature.
The smeared-out regions near the origin in Fig. 6 involve a stretch-bend trade-off, resulting in locally non-isometric deformations that are beyond the scope of this paper. The extent of the smeared region is doubtless p t (as explored for m = +1 cones in ref. 26), and would vanish in the vanishing thickness limit. In practice, the region extends a surprisingly large multiple of t in higher-order defects, which we attribute to the very high distributed curvatures near the origin in these samples. However, ultimately Gauss-Bonnet guarantees that all the curvature encoded at the origin must be accounted for, once the boundary of the integration disk is large enough that the surface is isometric at the boundary.

Gauss curvature induced by directional fields with varying magnitudes
Finally, we consider a programmed directional growth field in which l 8 and l > are also functions of position. Such programming Fig. 4 Computational results for the activated morphologies of constant-speed defects in nematic elastomer sheets. The unactivated sheets are thin disks of radius R and thickness t E R/100, as is typical for programmed nematic elastomer sheets. In all cases the programmed metric has l 8 = 0.75, l ? ¼ 1 ffiffiffiffi ffi l k p except for m = +5/2, where l 8 = 0.9 is depicted to avoid self-intersection. The m = +1 defect is programmed with a = p/4, to isolate the topological contribution to the Gauss curvature, resulting in a weak anticone and enabling comparison with the higher-order defects.   (23), forming a quantized ladder. Owing to finite-thickness effects, the curvature is not singularly concentrated at a point at the origin, but is spread over a finite core.

Discussion and conclusions
We have presented a simple result for the Gauss curvature developed by a flat sheet programmed with patterned directional growth. Our result is valid for all patterns of growth, even if they contain topological defects, and thus unifies previous (and superficially contradictory) results on continuously distributed Gauss curvature and singularly concentrated Gauss curvature at simple +1 defects, and extends these results to include all defect charges.
As seen most clearly in eqn (13), our result reveals a subtle interplay between the topology of defects and the precise structure of a given defect's realization, which in turn generates a subtle interplay between topological defects and intrinsically sharp points. The topological term in eqn (13) contributes a finite integrated curvature to the defect point, which depends only on the charge of the defect and the value of l 8 /l > . In particular, the topological term always gives negative contributions for positive charge defects and vice versa, and is invariant under taking the orthogonal dual of a pattern. In contrast, the structural term in eqn (13) can produce both distributed Gauss curvature, and contribute finite integrated curvature to points with discontinuous director. At topological defects, the director is discontinuous, and the structural term makes a finite contribution to the defect point. This structural contribution depends on the exact spatial form of the defect pattern, and is always inverted, K str -ÀK str , by taking the orthogonal dual of a pattern. As seen in eqn (20), the resultant integrated curvature at a point with a topological defect contains both structural and topological contributions, and can thus lie anywhere in a given range, centered on the topological term and spanned by the structural term. Generically the two contributions do not exactly cancel, so the defect point has a finite integrated curvature and encodes an intrinsically sharp point in the activated surface. However, it is possible to structure defects such that the two contributions exactly cancel. In this case the defect point carries zero integrated curvature and does not create an intrinsically sharp point in the activated surface.
Conversely, although all the intrinsically sharp points we have covered are generated at topological defects, it appears possible to program one via the structural curvature alone, at a point without a topological defect. A full treatment of the conditions required for such points lies beyond the scope of this paper. However, heuristically, the structural curvature within a patch is controlled by H b þ s ð ÞÁmdl, which may encode a finite value at any point with divergent b or s. Since bend and splay are the curvatures of director lines and orthogonal dual lines respectively, such a divergence certainly requires infinitely curved director/dual lines, and hence a discontinuity in director at the point. However, such a discontinuity would arise at any kink in a director curve, which need not carry a topological charge. An intrinsically sharp point thus requires a point of discontinuous director, but the discontinuity need not be a topological defect.
The key enabling approach to characterizing the curvature encoded at defects has been to use Gauss-Bonnet to evaluate the Gauss curvature within a patch in terms of its boundary properties, resulting in a formula for the Gauss curvature that needs integrating along the boundary but not over the (defectcontaining) area of the patch itself. In the particularly simple case where the patch boundary is constructed from u and v lines (integral curves of n and n*) we can go a step further and conduct the boundary integral along each segment to calculate the Gauss curvature in terms of properties at the corners. More precisely, we can conduct the integral in eqn (14) along each segment, which yields the change in angle of the path segment, so the Gauss curvature within is where Dc u and Dc v are the change in director angle (and hence change in path angle) along the u line and v line segments. The Gauss curvature within the patch is thus entirely determined from the director angles at the corners. For example, for the quadrangular path shown in Fig. 2, the Gauss curvature within is simply Going further, if the boundary is largely composed of u lines, the Gauss curvature within will be while if it is largely composed of just v lines it will be 1 À l ? l k 2pm: results that are familiar from radial and circumferential +1 defects, but now generalized to patches containing multiple higher-order defects with the correct shape of boundary. We anticipate that these geometric insights will be of use for designing patterns of growth to achieve particular surfaces. Current strategies for inverse design cannot produce patterns containing defects, 40 but our results highlight that defects not only offer a route to designing surfaces with sharp features, but can also occur without any concentrated Gauss curvature as part of the programming of a smoothly curved surface. Furthermore, our results show that the integrated Gauss curvature of a surface is intimately related to the total topological charge it contains, suggesting that the inclusion of topological defects may be crucial to the design of surfaces containing large integrated amounts of Gauss curvature, such as needed to wrap (or even double-wrap) a sphere.
Looking beyond our current results, we first note the possibility of encoding growth into an initially curved surface, so that it morphs into a different curved surface on activation. This situation is actually typical in biology, and, via ''4D'' printing, is increasingly accessible to engineers. 30 Although we defer a full investigation to a later date, we anticipate that our existing approach will generalize straightforwardly; indeed, as long as one uses nematic coordinates, the metric of an initially curved surface could be represented via eqn (3) without modification suggesting that many of our current results will be directly applicable. A second extension would be to compute Gauss curvature when the director is discontinuous along a seam, rather than a simple point defect. 45,49,50 Our current approach relies on the director being continuous on a patch boundary, then using Gauss-Bonnet to capture the curvature within. Our results thus apply directly to a seam that lies entirely within a patch (rather like a Kirigami cut) which can be circumnavigated without difficulty, and an appropriate m assigned from the winding of n around the boundary. However, if the seam passes through the patch boundary then director discontinuity also occurs on the boundary, and the concept of enclosed topological charge breaks down. As we will show in a forthcoming paper, in this case the Gauss-Bonnet approach still works, but yields rather different results.
Finally, we comment that many biological tissues also have nematic or polar order (stemming from the anisotropy of their cells), and there is growing appreciation of the physiological implications of defects: for example +1/2 and À1/2 defects in epithelial tissue are known to be drive cell shedding. 51 Similarly the patterns of directional elongation that generate both the exaggerated nectar spurs in Darwin's orchid 21 and the reproductive whorl in the green alga Acetabularia acetabulum 22,25 are centered on +1 defects, and these must be accounted for properly to compute the Gauss curvature and understand the shape programming. Unlike in LCEs, biological growth fields tend to be patterned in both direction and magnitude, but fortunately our results straightforwardly generalize to this case (eqn (29)). Furthermore, given the metric itself has quadrupolar symmetry, this formulation is not limited to tissues with nematic order, but encompasses the full spectrum of encodable metrics, and is likely to be useful in polar tissues with topological defects as well as their nematic counterparts.

Conflicts of interest
There are no conflicts of interest to declare.  Figure S1: Alternative views of the simulated surfaces in Fig. 4 of the main paper. Note that despite the 2(1 − m)-fold rotational symmetry of the m < 0 patterns, this symmetry is broken by the observed shapes, even though all the large 'lobes' popped the same way. The same is true for m = +5/2, and is expected to hold for greater m values too. See also the supplementary movies M1-10. Figure S2: Two other simulations of the +2 defect, showing that a variety of different shapes are possible if the sheet is allowed to 'pop' in different ways. These shapes had a higher residual energy than the +2 defect in Fig. S1, suggesting that they are local energy minima, while that in Fig. S1 is expected to be the global minimum. These might nonetheless be realised physically however.

S2 Computational details
Computations were carried out on disks of (reference-state) thickness t and radius R t, programmed with constant-speed defect metrics with λ ⊥ = 1/ λ . The disks were meshed with 2D triangles of edge length h t to resolve the largest curvatures. All meshes were created with the PyDistMesh implementation of the DistMesh algorithm [1]. The specific parameters for each disk are given in the table below; the finest mesh had > 1.8 million triangles. S1 Simulation(s) The initial states were perturbed slightly in the out-of-plane direction to encourage particular 'popping' patterns that are likely global energy minima. The shells are multistable and various other stable, differently-popped solutions were also encountered (see Fig. S2 above).
MorphoShell calculates the metric and second fundamental form for each triangle based on current node positions (see Section S3). The force on each node is then calculated as the gradient of the total energy ((24)+(25) in the main paper), and the nodes are moved accordingly via damped Newtonian dynamics or gradient descent. In our calculations equilibrium was deemed achieved when non-dimensionalized elastic force and speed both fell below 5 × 10 −5 . The MorphoShell code and associated files used to create Fig. 6  Radius R for each disk was chosen to be large enough that the integrated Gauss curvature attained its asymptotic far-field value away from the core (Fig. 6, main paper). Non-Euclidean shells are also known to have a non-isometric boundary layer within ∼ t of the outer perimeter, which would obscure the asymptotic behaviour, so the plots in Fig. 6 (main paper) are truncated before the boundary layer is reached. Furthermore the sheet was made thinner (by a factor ∼ 15) and stiffer within ∼ 10t of the boundary (without altering the bending modulus), to restrict the boundary layer to this region.

S2.1 Gauss curvature estimates for Figs 5, 6 of main paper
Figs. 5, 6 of the main paper required a robust estimate of Gauss curvature, on a triangulated mesh. We borrowed the well-known 'angle deficit' method from discrete differential geometry. This method is motivated by an exact discrete analogue of the Gauss-Bonnet theorem for a triangulated mesh surface S [2]: S2 where the sums are over the non-boundary and boundary nodes respectively, and is the sum of mesh-triangle interior angles {β ∆ } around the node i. Comparing with the continuum Gauss-Bonnet theorem, this provides a natural definition for the integrated Gauss curvature of an interior node in a triangulated mesh: the 'angle deficit', 2π − Ξ i . This is a widely used Gauss curvature measure [3], and is the measure used in Fig. 6 (main paper).
For Fig. 5 (main paper), the local Gauss curvature was also required. To obtain this, we simply divide each node's angle deficit by an area assigned to that node, taken to be A i = ∆ A ∆ /3, so that each triangle contributes 1/3 of its area to each vertex. Thus we obtain: This and similar estimates for K are not locally convergent in full generality [4]. However they are unproblematic for highly regular meshes such as ours, (hence their widespread use), and the estimate K i performed well in our case. MorphoShell evolves an unstructured triangulated mesh. Each node has a fixed 2D position vector (i) r in the initial (reference) state, and a dynamical 3D position vector (i) x in the deformed ('current') state. Figure S3: The reference and deformed states of a small piece of mesh.
The deformed metric for each triangle in the mesh is taken to be the unique 2x2 matrix mapping its initial squared side lengths to their deformed values.
We use a 'patch fitting' approach to estimate the second fundamental form, as is common in computer graphics. In contrast to some others [5], our estimate contains only node position S3 degrees of freedom.
For a given triangle T , six 'patch nodes' are determined: T 's three vertices, and three further nodes selected from nearby triangles. The deformed surface in the vicinity of the triangle is approximated by a 'patch' surface p(r) that is a quadratic function of Cartesian reference coordinates (r 1 , r 2 ): where ∆r 1 ≡ r 1 −r 1 and (r 1 ,r 2 ) are the reference coordinates of T 's centroid. This is nothing more than a truncated Taylor expansion for each component of p(r), valid as long as the surface's local radius of curvature is significantly larger than the element size.
Using B αβ = (∂ α ∂ β x) ·N , we combine Λ with the triangle face normalN to give an estimate S4 Figure S4: Sketch of the reference and deformed state for a triangle T and its six patch nodes, and the fitted patch surface p(r).
of the second fundamental form for triangle T : The three non-vertex patch nodes were simply chosen to be those closest to T 's centroid in the reference state, subject to M −T being suitably well conditioned.