Conformal Elasticity of Mechanism-Based Metamaterials

Deformations of conventional solids are described via elasticity, a classical field theory whose form is constrained by translational and rotational symmetries. However, flexible metamaterials often contain an additional approximate symmetry due to the presence of a designer soft strain pathway. Here we show that low energy deformations of designer dilational metamaterials will be governed by a novel field theory, conformal elasticity, in which the nonuniform, nonlinear deformations observed under generic loads correspond with the well-studied conformal maps. We validate this approach using experiments and finite element simulations and further show that such systems obey a holographic bulk-boundary principle, which enables an unprecedented analytic method to predict and control nonuniform, nonlinear deformations. This work both presents a novel method of precise deformation control and demonstrates a general principle in which mechanisms can generate special classes of soft deformations.


INTRODUCTION
Mechanical metamaterials use patterns of cuts, pores and folds to achieve nonlinear [1][2][3], programmable [4,5], polar [6][7][8] and other exotic behavior [9][10][11][12][13][14] in response to external forcing. Often these features rely on the careful arrangement of the cuts, pores and folds to emulate a mechanism, which is a special pathway of deformation that enables the metamaterial to change shape at very small (ideally zero) energy cost. For example, the canonical rotating square (RS) mechanism, which consists of perfectly hinged rigid squares, enables a uniform dilational motion (Fig. 1a) and has inspired the design of a range of metamaterials which collapse inward rather than bulge outward when compressed from the lateral direction [2,15,16]. While the dilational mechanism becomes a true zero energy motion in the limit of vanishing hinge size, for realistic metamaterials with finite hinges this uniform dilational motion is only observed for the particular case of a completely homogeneous loading condition. For generic, i.e. inhomogenous, loading conditions the response is more complex [17], and yet a general framework to describe the general nonuniform soft deformations of mechanism-based metamaterials is lacking.
Here, we aim to decode the nonuniform deformations of metamaterials based on a dilational mechanism. While elastic response of ordinary materials will locally be composed of some finite portion of shear (Fig. 1b), a dilational material will strive to expel shear everywhere in favor of the dramatically softer dilational strains. Hence, one may wonder whether the nonuniform deformations of dilational metamaterials may be locally composed of pure dilation with no shear? While not all spatial patterns of strain are realizable due to basic geometric restrictions from compatibility relations, conformal maps are well-known to constitute the full set of compatible smooth deformations which locally are composed of dila-tional strains only. These maps therefore provide a recipe to compose nonuniform soft modes from slow spatial variation of e.g. the RS mechanism (Fig. 1c). We therefore formulate our conformal hypothesis: under a generic and broad set of loading conditions, dilational metamaterials will respond with an angle-preserving conformal deformation (Fig. 1d), energetically much softer than conventional elastic strains.
In this work, we confirm this hypothesis using simulations, experiments, and a coarse grained elastic theory. We first show that the response of the RS metamaterial is indeed conformal at both global and local scales. We then present a reduced elasticity theory which facilitates analytic insight, including new methods of predicting nonlinear deformation. We then use the bulk-boundary correspondence principle obeyed by these metamaterials to introduce a recipe for on-demand activation of each soft configuration.

RESULTS
To test the hypothesis of conformal deformation, we investigate the elastic response of RS-based metamaterials at a range of hinge thicknesses using finite element simulations which preserve the intricate pore structure (see Methods). While the detailed microscopic strains may contain significant shear, we use the displacement of the square centers to extract coarse-grained shear strain and dilational strain fields, with the expectation that the latter should dominate as the hinges become small. Even under nonuniform strain, dilation indeed dominates over shear in simulations of three-point "bending" (Fig 2a), a local "dipole" dilation ( Fig 2b) and even when the system is subject to global "pure shear" via compression along one axis and expansion along the other (Fig 2c). Even for dilations that vary dramatically through space (Fig. 2d), c.

FIG. 1. Dilational Mechanism and Conformal Deformations. (a)
In the ideal rotating square mechanism, a structure of rigid squares (black) connected by frictionless hinges (red), may be dilated and contracted at zero energy cost when the squares are rotated opposite to their neighbors. (b) Applying fixed deformations (red arrows) at the boundary of a conventional elastic material leads to a spatially varying strain field that includes shear components that change the local angles between pieces of material, such as the initially perpendicular grid lines (blue, yellow). (c) In contrast, a pure dilational metamaterial designed around the RS mechanism, may accommodate the loading without shearing of the unit cells, so that even as grid lines rotate, the angles between them remain fixed. This angle-preserving behavior arises due to a local version of the RS mechanism, in which square elastic chunks of side length a rotate about small, flexible hinges of thickness to open and close pores according to α(x, y), the local linear dilation factor. (d), A fabricated sample of checkerboard material approximately preserves right angles (shown in yellow) under a generic nonlinear "foot" load, suggesting conformal deformation behavior.
the average amount of local shear compared to dilation is captured in the shear fraction, defined in Appendix 1, and varies from 10 −1.5 to 10 −3 , in contrast with values on the order of one in conventional structures (Fig 2f). The low amounts of shear also imply nearly preserved angles, which is the defining characteristic of conformal deformations. As confirmed in Fig. 2e,f, the average angle change and the shear fraction both decrease as the hinge thickness is reduced and the system approaches the ideal mechanism limit. Despite finite-size effects from the 16 × 16 unit cell lattice, the strong preservation of angles indicates that the response is locally conformal.  2. Deformations of dilational metamaterials are angle preserving. (a) A Finite-Element Simulation of a threepoint bending test (the "bending" load) locates force balanced (relaxed) states of the RS-based elastic metamaterial subject to displacements on the boundaries that are incompatible with uniform dilation. (b) The "dipole" load consists of closing a single central pore by displacement of points on opposite sides of the pore, thereby generating a localized dilation in the bulk of the material. (c) The "pure shear" load consists of compression along one axis and extension along the other, which is nevertheless compatible with a spatially varying pure dilation field. (d) The resulting dilation field under the "bending" load varies widely across the sample, including nonlinear compression and expansion. (e) Nevertheless, the angle changes |∆ψ| remain small, even in proportion to the average strain magnitude | |. Here, the histogram is shown on log-linear scale, with each square marker giving the fraction of unit cells with that angle, which tend to decrease as the hinge size decreases. The preservation of angles is the defining feature of conformal maps. (f ) For all loads (colored in correspondence with the box outlines in (a), (b), (c)) the dilation fraction (solid lines) is nearly unity, and the deviation from unity is captured by the shear fraction (dot dashed lines), which for linearly small shears is the complement of the dilation fraction, and at the smallest hinge size ranges between 0.001 − 0.05. Both the shear fraction and the normalized average angle change (dashed lines) decrease with decreasing hinge size. The plateau in the bending and pure shear loads can be attributed to non-continuum lattice effects.

Fitting Deformations with Conformal Maps
The observation of approximately conformal local behavior suggests that deformations of the elastic RS metamaterial may also correspond globally to a conformal map. This motivates an elegant formulation expressing positions in the plane as complex numbers r = (x, y) ↔ x + iy ≡ z which has previously been applied to elasticity in a variety of contexts [18]. In our case, it has the added utility that any exact conformal deformation takes the convenient form where f (z) is the deformed position of the material element initially located at z and the complex coefficients {C n } define the map. The ability to express such a map in this reduced series form, turns the search for the conformal map closest to our data into a linear algebra problem which is readily solved. Using only the first 20 C n coefficients, we use this procedure in Fig. 3a,b to identify conformal maps that are able to capture all but a small fraction ∆ conf (about 1%) of the observed displacements. Note here that, to appropriately assess the error in many different approximations of the observed square block displacements, we employ the fraction of variance unexplained ∆ 2 . This measure is equivalent to the mean square displacement error normalized by the mean square magnitude of displacements; an explicit formula for this is given in Appendix 2. Indeed, despite the nonzero local angle deviations, we find that RS-based metamaterials will respond to loading with a deformation that is closely matched with a conformal map even under nonlinear strains. This is in contrast with generic deformations, which fit poorly to equation (1) as shown in Fig. 3b and in general require terms with complex conjugation of z such as z 2z5 included in the expansion to be captured accurately.

Conformal Elasticity
Given this evidence that the RS metamaterial responds to loading with a near-conformal deformation, we present an elastic energy functional, which the system minimizes subject to the boundary conditions. Here, α is a (nonlinear) spatially varying field describing the dilation factor of the structure relative to its equilibrium, while s 1 and s 2 are the (linear) coarse pure shear and simple shear respectively, is the width of a hinge and a the width of a square piece. G 1,2 (α) are shear moduli, M (α) the dilation energy density and M (α) the modulus associated with spatial variations in the dilation. These fields may be defined more precisely in terms of the right Cauchy-Green deformation tensor C (the metric of deformation). When we choose to orient the axes of our reference space along the vectors connecting square centers to those of their neighbors they take the form where σ (1) and σ (3) are the first and third Pauli matrices.
While alternate analytic approaches to the nonlinear strain of anisotropic media have been established (e.g., [19]), this formalism is particularly convenient for the RS metamaterial, for which any additional terms such as couplings between dilation and shear are excluded by scaling arguments and by the the orthotropic symmetry (p4g) as presented in Appendix 5. This specific form in equation 2 may also be derived by a nonlinear coarse graining, shown in Appendix 4, which provides additional useful insight into the moduli. Although other continuum elasticity theories for the RS metamaterial have been developed in one [17] and two dimensions [11,12], and another theory for a bistable dilational material [13], this nonlinear analytic form, which centers the dilational mechanism and its gradients, has not been proposed previously. While this energy function holds a rare level of insight into nonlinear deformation, it is still quite difficult to solve analytically.
To gain more analytical traction, we employ perturbation theory in the limit in which the hinges are tiny relative to the unit cell size l/a 1 (i.e. the soft mechanism limit) and the material sample is composed of very many unit cells (i.e., the continuum limit). Here, the second two terms in equation 2 may be regarded as small perturbations, and the energy cost of deformations that include shear (first two terms of equation (2)) become prohibitively stiff by comparison. Therefore, in this limit, dilations will begin to act as a local symmetry and the conformal maps will constitute a degenerate space of ground states of the energy. Restricting our focus to this space of ground states, we recover a nonlinear notion of conformally invariant elastic theory, which was suggested by Sun. et al. [20] previously for the linear deformations of the Kagome lattice, but has otherwise been viewed as "unphysical" [21] until now.
With the nonperturbative part of the energy penalizing only shear, and the set of conformal maps forming a degenerate set of ground states, it would seem that an infinite number of conformal maps are equally likely to arise in response to generic loading such as that in the FEM simulations and experiments of Fig. 2. However, this is not the case, as the degeneracy is broken by the perturbative part of the energy functional which is now expressed in terms of the conformal map f : z → f (z) describing the local dilation and rotation, with f (z) ≡ αe iφ (equation (1)). While the usual problem of reducing to a constrained elasticity theory is typically done using Lagrange multipliers, conformal maps simply have s 1 = s 2 = 0 and α = |f |, allowing equation (6) to be obtained easily as the conformal limit of equation (2). The alternate route, utilizing Lagrange multipliers, is explored in Appendix 8, and yields useful information about stress. The energy in equation (6) arises purely from the last two terms in equation (2) and simultaneously breaks the conformal invariance and the ground state degeneracy, allowing predictions of specific conformal response to be generated; we refer to this procedure as "conformal elasticity". With this, we have reduced the difficult tensorial problem of conventional nonlinear elasticity of a material with pores down to a scalar theory which is much more analytically tractable. For small loading, minimizing the energy and thus predicting the deformation is reduced to a linear algebra problem, which is readily solved. As shown in Appendix 7, these predictions closely match the observed finite-element displacements with correlation coefficient R 2 ≈ 0.99−0.9999 for hinges with /a = 0.005. This result showcases both the accuracy of conformal elasticity and the mathematical convenience of conformal maps. The nonlinear deformation of a three-point "bridge" load can be well fit to a conformal map (magenta grid) chosen to minimize the average squared deviation between the mapped and observed displacements. (b) The fraction of variance unexplained ∆ 2 conf between observed deformations and conformal fits in both simulations and experiments is very small, demonstrating that the soft deformations correspond globally to conformal maps as well as locally to pure dilations. Experimentally, the fit is especially good in the unloading step, in which stresses arising from frictional interactions at the material boundary are able to relax. In contrast, the data for a conventional elastic material (pink) cannot be described as a conformal map. (c) According to the holographic property of conformal maps, the dilation field on the boundary of an arbitrarily chosen (gerrymandered) region, such as the one depicted in yellow for a simulation configuration, should uniquely determine the deformation within the region. (d) Indeed, the displacements inferred from the boundary dilations (orange arrows) match closely with those observed directly in simulations (blue arrows). (e) The observed displacements are now compared to the conformal maps inferred in this way, and the fraction of variance unexplained plotted as ∆ 2 inference . This demonstrates that the deformations are not merely conformal, but the precise conformal map is encoded on the boundary. Here, the predictions of experimental deformations are inferred from the full rectangular boundary, while predictions of FEM data are all inferred from the gerrymandered boundary points shown in (c) to show robustness to the choice of boundary. Note that the example chosen for display in (c,d) corresponds to the largest strain explored in the bending load, and according to (e) has the largest prediction error of the FEM data, and all other FEM data is predicted to even greater accuracy than depicted here.

Bulk-Boundary Correspondence Generates Accurate Nonlinear Analytic Predictions
Solving equation (6) analytically is only possible in the linear regime and relies on prior knowledge of the effective stiffnesses M andM . However, as shown in Fig. 3a,b, the conformal property itself extends to nonlinear deformations. We therefore devise a scheme to extend analytic predictions to the nonlinear regime, relying on a mathematical property of conformal maps: bulk-boundary correspondence. Within our formalism, this bulk-boundary correspondence works as follows: if we know the amount of dilation along the entire boundary of a section of RS material, we can predict analytically the local dilation everywhere in the material interior from the unique conformal map derivative f (z) consistent with the boundary conditions and can further integrate to infer the displacements. We illustrate the principle in Fig. 3c,d,e. We retrieve the dilation on an arbitrarily chosen yellow domain in Fig. 3c from the FEM data, use the bulkboundary correspondence to calculate the deformation inside this domain (See displacement field depicted by blue arrows in Fig. 3d), and show that these are in excellent agreement with the observed displacement data (See displacement field depicted by orange arrows in Fig. 3d). This method is robust to the arbitrary choice of bound-ary shape and captures > 99% of block displacements even at large strains (Fig 3e).

Bulk-Boundary Correspondence Generates a Method for Precise Deformation Control
The bulk-boundary correspondence can also be used for prescribing on-demand deformations. Given a target deformation field of some section of material, we can determine the suitable actuation pattern by simply examining the local target dilation at the material boundary. We illustrate this with three different actuation patterns in Fig. 4a,b,c. Within a simulation described in Appendix 10, we actuate the material at its boundary (the red bars represent actuators) and observe that the numerically force balanced deformation is in good agreement with the yellow dots that mark the target centers of the squares. Of course in addition to being conformal, this procedure is only able to actuate deformations that do not stretch or compress the material beyond the physical bounds of the mechanism itself. This task is aided by the maximum modulus principle of conformal maps, in which both the maximum and minimum dilations must occur on the material boundary. We can therefore guarantee that for any choice of actuation that varies slowly along the boundary, there is a physically valid soft conformal deformation that will be activated in the interior. We propose that such boundary-control may be achieved via the insertion of struts in the pores near the boundary, which, with remote control actuation, would allow these soft modes to be activated remotely on-demand. Exploration of this possibility, with applications in development of soft robotics, is reserved for future work.
While the deformations of Fig. 4 arise from relatively simple functions for the input boundary dilation patterns, each analytic form of the interior deformation turns out to be much more complicated. However, one limit is simple enough to offer intuitive insight. When the lengthscales at which the dilation varies are much shorter than those at which the boundary curves, the boundary may be treated as that of a semi-infinite plane. For conformal metamaterial occupying the infinite upper half plane (y > 0), and the dilations along the boundary (x-axis) constrained to satisfy α(x, y = 0) = α 0 (x), there is again only one allowed pattern of dilation, which takes the form where w(k) is the Fourier transform of the function w(x) = ln(α 0 (x)). In this form, we can see that the dilation (or, rather, the logarithm of the dilation) decays into the bulk exponentially according to the lengthscale of variation of the boundary dilation. For instance, when w(x) = w 0 cos(k 0 x), we will have α(x, y) = exp w 0 e −k0y cos(k 0 x) .

DISCUSSION
The intuitive concept of soft modes that are locally composed of a spatially varying mechanism is not confined to the RS mechanism explored here. The conformal modes explored here arise because they are the Nambu-Goldstone modes associated with the local dilational symmetry, rather than from the details of the specific microstructure. We therefore propose the analogy: as isometries are to thin elastic sheets [22], so are conformal deformations to dilational metamaterials. Consequently, dilational materials derived from a variety of different metastructures [20,[23][24][25][26] including fractal RS structures [25,26], disordered "pruned" networks [27], nanoscale [28], and three dimensional generalizations [3,29,30] will similarly have mechanical response controlled by a set of conformal soft modes. These alternate architectures provide the possibility for greater ranges of dilation [24], opening the door for even more dramatic soft deformations such as pictured for the kagome lattice in Fig 4d. Note, by example, that many possible loads applied to thin sheets (e.g. applying an in-plane stretch to a flat sheet) are incompatible with isometries, and an isometric soft mode theory must break down there. Similarly, a variety of overly strict loads will be incompatible with conformal deformation and our conformal elasticity theory must similarly break down. For instance, this will happen when attempting to constrain the displacements, rather than dilations all along a closed boundary. Conveniently, we have identified two particular scenarios with either a finite number of sufficiently spaced point displacements, or with dilations controlled along a closed boundary, each of which contain a broad variety of loading possibilities which are guaranteed to be compatible with a conformal map, and are therefore governed by our elasticity theory. Identifying theories that apply beyond the soft conformal response will provide an interesting avenue for future work.
We suggest that a broad class of generic mechanisms not confined to pure dilation will also generate families of soft modes that govern material response, as were indeed observed in so-called "kirigami" structures [31] and that this may become a fundamental principle for mechanism-based metamaterials, with potential applications from footwear to soft robotics. In addition, topologically polarized systems [6,32], which have additional mathematical structure controlling an additional set of exotic boundary modes, necessarily break the symmetry of the mechanism strain pathway that allows conformal modes and the notion of mechanism compatibility may be extended to their nondilational mechanisms while incorporating topological notions. Exploring how these new classes of soft modes may obey, e.g., a generalized bulkboundary correspondence may yield fruitful connections to exotic field theories [33] as well as black hole and string theories where holographic principles already play a vital role.

a.
b. c. d.

FIG. 4.
Bulk-Boundary correspondence allows for precise control of metamaterial deformation. (a, b, c) The bulk-boundary correspondence shown in Fig. 3 enables on-demand actuation of material deformation via the addition of boundary springs (red). Here, the locations of square centers observed in simulation match well with analytically obtained target positions (yellow dots) for force-balanced configurations of a simplified ball-and-spring model. (d) Because the principles presented here extend to any dilational metamaterial, the kagome lattice may likewise be placed in a low-energy conformal deformation. The kagome lattice is capable of undergoing a broader range of deformations without self-intersection, permitting a greater range of controllable deformations such as pictured here.

Generating deformations analytically
We analytically generate the deformations of twodimensional metamaterials that locally resemble the uniform mechanism and a rotation and hence eliminate the dominant contributions to the elastic energy, taking advantage of a map between the plane of deformation and the complex plane. In general, a deformation map x → X(x) has deformation tensor F ij ≡ ∂Xi ∂xj which deforms infinitesimal material vectors as dx i → F ij dx j . To correspond locally to a rotation of the mechanism field, the deformation tensor must be proportionate to a rotation matrix. We seek solutions in which the rotation angle φ and the dilation factor α vary over space, so that F = αR(φ). However, this generates nontrivial restrictions because of the requirement that the deformed state not tear, i.e. that dX = 0 along any path. By applying Stokes' theorem, this leads to the standard geometric or kinematic compatibility conditions: This generates two independent conditions on the two fields, suggesting a unique solution subject to certain boundary conditions. However, it becomes convenient to define new variables z = x + iy,z = x − iy, with i the imaginary unit. Enforcing the conditions ∂ z z = ∂zz = 1, ∂ zz = ∂zz = 0 then requires that By summing the second compatibility condition with (−i) times the first, we find, following algebra, that compatibility is equivalent to the complex condition ∂z [α exp(iφ)] = 0. Thus, the locally dilational deformations of the metamaterial are exactly the conformal maps of the plane, which may be expressed as analytic functions f (z) whose magnitude is the local dilation and whose argument is the local orientation. Displacements can be obtained via complex integration: f (z) = dzf (z). This method hints at a more general class of spatially varying deformations in unimode materials, but the results in this case follow from the fact that the shear-free deformations of the structure are precisely the angle-preserving ones well-known to be described by complex analytic functions.

Finite Element Simulation Protocol
For the 2D finite-elements simulations, we use the commercial software Abaqus/Standard. We use a neo-Hookean energy density as a material model, a shear modulus µ = 0.333 MPa*m, bulk modulus K 0 = 16.7 GPa*m and plane strain conditions with hybrid quadratic triangular elements (Abaqus type CPE6H). We construct the mesh so that the thinnest parts of the samples are at least two elements across. We perform two types of simulations (i) on the metamaterials unit cells with periodic boundary conditions, on which we perform a low strain static test (5 × 10 −6 )) and (ii) on the full structure, on which we perform static nonlinear analysis.
In the simulations (i), to implement periodic boundary conditions, we define constraints on the displacements of all of the nodes at the horizontal and vertical boundaries of the unit cell and impose displacement using virtual nodes [34]. We perform four types of simulations to extract the macroscopic linear bulk modulus , and shear moduli G 1 (α = 1), G 2 (α = 1). In order to extract K and G 1 , we apply biaxial compression ε comp . G 2 is extracted using the identical procedure as for G 1 applied to a unit cell that is rotated by π/4. We examine a unit cell with a = 12mm, pretwist θ 0 = 0.39 and hinge thickness = 0.1mm, finding (K, G 1 , G 2 ) ∼ (35, 5.2 * 10 4 , 2 * 10 4 )Pa*m. These modulus measurements are essential for obtaining the perturbative energy minimizing prediction of deformation, summarized briefly in the main text and given in more detail in Appendix 7.
In the simulations (ii), we use three kinds of boundary conditions, dipole, bending and pure shear, by imposing displacement of specific nodes, as described in Fig. 2. We extract the coordinates (position, angle) of all the squares and use these to compute the dilation and shear field (method described in Appendix 1) and to calculate the relative amount of shear, as well as the accuracy of conformal map fits and of the analytic bulk-boundary correspondence (Fig. 3). Simulations are performed for a series of strains at identical material geometry as (i) and separately at linearly small strain across a range of hinge thicknesses = {0.02, 0.033, 0.055, 0.092, 0.155, 0.258, 0.431, 0.719, 1.199, 2.0} mm.

Experimental Protocol
We 3D printed a elongated meta-beam of length 306mm, width 64mm and height 40mm, consisting of a checkerboard of lattice of 48 × 10 squares of side a = 4.8mm connected by hinges of thickness = 0.2mm using a Connex 500 stratasys 3D printed and the proprietary Stratatys Agilus 30 material (30 Shore A). This material is viscoelastic, its Young's modulus at short times is E= 3.3 MPa and E= 0.6 MPa at long times (see, e.g., Dykstra et al. JAM 2019 [35] for a calibration). For imaging purposes, we 3D printed black square-shaped pads (3.6mm) on one edge of the sample (proprietary Stratasys material: Vero Black).
The sample was tested under two sets of non-uniform boundary conditions, (i) by a foot-shaped indenter, which was 3D printed in ABS Materials using a Dimension 3D printer (Stratasys); (ii) a three points bending (i.e. bridge) test, using a universal uniaxial mechanical testing machine (Instron 5943). In parallel to the mechanical loading, we recorded frames using a high-resolution camera (6000 × 4000 pixels (Nikon D5600) with a telephoto lens (Nikon 200mm, F4), positioned at 4 meters from the sample. To ensure uniform lighting, the sample was illuminated by two led panels (Bresser, LG 900 54W) and white paper was used to ensure a uniform bright background. The 3D printed material tends to highly frictional and even adhesive, therefore, we used fine powder to lubricate the interaction with the bottom surface and the foot for (i). For (i), we performed a compression-decompression cycle at a rate of 0.1 mm/s up to a maximum 20 mm stroke, during which 400 frames were recorded (1 frame / sec). For (ii), we performed a compression-decompression cycle at a rate of 0.2 mm/s up to a maximum 40 mm stroke, during which 400 frames were recorded (1 frame / sec). The images were processed using standard image-tesselation and tracking techniques to extract the positions of the squares with subpixel detection accuracy of 0.02 pix (1 µm). We attribute the strongly different behavior observed between compression and decompression, to a combination of frictional, viscoelastic and self-adhesion effects.

Extracting Dilation and Shear Strains
To extract nonlinear measures of coarse dilation and shear strain magnitudes from deformation of RS-based structures, we track the square center displacements. Using the vectors connecting square centers to neighboring square centers as the infinitesimal material vectors in the reference (dx (1) , dx (2) ) and target spaces (dX (1) , dX (2) ), we may infer the deformation tensor F ij = ∂Xi ∂xj . As the reference dx (1) and dx (2) are orthogonal, the recipe is simply Then the induced metric of deformation is obtained via C = F T · F. Nonlinear dilation and shear strain magnitudes may be obtained from the invariants of the metric tensor via These recipes, which are derived in greater detail in Appendix 1, give the dilation as the local area change, while the shear comes from the nonconformal part of the deformation as a complex map h : z → h(z,z), satisfying Fitting displacement data to the closest conformal map To identify the the conformal map f (z) which most closely maps a set of n p points {z i } (expressed in complex form) to final positions {z i }, we minimize the error Utilizing the analytic expansion of f (equation (1)), we may find the set of coefficients {C n } which minimize the error by setting the partial derivatives of the error to zero, yielding equations where we have defined the matrix and the coefficients are now cut off at a maximum n c . This effectively reduces the least-squares error analysis to a straightforward linear algebra problem, which is readily solved. We note that the cutoff n c should be chosen to be much less than the n p data points, to avoid over-fitting.
To evaluate the accuracy of the fit, we define a similar function to the error in equation 12 where − z are the displacements proposed to fit the data. The functional ∆ 2 [ ], quantifies the "fraction of unexplained variance" in the displacements, and constitutes a general method to quantify the amount of the deformation captured by a candidate conformal function, used in Fig. 3b,e as well as in the Appendices.
Inferring the interior conformal deformation from boundary dilation data In our new bulk-boundary method for predicting deformation, we are able to infer the conformal deformation f : z → f (z) which will arise from a discrete set of M applied dilations {α k } at a corresponding set of locations {z k } along the boundary of a simply-connected planar domain. Here, α k = A k /A k is ratio of local area elements before A k and after A i deformation. The inference takes place in two steps: first an intermediate function g(z) = log(f (z)) is obtained, then this function is used to integrate for the displacements.
Because f (z) will be conformal, so is the function g(z) = log(f (z)) and therefore admits a similar expansion where the cutoff N should be significantly less than half the number of boundary points M in order to avoid overfitting. Note now that since g(z) = log(α) + iφ, the real part of g alone determines the dilation field. We may therefore infer the g that fits the boundary conditions by minimizing the error Inserting equation (16), we minimize this error with respect to the coefficients D n = A n + iB n yielding the equations where and we have expressed z k = r k e iθ k in a complex polar form. equation (18) reduces the inference of g to a linear algebra problem which may be readily solved using builtin tools in, e.g., Mathematica. Noted also in Appendix 9, the row and column corresponding to B 0 = φ 0 (the undetermined global rotation) are zero throughout and this row-column pair must be removed before numerically solving.
As described in Appendix 9, the function g(z), specified by the coefficients {D n }, is used to reconstruct the z-derivative of the conformal map via ∂ z f = Exp(g). Generating the predictions of block displacements u(z) = f (z) − z, shown in main text Fig. 3d and evaluated in Fig. 3e, is a matter of complex integration, which we accomplish using the built-in numerical integrators in Mathematica.

DATA AVAILABILITY
The experimental and finite element simulation data generated in this study have been deposited in the Zenodo database under accession code https://doi.org/ 10.5281/zenodo.4646672.

COMPETING INTERESTS STATEMENT
The authors declare no competing interests.

APPENDICES Appendix 1. Quantifying Local Dilations and Shears
To our knowledge, a parameter which quantifies what fraction of a nonlinear deformation is composed of local dilation, as opposed to local shear, has not been introduced previously. We therefore present a more detailed derivation in the context of Finite Strain Theory. This is an essential ingredient for supporting our hypothesis in the main text that nonuniform deformations of a dilational metamaterial are nonetheless composed primarily of local dilations.
The method described here applies to data describing deformation from a reference space X to a target space x(X), which are smooth enough to take derivatives. In the case of the RS metamaterial, the square centers provide a discrete approximation to a smooth deformation, which may then be used to define local material derivatives. More generally, one may use the lattice vectors of a material with repeating metastructure as material vectors to take these material derivatives. From these derivatives we assemble the deformation tensor For a smooth deformation, this tensor contains all required information to quantify dilation and shear at each point in space.
Using the polar decomposition theorem, we may write F = R · U, where R is a rotation and U is the symmetric "right stretch tensor". From this form, it is clear that an infinitesimal material vector in the reference space dX will deform first via multiplication with U, then will be rotated via R to reach a final state dx = R · U · dX. Since U is symmetric, there exists a local choice of coordinate system where U is diagonal and acts on the local elements dX as a combination of just dilation and pure shear. The eigenvalues (λ 1 , λ 2 ) of U are called the principal stretches and the local material dilation and shear magnitudes may be written intuitively in terms of them. An infinitesimal square sample of side lengths ( , ) oriented in the diagonal frame of U will transform into a rectangle with side lengths (λ 1 , λ 2 ) after deformation. Using this, the local dilation may be defined intuitively in terms of the area where ∆A is the change in area of our infinitesimal element and A 0 = 2 is the initial area in the reference space. For a similarly intuitive notion of shear, we turn to the anisotropy of the tensor U. We define the nonlinear shear magnitude in terms of the difference of the principal stretches as s 2 = (λ 1 − λ 2 ) 2 . We will see that this definition matches with conventional intuition for standard simple shear nonlinearly as well as general mixed linear shear. Importantly, this shear may be rewritten in terms of invariants of U as s 2 = Tr[U] 2 − 4Det[U], and is therefore also a rotationally invariant quantity.
Since the stretch tensor is generally less convenient to extract numerically, we show that this local dilation and shear may also be written in terms of the trace and determinant of the metric (equivalently the "right Cauchy-Green deformation tensor"). The metric is defined as The coordinate system which diagonalizes U also diagonalizes C, and it follows that the eigenvalues of the metric are the squared eigenvalues of U. Using this locally diagonalized frame, we may extract general relations between rotationally invariant quantities Using these relations, we may write the more useful formulas completing our recipe to extract the local dilation and shear magnitudes from the metric tensor. That these relations in equation (24) discriminate conformal from non-conformal deformation is not entirely obvious. We therefore present them in alternative form by connecting to the complex formulation where the vector mapping x → x(X) is captured by the complex function z → f (z,z). In this case, the components of the metric may be conveniently rewritten as where ∂ z = 1/2(∂ x − i∂ y ) and ∂z = 1/2(∂ x + i∂ y ) are standard complex derivatives [18]. With these relations, our metrics of conformal and non-conformal become where, remembering that the conformal condition is ∂zf = 0, the shear part now more transparently quantifies the non-conformal part of the deformation. While both dilation and shear are strain quantities, there is some ambiguity in comparing them nonlinearly. The expressions used here, composed entirely in terms of invariants of C, are constructed to be intuitive and geometric. In addition, rewriting in terms of the components of the displacement u = x − X reveals that our definition of the shear matches with the combined magnitude of pure and simple shears (first and second terms in equation 26 respectively) from linear strain theory. Note that, in the linear limit, these strain measures satisfy 1/2(d 2 + s 2 ) = ij ij . While alternate definitions of shear may be identified which also satisfy this in the linear limit, we have explored these and they do not qualitatively change the results presented here. Finally, for a conventional simple shear u = γyx, our formula reassuringly recovers s 2 = γ 2 to nonlinear order. As shown in Fig. 5, these formulae reveal that the typical magnitude of local dilation is much larger than that of local shear. While the shear may become large in small regions near the applied loads and where the dilation changes over short distances, dilation still dominates by nearly an order of magnitude at the worst. Nonuniform, nonlinear deformations of the RS-based mechanism in force-balance may then indeed be composed of local dilations with comparatively very little shear. We now identify a simple fitting method of locating the closest conformal map describing a discrete set of displacement data points. Per Appendix 1, the discrete data should be chosen carefully to approximate a smooth deformation. For lattice materials, this means tracking identical points in neighboring unit cells. For the RSbased metamaterial, we can choose the centers of the square elements.
The method consists of first noting that a general conformal map f (z) on a simple finite planar domain with no holes (i.e. simply-connected) may be expressed precisely as an analytic function with a series expansion In contrast, a generic non-conformal deformation must be expressed as a double sum over all products of z and the complex conjugatez. The complex expansion coefficients C n provide a complete description of the map and solving for them is therefore the goal of our analysis. Finding the nearest conformal map is a matter of minimizing the square error function where the sum runs over the n p discrete data points z i → f i . Inserting the form of f (z) (equation (27)) and minimizing with respect to the C n , we find  where we have defined the matrix and the coefficients are now cut off at a maximum n c . This effectively reduces the least-squares error analysis to a straightforward linear algebra problem, which is readily solved. We note that the cutoff n c should be chosen to be much less than the n p data points, to avoid over-fitting.
To evaluate the accuracy of the fit, we define a function where u i = f i − z i is the observed displacement data, u(z) ≡ f (z) − z are the displacements proposed to fit the data. The functional ∆ 2 [ ], known as the fraction of variance unexplained, is then a general method to quantify the amount of the deformation captured by a candidate conformal function, which we use in the main text and in other appendices. The method evaluates error using the displacements rather than positions in order to avoid a false inflation of the results at small deformation. Inserting the numerically solved closest fit into equation (31), we show in Fig. 6 that only around 20 coefficients are needed to fit the deformations of the RS metamaterial before the improvement starts to become negligible. This is good, as the more coefficients that are included, the more computationally expensive the fitting problem becomes. We remark that, in practice, application of this method works best if the coordinate system is first centered and rescaled so that the data points reside closer to unity away from the origin. This helps to avoid large and small number issues in the numerics, which includes higher and higher exponents of the positions as more coefficients are included.
This constitutes a fast and accurate method to extract the nearest conformal map from position data. Note that, while the method is reduced to linear algebra, no assumptions of linearity of the strains in either the deformation, or the fitting map, have been made. Indeed, this method shows the nonlinear deformations of the RS metamaterial are well-approximated as a conformal map to nonlinear order.

FIG. 6.
The accuracy of the conformal map with the increasing number of coefficients allowed in the expansion equation (27). Solid lines correspond to fits with FEM data at the smallest strains, while the dotted lines are for the largest strains explored in the main text analysis. In comparison to the conventional elastic material in the absence of pores, the fitting quickly reaches very small errors ∆ 2 conf which is the function ∆ 2 [ ] evaluated for the particular candidate fit u(z). For loads with discrete rotational symmetry, the errors go down in steps, since only every few coefficients may be non-zero.

Appendix 3. Energetics of Hinge Deformations
We now aim to derive a continuum elasticity theory for the RS metamaterial. However, to consider the RS elastic material at full detail requires allowance of all possible deformations of the porous elastic structure. Such a tensor field theory for deformations would be complicated and unwieldy, leaving very little room for analytic progress. Fortunately, due to the careful design of the pores, we may think of the material as composed of square elements of side length a connected by comparatively small "hinge" elements of thickness . This enables an effective description of the material as a collection of hinges, written in terms of simplified variables which prove to be very convenient for coarse-graining. We establish this effective "hinge-based" energetic formalism below.
In the limit that a, the hinge will become very flexible compared to the stiff square pieces. We therefore follow successful previous examples [11,17] in which the large square elements are approximated as rigid bodies and all strain deformation is assumed to take place at the hinges. This assumption is readily verified in the Finite Element deformation data, where the stress is generically found to be dramatically localized to the hinge region, as shown in Fig. 7.
To quantify the energetics in this effective description, we consider a single hinge connecting two elastic blocks as shown in Fig. 8-a. Knowledge of the deformation tensor F hinge ij (r) throughout the hinge completely determines the system configuration, as stress and strain are already assumed to be zero elsewhere. Here, r and R, respectively, are target and reference coordinate systems describing the fine structure of the hinge deformations, in contrast to x and X which will describe a coarsened picture of continuum deformations of the metamaterial at the mesoscale. Given the complete constitutive relation for the elastic material, the tensor field F hinge ij (r) provides sufficient information to compute the energy density at each point in space. For generality we write this simply as an unknown function e(r) = φ(F hinge ij (r)). We now confine our focus to hinges with aspect ratio ∼ 1, which is the case for the experiments and simulations in this manuscript. This condition generally removes opportunities for bistability and other pathdependent elastic behavior at the hinge. Therefore, prescribing the orientations and positions of the two squares attached to a particular hinge will leave only one possible force-balanced deformation F hinge ij (r) in the hinge. The energy of this hinge may therefore be rewritten as a function of the two square orientations and two 2d vector positions, six numbers altogether. Accounting for translational and rotational symmetry leaves only three scalar parameters necessary to describe the configuration of the hinge. We will focus on two distinct, yet convenient choices for these three parameters as described below.
Importantly, we note that arbitrary relative positions and orientations of a pair of squares will generally cost a great deal of energy and may stretch the hinge to the point of fracture. We wish to focus on low-energy configurations that do not stretch the hinge too much, which includes a particular set of nonlinear rearrangements due to the design based on an ideal-hinge mechanism. Therefore, in order to find an appropriate parametrization for these low-energy configurations, we must first quantify the finite-hinge (finite-energy) analogue of this mechanism. As shown in Fig. 8a,b, fixing the distance between the square centers (c) and allowing the square orientations to relax defines our twisting angle T (c) ≡ π/4 − ψ(c)/2. Matching with the common nomenclature for twisting in kagome and RS lattices, this represents the amount the squares must each be rotated away from the fully extended state to balance torque at fixed distance. To twist the squares oppositely according to T (c) while moving the square centers to distance c therefore constitutes our generalized mechanism. In the limit that the hinge size goes toward zero, and the strain is increasingly localized to an infinitesimal hinge, this motion must approach the ideal frictionless hinge version of the mechanism. We may therefore write the twisting function as where δ T (c) is an unknown dimensionless function capturing the first-order corrections to the twisting due to the finite size and stiffness of the hinge. Eventually we will proceed to ignore this first order correction, but must acknowledge that the function δ T (c) may diverge in the vicinity of c = √ 2a (the untwisted state). Both our hinge theory and subsequent coarse-graining of the Rotating Square (RS) lattice mechanics will break down near this untwisted state and we therefore leave analysis of this anomalous point in the configuration space for future work.
Alternatively, we may explore this mechanism by fixing the orientation ψ = π/2 − 2T and allowing the positions to relax to a square center spacing c(T ) which can be expressed as Again, the second term on the right-hand side is a correction term dependent on an unknown function δ c which will eventually be neglected in the small hinge limit. We now use these mechanism twisting functions to define two parametrizations of the hinge energy. In both cases, we first choose to move to a point along the mechanism pathway, setting the distance c and twisting angle to ψ = π/2 − 2T (c) as shown in Fig. 8 (left). We may then deviate from this configuration in two ways: In the first, we may choose to fix the square orientations and allow a small vector deviation δc in the relative positions as shown in the top right of Fig. 8. We will refer to this as the "vector parametrization". Alternatively, we could keep the square centers fixed in place and allow each to have a small reorientation away from the mechanistic twisting orientation as shown in Fig. 8. This defines deviation angles δθ A and δθ B and we will refer to this as the "angular parametrization". In the end, the choices will span the same configurations of the hinge structure.
Both parametrizations above are useful for isolating the "non-mechanistic" part of the hinge motion into small variables. It is clear for small hinges that comparatively large δc will put the hinge under large strains which scale roughly as ∼ |δc| and will put conventional materials at risk of tearing and fracture. We therefore work in the limit |δc| 1. Equivalently, (δθ A , δθ B ) 1 are the appropriate conditions in the angular parametrization. We will see below that the vector parametrization is useful for deriving scaling relations, while the angular parametrization is useful for the coarse-graining. Given some configuration of the hinge we may now write the hinge energy in terms of the vector parameters or equivalently in terms of the angular parameters where T 0 is the "pretwist" present in the reference configuration of the hinge ensemble. Finally, we note that it is useful to be able to convert between these descriptions. Neglecting second and higher order terms in the deviations δc and dimensionless hinge size /a, we may write Even with perfect knowledge of the hinge geometry and nonlinear material elasticity, finding the form of the functions g v , g a would require solving a difficult nonlinear mechanics problem. However, a lot can be discerned without knowing the exact form of the energy and instead relying only on scaling arguments. In the vector parametrization, our hinge deformation tensor is a function not only of the location within the hinge x, but also this tensor field must depend on the inputs T and δc. It is instructive to consider rescaling the hinge reference and final states by a factor γ while preserving shape. This means that points in the hinge r are mapped to pointsr = γr and the deformation tensor is preserved viã F hinge ij (r) = F hinge ij (r). Because the deformation tensor is preserved, so is the energy density. The only difference in computing the total energy of the rescaled hinge comes from integration over a rescaled domain, and we find E H = γ 2 E H . However, this rescaling also impacts the relative locations of the large square elements. Consider a. c. b.

d.
FIG. 8. In the absence of additional forces applied near the hinge, the elastic hinge energetics may be parametrized with three numbers describing the relative orientations and positions of the squares. First, the "mechanism" portion (a → c) may be traversed by fixing the center-to-center distance c0 → c of the squares and allowing the square orientations to relax. Then small deviations around this state are allowed by either (d) preserving the square distance and allowing small orientation changes (δθA, δθB) or (b) preserving the square orientations and allowing a small change (δc) in relative square positions. the internal vector e which points from the projected corner of square A to the projected corner of square B as shown in the zoom in Fig. 8-b. Due to the preserved shape of the hinge, this internal vector must also scale directly with γ. This internal vector also has an analytic form where δ c (T ) is again the correction to the mechanism from equation (33). Considering the configurations where δc = 0 in equation (37), we note that because e scales with γ, and because also scales directly with γ, then δ T must not scale with γ. Noting that δ T is independent of δd, it will remain scale-free when we also apply such a displacement δc. Because each side of equation (37) must scale the same way with γ, we have now established the result that the displacement vector δc indeed must scale linearly with γ.
Rescaling the hinge has no impact on the square relative orientations ψ nor on the "pre-twist" T 0 . This rescal-ing impacts only and δc and the hinge total energy E H as described above (in the vector parametrization). We may therefore rewrite the hinge energy function where g 1 is another unknown function, whose value may be thought of as a mean energy density in the hinge. Another way of arriving at equation (38) is the following: given we had started with a smaller hinge of the same initial shape, we must arrive at the same final force-balanced hinge shape by simply rescaling the applied displacement vector δc by the same factor as . Therefore, we must be able to write the hinge energy as in equation (38) as a function of the current twist, the pretwist, and the value of δc relative to .
We now apply our previous assumption that we will avoid extreme, high energy configurations of the hinge where material fracture becomes likely. Quantitatively, this is the condition δc 1. We may use this condition to expand the energy in equation (38) to second order, yielding The twist-dependent coefficients may be written in terms of g 1 via Using equations 36 & 40, we can convert this energy to the angular parametrization. This becomes wherec = c √ 2a is the dimensionless rescaling of the square center distance and the moduli are defined by g(c) = g T (T ( √ 2ac)) = g T (ArcCos(cCos(T 0 ))) While the function g 1 remains unknown, and the energies in equations 35 & 34 may not be evaluated directly, these expressions provide useful insight into the nonlinear macroscopic energetics of the RS metamaterial while excluding extremal configurations, as shown in Appendix 4. Importantly, as the function g 1 defined in equation (38) is scale-free, there is no remaining "hidden" dependance on the hinge size nor the square size a in the energy function and the scaling is directly apparent. We have therefore avoided a complicated tensor-field description of our structured elastic metamaterial, in favor of a hinge-based energetic formalism based on a greatly simplified set of degrees of freedom: the positions and orientations of the square elements. In this formalism, energetics and scaling permit analytic descriptions which prove helpful for coarse graining. We now derive the continuum elasticity theory for the rotating-square (RS) lattice, as presented in the main text. Following Appendix 3, we approximate the elastic RS metamaterial as a collection of rigid square elements connected by comparatively flexible hinges. The hinges and square elements may be made out of the same material, while the comparative flexibility arises due to the small size of the hinge.
We establish the local continuum energy penalty created by a mapping x → X(x). This mapping will control the location of the square centers, while the orentations are left free to realize torque balanced configurations. Following the main text result that soft deformations of the RS metamaterial are near-conformal, we restrict our focus to deformations composed of a large nonlinear dilation, a finite local rotation, and a small generic shear. This prescription is captured by the deformation tensor are Pauli matrices and R(φ) a standard rotation by φ.
To understand this tensor, we break it down piece by piece. First, recall that by definition this tensor describes the Jacobian of the deformation map X → x(X) via F ij = ∂xi ∂Xj . The infinitesimal material vectors are therefore understood to transform under deformation via dX → dx = F · dX. Consider, then, the action of each element of F going from right to left, in the order that they act on dX. The first (rightmost, in square brackets) term applies a linear strain composed solely of shears s 1 (X) and s 2 (X) defined as where u i is the i-th component of displacement, defined at this intermediate deformation for the purposes of understanding. Defined in this way, s 1 and s 2 correspond to the traditional notions of pure shear and simple shear, respectively. Following the application of the linear shear, the material vectors undergo an isotropic dilation α(X) and are then rotated by φ(X), bringing us to the final configuration. While this is a general form for such a (large dilation plus small shears) deformation, we note that these new variables φ, α, s 1 , s 2 cannot vary arbitrarily through space, as they must still obey a closure condition to guarantee geometric (mechanical) compatibility. This is because a generic, spatially varying candidate strain field may not correspond to a valid displacement field, just as a generic candidate electrical field does not always correspond to an electrostatic potential. In the vanishing shear limit the closure condition takes the form where ij is the 2-by-2 antisymmetric unit tensor. Reassuringly, equation (46) captures the Cauchy-Riemann condition for a deformation to be conformal. To see this, note that the derivative of a conformal function may be written in modulus-argument form as . Then, as conformal functions must be functions of z only, we must have ∂zf = 0. Writing out the real and imaginary parts of this condition leads to the xand y-components of equation (46). We now choose the vectors c ν connecting neighboring square centers (and where ν = {1, 2, 3, 4} indexes which neighbor) to represent the infinitesimal vector elements of the effective continuum material, and thus F ij prescribes all the relative square locations after deformation. If we can now find a similar recipe for writing down the square orientations in terms of continuum quantities, then our description of the metamaterial configuration is effectively complete. We will then be able to write down the hinge energies and sum these over the unit cell to get a continuum energy density. To identify the torquebalanced square orientations, we note that there are only two squares per unit cell. We therefore assume that square orientations are controlled by two angular fields which vary slowly through space. A unit cell contains squares A and B which after deformation are located at r A and r B and have orientations θ A and θ B . There is already some intuition for the expected orientations in terms of the counter-rotation mechanism of the squares arising from the local lattice dilation factor α(X). This is expressed by writing the orientations as where T () is the same mechanism function from Appendix 3, while ∆T, ∆φ are the two smooth fields required to capture the two orientation degrees of freedom in each unit cell. Note that while equation (47) relies on a well-motivated guess, there is no loss of generality due to the presence of the correction fields ∆T, ∆φ, and conveniently that if our guess is correct, then enforcing torque balance will simply set the value of these fields locally to zero. In the limit of vanishing shears and vanishing gradients of α and φ, we find that ∆T and ∆φ also vanish and therefore the magnitude of these angular correction fields must be on the order of the shears or ∇φ or smaller.
We are now equipped with a sufficient description to specify the energy (equation 41) of an arbitrary hinge connecting the square initially centered at X to the square initially centered at X + c 0 ν in terms of our continuum quantities (φ, α, s 1 , s 2 , ∆T, ∆φ). Constructing the required geometric quantities to evaluate the hinge energy, we first find via application of the deformation tensor that the length of c ν changes as whereê 1 ,ê 2 are the orthogonal basis vectors (x-and ydirections, respectively) in the reference space. The deformation also produces a reorientation of this local infinitesimal via This information, combined with equation (47), allows us to write the angles of deviation δθ (ν) (50) which will be input to determine the energy of hinge ν. For convenience, we have moved to the complex formulation rather than the vector formulation where z ≡ X 1 + iX 2 , and similar for vector quantities such as c ν . Our continuum energy density may now be constructed by summing the hinge energy equation (35) over the four hinges in the unit cell via where A cell = 4a 2 cos(T 0 ) is the area of the unit cell in the reference space and E H is the energy from equation (41) for hinge ν. Inserting equations 50 and performing the sum yields The energy in equation (52) includes all terms which are not higher order in either /a or the shears or |∇α|. To reduce to an elasticity theory in terms of spatial deformation variables, we minimize equation (52) with respect to the internal degrees of freedom ∆T and ∆φ. This leads to ∆T = 0 and ∆φ = 0, which indicates that our mechanism-based guess for the square orientations in equation (47) was correct. The continuum elastic theory for the RS lattice is then capture the energy cost of dilation, pure shear, simple shear, and dilation gradients, respectively. Note that in the small hinge limit, the strain gradient modulus may be inferred from the shear moduli via the simple relatioñ where we have used the analogue for the mechanism in equation (32) to replace excluding higher order terms proportional to /a. Here, we present an alternate derivation of an energy functional of the form equation (53) without recourse to the exact RS microstructure, based instead upon its fundamental symmetries.
We consider, then, a general elastic system that: • has a dilational mechanism • has standard elastic translational and rotational symmetries • has a four-fold rotational symmetry • has an x → −x mirror symmetry • has a y → −y mirror symmetry (as implied by the previous two symmetries) These last three properties arise from the p4g wallpaper symmetry group, which applies to all states of the lattice along the uniform dilational mechanism. A uniform application of the mechanism will generate an energy density that we again label ( /a) 2 M (α), having identified the scaling with hinge thickness in Appendix 3. Additionally, energy terms may arise from other components of the deformation tensor: pure shear (s 1 ), simple shear (s 2 ) and rotation (φ). However, rotation in particular can be easily eliminated: the system is invariant under a uniform rotation and equation (46) (Cauchy-Riemann) shows that for dilation-dominated deformations, rotation gradients can be expressed in terms of dilation gradients instead (up to a small correction from the shears). Overall, our energy density can therefore be expressed in terms of three fields: (s 1 , s 2 , α) and gradients thereof. In general, this allows for terms proportional to s 2 1 , s 1 s 2 , ∂ x α∂ y α etc. Terms such as ∂ x s 1 ∂ y s 1 and s 2 ∂ 2 x α are also permitted, but these will be higher order in gradients and shears which are small due to the comparative size of the unit cell and design based on a mechanism, respectively. Note that, for present convenience, we are using notation different from the previous appendix, where x → X 1 and y → X 2 .
Such energy terms are only permitted if they are invariant under the discrete symmetries of the undeformed or uniformly dilated lattice: a four-fold rotation about the center of a square and two mirror symmetries about the center of an open pore.
The action of the symmetry operations on the relevant fields are listed in Table. I. Enforcing that the energy remain the same before and after applying the s1 s2 ∂xα ∂yα C4 −s1 −s2 −∂yα ∂xα x-mirror s1 −s2 −∂xα ∂yα y-mirror s1 −s2 ∂xα −∂yα x-mirror symmetry eliminates (s 2 , ∂ x α, s 1 s 2 , ∂ x α∂ y α, s 2 ∂ y α, s 1 ∂ x α ) from the energy density. Similarly, the y-mirror symmetry eliminates the terms (∂ y α, s 2 ∂ x α, s 1 ∂ y α). Finally, the C 4 symmetry eliminates S 1 , and further enforces that the energetic coefficient of the (∂ x α) 2 term match the coefficient of (∂ y α) 2 . We do not consider terms involving second gradients of the dilation field, which are ruled out in the coarse-grained microscopic theory.
With this, the energy functional to lowest order in the small shears and dilation gradients must take the form (60) where the coefficients have been named to connect to the coarse grained theory in equation (53).
Importantly, while we have not invoked the specific nature of the RS mechanism or material properties to derive this elastic energy, is it specific to the point symmetry group of the lattice and a different energetic form is expected for other dilational metamaterials such as the kagome lattice.

Appendix 6. Stress of Near-Conformal Deformations of the RS Metamaterial
Having constructed the Energy in equation 53, we would like to know the corresponding stress tensor, as this is the more common quantity used in nonlinear elasticity. However, this theory is both nonlinear and includes strain gradient terms, and we require a relation which appropriately incorporates these effects. To do this, we start with the essential information that the functional derivative with respect to displacement gives the force density as In addition, we know that the force density is the divergence of a stress tensor and so the process of deriving the stress is a matter of taking the functional derivative in equation 61 and fitting it into the form of equation 62. Our energy functional is an integral over an energy density which is a function of our quantities α and s 1 and s 2 . These are strain quantities, and may be expressed in terms of the Lagrangian strain where C ij = F ki F kj is the right Cauchy-Green deformation tensor (i.e. the metric of deformation) and F ij = ∂ j u i is the deformation gradient tensor. The variables in our energy functional equation 53 may be written, to lowest order in s 1 , s 2 , in terms of this strain via Therefore, we should think of our energy functional as having the form With this, taking the functional derivative in equation 61 is a matter of using the chain rule for functional derivatives Here, jk may be thought of as a functional of u i and this derivative is taken using the definition of the functional derivative which becomes And using the form of equation 63, we find (71) Plugging this back into equation 68, and simplifying a bit, we may find where is the standard functional derivative. As the divergence of this stress in terms of the reference space coordinates gives the force density in reference space, this construction captures the Nominal stress (also known as Engineering stress and corresponding to the transpose of the first Piola-Kirchhoff stress) in accordance with the nonlinear elasticity literature [36]. Using the standard conversion formulae taking us between the different standard nonlinear stress definitions, we find that the second Piola-Kirchhoff stress is given by the simple expression Inserting the actual energy from equation 53, we will find While the 2nd Piola-Kirchhoff doesn't admit a direct physical interpretation in terms of real traction forces across material surfaces, it is mathematically convenient and can easily be converted to the more physically relevant Cauchy stress with the formula

Appendix 7. Analytic prediction of deformation
Having identified the effective continuum energy which will govern the deformation of the RS metamaterial, we would like to use it to predict material response to applied loads. Taking the limits of small hinge and small lattice spacing in equation (53) leads to an energy which is comparatively very stiff against deformations that include shears, and the shear terms dominate the energy functional. In this limit, the material will choose to expel shear whenever possible, leading directly to the conformal maps (which definitionally exclude local shear) as the lowest energy deformations. With only shear terms in the energy, the conformal deformations form a highly degenerate space of ground states. The point displacements applied in the Finite Element simulations from Appendix 1 are not sufficient to fully break the degeneracy -an infinite space of deformations satisfy the (point) boundary conditions without any shears, leaving no single prediction for the deformation.
To break this degeneracy, we turn to perturbation theory, adding back in the dilation M (α) and dilation gra-dientM (α) terms. Predicting the deformation then becomes a search for the conformal map which minimizes the perturbative energy Shown in Fig. 9-a, this reduced form is a reasonable approximation to the true energy, with the shear energy roughly an order of magnitude smaller at small hinge size. This energy gap will be even more significant, and the perturbative approximation better, for larger material samples with more unit cells which are closer to the continuum limit. As described in the main text, this process of obtaining and using this effective theory constitutes our conformal elasticity.
a. Analytic prediction of linear deformation using effective conformal theory The energy in equation (78) is highly nonlinear and minimization requires prior knowledge of the nonlinear stiffnesses k 1 (α), k 2 (α), g(α), which are not easily obtained. Even with knowledge of these moduli, analytic solutions to such nonlinear mechanics problems are hard to come by. However, taking the theory to lowest order in strains and rewriting in terms of a complex displacement u(z) = f (z) − z = u x (z) + iu y (z) yields where ∂ z = 0.5(∂ x − i∂ y ) is a standard complex derivative [18] and are the bulk and bulk gradient moduli, respectively. Note the unknown function g, defined in Appendix 3, which captures the mechanism energy. In order to incorporate the applied point motions (constraints), Lagrange multiplier terms have been added to the energy functional as b. a.  (53) to equation (78) is a good approximation by showing that the contribution to the energy from the shear terms is an order of magnitude smaller (or more) in the limit of small hinges. Energies are normalized by E0, which is the amount of energy required to create the same coarse strain magnitude in the absence of the pore structure, and which showcases also the comparative softness of the porous metamaterial. The energy in equation (78) may also be used to generate predictions of small deformations which have very small fractional mean squared error shown in (b).
where {u k } are the displacements prescribed at locations {z k }.
The energy in equation (79) will be used to select between conformal displacements u which may be written as an expansion u(z) = n C n z n . As higher coefficients C n generate sharper and sharper features, we can cut off the expansion and then solving the equations of energy minimization is reduced to a linear algebra problem.
The above is not yet sufficient to generate predictions, as the moduli K andK remain undetermined. Conveniently, the coarse-graining procedure from Appendix 4 offers a relation between these moduli which is exact in the limit of small hinges (82) With this, the strain level moduli (K, G 1 , G 2 ) obtained numerically in the methods section of the main text constitute all necessary information to make linear deformation predictions. For simulations with hinge thickness 0 = 0.1mm these moduli are (K,K) = (35Pa*m, 125Pa*m 3 ). To find predictions at different hinge thicknesses , we rely on the coarse graining results, which indicate that K ∼ 2 while all other moduli G 1 , G 2 ,K are independent of . Therefore, we write K( ) = 2 2 0 K( 0 ) in terms of the known K( 0 ) and use these values to solve the linear algebra problem for the coefficients {C n } and therefore our prediction for u(z). These predictions capture all but a small portion ∆ 2 pred of the error, as displayed in Fig. 9 In this Appendix, we employ a more conventional method of enforcing the shear-free conditions on the RS metamaterial by the introduction of Lagrange multipliers λ 1 , λ 2 . This produces insight into the shear stresses enforcing this shear-free condition. For simplicity, we consider a system without gradient terms, so that the energy takes the form where F ij is the deformation tensor, J = α 2 is the determinant of its matrix form and the general dependence on the reference coordinate has been suppressed.
The conditions for equilibrium are that the energy cannot be lowered by any movement of the material from its target position, r(R). This can be expressed as a functional derivative (See, e.g., Lazar and Kirchner [37] or Saremi and Rocklin [8]): Consequently, after some algebra, we obtain our equilibrium conditions: b (J) [F 11 ∂ 1 J + F 12 ∂ 2 J] + ∂ 1 λ 1 + ∂ 2 λ 2 = 0, (86) b (J) [F 11 ∂ 2 J − F 12 ∂ 1 J] − ∂ 2 λ 1 + ∂ 1 λ 2 = 0. (87) In obtaining this expression, we have already used that F 22 = F 11 , F 21 = −F 12 (the constraints). To this, we now add our compatibility condition, as discussed in the Methods section of the main text. We find that this compatibility condition is equivalent to or We now note that J = |f | 2 and we can take our two original real equilibrium conditions and add the first equation to the second equation multiplied by i, and obtain the equivalent complex condition b (J)|f 2 |f + ∂ z λ = 0, (93) λ(z,z) ≡ λ 1 + iλ 2 . (94) governing the constrained equilibrium states. Thus, there are stresses that are supported even in the absence of deformation (f (z) = 0): the structure supports shear stresses that are functions ofz only. That is, they are antiholomorphic, whereas the low-energy deformations are holomorphic, or conformal. In addition, there is an additional shear stress proportional to b (J) that forms when there is a nonzero bulk modulus and a spatially varying deformation tensor. Similar results arise more generally when dilation gradient terms are included, still allowing for an undetermined antiholomorphic space of shear stresses.
• integrate to find f (z) We remark that the analytic viability of this recipe depends on the ability to generate solutions to the Laplace equation on the chosen domain, and to integrate them. The examples of square and circular domains provide convenient solvable examples. However, the true strength of our method lies in the reverse application. Given a desired conformal deformation f (z), this may be prescribed by simply taking the derivative and evaluating the magnitude at the boundary, yielding the appropriate boundary dilation pattern to actuate. It is the uniqueness of the inverse recipe shown here (up to overall translations and rotations) which guarantees only the target map can satisfy these boundary conditions.
Finally, the particular choice of mechanism sets limitations to the amount the metamaterial can be dilated and contracted to avoid self intersection and material overextension. For instance, the RS mechanism requires α to be everywhere confined to the range 1/ √ 2 < α < 1, when starting from the fully dilated configuration. However, very conveniently, it turns out that any sufficiently smooth choice of boundary dilations which do not overor under-dilate the RS lattice on the boundary, will also generate a conformal deformation in the interior which does not go beyond the valid dilation range of the specific mechanism. This is due to the maximum modulus principle which states that both the maximum and minimum dilations of a conformal map will occur at the boundary.

a. Discrete Inference Method
For most material domain shapes and boundary dilation patterns, a closed-form solution for f (z) is not guaranteed. We therefore illustrate an alternative scheme that allows for quick and accurate inference of a map from discrete boundary data, and which generates the predictions in main text Fig. 4.
In this method, we first identify the boundary via a set of discrete points {z k } and the corresponding local dilations {α k }. Points should be densely spaced compared to the lengthscale of variation of α k . And, naturally, they should approximately trace a single continuous path enclosing a simply-connected region of the plane. Following this, we choose a cutoff N in the number of coefficients of N should be significantly less than half the number of boundary points M in order to avoid overfitting.
Next we demand that our function g(z) satisfy the boundary conditions by minimizing the error err = M k (Re[g(z k )] − ln(α k )) 2 . (100) Inserting equation (99), we minimize this error with respect to the coefficients C n = A n + iB n yielding the equations composed purely of Hookean springs of identical stiffness k connected at frictionless nodes. As shown in Fig. 10, each rigid square element (grey) is emulated by a grouping of six springs. While a single cross-spring is sufficient to render a square element rigid, the inclusion of both ensures that the spring ensemble will obey the same symmetry properties as the elastic RS metamaterial.
To actuate a particular soft conformal mode, stiff springs are added at the boundary, as shown in Fig. 10. These boundary springs are taken to have stiffness 10 4 k, so that they will act as rigid constraints, realizing their rest lengths quite accurately compared to the soft mechanics of the bulk material. While the rest lengths in the interior ( blue bonds) are chosen to leave the square shape at zero energy, the rest lengths of the boundary springs are varied to match some target local dilation. For three patterns of boundary dilation, we set these rest lengths and identify force balanced states using the conjugate gradient numerical minimization procedure "minimize" from the scipy.optimize toolkit in python. At the same time, using the methods from Appendix 9, and as illustrated in the main text Fig.3d&e, predictions for displacements may be generated from this input set of boundary dilations. As shown in main text Fig.4, these are in good qualitative agreement with the numerically identified force balanced configurations.