Topological elasticity of flexible structures

Flexible mechanical metamaterials possess repeating structural motifs that imbue them with novel, exciting properties including programmability, anomalous elastic moduli and nonlinear and robust response. We address such structures via micromorphic continuum elasticity, which allows highly nonuniform deformations (missed in conventional elasticity) within unit cells that nevertheless vary smoothly between cells. We show that the bulk microstructure gives rise to boundary elastic terms. Discrete lattice theories have shown that critically coordinated structures possess a topological invariant which determines the placement of low-energy modes on edges of such a system. We show that in continuum systems a new topological invariant emerges which relates the difference in the number of such modes between two opposing edges. Guided by the continuum limit of the lattice structures, we identify macroscopic experimental observables for these topological properties that may be observed independently on a new length scale above that of the microstructure.


I. INTRODUCTION
Mechanical metamaterials, defined as structures such as origami sheets or spring networks with repeating patterns of elements, possess properties not found in uniform slabs of material, such as superior strength to mass ratios and vanishing (pentamode) [1][2][3] or even negative (auxetic) [4][5][6][7] elastic moduli or Poisson ratios. Of particular interest are flexible mechanical metamaterials, which possess low-energy deformation modes that can be used to achieve shape-changing, programmed response and strong nonlinearities [8]. Some of these properties are now being demonstrated at microscopic length scales, via kirigami (cut) graphene ribbons [9], self-assembled patchy colloids [10], nanolithography [11] and DNA "origami" [12,13]. These cutting-edge techniques raise the question of what happens when such systems are manipulated on scales much larger than the unit cell. This is precisely the limit of conventional solid mechanics, which occurs well above atomic scales. However, conventional elastic theory assumes a smoothly varying strain field, and so cannot capture the short-distance rearrangements of the flexible unit cell.
To address these issues we develop a micromorphic continuum elasticity applicable to flexible structures. While conventional Cauchy elasticity depends only on the strain of an infinitesimal region, micromorphic elasticity considers regions with additional relevant structure [14]. While this is unnecessary for conventional atomic solids, it can lead to much richer mechanical response, somewhat as liquid-crystalline order modifies conventional fluid dynamics. We consider two situations in which such a theory arises: as the long-wavelength limit of microscopic lattice theories, and as an intrinsically continuum theory based solely on macroscopic observables without recourse to a particular microstructure. In either case, the repeating spatial structure gives rise to an energetic term based on the gradient of the elastic strain, which can be integrated to yield a surface term in the elastic energy.
Our particular focus is on topologically protected boundary modes. In the lattice theory, it has been shown that a topological invariant derived from the bulk structure controls the number of zero-energy modes on a given surface of mechanically critical lattices (having equal numbers of degrees of freedom and constraints, also called "Maxwell" or "isostatic") is determined by an integer-valued topological invariant derived from the bulk structure [15,16] and leading to a wide variety of physical implications [17][18][19][20][21]. This is derived from an integral around the topologically nontrivial Brillouin Zone, which is not present in the continuum limit. Instead, we present analysis resulting in a new invariant, which instead measures the difference between the numbers of zero modes on two opposing edges, analogous to the system's polarization rather than its surface charges. Because these systems are only marginally mechanically stable, it becomes important to consider how mechanical constraints generate energy costs based not only on strain but on strain gradients, and these latter terms break spatial inversion symmetry, necessary for the type of polarization we consider.
The remainder of the paper is organized as follows. In Sec. II, we show how a model class of microstructures generates mechanical constraints that couple both to macroscopic strains and to microscopic degrees of freedom. In Sec. III we show how this generates surface surface terms complementing bulk ones. In Sec. IV we show how a simple linear-algebraic procedure based on the self stresses (stressed states of equilibrium) captures the equilibration of the system and generates an effective theory in terms of smooth strain fields only. In Sec. V we introduce the topological invariant and relate it to the surface modes, relying on the continuum theory rather than a particular microstructure. In Sec. VI we characterize the length and energy scales associated with the boundary modes and how they can be experimentally observed. We conclude in Sec. VII.

II. EFFECTS OF SMOOTH STRAIN FIELD ON MICROSTRUCTURE
Here, we introduce the class of mechanical systems whose energy is stored in discrete, spring-like bonds and consider how these couple to external strain fields. Let u(r) denote the displacements that sites undergo at position r in the undeformed reference space. Consider a particular spring in a periodic crystal cell structure, and let r be the position of its cell, p the position of the center of the bond relative to the cell and b the vector from one end of the bond to the other, as shown in Fig. 1. The extension of the spring is, for sufficiently smooth displacement fields (repeated lower indices implying summation): In the preceding lines, we conduct an expansion assuming that the relative displacements of the ends of the bond remain small, yet considering a higher gradient in the displacement field. Because of the rotational invariance of the problem (reflected in the symmetric b i b j prefactor) even non-uniform rotations do not extend the spring and instead we obtain a result purely in terms of the symmetrized strain ij ≡ (1/2)(∂ i u j + ∂ j u i ) and its gradients. The strain gradient terms are usually ignored in the elasticity of rigid bodies, for which they are small for smooth strain fields. However, it is important in systems near the isostatic point, for which the contributions to the elastic energy of the conventional terms can vanish. Note also that this expression is evaluated at the position of the cell, not the bond, such that our object of interest is the strain field within the cell and the characteristics of the individual bond are encoded in p, b.
Notably, given that bonds repeat periodically throughout the structure, there is an ambiguity in which cell to assign which bond, the so-called discrete gauge symmetry of the lattice theory [15]. Choosing a different such assignment would shift p by some combination of lattice primitive vectors but shift the point at which the strain is evaluated by an equal and opposite amount, such that the physical observable of the bond extension at a particular point in space (as opposed to a particular cell) is unchanged. Keeping this in mind, we may use the above relationship to construct a tensor that relates the extensions e m (r) indexed by m to the strain field: Here, we introduce the initial rigidity map, R 0 , which determines the spring extensions resulting from a given set of strains. Although it is the continuum analog of the rigidity (or compatibility) matrix used in lattice theories [15,22,23], we term it "initial" here because it generally results in unbalanced forces on the sites connected via the bonds, a limitation that we resolve below. Note that since this object maps from dimensionless strains to spring extensions, its elements have units of length. An alternate approach would be to calculate the fractional extensions of the springs, such that the elements of the map were rescaled by spring lengths and thus mapped from the strain of the embedding space to the strains of the springs. FIG. 1. A periodic spring network has a periodic microstructure consisting of bonds connecting sites, as in the generalized kagome lattice shown here. A solid bond is located at p relative to the center r of the cell in which it lies. The bond connects two sites with relative position b, undergoing displacements u1(r), u2(r) that cause extension/compression of the bond. For continuum fields, the displacements may differ greatly for each site within the repeating cell, but vary smoothly across the cells.
In the preceding Section, we described a purely structural, geometrical relationship between strain fields and mechanical constraints. Here, we relate the violations of these constraints (in particular, the stretching of springs) to elastic energy. For simplicity, we consider a system with bonds of a single spring constant, and choose units such that this constant (per unit volume) is unity. The resultant energy comes from the sum of the squared spring extensions, integrated across the bulk of the system: Given the dependence of the spring extensions on both strains and gradients [(Eq. (1a)] each spring results in four energetic terms. Following a number of manipulations (see Supplementary Material) the two terms that break spatial inversion symmetry may be expressed as a total divergence term and hence reduced purely to a surface contribution E s to the energy expressed in terms of the outward-facing surface normaln: These terms therefore do not affect the bulk physics, but do modify the response of boundaries and interfaces, analogous to total gradient terms which are not present in nematics but do occur in cholesteric liquid crystals [24]. In addition to this surface term, we have spatial inversion symmetric terms due to the conventional strains and their gradients: Hence, we recognize that the energy consists of three types of terms: conventional bulk strain energies leading to elasticities of the sort described in [22], surface strain terms, and higher-order bulk terms consisting of strain gradients. Expanding the initial rigidity map further would generate higher-order bulk and surface terms. The bulk terms are even under spatial inversion and the boundary terms odd. Note that we have not yet made any assumptions about either mechanical equilibrium or mechanical criticality. As such, even when the conventional elastic force balance is achieved (i.e., the divergence of the stress vanishes in the bulk) there are unbalanced forces on the individual degrees of freedom shown in Fig. 2(a). Before we rectify this, let us illustrate the bulk and boundary energies with a minimal example, occurring in one dimension and with a single strain component (x). Let the extension of a bond be given by e = (a − b∂) , resulting in an energy Requiring force balance in the bulk (namely, that the differential of the energy with respect to the strain function vanish) generates strain profiles exponentially localized to the two edges: which are independent of the surface term −ab 2 , as indeed force-balanced terms must be. However, the surface term does determine the energy of these modes, such that depending on the sign of a/b either the left-hand or the right-hand term costs zero energy. This is the linear limit of the nonlinear soliton chain considered in [25]. In this way, we see that the surface term determines the existence and position of zero-energy boundary modes. The above equations also reveal that, in general, this surface term and hence the zero-energy modes are independent of the bulk energetics. Therefore, in the following sections we return to a general form of the bulk constraints in order to consider the effects of mechanical equilibrium, mechanical criticality and topological protection. Note finally that the energy of Eq. (7) is characteristic of a mechanically critical lattice; more generally a system could have additional energetic terms such that neither surface mode cost zero energy.

IV. MECHANICAL RELAXATION AND EQUILIBRIUM IN THE MICROSTRUCTURE
From the energy functional of Eq. (4), we may obtain the stress tensor σ ij (r) = δE/δ ij (r) with the details of the functional differentiation described in the Supplementary Material. Just as each bond extension was obtained as a linear operator on the strain, now the stress is obtained as a linear operator on the bond extensions, which we denote as the initial equilibrium map Q 0 ij,m , again in analogy to the equilibrium matrix of lattice theories: The above equations describe the elastic relationships between stress, strain, bond extensions and energy. Note that these results mirror those of the lattice theories [22]. In particular, if we consider modes with spatial variation exp(iq · r), the equilibrium map is simply the transpose of the rigidity map at the opposing wavevector: However, the smooth strain fields assumed here are unrealistic for isostatic systems, which undergo shortwavelength non-affine relaxation events to dramatically lower their energy in response to unbalanced forces within the unit cell. Consider, as we have above, the extensions of bonds that occur when a strain field is imposed externally. These extensions generate tensions, which in turn generate unbalanced forces on the sites. These result in additional displacements of the sites, achieving force balance and energy minimization. These displacements may of course be obtained by solving whatever microscopic force-balance equations are appropriate to a given microstructure, but a more elegant and efficient method exists for systems close to the isostatic point. For such systems, there are only a few states of self stress [15,23], which are sets of tensions that generate no force on any site. As has been shown [16], the post-relaxation bond extensions e are precisely the projection of the prerelaxation bond extensions e 0 into the space spanned by an orthonormal set of states of self stress, {s i }.
For Maxwell lattices, which have equal numbers of site degrees of freedom and bond constraints, this technique has proven useful to obtain the uniform elastic response. Such a lattice has d modes of translation that result, via the Maxwell-Calladine index theorem, in d states of self stress. The d(d + 1)/2 components of the elasticity tensor are all obtained from these, implying that d(d − 1)/2 strains cost zero energy. These are known as Guest modes [26]. This method, which accounts naturally for the large intra-cell relaxation, yields the correct moduli, whereas assuming a uniform strain field overestimates them substantially.
In the present analysis, we wish to consider external fields that are applied over length scales large compared to the unit cell, such that displacement fields necessarily vary smoothly over these scales. At the same time, the fields may vary dramatically within the cell, as when, e.g., a rigid triangular unit undergoes a rotation. To that end, we wish to consider the behavior of structures such that the center of mass of each cell is coupled to external fields (such as induced by smooth, rigid mechanical barriers) while the other, internal degrees of freedom (rotations and shears) are allowed to relax and achieve mechanical equilibrium. In this way, we retain a theory expressing the energetics in terms of the slowly-varying strain field across the periodic structure, while still permitting deformations that relax energy within the crystal cell.
To that end, we identify the smooth displacement (alternately strain) fields u(r) = u exp (iq · r). These are the fields we wish to retain, while the short-range relaxations within the unit cell are allowed to occur to lower the system energy. Note that although u has only d components, it captures all of the allowed strains. For q = 0, there exist d(d − 1)/2 mechanical compatibility conditions [26] that are satisfied by the d(d + 1)/2 components of the strain, such that they do indeed derive from a valid displacement field defined by d components.
How, then, do we allow relaxation of the remaining degrees of freedom? Previously, we had required that all the components of the force f = Qe on the sites resulting from the extensions must vanish. Now, though, we wish to allow forces on the smooth strain fields to remain unbalanced. For example, if we consider a set of strain fields that generates a force on the center of mass of a cell and a torque on one of its rigid elements, we would allow the cell to relax to alleviate the torque while retaining the center of mass force. That is, if {|v i } are the set of externally imposed center-of-mass displacements on the cells we wish to allow no center-of-mass forces, such that we allow v i |f = 0 but require that the other components of the forces vanish. To that end, our self stresses are now defined as an orthonormal basis for the nullspace of the equilibrium matrix with these modes projected out: For the periodic systems currently under consideration, modes at different wavevectors are orthogonal. Consequently, we can consider rigidity and equilibrium maps at a particular wavevector, such that ∂ j → iq j and a strain at a given wavevector results in spring extensions only at the same wavevector. Consequently, there are only d modes at any given wavevector that need be projected out in Eq. (12). Each such mode consists of displacements of sites along one of the Cartesian directions, with spatial dependence exp(iq · r). The result is d states of self stress.
We thus obtain the primary object under consideration in the present work, the rigidity map R m,ij , which describes the extensions which result from a smooth strain field following the short-distance relaxations within the unit cell. From this, the equilibrium map may be derived in precisely the same way as the initial equilibrium map was derived from the initial rigidity map. In terms of the initial rigidity map of Eq. (3) and s m,n , the n th component of the m th state of self stress, The above analysis implies that, if we think of the rigidity map as a matrix acting on the d independent components of the strain, it similarly results in d independent self stresses, and is thus a square matrix. That is precisely the Maxwell condition for lattice theories, that the energetic constraints and the degrees of freedom are equal in number. Thus, we arrive at a natural condition to extend the Maxwell condition to theories of continuum fields: that the configuration spaces and constraints be equal in number.
A final property proves crucial to the low-energy response of the lattice. The aforementioned Guest modes ensure that the rigidity maps have zero modes at low wavevectors. That is, generically the eigenvalues of the rigidity map would have O(q 0 ) components, yet for the Guest mode the leading contribution would be O(q 1 ).

V. TOPOLOGICAL POLARIZATION IN THE CONTINUUM
We have now generated a continuum map that describes the relationship between local strain degrees of freedom and a like number of energetic constraints, retaining the Maxwell criterion of the lattice theories. However, it is not a priori clear that this formulation captures the phenomena of the lattice theory, including topological protection of modes at interfaces, edges and defects. Such modes are generated and protected by topological invariants defined by the homotopy class of loops across 2. (a) On a periodic system, we apply a particular strain ij (r) = ( xx, yy )e i(q·r) . This causes bonds to stretch (green) or compress (red). (b) The system then relaxes by projecting onto space of self-stresses that mostly capture particles displacing over short distances but reducing the energy cost overall.
the Brillouin Zone, which no longer exists in our continuum formulation. Do the zero modes have continuum descriptions, and do they retain topological protection? We resolve these questions in the affirmative.
Since we mean to include edge modes, we allow q to be complex, with the imaginary part representing the rate at which the mode decays away from the edge. In particular, suppose we have an edge in the second spatial direction, on the line (0, r y ), with the system extending to its right, with r x > 0. A mode extending along and exponentially localized to this edge would have Im(q y ) = 0 real and Im(q x ) > 0, while one on the opposite edge would have Im(q x ) < 0.
Because any extensions of bonds cost energy, zeroenergy modes are precisely those that lie in the nullspace of the rigidity map of Eq. (13). Note that the use of this map means that we are allowing the local structure to relax to minimize energy, in this case to zero. Thinking of this map as a square matrix, we may take its determinant, det(q), which vanishes if and only if a zero mode exists. Because of the two modes of translation, det(q) has a double zero at q = 0 and hence takes the following form: with A n1n2 real coefficients set by the microstructure. Because our rigidity map has units, A ij q n1 x q n2 y has units of the volume of d-dimensional space (two, here). Terminating the expansion to order n in q indicates the presence of n zero modes. Of these, two take the form while the others are short-wavelength modes of order q 0 y . α is determined by the coefficients of order O(q 2 ), while the signed inverse decay length κ ± ≡ β ± q 2 y depends on the second and third order terms of Eq. (15). While a lattice may have short-wavelength modes, our continuum formulation deliberately excludes them and hence the short-wavelength modes here are non-physical and dictated by the order at which we terminate our expansion. Our long-wavelength zero edge modes are then restricted to wavevectors obeying the above equation. The trivial state is the one in which κ ± take opposing signs, whereas in a polarized state both κ ± would take the same sign. Thus, to capture the edge mode polarization, a topological invariant must count the numbers of positive and/or negative κ ± .
Because the loss of periodicity destroys the Brillouin Zone, a natural approach would be to perform an integral instead over all real values of q y . This could lead, via the extended real number line, to a non-contractible loop such as one has in the periodic case, analogous to an approach that was applied in [27,28] to determine Chern numbers on the surfaces of spheres rather than tori. However, this procedure would not be appropriate here since it would include the fictitious short-wavelength modes.
Instead, to capture the long-wavelength behavior, we introduce |l 1 |, the length of the first lattice primitive vector, and consider edge modes for which the transverse component of the wavevector is small: q y |l 1 | 1. In selecting modes over which to integrate, we wish to ensure that we consider values of the remaining component of the wavevector that extend far beyond those while also ignoring fictitious short-wavelength modes. That is, we wish to consider values of q x which satisfy q y |l 1 | q x |l 1 | 1. We choose then for ≡ q y |l 1 | to integrate over values of q x |l 1 | extending to 1/2 . As described in detail in the Supplementary Material, there are alternate methods which may capture the transition more sharply. The resulting topological invariant describing the numbers of modes on the left and right edges is:

Im[q x ]
FIG. 3. The contour used to establish the relationship between the bulk structure and the locations of the zero modes. Via complex analysis, the full contour counts the difference in numbers of long-wavelength zero modes on the left and right edges while remaining insensitive to the short-wavelength modes. As described in the text, for this particular contour the dashed portions may be neglected, such that the locations of the zero modes are determined solely by the bulk physics captured in the solid lines.
This relationship follows from the contour in the complex plane shown in Fig. 3. In that figure, we search for a relationship between the behavior of the bulk modes, shown on the real axis, and the zero-energy surface modes, represented as black points. Applying the Argument Principle means that the long-wavelength modes will cause the phase of the determinant of the rigidity map to wind by ±2π for modes on the left (right) edges. The curved portions of the contour themselves have imaginary components of the wavevector and hence are boundary modes that cannot be considered in establishing a bulk-boundary correspondence, a key signature of topological protection. However, while the contributions of the curved portions are never small, their difference does vanish in the given limit. As such, we may determine the presence of zero modes solely through the portions of the contour corresponding to bulk modes (solid lines). In this way, the topological invariant becomes quantized (integer) precisely in the long-wavelength limit, as shown in detail in the Supplementary Material. Indeed, this notion of the topological invariant emerging in the longwavelength limit is implicit in other continuum treatments, insofar as they likewise consider length scales over which atomic details may be neglected. This continuum object in fact bears a closer relationship to conventional polarization than does the discrete analog, which counts the number of zero modes on a particular edge rather than the difference between the two edges.
To demonstrate that nonzero topological polarizations do occur, we consider a class of generalized kagome lattices, consisting of a periodic cell with three sites in two dimensions, joined by three intracellular and three intercellular bonds as shown in Fig. 4. In that figure, we consider a 1D family of such lattices that undergoes a topological transition at which deformation of the lattice shifts a boundary mode from one edge to the opposing edge, polarizing the system. As shown in Fig. 4(a), our numerical technique generates noticeable error as described in the Supplementary Material very close to the transition point, which can also be identified by direct geometrical means.

A. Topological polarization as a vector
In the continuum, it becomes natural to ask about the polarization not only at a certain interface -for now, we've focused on a vertical interface, with decay along the horizontal direction (x-direction) -but for every possible interface. Specifically, we look for ∆N (θ n ), the difference between the number of zero modes on the "left" and "right" edges when the vector pointing from the left to the right edge isq n (θ n ) ≡ cos θ nx + sin θ nŷ . In this way, we can recast the wavevector in terms of the components along this normal direction, q n ≡q n · q and the component q t along the tangent directionq t (θ n ) ≡ − sin θ nx + cos θ nŷ . This allows us to re-write the determinant in the new (q n , q t ) basis and find the longwavelength modes of present interest: det(q n , q t ) = det q x = cos θ n q n − sin θ n q t , q y = sin θ n q n + cos θ n q t = A 2,0 q 2 n + A 1,1 q n q t + A 0,2 q 2 t + iA 3,0 q 3 n + . . .
where once again the sign of κ ± (θ n ) ≡ β(θ n )q 2 t determines the topology in the particularr n (θ n )-direction. This new expression allows us to describe how the topological polarization changes at different angles. What we observe is that the modes flip from one side to the other as the directionr n (θ n ) crosses a particular set of directions (indicated by orange lines on Fig. 5). Those regions are called soft directions and we address their nature in the following section.

B. Soft Directions
Focusing on the long-wavelength zero modes established in Eq. (16), we find the existence of soft directions (indicated by orange lines on Fig. 5(b)), characteristic of the lattice [29] and determined by the value of the α ± coefficients: We can then re-write the determinant of a generic wavevector q in this new basis q = q +q+ + q −q− : det(q) = A 1,1 q + q − + iA 3,0 q 3 + + iA 2,1 q 2 + q − + . . . , (21) so that if a wavevector lies exactly on either soft direction (i.e q = q +q+ or = q −q− ), it is by definition a zero of our rigidity map at least to 2nd order.

VI. EXPERIMENTAL LENGTH AND ENERGY SCALES
The triumph of conventional elasticity is its ability to capture mechanical response at length scales extending to system size, far beyond those of the underlying, unobserved atomic interactions. It is not clear a priori whether the boundary modes which we consider extend to macroscale systems or whether they are valid only on the "atomic" length scales of the unit cell, those already captured by the lattice theory. Indeed, the dramatic change in edge stiffness predicted by simple central-force models [30] have resulted in relatively modest differential stiffnesses in 3D-printed systems [31]. In fact, as we show FIG. 6. Shape of the edge modes in a system with polarization ∆N = 2. We impose a mode with wavelength λ = 15.6 |l1| on the boundaryrt (where l1 is the first lattice primitive vector). This mode then propagates through the bulk and decays along the directionrn. Each mode (red, blue) varies sinusoidally along its corresponding soft direction and therefore also varies along the boundary with a longer wavelength. The modes decay into the bulk over length scales much longer than this wavelength.
here, the boundary modes extend to wavelengths intermediate between the unit cell and the system size, establishing new criteria for experimentally realizing strong topological effects.
Let us consider in particular the imposition of a distortion on a boundary with wavenumber q t on an interface with tangent directionr t (θ n ). In order to minimize energy, the system will undergo a distortion associated with one (or both) of the two soft directions described in the preceding section. We thus search the space of zero modes described in Eq. (21) for a zero mode with the appropriate component along the surface tangent direction, resulting in a wavevector of the form where |q + ,q − | denotes the determinant of the matrix of the given columns. Expressed in terms of the angles θ ± , θ n of the soft directions and the normal direction, this is We note that the decay part of this expression includes parameters of the systems that come from the determi-nant of our rigidity mapping, i.e. the A terms of equation (23), and hence can't be measured by the bulk response.
We now wish to express this relationship in terms of unitless ratios (in brackets) between the quantities with units of the lattice length scale, the decay length and the wavelength of the surface distortion, λ = 2π/q t : From the above expression, we see that the number of cells over which the topological mode decays is generically on the order of the square of the number of cells over which it extends on the surface. For example, if the wavelength extends over thousands of unit cells, the boundary mode generically decays on the scale of millions of cells. Hence, the surface theory extends from microscopic wavelengths up to the order of the geometric mean of the system size and unit cell size. However, this geometric relationship is strongly modified by the direction of the interface. When the soft direction is nearly aligned with the normal direction, the wavelength of the modes in the bulk becomes much smaller than that on the boundary, and the decay length is correspondingly reduced, meaning that the effect may be observed in smaller systems.
Beyond these concerns of length scale, the surface theory also relies on appropriate energy scales. While a full energetic analysis is beyond the scope of the present work, consideration of microscopic interactions and length scales permits us to estimate additional criteria to observe continuum topological polarization. Real systems contain additional constraints beyond those included in the Maxwell rigidity map. In a 3D-printed frame, in additional to strong energetic costs to centralforce compressions and extensions, beams possess finite bending stiffness not found in central-force springs. In our topological modes, these bending costs are imposed on a length scale proportional to the square of the wavelength, whereas a conventional Rayleigh-type surface mode would impose greater central-force energy costs over a smaller volume, extending to a depth proportional only to wavelength. Hence, the topological surface modes should only be strictly observable on wavelengths sufficiently short that they lie below the conventional modes in energy.

VII. CONCLUSIONS
We have presented a new method for describing the energetic effects of imposed macroscopic strain fields and their gradients on a microstructure that undergoes microscopic relaxation. We arrive at separate expressions for bulk and surface energies, and establish for criticallycoordinated systems an elastic bulk-boundary correspon-dence between the bulk structure and topologically protected zero-energy modes. The underlying topological invariant establishes the relative numbers of zero modes on two opposing boundaries, as either number separately would rely on short-wavelength physics. This establishes a sort of mesoscale elasticity, with topological surface modes existing on length scales far beyond the unit cell but necessarily far below the system size. Notably, this may lead to experimental demonstrations of topological polarization based on macroscale strains even when microscopic structure remains unobserved. The reader is encouraged to compare our approach with the illuminating [32], which was posted independently during the preparation of this manuscript.
These extensions of the lattice theory to the continuum present several avenues for further study. Lattice theories have considered bulk response [20], topological defects [17], buckling failure [18] and fracture [21]. All of these phenomena may extend to the continuum in intriguing ways. Following the realization of topological states in kirigami sheets [19], incorporating curvature into our continuum model may yield new physics.
Intriguingly, recent work has also shown that topological boundary modes can arise even in amorphous systems due to mesoscale structure [33]. Similarly, our surface theory relies on the short-length structure, suggesting that it may extend even to non-periodic systems, though topological polarization still relies on breaking spatial inversion symmetry. If indeed such local structure can describe directional response, it may underlie artificial structures which program intricate mechanical responses [34][35][36], themselves inspired by biological allostery.
Nonlinearities too may prove tractable in the general continuum theory, as when they lead to topologically protected solitons in one-dimensional systems [25].
In the present work, we considered a critical balance between constraints and independent components of the strain field that was derived by the critical coordination of the microscopic system. It is an open question as to whether such theories may emerge in the continuum without being present in the microscopic system. Structures such as those consisting of rigid square pieces joined at corners have a nonlinear zero-energy dilation mode [4] and are well-described by continuum theories [37,38]. However, such systems rely on symmetry to achieve their deformation mode and do not have zero-energy boundary modes.
we apply the well-known result that such functional differentiation flips the sign of the gradients. This may be seen by explicitly writing out gradients of the delta functionals that are used in functional differentiation and integrating by parts. The result, in this case, is that leading to the desired result in real space: We now return to incorporating our system's translational invariance and (semi-)local interactions, such that In reciprocal space, such that ∂ → iq, we may then write the relationship between equilibrium and rigidity maps as Our linear relationships then have the form:

C. ENERGY AND SURFACE TERMS
In this section, we derive how the total energy of our system breaks down between surface and bulk terms, using the expression of our rigidity map.
The second term of Eq. (38b) is a bulk term and can't be simplified any further. However the first term can nicely be expressed on the surface only using the divergence theorem: (39b)

D. WINDING NUMBER
As described in the main text, the task of relating the number of modes on edges of the systems to their bulk systems reduced mathematically to counting the numbers of zeroes in the complex plane. Here, we supply mathematical details and quantify the error associated with the long-wavelength approximation. This section implicitly uses a length scale so that we may treat physical quantities with units of length as pure numbers.
Physically, the points in the complex plane correspond to complex wavevectors. We draw a contour as shown in Fig. 3 of the main text which encloses all zeroes close to the origin (of order of the wavevector, which we term here O( )) and avoids those far from the origin (of order 1). The latter are non-physical, since our theory only describes the long-wavelength limit. Due to the shape of the contour, zeroes in the upper half-plane, corresponding to modes on the left edge, are enclosed in a positive orientation, while those in the lower half-plane are enclosed in a negative one. Note also that the contour in question allows a branch cut to be defined stretching along the negative real axis, important for evaluating the phase of complex numbers.
Were we to retain all components of the contour, the Residue Theorem would ensure that we could exactly evaluate the number of zeroes enclosed by the contour (we choose a gauge in which no poles are present). However, in order to achieve bulk-boundary correspondence, we must identify the complex zeroes by only considering the bulk modes, which lie on the real axis. We thus neglect the curved components of the contour. We now show that this approximation nevertheless recovers the result given in the main text, up to a small, controlled error.
Without loss of generality, we may consider the phase added to our contour integral from each zero separately. As such, we wish to obtain the change in phase of a complex function of the form as z winds along a contour re iθ with θ going from 0 to either π or −π. Here, we will choose r to be much larger than long-wavelength zeroes (O( 1 )) but much smaller than the short-wavelength ones (O( 0 )).
Because the argument of the complex number is simply the imaginary part of its logarithm, we immediately obtain two results. First, the contribution of the straight sections of the contour, which we retain, is the expression given in the main text. Second, that the contribution from the curved parts which we neglect is Noting here that although the long-wavelength zeroes have real part O( 1 ) they have imaginary parts O( 2 ), we find that for |r| |z 0 | the error is of order O( 2 /r). In contrast, the error from the short-wavelength, spurious modes is of order O(r), assuming that our original contour actually encloses all of the long-wavelength zero modes. However, if this condition is not met, the error in the calculation of the change in phase increases abruptly to O(1). To minimize error, we should then select the bounds of our contour to minimize their size while ensuring that they capture the essential physics. In the main text, we make the natural choice for r to lie intermediate between the two regimes, r = 1/2 in units in which unit cell size is of order one, generating error of O( 1/2 ). More aggressive schemes, such as r = 3/4 , . 99 could more closely characterize the topological transition, but would require greater knowledge of the microscopic details or willingness to tolerate occasional large errors.
As discussed above, the error from the true determinant function may be obtained simply by summing over a few cases of these zeroes, and is thus of the same order. Thus, we have shown using complex analysis that in the long-wavelength limit the correspondence between the bulk structure and the imbalance of topological edge modes given in the main text is mathematically justified. We emphasize here that this is not because the bulk modes along the curved parts of the contour are negligible but because the choice of contour causes them to cancel out. Thus, bulk-boundary correspondence in the long-wavelength limit links together two boundaries at once, in contrast to the lattice theory.

E. SOFT DIRECTIONS, LENGTH SCALES AND CHOICE OF MATHEMATICAL BASIS
In the main text, we introduce a number of different bases for the wavevectors. First, we have a simple Cartesian basis (q x ,q y ). Second, we present a rotated version thereof, in order to accommodate boundaries with varying orientations, in terms of inward-facing normal directionq n and the tangent directionq t . Finally, we consider the two soft directions, along which there lie modes which satisfy the structural constraints to linear order, pointing along (non-orthogonal directions) (q + ,q − ). We use the symbols α ± , α ± (θ n ) to relate these: (α ±qx +q y ) = 1 1 + α 2 + (θ n ) (α + (θ n )q n +q t ) . (42) As discussed in the main text, the bulk structure requires that zero modes occur at wavevectors satisfying the following constraint: We now wish to consider a boundary mode which has component q t along the tangent direction and thus takes the form q nqn + q tqt . By requiring that this mode satisfies the above constraints to second order, we obtain a result for the coefficients α ± (θ n ) which define the soft directions: q n = α ± (θ n )q t + iβ ± (θ n )q 2 t .
In theq n ,q t basis, we know the expression of q. Let's assume that q is a zero of our map, therefore it has the following form: q = q t (α + (θ n )q n +q t ) + iκ + (θ n )q 2 tqn , We recognize that the leading-order contribution to this zero-energy edge mode is simply the bulk soft mode. Expressing it in this form explicitly, we may write: where q + represents the value of the mode along the soft direction, while q t is the value of the mode along the transverse direction of the normal. Consider now such mode which, let's say, lies to first order along theq + direction. For it to be a zero mode, it must contain a second order correction along theq − direction that satisfies the zero-energy condition, which in this basis has a simplified expression: det q = A 1,1 q + q − + i(A 3,0 q 3 + + A 2,1 q 2 + q − + A 1,2 q − q 2 + + A 0,3 q 3 − ) = 0. (47) Based on a method relying on perturbation theory and using trigonometric identities, we're able to express the value of such zero mode as a function of the system's parameters only: q = q t sin(θ + − θ n )q + + iq 2 t A 3,0 A 1,1 sin(θ + − θ − ) sin 3 (θ n − θ + )q n (48a) Finally, as we explain in the main, we find that some of the system's parameters are dependent of the size of the cell |l 1 |, particularly that the quantity A 3,0 /(A 1,1 |l 1 |) is dimensionless. In addition to expressing the wavevector q t in terms of its wavelength λ (q t = 2π/λ), we can then derive a expression for a zero mode as a function of the properties on the boundary and the associated decay length through the bulk ζ + :