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

Abstract

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 of 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.

---

1 Introduction

Differential growth and differential muscular contraction underpin many impressive and important shape transformations found in biology,¹ 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 moves – isotropic 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–6 and for use in soft machines^7,8 and deployable structures.⁹

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–13 patterned contraction/elongation in liquid crystaline elastomers/glasses (LCE/Gs) subject to heat or light,^14–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.

Fig. 1 A sheet of shape-shifting material contracts along a programmed direction (blue lines) on activation. If the direction is constant the sheet remains flat (top). Concentric rings of contraction morph a disk into a cone, while a radial contraction pattern makes an anticone.

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₁R₂)) 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²⁷ 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–33} Such sheets morph dramatically into cones (Fig. 1) demonstrating that topological defects offer a way to produce sharp features with divergent Gauss curvature.²⁶

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.³⁴ This curvature-induced rigidity, which is familiar from corrugated cardboard and pizza consumption,³⁵ leads to the general result that, whilst flat sheets buckle in a pure bending mode at a characteristic compressive force F ∝ t³, curved shells require both stretch and bend to buckle, leading to a much higher buckling force F ∝ t² dictated by the harmonic mean of the bending and stretching moduli.^36–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³² 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–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.⁴² 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 ψ,sin ψ) along which the material will contract/elongate by a factor of λ_‖ upon activation, accompanied by, a sympathetic change of λ_⊥ 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 < λ_‖ < 1, while the lateral expansion, λ_⊥ = λ_‖^−ν, is determined by the opto-thermal Poisson ratio, which is strictly ν = 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, λ_‖ = 1, but ideally contracts laterally by λ_⊥ = 2/π = 0.63….

As shown in Fig. 2, we may also define an orthogonal dual director on the sheet with ψ → ψ + π/2, such that n* = (−sin ψ,cos ψ). Then, in the undistorted sheet, the infinitesimal length element dl = (dx,dy) has length dl² = dl^TIdl, which, on activation, becomes

dl_A² = dl^T(λ_‖²nn + λ_⊥²n*n*)dl ≡ dl^Tādl.

Fig. 2 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*.

The activated sheet (indicated via subscript A) must deform into a surface following the programmed metric, ā. Suitable 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.⁴⁰

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 rλ_⊥ > r, and new circumference 2πrλ_‖ < 2πr. This geometric contradiction is resolved by the disk morphing into a cone, with half-angle sin⁻¹(λ_‖/λ_⊥). 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, , 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 λ_‖ and λ_⊥, gives a finite integrated curvature .^26,48

With suitable care, differentiation can yield delta functions; for example in 3D electrostatics we routinely write ∇²r⁻¹ = −4πδ³(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 (ψ → ψ + π/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,²⁴ 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 χ(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 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 at intrinsically sharp points where K_A itself is infinite.

2.1 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 the unit-speed curve (u(l),v(l)) takes the simple form (²⁴p. 296) sometimes known as Liouville's formula:Above, ϕ 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,⁴² 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 η(u,v) and β(u,v) are computed from the geometry of the pattern in question.⁴² 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⁴⁰

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 ϕ between the path and a u-line in the unactivated flat (aka reference) state is given by which, upon activation becomes . This allows us to compute the ϕ_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 = η du n + β dv n*:

Along a u line or a v line we have ϕ_A′ = 0, as the angle between the curve tangent and the coordinate line is fixed. Setting λ_‖ = λ_⊥ = 1 to interrogate the non-activated flat sheet, we see that the curvature of a u-line is −η_v/(ηβ), and the curvature of the v-line is β_u/(ηβ). However, these curvatures are simply the 2D bend and splay (i.e., 2D curl and divergence) of the flat nematic pattern,⁴²

which allows us to eliminate the unknown fields η and β 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 arc length, to get

2.2 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 dϕ_A = dϕ = 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 , 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 +π/2 contributions in the Gauss–Bonnet formula, which, given the region is topologically a disk with χ = 1, exactly cancel the 2πχ.

The value of this integral is unchanged by an equal rotation of each vector. Taking a clockwise π/2 rotation, the line element rotates to point along the loop's outward normal, dl → v̂dl, 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 = λ_‖λ_⊥dA, we find that the local value of Gauss curvature is K_A = λ_⊥²(−b² + (n*·∇)b) + λ_‖⁻²(−s² − (n·∇)s), in exact agreement with⁴² and.⁴⁰ Given the dependence on and λ_‖ and λ_⊥, this form doesn’t immediately resemble eqn (1), although direct substitution of (smooth) b, s, n and n* reveals the two do indeed agree.

2.3 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 2π around the closed loop. Of these, it will wind by 2πm in the u and v line segments, and hence the contribution from the corners must be 2π(1 − m). This introduces an additional 2πm into eqn (11) leading to
where the m_i are the topological charges enclosed. This form can be simplified by noting that ∇ψ = bn + sn*, and henceso we can integrate the b − s component of our previous result, to get

Fig. 3 Examples of defect patterns with topological charge ±1/2. Note that a path around the defect consisting of u and v segments must have different numbers of segments for each charge.

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 = λ_‖λ_⊥dA, this gives

which manifestly has the same dependence on λ_‖ and λ_⊥ 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.

3 Constant-rotation m = 1 defects

We first test our new form on the familiar constant-rotation +1 defect patterns²⁶ which, in plane polar coordinates (r,θ) have the form ψ = θ + α and hence ∇ψ = e_θ/r, n = cos(α)e_r + sin(α)e_θ, and n* = −sin(α)e_r + cos(α)e_θ. The director curves thus form log-spirals with constant angle α between the director and the radial direction, such that α = 0 is a radial defect and α = π/2 produces circles. In this case the bend and splay are given by

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 ∇·(e_r/r) = 2πδ²(r) and ∇·(e_θ/r) = 0, so ∇·b = −2πδ²(r)sin² α and ∇·s = 2πδ²(r)cos² α.

We can compute the same result from eqn (13) by noting that to get:

Either way, setting α = 0, π/2 we recapitulate the familiar results for anticones/cones formed from radial/circumferential patterns respectively. However, for an intermediate angle, cos(2α_c) = (λ_‖ − λ_⊥)/(λ_‖ + λ_⊥), 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.²⁶

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 ∝ cos(2α) 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.

4 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 ψ = θ + α(θ,r), where α is now a varying function that will wind m − 1 times on a path around the origin. Computing the bend and splay, we now getso, defining α₀(θ) = α(θ,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 is strictly negative, we conclude that positive charge defects are inclined to produce negative Gauss curvature, and vice versa. Furthermore, given −1 < 〈cos(2α₀)〉 < 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 λ_‖/λ_⊥, will be possible to satisfy over a range of m centered on zero.

4.1 Constant-speed higher-order defects

The simplest manifestations of higher-order defects are those with constant rotational speed, which, in plane polar coordinates (r,θ) have the form ψ = mθ + γ, and hence α(θ,r) = α₀(θ) = (m − 1)θ + γ. These patterns are also of particular interest because they have been experimentally realized.²⁹ 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 ∝ cos(2α) = cos(2(m − 1)θ + 2γ), in agreement with the calculations in ref. 43 and 45.

The integrated curvature at the defect point itself is characterized by eqn (20). Along an infinitesimal circle around the defect we have 〈cos(2α₀)〉 = 0 for all m ≠ 1, so the structural term vanishes leaving only the topological term. For m = 1 the angle α is constant, so 〈cos(2α₀)〉 = cos(2α) and both terms contribute. Within the circle we thus have

where the Kronecker δ_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 α_c discussed previously.

This result highlights a profound difference between m = 1 and m ≠ 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 ≠ 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, (λ_⊥⁻² − λ_‖⁻²)πm cos(2α)δ²(r). Then, integrating the combined (singular + distributed) divergence over an infinitesimal disk, the cos(2α) variation causes the resulting structural curvature to vanish for all m ≠ 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 α everywhere fundamentally changes the structure of an m = 1 pattern, but simply rotates an m ≠ 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 ≠ 1.

4.2 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,⁵ 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_iα = ∂_αx_i, and hence the local metric of the deformed sheet as dx^Tdx = dr^TF^T·Fdr ≡ dr^Tadr. If this current metric differs from the activated target metric ā = (λ_‖²nn + λ_⊥²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 μ and (unactivated) thickness t. This energy is minimized (and vanishes) if the sheet achieves an isometry, 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_αβ = (∂_α∂_βx)·N̂, where N̂ is the current surface normal. Following,⁵ 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 constant-speed 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,¹⁵ 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†).

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 ≈ R/100, as is typical for programmed nematic elastomer sheets. In all cases the programmed metric has λ_‖ = 0.75, except for m = +5/2, where λ_‖ = 0.9 is depicted to avoid self-intersection. The m = +1 defect is programmed with α = π/4, to isolate the topological contribution to the Gauss curvature, resulting in a weak anticone and enabling comparison with the higher-order defects.

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,θ), 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,θ), 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).

Fig. 5 Activated Gauss curvature as a function of θ at r ≈ 194t (far from the defect core), for an m = − 2 constant-speed defect with γ = π/4 and thus cos(2α) = sin(6θ). The values estimated numerically on a simulated surface show good agreement with eqn (22).

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, , 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, , 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 λ_‖ and λ_⊥ of the curvatures in Fig. 5 and 6, which highlight the distinction between the topological and structural contributions to the curvature.

Fig. 6 Integrated Gauss curvature (up to radius r) for 10 simulated defects. Each defect generates the topological curvature predicted by eqn (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.

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 ∝ 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.

5 Gauss curvature induced by directional fields with varying magnitudes

Finally, we consider a programmed directional growth field in which λ_‖ and λ_⊥ are also functions of position. Such programming has not yet been demonstrated in LCEs (although it could be achieved by additionally patterning the crosslink density or the strength of the order parameter), but is commonly found in biology: for example the exaggerated nectar spurs in Darwin's orchid develop via a process of cell elongation,²¹ but the degree of elongation varies with distance from the tip of the spur. Given the local metric now has three degrees of freedom, ψ, λ_‖ and λ_⊥, it is possible to represent any programmed metric in this form. However, this is likely to be a natural representation in cases where the direction n is physically meaningful, such as being the direction of cell polarization. In this case we can again adopt the nematic coordinate system, (u,v); the metric is still diagonal, but now λ_‖ and λ_⊥ are functions of u and v. Applying eqn (2), the geodesic curvature is thus

To bring this into a coordinate-independent form, we again use dl = η du n + β dv n* and recognize the (scalar) bend and splay, but must also recognize that and , which allows us to eliminate the unknown fields η and β and write

As an aside, we note that this can be neatened up to give

However, returning to eqn (27) and again applying Gauss–Bonnet to a patch built from u lines and v lines, and containing defects with charge m_i, we can compute the contained Gauss curvature as

If desired, eqn (12) can again be deployed to either bring the topological term inside the integral, or to evaluate the b − s component of the integral. The resulting expressions are identical to eqn (13) and (14), except with an additional contribution, , to the structural curvature arising from the varying magnitudes.

6 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 λ_‖/λ_⊥. 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 , 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 (defect-containing) 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 Δψ_u and Δψ_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 : 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,⁴⁰ 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.³⁰ 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 viaeqn (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.⁵¹ Similarly the patterns of directional elongation that generate both the exaggerated nectar spurs in Darwin's orchid²¹ 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.

Acknowledgements

J. S. B. is supported by a UKRI Future Leader Fellowship, grant no. MR/S017186/1. D. D. is supported by the EPSRC Centre for Doctoral Training in Computational Methods for Materials Science, grant no. EP/L015552/1, and thanks the creators of Distmesh.⁵²

Notes and references

1. D. W. Thompson, On growth and form, Cambridge Univ. Press, 1942.
2. M. A. McEvoy and N. Correll, Science, 2015, 347, 1261689.
3. M. Warner, Annu. Rev. Condens. Matter Phys., 2020, 11, 125–145.
4. M. B. Amar and P. Ciarletta, J. Mech. Phys. Solids, 2010, 58, 935–954.
5. E. Sharon and E. Efrati, Soft Matter, 2010, 6, 5693–5704.
6. J. Dervaux, Y. Couder, M.-A. Guedeau-Boudeville and M. B. Amar, Phys. Rev. Lett., 2011, 107, 018103.
7. B. Gao, Q. Yang, X. Zhao, G. Jin, Y. Ma and F. Xu, Trends Biotechnol., 2016, 34, 746–756.
8. M. Cianchetti, T. Ranzani, G. Gerboni, T. Nanayakkara, K. Althoefer, P. Dasgupta and A. Menciassi, Soft Robot., 2014, 1, 122–131.
9. W. Wang, H. Rodrigue and S.-H. Ahn, Sci. Rep., 2016, 6, 1–10.
10. Y. Klein, E. Efrati and E. Sharon, Science, 2007, 315, 1116–1120.
11. J. Kim, J. A. Hanna, M. Byun, C. D. Santangelo and R. C. Hayward, Science, 2012, 335, 1201–1205.
12. J.-H. Na, N. P. Bende, J. Bae, C. D. Santangelo and R. C. Hayward, Soft Matter, 2016, 12, 4985–4990.
13. A. S. Gladman, E. A. Matsumoto, R. G. Nuzzo, L. Mahadevan and J. A. Lewis, Nat. Mater., 2016, 15, 413–418.
14. L. T. de Haan, C. Sánchez-Somolinos, C. M. Bastiaansen, A. P. Schenning and D. J. Broer, Angew. Chem., Int. Ed., 2012, 51, 12469–12472.
15. T. H. Ware, M. E. McConney, J. J. Wie, V. P. Tondiglia and T. J. White, Science, 2015, 347, 982–984.
16. H. Aharoni, Y. Xia, X. Zhang, R. D. Kamien and S. Yang, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 7206–7211.
17. M. Barnes and R. Verduzco, Soft Matter, 2019, 15, 870–879.
18. E. Siéfert, E. Reyssat, J. Bico and B. Roman, Nat. Mater., 2019, 18, 24–28.
19. M. Warner and E. Siéfert, Proc. R. Soc. A, 2020, 476, 20200047.
20. G. Mitchison, J. Theor. Biol., 2016, 408, 155–166.
21. J. R. Puzey, S. J. Gerbode, S. A. Hodges, E. M. Kramer and L. Mahadevan, Proc. R. Soc. B, 2012, 279, 1640–1645.
22. K. A. Serikawa and D. F. Mandoli, Planta, 1998, 207, 96–104.
23. C. F. Gauss, Disquisitiones generales circa superficies curvas, Typis Dieterichianis, 1828, vol. 1.
24. B. O’neill, Elementary differential geometry, Academic Press, 2014.
25. J. Dervaux and M. B. Amar, Phys. Rev. Lett., 2008, 101, 068101.
26. C. Modes, K. Bhattacharya and M. Warner, Proc. R. Soc. A, 2011, 467, 1121–1140.
27. P.-G. De Gennes and J. Prost, The physics of liquid crystals, Oxford University Press, 1993, vol. 83.
28. M. A. Dias, J. A. Hanna and C. D. Santangelo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 036603.
29. M. E. McConney, A. Martinez, V. P. Tondiglia, K. M. Lee, D. Langley, I. I. Smalyukh and T. J. White, Adv. Mater., 2013, 25, 5880–5885.
30. C. P. Ambulo, J. J. Burroughs, J. M. Boothby, H. Kim, M. R. Shankar and T. H. Ware, ACS Appl. Mater. Interfaces, 2017, 9, 37332–37339.
31. A. Kotikian, R. L. Truby, J. W. Boley, T. J. White and J. A. Lewis, Adv. Mater., 2018, 30, 1706164.
32. T. Guin, M. J. Settle, B. A. Kowalski, A. D. Auguste, R. V. Beblo, G. W. Reich and T. J. White, Nat. Commun., 2018, 9, 1–7.
33. M. López-Valdeolivas, D. Liu, D. J. Broer and C. Sánchez-Somolinos, Macromol. Rapid Commun., 2018, 39, 1700710.
34. V. Pini, J. Ruz, P. M. Kosaka, O. Malvar, M. Calleja and J. Tamayo, Sci. Rep., 2016, 6, 1–6.
35. M. Taffetani, F. Box, A. Neveu and D. Vella, EPL, 2019, 127, 14001.
36. R. Zoelly, PhD thesis, ETH Zurich, 1915.
37. W. Koiter, PhD thesis, Delft University of Technology, 1945.
38. B. Audoly and Y. Pomeau, Peyresq Lectures On Nonlinear Phenomena, World Scientific, 2000, pp. 1–35.
39. C. Mostajeran, M. Warner, T. H. Ware and T. J. White, Proc. R. Soc. A, 2016, 472, 20160112.
40. I. Griniasty, H. Aharoni and E. Efrati, Phys. Rev. Lett., 2019, 123, 127801.
41. A. Lucantonio and A. DeSimone, Mech. Mater., 2020, 103313.
42. I. Niv and E. Efrati, Soft Matter, 2018, 14, 424–431.
43. H. Aharoni, E. Sharon and R. Kupferman, Phys. Rev. Lett., 2014, 113, 257801.
44. C. Mostajeran, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 062405.
45. C. D. Modes and M. Warner, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 021711.
46. M. Warner and C. Mostajeran, Proc. R. Soc. A, 2018, 474, 20170566.
47. B. A. Kowalski, C. Mostajeran, N. P. Godman, M. Warner and T. J. White, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2018, 97, 012504.
48. M. M. Müller, M. B. Amar and J. Guven, Phys. Rev. Lett., 2008, 101, 156104.
49. P. Plucinsky, B. A. Kowalski, T. J. White and K. Bhattacharya, Soft Matter, 2018, 14, 3127–3134.
50. F. Feng, J. S. Biggins and M. Warner, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2020, 102, 013003.
51. T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, Nature, 2017, 544, 212–216.
52. P.-O. Persson and G. Strang, SIAM Rev., 2004, 46, 329–345.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0sm01192d

---