Continuum mechanics of moving defects in growing bodies

Growth processes in many living organisms create thin, soft materials with an intrinsically hyperbolic geometry. These objects support novel types of mesoscopic defects - discontinuity lines for the second derivative and branch points - terminating defects for these line discontinuities. These higher-order defects move"easily", and thus confer a great degree of flexibility to thin hyperbolic elastic sheets. We develop a general, higher-order, continuum mechanical framework from which we can derive the dynamics of higher order defects in a thermodynamically consistent manner. We illustrate our framework by obtaining the explicit equations for the dynamics of branch points in an elastic body.


Introduction
Hyperbolic sheets abound in nature (see Fig. 1). As Margaret Wertheim writes in her delightful essay "Corals, crochet and the cosmos: how hyperbolic geometry pervades the universe" [Wer16] -We have built a world of largely straight lines -the houses we live in, the skyscrapers we work in and the streets we drive on our daily commutes. Yet outside our boxes, nature teems with frilly, crenellated forms, from the fluted surfaces of lettuces and fungi to the frilled skirts of sea slugs and the gorgeous undulations of corals. A natural question is -Why these shapes? One suggestion is that cells in living organisms proliferate to "maximize" their number (area) subject to any applicable constraints [Wer16] and this naturally results in hyperbolic geometries. This is a "static" argument the relates the mechanisms of growth to the resulting (quasi-2D) intrinsic geometry of living organisms. In this paper, we attempt to go beyond this "static" argument and develop models, based on thermodynamic considerations, to gain a quantitative understanding of the interplay between growth, mechanics and dynamics in soft objects. These models have the potential to describe the dynamical processes that result in the observed intricate three-dimensional (i.e. extrinsic) morphologies in nature. Figure 2: A free swimming sea slug Hexabranchus Sanguineus. The frames are 2s apart. Images used with permission from the copyright holders of the original video [Jon10].
Particularly striking examples dynamical behaviors in organisms with differential growth (hyperbolic geometries) occur in sea slugs (Nudibranchia) and marine flatworms (Polycladida). These marine invertebrates are found in many environments, particularly in coral reefs. While most of them crawl on the sea floor, a few are capable of free swimming [New03]. They move/swim by sending waves of undulations from the front to the back along their skirts (for sea slugs) or across their entire body (for flatworms). Fig. 2 shows 4 frames from a video of a free-swimming sea slug [Jon10]. The geometry of the slug is clearly hyperbolic. It has multiple undulations and undergoes significant bending deformations in the course of one swim cycle. While it is hard to quantify the strains within the organism it is not unreasonable to consider them small in comparison to the obvious large rotations/twist of the body.
In a different context, the interplay between growth and dynamics is also relevant to the development of leaves, flowers and other plant tissues that can be modeled as thin laminae [LM09,Bou10,LM11,Gor17,SS18]. Laboratory experiments using hydrogels [KES07, KHB + 12] have led to a semi-quantitative understanding of time-dependent, dissipative deformations of thin soft materials with a prescribed prestrain. In living organisms, however, the prestrain is not prescribed a priori, and how the prestrain development may be related to mechanics is not clear. Complicated physico-chemical processes are involved that need to be incorporated into mathematical models. It therefore seems reasonable to derive systematic constraints on the mathematical description based on a careful consideration of the non-standard kinematics involved and the general principles of continuum thermomechanics.
Earlier work, reviewed briefly in Sec. 2, implicates higher-order defects, in contrast to disclinations and dislocations, as playing a key role in the mechanics of intrinsically hyperbolic elastic sheets [GV11,GV13,GSSV16,SV]. This points to the need for tools to describe the evolution of (terminating) discontinuities of the second-gradient of the displacement field -when viewed at the macroscopic scale -for a proper description of the soft material deformations involved. It turns out that, within a continuum mechanical perspective, this fits in nicely within the question of describing the coupled mechanics of discontinuities and singularities of the elastic displacement field and its higher derivatives up to order three. This is the question that is addressed in this paper.
While, as evident from Fig. 2, it is natural, and necessary, to consider unrestricted finite deformations when dealing with soft materials, we restrict attention to 'small deformation' kinematics in this first effort due to the extra subtleties involved with higher order defect kinematics. Thus, we consider deformations of a fixed reference configuration that may or may not be stress-free. When the configurations attained by the system remain in close proximity to this fixed configuration, this is an adequate assumption. We consistently invoke Occam's razor as a guiding principle in our development -for instance, we restrict to the use of only ordinary stresses and couple stresses since forces and moments are the only agents of mechanical stimuli that we have some intuition for. Similarly, if branch point and surface defect velocities are to be the only dissipative mechanisms requiring constitutive specification without involving their spatial derivatives, then it turns out that the appropriate variables for the analysis of thermodynamics is in terms of the 'singular parts' of the first and higher order displacement gradients, instead of the more natural singular parts of the corresponding elastic distortion gradients that naturally arise in the analysis of defect kinematics. This is in sharp contrast to dislocation and g.disclination mechanics [AF15] where this distinction does not arise because of the relatively lower order kinematics involved. We develop the relationship between the two types of entities in this paper.

Statics and equilibria of non-Euclidean elastic sheets
One approach to modeling the mechanics of a growing hyperelastic body, borrowed from the literature of finite elastoplasticity, is to assume a reference configuration S and a deformation y : S → R 3 along with a multiplicative decomposition of the deformation gradient F = ∇y as F = EG (or F = F e F p in the plasticity literature) where the two-point tensor G models the effect of the growth processes in the material and E is the "residual" elastic deformation [Gor17]. The energy of the configuration defined by y is then given by W (E) = W (F G −1 ), where W denotes a hyperelastic energy density, vanishing on SO(3) [Gor17,LMP14]. In particular, the material is "stress-free" if F G −1 ∈ SO(3), although, for a general G, there might not be any deformation of the body y : S → R 3 whose gradient is a rotation times G. Such objects, with no stress-free configurations in R 3 , lead to incompatible elasticity.
The non-Euclidean formalism of thin sheet elasticity [ESK09] is a reduced dimensional description of thin elastically incompatible objects. The reference manifold S = Ω × [− t 2 , t 2 ], where t, the thickness, is "small" compared to the "in-plane" dimensions of the center surface Ω ⊂ R 2 . In this setting, the effect of the growth has a reduced dimensional description as a 2-manifold (Ω, g, b) where g, b are symmetric (0, 2) tensors. These tensors denote, respectively, the 'target' 1st and 2nd fundamental forms of the stress-free state of the sheet, pulled back to the reference manifold [ESK09]. This framework also applies to incompatible elasticity [BAG05,LMP14,BLS16] where, in general, there exists no deformation f : Ω → R 3 realizing a surface in ambient three-dimensional space whose first and second fundamental forms match (the push-forward of) the targets g, b (by f ), i.e., incompatible sheets have no stress-free configurations in our three dimensional space.
Assuming the Kirchhoff-Love hypothesis [Fun65], so that the (3D) deformation of a thin sheet is determined by the (2D) mapping y : Ω → R 3 on the center surface. This allows for and asymptotic expansion of the elastic energy as a sum of stretching and bending contributions [ESK09,ESK13] : where the oriented normal field N : Ω → S 2 , also called the Gauss Normal map [Sto89], is obtained from ∇y T · N = 0. Q 3 is a non-degenerate quadratic form, on symmetric 2 × 2 matrices, that depends on the Poisson's ratio ν of the material [ESK09], and dA is the area element on (Ω, g).
For various choices of g and b and boundary conditions, the energy functional (1) describes a variety of phenomena in thin sheets, including multiple-scale buckling in free sheets with 'excess length' near an edge, e.g. torn plastic or flat leaves treated with an Auxin near the edge [SRM + 02, SRS07,SMS04]. The excess length near the edge is modeled by a metric g with negative intrinsic curvature [ESK13].
Starting with a fully 3D elastic energy, Lewicka and Pakzad [LRP11] have obtained a reduced dimensional model for the limit t → 0 using Γ -convergence. In particular, they showed that for an appropriate quadratic form Q 2 . This energy has clear similarities with the energy in (1), although the details are somewhat different. Nonetheless, in either framework, the elastic energy scales like t 3 in the thin limit t → 0 if and only if there exist finite bending energy (mathematically y ∈ W 2,2 ) isometric immersions y : (Ω, g) → R 3 . To illustrate the physical import of this theorem, we remark that W 2,2 surfaces necessarily have a continuous tangent plane and normal [Eva98]  This theorem is central to the work in this paper. We are mainly interested in the mechanics of thin objects whose energy scale is O(t 3 ) as t → 0 (which includes smaller energies scaling with powers of t greater than 3). Such objects are naturally "floppy" since they are governed by "bending" and weaker forces, and the relevant defects are higher order in contrast to dislocations or (g.)disclinations.

Branch points and lines of inflection
The preceding remark highlights the role of the regularity of isometries. Beyond the existence/nonexistence of isometries, it is crucial whether a candidate isometry is in W 2,2 . This motivates the problem: Find y : Ω → R 3 such that ∇y T · ∇y = g and We have rigorous results showing that the problem (2) is flexible and solutions are plentiful [GSSV16,SV] (with prescribed zero-traction and moment boundary conditions, i.e. for free sheets). The proof is constructive, and uses ideas from Discrete Differential Geometry DDG [BSSZ08,GSSV16]. This lack of uniqueness in admissible static configurations with prescribed boundary conditions underscores the necessity of a dynamical model to 'choose' between acceptable configurations and/or describe the transitions between multiple admissible states [GV13].
If y : Ω → R 3 is C 1 , the Gauss Normal map is given by N = ∂ 1 y × ∂ 2 y ∂ 1 y × ∂ 2 y , where ∂ i = ∂ ∂x i for (arbitrary) coordinates (x 1 , x 2 ) on Ω. Further, if y and g are C 2 , Gauss' Theorema Egregium implies that (2) is equivalent to the Monge-Ampere Exterior differential system (EDS) [BCG + 13, IL03]: where dΩ is the area form on the sphere S 2 and κ is the Gauss curvature. Classical results in differential geometry imply that smooth solutions of (3) with κ < 0 are hyperbolic surfaces and locally saddle shaped. In contrast, the curly mustard leaf in Fig. 3b is Figure 3: Hyperbolic surfaces in R 3 . The inscribed (geodesic) triangle in the smooth saddle has angles that sum up to less than π, illustrating the connection between the extrinsic and intrinsic geometries -Gauss' Theorema Egregium.
"frilly", i.e buckled on multiple scales with a wavelength that refines ("sub-wrinkles") near the edge [SMS04]. This "looks" very unlike the smooth saddle in Fig. 3a. Why do we see frilly shapes in natural surfaces, as in Fig. 3b, rather than the smooth saddles of Fig. 3a? Indeed, any finite piece of a smooth hyperbolic surface can always be smoothly and isometrically embedded in R 3 [HH06] as "non-frilly" surfaces.
We have addressed this puzzle in recent work [GV11, GV12, GV13, GSSV16,SV] and the short answer is that, for a given metric g, the frilly surfaces, somewhat counterintuitively, can have smaller bending energy than the smooth saddle. It is true that C 2 (twice continuously differentiable) hyperbolic surfaces are saddle-like near every point. We find a topological invariant [SV], the index of a branch point -intimately related to the quantity Σ α (3) n da that emerges in Sec. 4 and the quantity Γ of Sec. 3.1 -that distinguishes sub-wrinkled surfaces from saddles locally. With branch points, the surfaces are only C 1,1 , but gain the additional flexibility to refine their buckling pattern, while lowering their energy [GSSV16]. This flexibility is not available to smooth saddles, and constitutes a key property of branched (sub-wrinkled) surfaces [GSSV16,SV]. Figure 4 shows a non-C 2 monkey saddle -a piecewise quadratic surface made from 6 sectors congruent to the wedge w = x 2 −3y 2 , x ≥ √ 3|y|, patched together by odd reflections about the lines x = ± √ 3y, x = 0 [GV11]. Although this surface is not C 2 , it is indeed W 2,2 and has a continuous normal vector and bounded curvature everywhere. Its "defects" include the point in the middle -a branch point and the 6 rays through this point -lines of inflection, which together constitute the asymptotic skeleton of the surface [SV]. This construction can be extended to generate C 1,1 hyperbolic surfaces with multiple distinct branch points, and an interesting question is how these defects interact with and influence each other [GSSV16].
Defects are of course ubiquitous in condensed matter systems. A key feature of defects in systems driven by a free energy is that the energy density typically diverges in the vicinity of a "bare" defect (and in some cases even the total energy diverges), and as a consequence, defects are always regularized, i.e. "cored" in physical systems. This is true for dislocations and disclinations in elastic objects, for creases in crumpled sheets, for defects in liquid crystals and many other types of defects. Uniquely, branch points and lines of inflection do not carry a singular energy density [GV11,GSSV16], and thus do not "need" a core for energetic reasons. Nonetheless, force and moment balance implies that these defects are indeed regularized into boundary layers, of width t 1/3 , mediating jumps in the normal curvature [GV12].
Like other defects in condensed matter, branch points and lines of inflection are thus mesoscopic. They contain large numbers of atoms (microscopic units) and are amenable to a continuum description, but are yet much smaller than the typical size of the sheet. Arguments from energy minimization, while implying the existence of these higher order defects, do not address the question of their evolution. One has to necessarily go beyond the elastic energy (1) and incorporate dissipative effects that are crucial in determining a thermo-mechanically consistent description of the coupled evolution of the shape y : Ω → R 3 and the internal geometry, given by the tensors g and b.

Notation
We define the notation employed in the paper in one place for convenience.
When a function on a domain is discontinuous across a (non-planar) surface S, we assume that its values along any sequence of points from either side of the surface approaching any fixed point on the surface take on a unique pair of limiting values, each element in the pair corresponding to the limit from one side. The difference of these limiting values, one for each point on the surface, is defined as the jump (denoted by · ) of the function on the surface. If ν(x) is the unit normal to S at x ∈ S, we say that x ± is a point on the ± side of S at x depending on ( We think of an n th order tensor as a linear transformation between the space of vectors (in the translation space of three-dimensional Euclidean space, also 1 st -order tensors) to the space of (n−1) th -order tensors, with its transpose defined in the natural way as being a linear transformation from the space of (n − 1) th -tensors to the space of vectors. All tensors components will be written w.r.t. the basis, (e 1 , e 2 , e 3 ) of a fixed Rectangular Cartesian coordinate system and all partial derivatives, denoted often by a subscript comma, will be w.r.t coordinates of this system. The Einstein summation convention will be used unless otherwise stated. Superposed dots will represent partial derivatives w.r.t. time. If A is a p th -order tensor then the operators ∇, div, curl may be (with invariant meaning independent of the choice of coordinate system and its basis, of course). The range of all indices above is 1 to 3 and e ijk represents a component of the third-order alternating tensor.
The symbol · i represents a contraction on i indices between two tensors. For any tensor A, we define the tensor obtained by symmetrizing in the first two indices as A (s) and the one obtained by antisymmetrizing in the first two indices from the left as A (a) . We denote the deviatoric part of a second-order tensor by the superscript dev.

Motivation for kinematics of the theory
In this section we provide some intuition on the defect kinematics we adopt for our theory of branch point singularities. This is first done by explicitly constructing a continuously differentiable deformation of a non-simply connected domain whose second derivative has a prescribed, constant jump across a planar surface in the body.
With reference to Fig. 5, we think of Ω occupying a simply connected d = 2 or 3-dimensional domain of ambient Euclidean space. Here, it may be viewed either as a right-cylinder (d = 3) or a cross-section perpendicular to its axis (d = 2). We choose a rectangular Cartesian coordinate system with the z-axis as the axis of the cylinder; e i , i = 1, 2, 3 are the unit vectors along the x, y, z directions, respectively. Ω c is a cylindrical subset of Ω with rectangular cross-section centered on the z-axis. The region Ω h := Ω\Ω c is not simply-connected. S is a surface in Ω h such that D := Ω h \S is simply connected. The layer L is defined as L = (x, y, z) ∈ Ω h | x < 0, − l 2 ≤ y ≤ l 2 and the surface S = {(x, y, z) ∈ Ω h | x < 0, y = 0}. We will refer to Ω c as a core.
Our goal in this section is to construct a vector fieldỹ (l) : Ω h → R n , n ∈ N, n ≥ d, 0 < l ∈ R, withỹ (l) ∈ C 1 (Ω h ) and the jump in ∇ 2ỹ(l) across S a specified constant, with the jump blowing up as l → 0 maintaining lim l→0 ∇ỹ (l) ∈ C 0 (Ω h ).
A necessary condition forỹ (l) ∈ C 1 (Ω h ) is that the jump in its second derivative across S be of the form ∇ 2ỹ(l) = A⊗ν, where ν is the unit normal field on S (with arbitrarily chosen orientation) and A is a R n×d valued matrix field on S. Noting that l −1 ∇ 2ỹ(l) may be formally considered an approximate discrete directional derivative of ∇ 2ỹ(l) in the direction ν (if the discontinuity were ignored), we define the field with A and ν = e 2 constants, and seek to construct solutions to the equations The restriction of Z to D is used since, even though Z is (distributionally) curl-free in Ω h (we interpret the curl of a matrix field as row-wise curls), Ω h is not simply connected but D is and hence we are guaranteed a solution Y in D, unique up to a constant. For any such Y field, we note that curl Y = 0 in D by the symmetry in the last two entries of Z, i.e. (∇Y e l )e k − (∇Y e k )e l = (Ze l )e k − (Ze k )e l = 0. Thus W satisfying (5) can be constructed, and W is also unique in D up to a constant for a given Y .
Arbitrarily fix one of the available Y fields. Such a Y has the explicit representation fopr x 0 being any point on the surface S, and the line integral is along any path from x − 0 to x contained in D. Now choose any path going from x − 0 to x + 0 (see Fig. 5) contained in D with the stipulation that it go through the points x 0 ± l 2 e 2 and the segments between x 0 and x 0 ± l 2 e 2 , respectively, are parallel to e 2 . Then, along the segment Y (x(s)) remains constant along the path between x 0 − l 2 e 2 and x 0 + l 2 e 2 . Therefore, using the segment Hence, the jump in Y at x 0 is given by Since x 0 ∈ S and Y such that ∇Y = Z in D were chosen arbitrarily, (7) holds for all x 0 ∈ S and admissible Y in the specified class. Thus Y is unique in that class, independent of position on S, and given by the constant −A ⊗ ν. We note that Y * : Ω h → R n×d×d may be viewed as a discontinuous function with the specification where the points x − ∈ D belong to the − side of S at x. We now evaluate the jump in the field W on S.
As already observed, for any Y satisfying (5) a W field in D can also be constructed and this will have the representation for any path linking y to x in D. We now arbitrarily fix an admissible field Y and choose the same path from Since Y remains constant at the value given by (6) along the chosen path from x 0 − l 2 e 2 to x 0 + l 2 e 2 , using (6) and (8). Now and Y (s) along the segment x(s) = x 0 + l 2 e 2 − se 2 , 0 ≤ s ≤ l 2 − s + is given by so that and therefore (9), (10), and (6), Thus, We now define the function W * : Ω h → R n×d as where the points x ± belong to the ± sides of S at x, respectively. W * is a continuous function on Ω h . We now assume that the constant A is of the form A = a ⊗ ν for a ∈ R n . Then, in D, This constant is free to choose, without loss of generality (related to the choice of Y (x − 0 ), for instance), and we assume that it satisfies (Ce i )e j = (Ce j )e i for i, j = 1, . . . , d. Then curl W = curl W * = 0 in D. This further implies that the line integral W * dx = b ∈ R n is a constant for any closed contour encircling Ω c .
If b = 0, then we define If not, we explicitly solve the system Solutions exist to this system (e.g. an explicit solution on star-shaped domains can be written down by using the Riemann-Graves integral operator [Ede85]) that belong to C 1 (Ω h ). Forcing by the Dirac distribution is not necessary; functions of (x, y) with support in a cylinder contained in Ω c satisfying Aα e 3 da = −b for any area patch A threaded by the cylinder also suffice for generating such solutions [Ach01]). Then defining we note that W dx = 0 for any closed contour encircling Ω c and that W ∈ C 0 (Ω h ). Then we defineỹ : for arbitrarily fixed z ∈ Ω h and a constant p ∈ R n . Clearly,ỹ satisfies ∇ỹ = W on Ω h andỹ ∈ C 1 (Ω h ).
Consider the constant vector a ∈ R n to be parametrized by the layer width l as All fields constructed with the use of A = a (l) ⊗ ν are denoted by a superscript (l). We have the following properties: This conclusion also holds for any value of β ≥ 0 when l > 0 is held fixed.
Remark 3.1 While the above considerations have dealt with one singular surface, the linearity of the construct on the prescribed field Z makes it clear that exactly similar arguments hold for the superposition of a set of deformations, each element of which contains a single planar surface of discontinuity of arbitrary orientation in Ω h terminating on Ω c . Considering y i , i = 1 to n ∈ Z + , each corresponding to a specified Z i field, the composite, superposed deformation n i=1 y i is C 1 (Ω h ), with generally discontinuous second derivatives on each of the S i corresponding to the specified Z i field. This corresponds to situations with a single branch point [GV11,GV12,GV13] as exemplified by the piecewise quadratic monkey-saddle that we discussed in Sec. 2.1.
Furthermore, given a fixed, simply connected domain Ω, let Ω i c ⊂ Ω, i = 1 to n, be a set of non-intersecting cores with Ω i h := Ω\Ω i c . Let each Z i now be specified on the domain Ω i h . Then each This corresponds to configurations with multiple branch-points [GSSV16,SV].
Remark 3.2 For thin objects modeled by d = 2, the construction above is a representation of folds without ridges. In Sections 4 and 5 we develop a continuum mechanical theory that encompasses the mechanics of such folds in simply connected domains within a setting that allows for deformations with less smoothness.
Remark 3.3 Consider d = 2, n = 2 and b = 0, and assume that W * (x), x ∈ Ω h , is invertible. A field y * : D → R 2 satisfying ∇y * = W * in D can be constructed that may be interpreted as a discontinuous deformation of Ω h . Now consider the metric g := W * T W * on Ω h . By the Nash C 1 embedding theorem, there exists a C 1 deformation z : Ω h → R 3 with (∇z) T ∇z = g = (∇y * ) T ∇y * .
For a mechanistic interpretation, consider the configuration in R 3 defined by z(Ω h ) as the stressfree, global reference configuration in a higher dimensional space (R 3 ) corresponding to a stressed body with a dislocation (with excluded core) in R 2 represented by Ω h . The stress-free reference cannot be represented by a compatible mapping of Ω h in the lower-dimensional space R 2 ; instead, one of its stress-free representations in R 2 is defined by the configuration y * (Ω h ). The stress-producing elastic Right-Cauchy Green tensor field is given by (W * −1 ) T W * −1 on Ω h .

3.1
The discontinuity of the deformation of a non simply connected domain with prescribed third 'deformation gradient' Consider the domain Ω h of Fig. 5 which is rendered simply connected by a single cut-surface S which is not necessarily planar. As before, we refer to Ω h \S =: D. We consider Z : Ω h → R n×d×d×d as a given field for which ((Z(x)e l )e k )e j is invariant w.r.t interchanges of e j , e k , e l for any values of j, k, l ∈ {1, . . . , d}. Furthermore, we assume that Z ∈ C 0 (Ω h ), and curl Z = 0 in Ω h . We are now interested in the construction of a field y : D → R n that satisfies which is possible since curlZ = 0 and D being simply connected. We note that Y Ijk (x) − Y Ikj (x) = Y Ijk (y) − Y Ikj (y) for x, y ∈ D, due to the symmetry of Z Ijkl in j, k and the connectedness of D.
Since the construction of Y allows the free specification of its value at one point of D, it can be assumed without loss of generality that Y Ijk = Y Ikj in D. Equation (15) and the symmetry of Z in the last two indices imply curl Y = 0 in D. Thus it is also possible to construct W : D → R n×d satisfying Furthermore, (16) and the symmetry of Y in its last two indices imply that a function y : D → R n can be constructed satisfying Now, because Z is curl-free in Ω h , we have by Stokes' theorem that Z dx =: Γ ∈ R n×d×d a constant, for the line integral over any closed loop encircling Ω c . (18) By (15) and (16), this further implies that Let x 0 , x ∈ S be connected by a curve c contained in S. Consider curves c + and c − on the ± sides of S connecting x ± 0 to x ± . Then as c ± → c. Similarly, (17) implies Now, due to the symmetry of Γ in its last two indices, Γ Ijk c k dc j ds = 1 2 d ds (Γ Ijk c k c j ) and the last line integral in (21) evaluates to 1 Remark 3.4 The jump in the deformation y across the cut-surface S is not arbitrary, being characterized by a finite set of parameters. One choice for this parameter set is the jump of the deformation at an arbitrarily fixed point on S, the jump of W at the same point, and Γ , the latter being a constant decided by the given field Z.
Remark 3.5 W is not constant on S even though curl Y = 0 in D unless the vector joining any two points on S lies in the null-space of Γ by (20). For S a planar surface with unit normal ν and Γ of the form a ⊗ ν ⊗ ν, a ∈ R n , ν ∈ R d constants, (20) implies that W is constant on S. If, moreover W − ( W ν) ⊗ ν = 0, then y is also a constant on S. These are all conditions satisfied by the example worked out in the preamble of this Section.
Remark 3.6 The argument remains unchanged for the case Ω h is just a punctured domain, i.e. Ω c shrinks to a point (a curve).

Kinematics
In this section we propose the kinematics for a model of the type of discontinuities treated in Sec. 3, to be broadly applied to the mechanics of materials. For that purpose, it is essential to deal with simply connected, compact domains containing the said discontinuities. The excluded core regions are now included in the domain as are the excluded surfaces of discontinuity. Roughly speaking, we consider an additive split of fields into 'regular' and 'singular' parts whenever the field in question contains high magnitudes concentrated in 'thin' regions approximating smooth lowerdimensional (< d) sets; the support of the singular part of the field contains these regions of high concentration and that of the regular part contains the support of the rest of the field, including regions supporting approximate discontinuities. Importantly, both the singular and regular parts are assumed to be at least integrable functions as we want to write governing equations for these fields in the form of pde that can at least be made sense of in some weak manner. Thus, we take a somewhat microscopic point of view, assuming that discontinuities and singularities of certain fields when viewed from a macroscopic scale have a smoother definition at a microscopic scale that we describe by additional 'eigenwall' fields. We also adopt the point of view that once macroscopic theories generate discontinuities and singularities, in most circumstances additional physical insight beyond the constraints placed by the governing equations of the macroscopic theory are required to define evolution with a modicum of uniqueness. We develop such a model in the rest of the paper. We refer to a fixed reference configuration, a simply connected compact region as B. In terms of the displacement field u and the i-eigenwall fields S (i) , i ∈ {1, 2, 3}, we define the i-elastic distortions Y (i) , i ∈ {0, . . . , 4}, as (Y (0) is analogous to the field y of Sec. 3, Y (1) to W , Y (2) to Y , and Y (3) to Z). Thus Y (0) = u and Y (4) , the gradient of the regular part of the gradient of the 3-elastic distortion, are assumed to have no 'singular' parts. We now define the 'composite' eigenwall fields S (i) , i = 1, 2, 3, as and we note that Physical considerations related to predicting stress fields of terminating twin boundaries and the stress-free, compatible, elastic, twinning shear distortions of through-twin boundaries [ZAP18] motivate the introduction of the following Stokes-Helmholtz (SH) decompositions: i ∈ {1, 2, 3}, We will also consider exactly analogous SH decompositions for the fields Combining (25) and (27) and noting the uniqueness of the SH decomposition we have up to at most a spatially constant function of time which we will assume to be a time-independent constant. Defining (noting that H (4) = 0), we define the i-defect density tensors for i ∈ {1, 2, 3} from (23) and (29) as using (25) and S (4) = χ (4) = 0. Since α (i) are defined locally as a curl, the local forms of the conservation laws for topological charge content, Σ α (i) n da, of an arbitrary area patch Σ is given bẏ where V (i) , for each i, is a vector field. V (i) is the velocity field of the i-defect density field. Combining (30) and (31), we have that for some F (i) that can be prescribed. Equations (24) and (32) implẏ with the last sum vanishing for i = 1. By kinematical arguments related to allowing for transverse motion of walls characterized by localized S (i) fields on surfaces, a part of F (i) is of the form F (i) = S (i) V ⊥(i) , where V ⊥(i) is the velocity of the i-eigenwall field. Guided by simplicity in thermodynamic arguments that precludes the appearance of (unremovable) gradients of dislocation and eigenwall velocity fields in the expression for dissipation of the body (see Sec. 5), we make the following choice In (32) and (33), incorporating (34), V (i) and V ⊥(i) are to be constitutively specified, minimally consistent with the second law of thermodynamics to be globally satisfied for all processes of any body modeled by this theory. Surfaces of displacement discontinuity (e.g. stacking faults) are not known to move transverse to themselves; moreover, such discontinuitites are often not identifiable based on knowledge of only the current state (and not of the distinguished coherent reference from which displacements are measured). Hence, we will assume V ⊥(1) ≡ 0. Elastic phase boundaries, i.e. localizations of the S (1) field along surfaces are known to move transverse to themselves, and not much is known about transverse motions of surfaces of discontinuity of the second gradient of elastic distortion, i.e. surfaces of inflection. Thus, we allow V ⊥(i) , i = 2, 3 to be nonvanishing fields in general. Hence, we have the following evolution equations for the eigenwall fields:

Thermodynamics
We assume a free-energy density function of the body with the following dependencies: using (29), (24), (28), and noting that H (4) = 0 (where the argument fields of each of the functions are evaluated at (x, t) to give the value of ψ(x, t)). Roughly speaking, the dependencies of ψ * on Y (i) , α (i) , i = 1, 2, 3 are expected to be convex and those on S (i) , i = 1, 2, 3 to be multi-well, nonconvex.
The balances of linear and angular momentum are given by where ρ is the mass density, v is the material velocity vector, T is the stress, Λ is the couple stress, and b, K are the body force and body-couple densities per unit volume, respectively. As usual in solid mechanics, we assume balance of mass is satisfied once the deformation map at any instant is determined by evaluating the density field on the deforming body from the formula ρ = ρ 0 det(I+∇u) , where ρ 0 is the density field on the reference configuration.
The mechanical power supplied to the body is defined as [MT62] P : using the balances of linear and angular momentum, where n is the outward unit normal to the boundary of the body, ω := 1 2 curlv = − 1 2 X · 2 Ω is the rotation vector where Ω := 1 2 ∇v − (∇v) T is the rotation-rate tensor, D := 1 2 ∇v + (∇v) T is the strain-rate tensor, and M := ∇ω. Denoting the mechanical dissipation, D, or the difference between the power supplied to the body and that stored in it, is given by In the following, we deduce guidelines for constitutive specification in our model that ensure that the mechanical dissipation vanishes in the absence of eigenwall and defect field evolution in any process and is positive otherwise, a minimal necessary condition for the mathematical model to be well-posed.
To facilitate the derivation of the thermodynamic driving forces for the various defect density and eigenwall fields, we will need the following auxiliary fields P (i) , i ∈ {2, 3} defined by the solutions of the following Poisson equations: which requires that the free-energy density function should satisfy the constraint (This is formally easily arranged by taking any arbitraryψ with the dependencies of (36) 3 , and defining ψ =ψ − 3 i=2 |Ω| −1 Ω ∂ H (i)ψ dv · i H (i) , but its physical and rigorous mathematical implications need to be understood). Defining (that exist by a unique Stokes-Helmholtz resolution of R (i) ), will aso be required in the sequel for deriving the thermodynamic driving forces. A long computation involving (36) 3 and the kinematics of the model defined in Sec. 4 reveals that the mechanical dissipation may be expressed in the suggestive form Thus, a set of constitutive equations, driving forces for dissipative mechanisms (denoted below by the symbol ), and some boundary conditions for the model are [−∂ ∇ 2 u ψ n + (div ∂ ∇ 3 u ψ) n] (s) (it can be checked that the rhs of (52) is deviatoric). Equations (51)-(54) along with the constitutive choices for the defect and eigenwall velocities to be in the direction of their respective driving forces, mediated by a positive, mobility/drag scalar required on dimensional grounds, ensures non-negative dissipation. Of course, other choices consistent with positive dissipation are possible as well. The boundary conditions (53)-(54) are not the most general, but a compromise between including higher order stress tensors with dubious physical meaning beyond couple stresses and simplicity in an already involved higher order theory of defects. It is clear from (51)-(52) and (37) that the governing equations lead to sixth-order pde in the displacement field u (see Sec. 6 below).
Remark 5.2 The composite eigenwall fields are coupled to each other through (35) and through the displacement field, appearing in the driving forces for the defect and eigenwall velocity fields, governed by (37). The results of Sec. 3 shows how the presence of a higher order defect (characterized by Γ = 0 in a non-simply connected domain) induces a lower order defect ( y = 0) that, in general, induces stress in the body (Remark 3.3).

Example: a model of branch-point defects in an elastic body
We specialize the general formalism to a specific case by making the simplest possible choice for the free energy density (36): with the ansatz that S (1) = 0, S (2) = 0, so that S (3) = S (3) =: S. Here, C is the 4 th -order tensor of elastic moduli with major and minor symmetries, c 2 , c 3 are non-negative scalars (in place of sixth and eighth order tensors!), d 3 is a positive scalar (that could also be a positive scalar-valued function of |curlS|), and l, 3 are positive scalars. The physical dimensions of c 2 , c 3 , d 3 , l, 3 are stress.(length) 2 , stress.(length) 4 , stress, length, and stress.(length) 6 , respectively. Since the equilibria we envisage are of nominally elastic bodies that show non-trivial shapes under no applied loads, f is generally expected to be a multi-well nonconvex function with the bottom of one well at the argument 0. Thus we are looking for the mechanics of surfaces of inflection and branch line defects in bodies with an evolving stress-free reference characterized by the choices S (1) = 0, S (2) = 0, S (3) = S, and since it is impossible to construct a displacement field of a 3-d body with vanishing strain, i.e., ∇u) (s) = 0 , whose third gradient is non-vanishing, the energy/stress-free reference for our body is never immersible in three-dimensional Euclidean space whenever S = 0, i.e. the stress-free state is necessarily incompatible or non-realizable.
The balances of linear and angular momentum (37) are solved by taking a curl of (37) 2 to obtain div T (a) = 1 2 curl div Λ dev + 1 2 curl K, that on substitution in (37) 1 leads to Constitutive equations (51)-(52) are used to solve for a displacement field from (58) (when the defect fields are assumed given), thus satisfying (37) 1 , and (37) 2 is then satisfied, in terms of this displacement field, by simply evaluating T a from the equation making the assumption that the constitutively undetermined trΛ = 0, without loss of generality. For the constitutive choice (57) and T (s) im = C imkl u k,l − c 2 u (i,m)ll + c 3 u (i,m)lppl − c 3 S (im)lp,pl ; div T (s) i = C imkl u k,lm − c 2 u (i,m)llm + c 3 u (i,m)lpplm − c 3 S (im)lp,plm so that the governing equation for the displacement field ((58)) may be written as where ∆ 3 (∆ 3 (·) = (·) ,iijjkk ) and ∆ 2 (∆ 2 (·) = (·) ,iijj ) are the triharmonic and the biharmonic operators, respectively. To develop the evolution equation for the field S we assume V ⊥(3) = 0 for simplicity. Since ψ in (57) does not depend on H (3) , we have P (3) = 0 in (55) 3 . The governing equation for the evolution of S therefore is given bẏ where B is a drag coefficient with physical dimensions of stress.(length) −2 .time.
Remark 6.1 Spatial derivatives of the 3-eigenwall field serve as a source term in (62) where ν is the unit normal to a planar surface, g is a scalar-valued function of the spatial coordinate along ν given by ζ = ν · x (say a Gaussian centered at ζ = 0), and b is a constant vector, this forcing is of the form d 3 g dζ 3 b. Equation (63) implies that there is no evolution of the eigenwall field at locations where curl S = 0, regardless of the energetic driving force there. For example, the field S(x) = g(ν · x) b ⊗ ν ⊗ ν ⊗ ν has no 'longitudinal' variation and does not evolve according to (63). However, where t is orthogonal to ν does evolve. Physically, the eigenwall field is 'dragged' by the evolution of its core.
Continuous dependence w.r.t initial data of the Cauchy problem for the evolution of displacement requires c 3 ≥ 0. When c 3 = 0, c 2 must be non-negative, with the requirement that µ ≥ 0 and λ + 2µ ≥ 0 if c 2 = 0. Within these parameter regimes, linear instabilities can arise for wavenumber and parameter combinations resulting in c d or c s taking complex values.

Uniqueness of the displacement field and boundary conditions
Our model encompasses a model of third-order elasticity in the absence of dissipative defect evolution, and involves the thermodynamically motivated higher-order boundary conditions (53)-(54).
Here, we use a uniqueness argument (in a putative smooth class of solutions) to deduce a full set of boundary conditions for the problem (62) when the S field is assumed specified. We abstract the results of the exercise in this special case related to the 'quadratic' energy (57) to identify a likely set of sufficiently general boundary conditions for the determination of the displacement field for processes consistent with the general constitutive statement (36). Consider two solutions u (1) and u (2) of (62) corresponding to identical S, K, b fields. Denote the difference displacement as u := u (1) − u (2) and its velocity v =u. Then u satisfies ρü = c 3 ∆ 3 u − c 2 ∆ 2 u + div (C∇u), and taking the inner-product of the difference velocity with the equation and integrating in space, we have Let us now assume that both u (1) and u (2) satisfy (53)-(54) consistent with (57). Then the last line of (64) vanishes due to the boundary condition (54) and the line before that due to (53). Let the stress field arising from (u (i) , S), i = 1 2, be T (i) = T (s)(i) + T (a)(i) , in accord with (59), (60), and (61). Then the third line from the bottom of (64) may be interpreted as and if we now additionally require that solutions satisfy specified tractions and velocities (or displacements) on mutually complementary parts of the boundary of the body for all times, then this term vanishes.
Consequently, we are left with d dt and if u (1) and u (2) both satisfy specified initial conditions on the displacement and velocity fields, then the bracketed quantity, an integral of sums of squares (in fact, the potential and kinetic energies of the body subjected to the difference displacement) vanishes at all times. This proves that the difference velocity vanishes point-wise, and the initial condition on the difference displacement implies that u (1) = u (2) for all (x, t). Obviously, the dynamic problem allows the prediction of unique rigid motions. In statics, i.e., when the inertia term is absent, one takes the inner product of the governing equation for the difference displacement with the difference displacement, and obtains, for the same boundary conditions (except only the displacement can now be specified on the part of the boundary complementary to where tractions are specified), All integrands are non-negative implying that the strain, or the symmetrized displacement gradient, vanishes (recall the minor symmetries of C) which, by compatibility, further implies that the displacement field is unique if a displacement boundary condition is specified and otherwise it is unique up to an infinitesimally rigid deformation. Thus, the higher order boundary conditions (53)-(54), along with classical displacement and traction boundary conditions may be expected to define a well-set problem (for the displacement field) in the case of the general constitutive equation (36). Of course, the traction now involves a stress tensor that has an antisymmetric part, and is constitutively dependent on higher order displacement gradients.

A 'plate' idealization
Let the reference B be a plate of thickness 2t, i.e., where B 2 is a flat 2-dimensional simply connected domain. Defining the through-the-thickness average of a function as , we now seek the governing equations for u and S, under the ansatz that S = S and ρ = ρ, i.e., S and ρ do not vary through the thickness of the plate, and K = b = 0. It is also assumed that a component of S vanishes if any of its last three indices takes the value 3. We use the notation that all lowercase Greek indices vary from 1 to 2 while lowercase Latin indices span from 1 to 3.
While not essential, the assumptions l = 2t, c 2 = Et 2 and c 3 = Et 4 , where E is the Young's modulus of the material can be made to draw an analogy with classical plate theory (the curvaturerelated elastic energy term in the thickness-integrated expression of (57) would then be proportional to t 3 ). For 0 < t 1, whenever S = 0, there is energy and stress in the body, possibly small, with the corresponding thickness-integrated 'elastic' energy of the plate (arising from the first three terms in (57)), alternatively the 'plate elastic energy', scales as ∼ t 5 , assuming energy is minimized, there are no external forcing or constraints, and 3 > 0 to rule out any possibility of a singular energy. Our governing equations (62) or (65) do not require that energy be minimized, so that scaling of the thickness-integrated elastic energy w.r.t t as t → 0 in the model can well contain lower order bending O(t 3 ) , and even stretching (O(t)), contributions.
Applying the averaging operator to (62) and noting that C ijkl u k,lj = C iβkα u k,αβ + C iβk3 u k,β3 + C i3kα u k,α3 + C i3kα u k,33 S ijkl,jkl = S iαβγ,αβγ + (S iαβ3 + S i3αβ + S iα3β ) ,αβ3 + (S i333γ + S i3γ3 + S iγ33 ) ,33γ + S i333,333 , we obtain Similarly, In equation (65), the terms beyond the first line represent forcings in the transverse direction to the plate and need to be specified (it would be physically legitimate to assume many of these terms to vanish); the third line of (66) has similar meaning and needs specification. The functions u, S represent the fundamental fields of the plate theory, governed by (65)-(66). Evaluating T (a) from (59) in terms of (u, S) solving (65)-(66) and K would imply the satisfaction of balance of angular momentum (i.e., moment balance) in the through-the-thickness averaged sense.
Remark 6.3 We note that non-evolving and non-vanishing S (1) , S (2) 'target' composite eigenwall fields can be included in the considerations of this Section (Sec. 6), with only slight increase of tedium in bookkeeping.
Within the context of energy minimization and for t > 0, if curl curl S (1)(s) T = 0, i.e.
S (1)(s) satisfies the St.-Venant compatibility condition, then an infinitesimal isometry exists (the reference configuration is assumed to be simply-connected) and the plate elastic energy scales as ∼ t 3 or of smaller magnitude; if S (1)(s) is not compatible, then the energy has to scale as ∼ t. We note that when S (1)(s) is compatible, unless S (2) = ∇ 2 v, where v is s.t. (∇v) (s) = S (1)(s) so that ∇ 2 v is unique, the plate elastic energy is going to scale as ∼ t 3 . The requirement S (2) = ∇ 2 v is nongeneric for a freely-specifiable S (2) field that, however, is satisfied by the choice S (1) = 0, S (2) = 0. Thus, in most circumstances the plate energy is expected to scale as ∼ t 3 , if the plate energy is minimized.

Discussion
Starting from the work of the brothers Cosserat [CC09, as presented in [TT60]], through those of Toupin [Tou64], Green and Rivlin [GR64], Mindlin [MT62,Min64], on to that of Fleck and Hutchinson [FMAH94,FH01,Hut12] and of Gurtin [Gur02,GA09], higher order theories of continuum mechanics have made an appearance off and on and have been noted for their intricacy and elegance, but always, arguably, with the nagging question of the physical justification (in their details 1 ) in view of their added complexity. Our work aims to provides a concrete, tangible, and compelling justification -that the precise treatment of defects in the deformation and its higher order gradients is the raison d'être for higher order theory in continuum mechanics. Our work is in the context of non-Euclidean elastic sheets with negative in-plane Gauss curvature. These objects are ubiquitous in nature and they display varied and intricate multi-scale behaviors [SRM + 02, AB03, KES07, KHB + 12, GV13]. Their elastic behavior is significantly different from that of elastic plates or spherical shells [GSSV16,SV]. In particular, they have "large" continuous families of low-energy states obtained from piecewise isometries, with each piece possessing additional "bending" degrees of freedom. Thin hyperbolic free sheets are thus easily deformed by weak stresses and their morphology is strongly dependent on the dynamics of the growth/swelling processes, material imperfections, or other weak external forces. This naturally motivates the need for tools to describe singularities/defects in these sheets, their interactions and the resulting dynamics.
Mesoscopic defects in hyperbolic sheets, associated with their "soft" modes of deformation, include lines of inflection that terminate at branch points [GV12,SV]. These are higher-order defects (termination of jumps in curvature) unlike the more common types of defects, disclinations and dislocations. Irreversible effects in the dynamics of disclinations and dislocations are associated with (macroscopic) plastic behaviors -stress-free large deformations, internal stresses, and microstructure -in solids. A natural question therefore is -what are the macroscopic manifestations of moving lines/surfaces of inflection and branch points/lines?
In this work we have begun to address this question in the context of 'small deformations' from a (potentially stressed, when occupied) reference configuration. A detailed analysis and characterization of the kinematics of branch point defects and the discontinuities in the deformation that they induce is achieved. This analysis, in its essence, is a non-trivial adaptation and extension of the ideas of Weingarten [Wei01] and Volterra [Vol07], from the dawn of elastic defect theory, to a context not restricted within the kinematics of only strain (the symmetrized gradient of the displacement, as well as its nonlinear analog) and its incompatibilities, and shows the natural way forward for deducing the constraints on possible jumps in deformation, i.e. global constraints, for locally compatible higher order deformation gradients, albeit on domains with the simplest non-trivial topology 2 . We then develop a thermodynamically consistent theory for the dissipative dynamics of such defects in a nominally elastic solid, allowing for their interaction with dislocation, g.disclination, grain, and phase boundary defects. The constitutive guidance provided by this thermodynamic argument ensures that the model is equipped with an energy (in)equality, a crucial necessary condition for its physical and mathematical well-posedness. The analysis uncovers the non-Newtonian, energetic driving forces on these defects that couple their dynamics and mutual interactions to applied loads and the deformation of the body 3 . Evolution of the defect fields subject to such driving forces necessarily reduces the system free-energy by design, within an overall dynamics that accounts for material inertia and is not restricted to its free-energy decreasing with time (depending on the external driving). As an example, we explicitly demonstrate the full set of governing equations for the case of branch point defects in an elastic material and develop a 'plate' theory idealization for it. The development of the finite deformation version of the model poses no conceptual or technical barriers 4 based on our prior work in g.disclination mechanics [AF15], but this same work makes it clear that the bookkeeping tasks in pushing through the analysis are going to be formidable.
We observe in passing that while we have been interested in developing a theory for branch point/line defects and lines/surfaces of inflection, i.e. a theory for the discontinuities and singularities of the deformation and its gradients up to order three, the analysis makes it clear that the mathematical/continuum mechanical formalism extends to describing the discontinuities and singularities of any finite integer order gradient of the deformation, while including only stresses and couple stresses. As already observed in [AF15], using the Second Law in global form is crucial for this, albeit at the expense of the application of limited (but adequate, as we show in Sec. 6.1) higher-order boundary conditions about which not much is physically known anyway.
As a final comment, we note that a geometric model of growth mechanics, based on Riemannian geometry and including evolution, has been proposed in [Yav10]. The viewpoint is different from ours and, in particular, the mechanics of incompatibility based on a Riemannian metric cannot describe (without non-trivial extension) the 'softer' branch point defects we focus on. We expect that one can recast our continuum mechanical kinematic constructs within a differential geometric structure involving the specification of a moving frame, and higher-order constructs based on such a field, thereby making connections with the "geometric" viewpoint of growth mechanics.