# Orisometry formalism reveals duality and exotic nonuniform response in
origami sheets

Michael Czajkowski  Dept. of Physics, Georgia Institute of Technology  James
McInerney  Dept. of Physics, University of Michigan, Ann Arbor, MI 48109, USA
Andrew M. Wu  Dept. of Physics, Georgia Institute of Technology  D. Zeb
Rocklin  Dept. of Physics, Georgia Institute of Technology

###### Abstract

Origami metamaterial design enables drastic qualitative changes in the
response properties of a thin sheet via the addition of a repeating pattern of
folds based around a rigid folding motion. Known also as a mechanism, this
folding motion will have a very small energy cost when applied uniformly; and
yet uniform activation of such remains highly difficult to observe, these
sheets instead generically displaying nonuniform response patterns which are
not yet well understood. Here, we present a purely geometric continuum theory
which captures the nonuniform, nonlinear response to generic loading as
composed locally of the planar mechanism, as well as previously identified
“twist” and “bend” modes which enable the patterned sheet to curve out of the
plane across long distances. Our numerical analysis confirms that these three
modes govern the observed nonuniform response, varying smoothly across the
sheet according to three PDEs which guarantee compatibility. In analogy with
the recently solved case of planar mechanism metamaterials, these
“Orisometries” (origami + isometries, so named by us) are subextensive but
infinite in number, with each mode displaying, in the linear limit, “sheared
analytic” spatial patterns which are controlled by the Poisson’s ratio of the
uniform folding mechanism. Furthermore, the “planar” mechanism-based
deformation patterns superimpose with a mathematically dual space of “non-
planar” twist/bend deformations to span the available soft linear response.
Together, our findings furnish the first quantification of the number of soft
response modes available, as well as the first intuitive quantification of
their spatial distribution.

##  I Introduction

A conventional uncreased sheet of paper may be readily deformed into
cylindrical and conical shapes with minimal applied force, but is not so
easily stretched or sheared along its plane. This familiar behavior does not
arise due to the specific composition of the paper, but rather is a geometric
property which extends beyond paper: the thinner an elastic sheet, the more
difficult it quickly becomes to stretch in comparison to bending Jellett
(1849). Deformations composed of bending without (in-plane) strain dominate
the response properties of thin sheets; these are known as _isometries_ , the
study of which extends over a century and a half Jellett (1849); Cerda and
Mahadevan (1998); Witten (2007) illuminating the otherwise cryptic nonlinear
deformations that are observed.

More recently, precise patterns of folding (i.e. Fig. 1a) have been revealed
which endow a conventional sheet with unconventional stretching capabilities
MIURA and Koryo (1985); Tachi (2009); Schenk and Guest (2013); Pratapa  _et
al._ (2019); Dieleman  _et al._ (2019); Feng  _et al._ (2020); McInerney  _et
al._ (2022). This is achieved via the design based on a rigid folding
_mechanism_ : an elaborate pathway of deformation which costs zero energy in
the limit that the creases act as perfect frictionless hinges. When one
ignores the fine details of the folding pattern in favor of the coarse shape
changes of the sheet, such carefully corrugated sheets appear to subvert the
inextensibility condition of isometric deformations: the sheets readily extend
and change shape in the plane that they approximate (see Fig. 1b). It is
important to note that, while this mechanism folding motion is in fact
isometric in the details of the crease microstructure, this is not a smoothly
differentiable isometry and the motion leads to highly non-isometric behavior
in the coarse macroscopic picture. As a result, the static response of a
mechanism-based creased sheet is in sharp qualitative contrast with the
response of the uncreased counterpart under identical loading (see Fig. 1c).
It then seems, in this coarse description, that the mathematics of the smooth
isometries will no longer form a useful description of the lowest energy
response available to a sheet with such a crease pattern. What then is the
spatial nature of these origami-based soft deformations? Do they admit a
differentiable continuum mechanics description?

Answers to these questions reach considerably beyond individual folded sheets.
Similar to the case with isometries, the exotic response properties of
patterned origami sheets tend to arise from geometry rather than from the
specifics of scale or material composition Mahadevan and Rica (2005); Audoly
and Boudaoud (2008); Schenk and Guest (2013); Christensen  _et al._ (2015);
Bertoldi  _et al._ (2017). Such geometric insights into these materials are
therefore particularly powerful because they will apply to origami at the
submicroscopic scale of molecule just as well as the scale of a skyscraper
Brittain  _et al._ (2001); Temmerman  _et al._ (2007); Randhawa  _et al._
(2011); Rogers  _et al._ (2016); Li  _et al._ (2019). As a result, these
geometric approaches have enabled a broad variety of applications including in
roboticsOkuzaki  _et al._ (2008); Hawkes  _et al._ (2010); Randhawa  _et al._
(2011); Jeong and Lee (2018); Lazarus  _et al._ (2019); Son  _et al._ (2022),
surgery An and Rus (2012); Felton  _et al._ (2014); Miyashita  _et al._
(2016), and even in the development of deployable shelters for the unhoused
Temmerman  _et al._ (2007).

Despite all the advances they have enabled, the mechanisms themselves do not
tell the complete story of geometry’s role in origami response. For instance,
although the quasi-planar mechanism motion programmed into these origami
sheets is ideally zero energy, it remains quite difficult to uniformly
activate Misseroni  _et al._ (2022). Instead, a variety of nonuniform
responses are observed Misseroni  _et al._ (2022), and unreliable folding is
pervasive even beyond these mechanisms Waitukaitis  _et al._ (2015); Stern
_et al._ (2017); Chen and Santangelo (2018); Misseroni  _et al._ (2022); Lee-
Trimble  _et al._ (2022).

Meanwhile, a recent series of investigations has begun to establish that for
any uniform planar mechanism there must exist an associated space of
nonuniform deformation modes which also cost zero energy in the continuum
limit (i.e. with infinitesimally fine folding structure) Czajkowski and
Rocklin (2021); Czajkowski  _et al._ (2022); Czajkowski and Rocklin (2022);
Zheng  _et al._ (2022). This suggests that paired with the quasi-planar
mechanisms, such as those of the the canonical Miura-ori and eggbox patterns,
there will atleast be an affiliated space of quasi-planar nonuniform soft
modes. In addition to these modes, it is well established that ubiquitous
“twist” and “bend” modes enable such mechanism-based folded sheets to curl
over long distances into nontrivial three dimensional shapes, again at zero
energy cost in the continuum limit. These additional geometric effects have
been utilized by Nassar et al. Nassar  _et al._ (2017, 2022) to identify a
striking variety of nontrivial, nonuniform soft deformations, and Xu et al. Xu
_et al._ (2023) have contributed further analytic configurations along with
proof of softness and a continuum mechanics theory. However, a complete theory
revealing the full spectrum of soft response in mechanism-based origami sheets
is lacking. And so we refine the question addressed in this work: Does there
exist a mathematical formulation analogous to the isometries which
characterizes the coarse response of these sheets?

Figure 1:  (a) The “Morph” family of crease patterns are defined by two
parallelograms, at angles $\alpha$ and $\beta$, and admit a quasi-planar
folding motion (b), parameterized by the angle $\theta$ which, compared to the
initial reference folding state at $\theta_{0}$, enables apparent stretching
of the surface the sheet approximates. (c) Two sheets under identical
displacement loading (red arrows) and with the opposite corners held in place
(red dots). In comparison to an uncreased sheet (top), a sheet which is
creased according to a folding mechanism, (bottom) will display drastically
different nonlinear nonuniform response, whose coarse shape change does not
resemble any smooth isometry.

Here, we answer this question by assembling these varied perspectives on
mechanism-based origami sheets to reveal a unified theory controlling the
available soft deformations therein. Just as isometries of uncreased sheets
take control as the thickness becomes small, the geometric limit of finely-
detailed folding patterns (i.e. the continuum limit) reveals the role of our
_orisometries_ in the response of mechanism-based origami sheets. In this
theory, the three previously identified low-energy modes available to the
sheet control deformation at the local level of just a few unit cells, rather
than uniformly across the entire sheet. This confines the sheet to a
“hyperribbon” (so named by our use here) in the space of the sheet fundamental
forms. Following related previous work Nassar  _et al._ (2017, 2022), we allow
these modes to vary slowly across the sheet, subject to spatial constraints,
which reveals an infinite but subextensive space of available soft response.
Furthermore, we find that the small linear deformations in this theory divide
into the (quasi-) planar motions and the non-planar motions which are not only
independent, but are mathematically dual to one another. These deformations
are then controlled by an exceptional point in the system Poisson’s ratio
which divides bulk versus boundary soft mode patterns according to the auxetic
versus anauxetic folding motion character. Our analytic results are confirmed
using a variety of numerical force-balancing simulations within a particularly
useful yet still broad class of “Morph”Pratapa  _et al._ (2019) fold patterns,
altogether furnishing a governing framework which has remarkable similarities
to the recent investigations in planar unimode metamaterials Czajkowski and
Rocklin (2022); Zheng  _et al._ (2022).

##  II Local soft motions of mechanism-based origami sheets

The investigation presented here concerns the “Morph” class of folding
patterns Pratapa  _et al._ (2019), illustrated diagrammatically in Fig. 1a,
and discussed in Appendix section A which includes the canonically studied
examples of the Miura-ori and eggbox folds. We first enumerate the essential
local ingredients of soft deformation which inform the subsequent analysis.
These folding patterns fall into the category of lattice origami, in which the
undeformed folding pattern repeats itself upon translation by any integer
combination of two lattice vectors
$\mathbf{l}^{(1)}_{0},\mathbf{l}^{(2)}_{0}$. For loading applied at
lengthscales much longer than the unit cell size (the focus of this work), a
coarse continuum description becomes useful. This is achieved by the origami
sheet “midplane”, depicted in Fig. 2 (bottom), which smoothly connects the
unit cell centers. In this picture, the lattice vectors locally approximate
the tangent vectors of the sheet, and travel along these vectors by continuous
amounts $X_{1},X_{2}$ allows the midplane surface $\mathbf{R}(X_{1},X_{2})$ of
the undeformed and corresponding deformed sheet to be parameterized.
Specifically, this parameterization is defined so that
$\mathbf{l}^{(i)}=\partial_{i}\mathbf{R}$ and differential geometric
quantities such as the metric
$g_{ij}\equiv\partial_{i}\mathbf{R}\cdot\partial_{j}\mathbf{R}$ and curvature
$h_{ij}\equiv-\partial_{i}\hat{\mathbf{N}}\cdot\partial_{j}\mathbf{R}$
(arising from gradients in the surface normal $\hat{\mathbf{N}}(X_{1},X_{2})$)
induced by the deformation may be estimated both numerically and analytically
in the continuum limit approximation.

Figure 2:  The origami fold patterns locally have access to a planar mode (the
mechanism), as well as a “bend” and a “twist” mode. Nonuniform patterns of
these three modes comprise the nonuniform response to generic nonuniform
loading, as approximated by the sheet midplane.

As the Morph patterns are contained within the recently studied “four-
parallelogram” class McInerney  _et al._ (2022), they each exhibit a nonlinear
quasi-planar folding mechanism as displayed in Fig. 1a. Quasi-planar means
that this motion does not generate any curvature of the midplane surface.
Instead, the mechanism motion induces a change in the surface metric away from
the identity

$g_{ij}=\frac{\mathbf{l}^{(i)}(\theta)\cdot\mathbf{l}^{(j)}(\theta)}{|\mathbf{l}^{(i)}(\theta_{0})||\mathbf{l}^{(j)}(\theta_{0})|}\rightarrow\begin{bmatrix}\lambda_{1}^{2}(\theta)&0\\\
0&\lambda_{2}^{2}(\theta)\end{bmatrix}$ (1)

according to the change in the folding parameter $\theta$ away from an initial
reference (rest) value $\theta_{0}$ shown in Fig. 1a,b. Note that the surface
metric and curvature tensors capture what are known as the first and second
fundamental forms, and we use the language interchangeably in this work. Here,
we have introduced the principal stretches $\lambda_{1},\lambda_{2}$, known
from finite strain theory, which here capture the relative extension of the
lattice vectors so that
$\mathbf{l}^{(1)}(\theta)=\lambda_{1}(\theta)*\mathbf{l}^{(1)}(\theta_{0})$;
this description accommodates the nonlinear changes such folded sheets are
capable of. Note that the expression in this form is specific to the
symmetries of the Morph patterns, and the metric need not remain diagonal
along the mechanism motion in all cases. Furthermore, considering a small
increment of the mechanism folding allows us to identify the planar Poisson’s
ratio

$\nu_{\text{in}}(\theta_{0})\equiv-\frac{\lambda_{2}^{\prime}(\theta_{0})}{\lambda_{1}^{\prime}(\theta_{0})}$
(2)

at the chosen reference state of this mechanism.

Figure 3:  (a) In the space of fundamental forms that describe the sheet
deformation, response is confined to the “hyperribbon” we introduce here, with
one long direction along the nonlinear planar mode and two short (linear)
directions along the bend and twist modes. This is illustrated in two slices
of the six dimensional fundamental form space for the numerically deformed
configuration above. (d) Quantified root-mean-square distance from the
hyperribbon (normalized by the root-mean-square spread along the hyperribbon)
in both the strain (solid lines) and curvature (dashed lines) subspaces across
a variety of fold patterns with different rest state Poisson’s ratios
$\nu_{\text{in}}$. Due to the normalization, the smallness of these quantities
on an absolute scale indicates a high quality of adherence to the hyperribbon,
which continues as the average strain is amplified to nonlinear values.

When the origami sheet panels are perfectly rigid, The quasi-planar motion is
the only possible deformation. However, for typical origami sheet material
which is thin and flexible, additional isometric deformations allow each panel
to bend along the diagonals. This enables two additional modes of deformation,
referred to as the “twist” and “bend” modes, which are known to only influence
the second fundamental form of the sheet McInerney  _et al._ (2022). These
modes, described in detail in Appendix. B only extend to linear order before
generating energy-costly stretching of panels and may be superimposed together
to generate the coarse curvature of the sheet

$h_{ij}\rightarrow
b\begin{bmatrix}\lambda_{1}(\theta)\lambda_{1}^{\prime}(\theta)&0\\\
0&-\lambda_{2}(\theta)\lambda_{2}^{\prime}(\theta)\end{bmatrix}+t\begin{bmatrix}0&1\\\
1&0\end{bmatrix}\,,$ (3)

where $b$ and $t$ are amplitudes of the bend and twist modes respectively. The
expression for the bending curvature in Eq. 3 is determined from the
established result that the bending mode generates an out-of-plane Poisson’s
ratio $\nu_{\text{out}}\equiv-\frac{h_{22}}{h_{11}}$ which is equal and
opposite to the in-plane ratio $\nu_{\text{out}}=-\nu_{\text{in}}$ Schenk and
Guest (2013); Wei  _et al._ (2013); Pratapa  _et al._ (2019); McInerney  _et
al._ (2022).

##  III Origami sheet deformations adhere to a shape “hyperribbon”

When applied uniformly, each of the three modes established in the previous
section incur only a small energy cost proportional to the crease folding
stiffness $k_{f}$ and panel bending stiffness $k_{p}$ without panel stretching
$k_{s}$. As discussed in greater detail in Appendix Section. C.2, the geometry
of the thin sheet dictates a qualitative difference between these moduli
$k_{s}\gg k_{f},k_{p}$, and hence a qualitative difference between these soft
modes and other arbitrary uniform deformations which include panel stretching.
However, promoting these soft mode amplitudes $\theta,t,b$ to scalar fields
$\theta(X_{1},X_{2}),t(X_{1},X_{2}),b(X_{1},X_{2})$ and allowing these modes
to vary slowly across the sheet may impose an additional energy cost (see Fig.
2). For very slow gradients, the additional cost must be proportional to the
ratio of the unit cell size $\sim a$ to the lengthscale of variation of the
mode $\zeta$. It follows that in the continuum limit where $\zeta\gg a$, this
additional energy from spatial variations must go to zero and such nonuniform
soft mode deformations will again have (low) energy proportional to the soft
crease and panel stiffnesses. Whenever possible, the sheet will therefore use
these slowly varying soft modes to respond to loading, rather than more
energy-costly deformations which stretch the constitutive sheet.

According to the rigidity theorem of Bonnet Toponogov and Rovenski (2006), the
shape of any smooth surface is completely described by the first and second
fundamental forms (captured by the components of the metric and curvature
tensor, respectively) across the sheet. Each point on the surface may be
assigned a 6-dimensional coordinate (3 coordinates for each of the two
symmetric tensors), and the entire sheet may be interpreted as a volume in
this 6-dimensional space of these forms. Therefore, a smooth differentiable
surface may be viewed as a continuous volume extending in the 6-dimensional
space of the components of these forms, each point on the surface having
coordinates ($g_{11},g_{12},g_{22},h_{11},h_{12},h_{22}$). However, for the
response of an origami sheet, the lowest energy possibilities will be
described by the amplitudes of planar($\theta$), twist($t$), and bend($b$)
modes. Therefore, such response must be confined to a 3 dimensional manifold
within the 6-dimensional fundamental form space. As this surface is
parameterized by one nonlinear number ($\theta$) and two small linear numbers
($t,b$), we refer to such as the hyperribbon.

To evaluate the adherence of real origami sheets to their shape hyperribbon,
as well as the subsequent analyses herein, we numerically investigate force-
balanced deformations using the MERLIN2 simulation framework for an “enhanced”
bar-hinge model of the origami sheet Filipov  _et al._ (2017); Liu and Paulino
(2017); liu2018highly. To facilitate algorithmic convergence, the data
displayed in Figs. 3,4&5 are obtained under arbitrary point-loading conditions
and using “rigid” parameters in which the crease folding incurs no energy
penalty ($k_{f}\rightarrow 0$). In contrast, to better represent real world
response the data illustrated in Figs. 1&2 are obtained using “realistic”
parameters emulating a constitutive sheet having thickness $1\%$ of the panel
dimension (i.e. roughly capturing printer paper), and thereby having a small,
but finite crease folding modulus $k_{f}$. All of these numerical results, as
well as the applied loading conditions, and simulation protocol are described
in Appendix. C.

As shown in Fig. 3, the numerically force-balanced _nonlinear_ deformations of
a variety of these origami sheets exhibit very close adherence to this
hyperribbon ($<1\%$ error) out to strains of $10-20\%$. As suggested in the
energetic argument above, the adherence to the hyperribbon ubiquitously
improves as the continuum limit is approached in Fig. 4. This non-hyperribbon
error follows a suggestive power-law scaling, in which the curvature error
scales inversely with the system size, while the metric error improves faster,
with the scaling $\sim N^{-2}$.

##  IV Field theory for origami sheets from compatibility

Figure 4:  The averaged fractional deviation from the hyperribbon is very
small and ubiquitously improves (becomes smaller) as the system size is
increased and the continuum limit is approached. The error in the curvature
(dashed lines) scales approximately with the lattice coarseness ($1/N$), while
the error in the metric (solid lines) scales approximately with the square
coarseness ($1/N^{2}$).

The nonuniform soft response of the origami sheet will be composed from these
three spatially varying modes
$\theta(X_{1},X_{2}),t(X_{1},X_{2}),b(X_{1},X_{2})$. However, it is well-known
in the study of differential geometry that not all spatial patterns of the
metric and curvature tensors correspond to any realizable (i.e. immersible)
surface: regardless of deformation mechanics and energy, geometric
compatibility restricts the spatial evolution of these tensors across the
surface. Therefore, as the soft modes
$\theta(X_{1},X_{2}),t(X_{1},X_{2}),b(X_{1},X_{2})$ control the metric and
curvature tensors in our origami sheet, they are in turn governed by geometric
compatibility conditions. Again relying on the rigidity theorem of Bonnet, it
is known that the Gauss-Codazzi-Mainardi-Peterson equations are both necessary
and sufficient to guarantee such a geometrically valid sheet Toponogov and
Rovenski (2006). These equations, discussed in detail in Appendix B.3
constitute a set of three PDEs, and with Eqs. 1&3 the general form of these
PDEs may be reduced to a form specific to this category of origami sheet. In
full form, these equations are:

$\displaystyle\frac{1}{\lambda_{1}\lambda_{2}}$
$\displaystyle\epsilon^{ij}\partial_{i}\left(B_{jk}(\theta)\partial_{k}\theta\right)$
$\displaystyle=-\frac{1}{\lambda_{1}^{2}\lambda_{2}^{2}}\left[b^{2}\lambda_{1}^{\prime}\lambda_{2}^{\prime}+t^{2}\right]$
(4) $\displaystyle 0=$
$\displaystyle[\partial_{1}t-\lambda_{1}^{\prime}\partial_{2}b-b\lambda_{1}^{\prime\prime}\partial_{2}\theta]+t(\frac{\lambda_{2}^{\prime}}{\lambda_{2}}-\frac{\lambda_{1}^{\prime}}{\lambda_{1}})\partial_{1}\theta$
$\displaystyle+b\left(\frac{\lambda_{1}^{\prime}}{\lambda_{1}^{2}}-\frac{\lambda_{2}^{\prime}}{\lambda_{2}^{2}}\right)*\lambda_{1}\lambda_{1}^{\prime}\partial_{2}\theta$
(5) $\displaystyle 0=$
$\displaystyle-[\partial_{2}t+\lambda_{2}^{\prime}\partial_{1}b+b\lambda_{2}^{\prime\prime}\partial_{1}\theta]+t(\frac{\lambda_{2}^{\prime}}{\lambda_{2}}-\frac{\lambda_{1}^{\prime}}{\lambda_{1}})\partial_{2}\theta$
$\displaystyle-b\left(\frac{\lambda_{1}^{\prime}}{\lambda_{1}^{2}}-\frac{\lambda_{2}^{\prime}}{\lambda_{2}^{2}}\right)*\lambda_{2}\lambda_{2}^{\prime}\partial_{1}\theta\,,$
(6)

where

$\mathrm{B}=\begin{bmatrix}0&-\frac{\lambda_{1}^{\prime}}{\lambda_{2}}\\\
\frac{\lambda_{2}^{\prime}}{\lambda_{1}}&0\end{bmatrix}\,.$ (7)

The first (Eq. 4) of these arises from the Gauss equation and is second order,
while Eqs. 5&6 are the first order Codazzi equations. While this set of
nonlinear PDEs is highly nontrivial to solve, the balance of the number of
fields and equations suggests that boundary conditions will be sufficient to
determine the soft mode patterning through the bulk of an origami sheet.
Despite the natural exceptions to this, it therefore appears that the number
of soft nonuniform deformations available to the origami sheet globally will
be infinite but subextensive, similar to the planar unimode case Czajkowski
_et al._ (2022); Zheng  _et al._ (2022); Czajkowski and Rocklin (2022). Note
that in the absence of the twist and bend fields, this reduces as expected to
the planar unimode compatibility theory established in Refs. Czajkowski  _et
al._ (2022); Zheng  _et al._ (2022); Czajkowski and Rocklin (2022).

###  IV.1 Dual spaces of sheared analytic functions capture linear
deformations in the origami sheet

Figure 5:  Under linear deformations, the spatial distribution of soft mode
amplitudes may be described using sheared analytic functions. (a) The
nontrivial spatial patterns of the three modes are reproduced to high
qualitative accuracy by closest fit sheared analytic function. (b) The
accuracy of these fittings is a proxy for the degree of success of the soft
mode theory. As expected, the error in fitting the metric (solid lines) and
the curvature (dashed lines) becomes small and the theory becomes accurate
only in the limit that the ratio of the panel stretching stiffness to the
panel bending stiffness becomes very large ($\sim 10^{5}$).

The theory in Eqs. 4-6 governs the nonuniform soft modes, yet analytic
solutions are generally challenging to obtain. We therefore consider the limit
of linearly small deviatoric mode amplitudes $\delta\theta,\delta t,\delta b$,
where these equations may be written in the compact form

$\displaystyle 0$ $\displaystyle=\partial_{w}\partial_{\bar{w}}\delta\theta$
(8) $\displaystyle 0$ $\displaystyle=\partial_{\bar{w}}H$ (9) $\displaystyle
0$ $\displaystyle=\partial_{w}\bar{H}$ (10)

where

$\displaystyle w$ $\displaystyle=X_{1}+\sqrt{\nu_{\text{in}}}X_{2}$ (11)
$\displaystyle\bar{w}$ $\displaystyle=X_{1}-\sqrt{\nu_{\text{in}}}X_{2}$ (12)

are a transformed “sheared” set of spatial variables established in Ref.
Czajkowski and Rocklin (2022) and

$\displaystyle H$
$\displaystyle=t-\sqrt{\frac{1}{\nu_{\text{in}}}}\lambda_{2}^{\prime}b$ (13)
$\displaystyle\bar{H}$
$\displaystyle=t+\sqrt{\frac{1}{\nu_{\text{in}}}}\lambda_{2}^{\prime}b$ (14)

describes the curvature. These new variables $H,\bar{H}$ determine the twist
and bend fields which, in the hyperribbon description of the origami sheet,
thereby determine the entire curvature tensor according to Eq. 3. Although
energetic couplings between the metric and the curvature tensors are not ruled
out, there is no coupling from geometric compatiblity in this linear limit,
and the metric and curvature may be determined independently. Therefore, we
may think of these deformations as composed of two steps: first a quasi planar
soft deformation is applied, followed by a “non-planar” deformation, which
does not contribute to the planar strains.

Consider first these non-planar deformations. In the limit of pure dilational
planar folding mode ($\nu_{\text{in}}=-1$), the sheared coordinate
transformation is the conventional complex-variable representation of the
plane $w,\bar{w}\rightarrow z,z^{*}$, and according to Eqs. 9&10, the modes
will be controlled by a conventional holomorphic $\partial_{\bar{z}}H=0$ (i.e.
conformal, analytic) and anti-holomorphic $\partial_{z}\bar{H}=0$ mapping.
More generally (i.e. for arbitrary $\nu_{\text{in}}$) the solutions for the
soft curvatures may be written down using “sheared analytic functions” as
introduced in Ref. Czajkowski and Rocklin (2022). Similar to the dilational
case, such functions are determined solely by one or the other of the sheared
analytic variables $w,\bar{w}$ but not both, and in the present case we have
$H(w,\bar{w})\rightarrow H(w)$ and
$\bar{H}(w,\bar{w})\rightarrow\bar{H}(\bar{w})$. This formalism, which
determines the spatial distribution of the soft modes, is described in greater
detail in Appendix B.1. These spatial patterns are controlled by an
exceptional point at $\nu_{\text{in}}=0$, where the sheared spatial variables
transition from being real-valued when $\nu_{\text{in}}>0$ (anauxetic folding
mechanism) to being complex-valued when $\nu_{\text{in}}<0$ (auxetic folding
mechanism). Tuning the system across this transition, the bulk oscillatory
modes of the anauxetic systems transform into evanescent modes which decay
away from the boundary in the auxetic systems.

The remaining question of planar folding mode compatibility may be again
addressed directly relying on previous work. By virtue of decoupling from non-
planar effects, the planar mode compatibility equation in Eq. 8 must be
equivalent to the linear analysis of planar unimode metamaterials addressed in
Ref. Czajkowski and Rocklin (2022). As a result of this, the second Eq. 8 may
be exactly rewritten as two first order PDE equations

$\displaystyle 0$ $\displaystyle=\partial_{\bar{w}}\varepsilon$ (15)
$\displaystyle 0$ $\displaystyle=\partial_{w}\bar{\varepsilon}$ (16)

where we will define $\phi$ as the local coarse rotation of each unit cell in
the sheet reference plane and
$\varepsilon\equiv\delta\theta{\lambda_{1}^{\prime}}-\frac{1}{\sqrt{\nu_{\text{in}}}}\phi$
and
$\bar{\varepsilon}\equiv\delta\theta{\lambda_{1}^{\prime}}+\frac{1}{\sqrt{\nu_{\text{in}}}}\phi$
are special defined variables capturing the planar components of deformation.
Eqs. 15&16 are clearly identical equations to those controlling the non-planar
modes (Eqs. 9&10), simply replacing the variables
$H,\bar{H}\rightarrow\varepsilon,\bar{\varepsilon}$. Therefore, the same
sheared analytic functions used above to generate a pattern of the non-planar
modes $t,b$, can be instead used to generate a pattern of the planar modes
$\delta\theta,\phi$. Each non-planar deformation then corresponds precisely
with one of the planar deformations, and vice versa, thereby establishing a
mathematical duality, similar to that revealed in Ref. Czajkowski and Rocklin
(2022). Any arbitrary linear deformation of the origami sheet is now generated
by choosing one set of sheared analytic functions controlling the planar modes
and a second set of sheared analytic functions controlling the non-planar
modes, and superimposing these deformations in the origami sheet. It is
important to note that in a given deformation the nonplanar mode need not be
the corresponding dual pattern to the superimposed planar mode and vice versa:
this represents a special case of a “self-dual” deformation, the physical
implications of which will require further investigation.

In Fig 5 the distribution of the fundamental forms along the sheet is
described to high accuracy (fractional error only reaching as high as $5\%$)
using these sheared analytic modes when the sheet bending $k_{p},k_{f}$
becomes much more energetically costly than sheet stretching $k_{s}$. This
confirms the notion that mechanism-based origami composed of sufficiently thin
sheets will have response governed by the orisometries presented herein.

##  V Discussion

In contrast to related previous investigations into nonuniform patterns of
folding in lattice origami sheets Nassar  _et al._ (2017, 2022); Xu  _et al._
(2023), the formalism presented here in terms of these spatially varying mode
amplitudes, is unique. Our approach to the compatibility requirements which
constrain these modes brings new qualitative insight to the breadth and
spatial character of response.

This orisometries formalism reveals also that the soft deformations correspond
in the linear limit to the sheared analytic modes introduced previously
Czajkowski and Rocklin (2022). These are governed by an exceptional point in
the Poisson’s ratio at $\nu_{\text{in}}=0$, which qualitatively divides the
character of the spatial response. Further, the in-plane and out-of-plane
modes which superimpose to span the space of soft response are mathematically
dual spaces, suggesting connections to other mechanical dualities Pretko and
Radzihovsky (2018); Fruchart  _et al._ (2020); Czajkowski and Rocklin (2022).
While this gives great insight into the character of the soft modes available
to these origami sheets, it also adds further emphasis to the question: what
is the deeper principle that these dualities in mechanism-based metamaterials
arise from?

While the analysis herein focuses on a particular mode-based formalism which
serves our purposes, many long-established alternatives exist for the
description of manifold deformation in three dimensions. Our analysis
nonetheless indicates that any such formalism, in describing the low energy-
deformations of such a sheet, must still exhibit this duality therein. This
new duality compounds with the duality identified previously between the
spatial patterns of soft strains and the stresses in planar unimode systems
Czajkowski and Rocklin (2022); it is expected that they will coexist in the
origami sheet, yet an in-depth analysis of the constitutive laws that may
control both the sheet strains and curvatures should be established first. The
formalism of spatially varying soft mode variables introduced here will lend a
useful framework for development of an energy functional to predict specific
deformation response.

Interestingly, according to constraint counting, the presence of three local
deformation modes is not unique to the patterns explored here. All fold
patterns may be considered to be effectively triangulated with soft panel
bending in this manner and, starting from a quasi-planar state, three modes
must always be present and governed by the equations of compatibility
McInerney  _et al._ (2020). Therefore, while the Morph family of fold patterns
explored here are chosen for analytic convenience, a subextensive but infinite
space of soft modes is expected to be the norm rather than the exception for
periodic fold patterns. However, in such origami sheets which approximate a
cone or a cylinder in their resting state, the mathematics of soft deformation
must change. In particular, such screw-periodic geometries exhibit two linear
isometries rather than the three for the spatially-periodic sheets considered
herein. It will be interesting to explore the possible soft modes and the
challenges that may arise in discretized descriptions of e.g. curvature in
this setting.

The orisometry approach presented here also presents new challenges and
questions. For instance, we pose the question of what governs the suggestive
and distinct continuum limit scaling exponents in the hyperribbon fitting
error shown in Fig. 4. Beyond this, our theory presents a series of nonlinear
equations, and the solutions to these equations are largely unexplored and
potentially rich with possibilities. Finally, in analogy with the phenomena of
stress focusing and developable cones in thin sheet response, these is a
question of where and when this theory here may break down, requiring more
singular and localized phenomena in place of smoothness.

##  Appendix A Morph-class of fold patterns

In this work we focus on a class of fold patterns dubbed the “Morph” Pratapa
_et al._ (2019). Each pattern in this class is composed of a repeating pattern
of four parallelogram panels (in two identical pairs) as shown in Fig. 6. To
generate an arbitrary member of this class of fold patterns, two “sector
angles” $\alpha$ and $\beta$, and three edge lengths $a_{1}$, $a_{2}$ and $c$
must be provided, as shown in Fig. 6a. For rigid origami (i.e. panels can
neither bend nor stretch) this information is only sufficient to determine the
crease pattern of the origami sheet, but not to determine the stress free
reference state. This is because the origami sheet has a single _mechanism_
degree of freedom, which extends nonlinearly, as shown in main text Fig. 1.
Such a mechanism, in the absence of any energy cost of folding along creases,
provides a degenerate space of ground states which the sheet may realize.

To understand this rigid mechanism folding motion, consider a single isolated
vertex at the intersection of the four parallelogram panels (again, two pairs
of identical parallelograms). We may consider the four dihedral angles of the
creases between the panels to form a sufficient set of nontrivial degrees of
freedom describing the vertex shape. These dihedral angles quantify folding,
and for a particular folding configuration around a single vertex to be
geometrically valid (i.e. to permit flat, unstretched panels between them),
these angles must satisfy the well-known Belcastro-Hull compatibility
conditions Belcastro and Hull (2002). This compatibility requirement
generically imposes three distinct nontrivial constraints on the dihedral
angles connected to a single vertex. For the Morph folding patterns, four
creases meet at the vertex and, accounting for the three constraints, one
nonlinear degree of freedom is left over to vary the vertex configuration.
This captures the rigid folding motion available to the sheet. We parameterize
this folding motion using a single scalar parameter $\theta$, as shown in Fig.
6b. This vertex, along with the four panels it connects, defines a unit cell
of the repeating fold pattern. It is, however, neither generic nor trivial
that repeating copies of such a unit cell can be assembled together to form a
single coherent origami sheet. The key restriction, as shown in Ref.
Waitukaitis  _et al._ (2015); McInerney  _et al._ (2022), is that each panel
must be chosen to be parallelogram shaped. This is, by definition, true of the
family of Morph patterns considered here. With a bit of thought it is clear
that the edges on opposite sides of the unit cell must match in both length
and orientation, and the repeating crease pattern is composed of such a unit
cell, translated by the lattice vectors and stitched together along crease
lines. And, as this is possible for arbitrary choices of the folding angle
$\theta$ of the single unit cell, this rigid folding motion is preserved in
the infinite sheet.

Given the defining parameters $a_{1},a_{2},c,\alpha,\beta$ and the rigid
folding motion parameter $\theta$, the procedure for generating a sheet
configuration in 3d space is straightforward but analytically cumbersome
McInerney  _et al._ (2022). Conveniently, in the case of the Morph patterns,
symmetry simplifies this analysis. For simplicity, we restrict our analysis to
the case $a_{1}=a_{2}\equiv a$. Consider the four vectors (orange
$\mathbf{a}_{1},\mathbf{a}_{2}$ and pink $\mathbf{c}_{1},\mathbf{c}_{2}$)
emanating from a vertex, as shown in Fig. 6b. These edge vectors may be used
to assemble the lattice vectors
$\mathbf{l}^{(1)}=\mathbf{c}_{1}-\mathbf{c}_{2}$ and
$\mathbf{l}^{(2)}=\mathbf{a}_{1}-\mathbf{a}_{2}$ connecting identical vertices
in neighboring unit cells. Because the Morph patterns are symmetric on
mirroring through the $\mathbf{l}^{(1)}$ direction, these lattice vectors are
orthogonal for all choices of the folding angle. These vectors may then be
used to establish a local orthonormal basis
$\hat{l}^{(1)},\hat{l}^{(2)},\hat{\mathbf{N}}=\hat{l}^{(1)}\times\hat{l}^{(2)}$.
The edge vectors may then be written in terms of their components along these
basis elements via

$\displaystyle\mathbf{a}_{1}$
$\displaystyle=a\left(\cos(\theta)\hat{l}^{(2)}-\sin(\theta)\hat{\mathbf{N}}\right)$
(17) $\displaystyle\mathbf{a}_{2}$
$\displaystyle=-a\left(\cos(\theta)\hat{l}^{(2)}+\sin(\theta)\hat{\mathbf{N}}\right)$
(18) $\displaystyle\mathbf{c}_{1}$
$\displaystyle=c\left(\sin(\psi)\cos(\phi)\hat{l}^{(1)}+\sin(\psi)\sin(\phi)\hat{\mathbf{l}^{(2)}}+\cos(\psi)\hat{\mathbf{N}}\right)$
(19) $\displaystyle\mathbf{c}_{2}$
$\displaystyle=c\left(-\sin(\psi)\cos(\phi)\hat{l}^{(1)}+\sin(\psi)\sin(\phi)\hat{\mathbf{l}^{(2)}}+\cos(\psi)\hat{\mathbf{N}}\right)\,,$
(20)

where we have implicitly defined a polar $\psi$ and an azimuthal $\phi$ angle
to describe the $c$-vector orientations in spherical coordinates. Here, we
have used the symmetry of the fold pattern to relate $\mathbf{c}_{1}$ to
$\mathbf{c}_{2}$ as the mirror across the plane with normal $\hat{l}^{(1)}$.
The goal is to determine the changes to the edge orientations which will
preserve the sector angles $\alpha,\beta$ between adjacent edges. Again this
is simplified into two constraining conditions by the symmetry across the
$\hat{l}^{(1)}$ direction. We require

$\displaystyle\mathbf{a}_{1}\cdot\mathbf{c}_{1}$
$\displaystyle=ac\cos(\alpha)$
$\displaystyle\mathbf{a}_{2}\cdot\mathbf{c}_{1}$
$\displaystyle=ac\cos(\beta)\,,$

and, plugging in, these conditions becomes

$\displaystyle ac\cos(\alpha)$
$\displaystyle=ac\left(\cos(\theta)\sin(\psi)\sin(\phi)-\sin(\theta)\cos(\psi)\right)$
$\displaystyle ac\cos(\beta)$
$\displaystyle=-ac\left(\cos(\theta)\sin(\psi)\sin(\phi)+\sin(\theta)\cos(\psi)\right)\,.$

These are solved by the relations

$\displaystyle\cos(\psi(\theta,\alpha,\beta))$
$\displaystyle=-\frac{\cos(\alpha)+\cos(\beta)}{2\sin(\theta)}$ (21)
$\displaystyle\sin(\phi(\theta,\alpha,\beta))$
$\displaystyle=-\frac{\cos(\alpha)-\cos(\beta)}{2\cos(\theta)\sin(\psi(\theta,\alpha,\beta))}\,.$
(22)

Here, the geometric angle $\theta$ is the chosen parameter of the rigid
folding motion. Therefore, we view $\psi(\theta,\alpha,\beta)$ as a function
of this angle and the sector angles that define the vertex, whose value is
obtained by solving Eq. 21. Similarly, $\phi(\theta,\alpha,\beta)$ may be
viewed as a function determined by the vertex geometry and folding, with a
value that can be determined by using $\psi$ to solving Eq. 22.

Importantly, these geometric relations allow us to write the lattice vectors
analytically in terms of the fold angle and sector angles

$\displaystyle\mathbf{l}^{(2)}$ $\displaystyle=2a\cos(\theta)\hat{l}^{(2)}$
(23) $\displaystyle\mathbf{l}^{(1)}$
$\displaystyle=2c\sin(\psi(\theta,\alpha,\beta))\cos(\phi(\theta,\alpha,\beta))\hat{l}^{(1)}$
(24)
$\displaystyle=2c\sqrt{1-\left(\frac{\cos(\alpha)+\cos(\beta)}{2\sin(\theta)}\right)^{2}}\times$
$\displaystyle\,\,\sqrt{1-\left(\frac{\cos(\alpha)-\cos(\beta)}{2\cos(\theta)\sqrt{1-\left[\frac{\cos(\alpha)+\cos(\beta)}{2\sin(\theta)}\right]^{2}}}\right)^{2}}\hat{l}^{(1)}$

With this information, all of the vertices of the unit cell may be defined by
transporting these vectors along one another, and the vertices of the infinite
sheet may be defined from the unit cell by translation along the lattice
vectors.

Figure 6:  Construction and simulation of Morph patterns in a bar-hinge model.
(a) The Morph unit cell is composed of four parallelogram panels, which come
in indistinguishable pairs. Parameters to generate these panels, as well as
the recipe to connect them into a unit cell illustrated here. (b) Once
constructed, such a Morph pattern sheet is still able to change shape via the
folding motion, parameterized by the folding $\theta$ shown here for a single
unit cell. (c) A representative model of such a sheet is achieved by adding
diagonal crease lines across the parallelogram panels and adding energy
penalties for changing the dihedral angles across these panel diagonals
$\rho_{p}$, changing the dihedral angles across the creases at the panel edges
$\rho_{q}$ as well as spring energy penalties for stretching the crease lines
along their axis $\mathbf{r}_{u}v$ stretching from node $u$ to $v$.

##  Appendix B Coarse formalism for origami deformations

Here we present a formalism to describe the coarse motion of an origami sheet,
patterned with a repeating microstructure of folds. The Morph patterns
constitute a representative set of examples, and the methods may be readily
generalized to other fold patterns, such as those explored in Ref. McInerney
_et al._ (2022).

To describe the coarse macroscopic shape changes of an origami sheet, we
obscure our view of the fine details of deformation within each unit cell, and
instead focus on how each unit cell is positioned in space relative to its
neighbors. This is accomplished using the smoothed midplane of the sheet,
which intersects with the geometric center of each unit cell. In the case of
the Morph, the unit cell geometry is represented by four vertices contained
therein (see color-coded dots in main text Fig. 1a,b). In this work we choose
the unit cell center as the geometric average of the four vertex positions.

The smooth sheet which the unit cell centers approximate, taken as a
differentiable manifold, is governed by familiar quantities from differential
geometry. For instance, it is well-known from the theorem of Bonnet Toponogov
and Rovenski (2006) that knowledge of the first and second fundamental forms
is completely sufficient to describe the shape of this sheet, up to overall
rotations and translations. These fundamental forms are defined from a
parameterization of the smooth surface using two scalar parameters
$(X_{1},X_{2})$. Within this picture, the position of a point in the sheet is
represented by a vector function $\mathbf{R}(X_{1},X_{2})$. The
parameterization here is chosen to be the distance along each lattice
direction in the undeformed lattice, such that the first fundamental form
captures the metric of deformation, and may be used to quantify the strain
induced in the sheet relative to the quasi-planar reference state. The tensor
components of the first fundamental form are defined from the tangent vectors
$\partial_{i}\mathbf{R}$, where $\partial_{i}$ denotes the usual partial
derivative with respect to the $i$th parameterization coordinate. Given the
choice of parameterization, the tangent vectors then take the form
$\partial_{i}\mathbf{R}=\mathbf{l}^{(i)}/|\mathbf{l}^{(i)}_{0}|$, where
$\mathbf{l}^{(i)}$ is the lattice vector connecting the unit cell position to
its neighbor in the $i$th direction and $\mathbf{l}^{(i)}_{0}$ is the same
lattice vector prior to deformation, in the quasi-planar reference state. The
tensor components of the first fundamental form are then written

$g_{ij}\equiv(\partial_{i}\mathbf{R})\cdot(\partial_{j}\mathbf{R})=\frac{\mathbf{l}^{(i)}\cdot\mathbf{l}^{(j)}}{|\mathbf{l}^{(i)}_{0}||\mathbf{l}^{(j)}_{0}|}\,.$
(25)

This tensor quantity may be used to define the Green-Lagrange nonlinear strain
of the surface in its local plane via

$\varepsilon_{ij}=\frac{1}{2}(g_{ij}-g_{ij}^{(0)})$ (26)

in which $g_{ij}^{(0)}$ is the reference metric of the material, typically
represented by the identity. Note that this expression for the metric relies
on the assumption that these lattice vectors are much smaller than all other
lengthscales of spatial variation of the midplane, and therefore that they
approximate the sheet infinitesimals.

To capture the second fundamental form, we use these midplane surface tangent
vectors to define the surface normal via:
$\hat{\mathbf{N}}\equiv\frac{\mathbf{l}^{(1)}\times\mathbf{l}^{(2)}}{|\mathbf{l}^{(1)}\times\mathbf{l}^{(2)}|}$.
From this, the components of the second fundamental form are obtained with the
formula $h_{ij}=-\partial_{i}\mathbf{R}\cdot\partial_{j}\hat{\mathbf{N}}$.
Note that due to the choice of sheet parameterization with units of length,
that components of the second fundamental form have units of inverse length.
In particular, the contraction with the inverse metric $g^{ik}h_{kj}$
evaluated in a diagonal frame captures the _principal curvatures_ , which are
the eigenvalues, and which reflect the inverse of the radii of curvature of
the midplane of the sheet. Note that these curvatures vanish for the
spatially-periodic rigid ground states of the sheets, as the sheet normal is
then uniform.

###  B.1 Triangulation and linear modes

As discussed in Sec. A, a fold pattern in the Morph family has only a single
rigid folding motion. However, each constitutive parallelogram panel is itself
a thin elastic sheet. As explored in greater detail in Sec. C.2, thin sheets
are qualitatively much softer to bending deformations than to stretching. In
general cases of soft deformations of a thin sheet one might typically
consider the highly dimensional space of isometries, in addition to the
possibilities of singular structures such as d-cones and ridges. However, each
of the parallelogram panels of an origami pattern such as the Morph are
constrained to match to neighboring panels along straight edges, which greatly
reduces the isometries available. As has become a standard protocol in origami
modeling and analysis Schenk and Guest (2013); Wei  _et al._ (2013); Filipov
_et al._ (2017); Liu and Paulino (2017); liu2018highly; McInerney  _et al._
(2020) the softest deformations of the panels maybe approximated by the
addition of a crease across the diagonal of the quadrilateral face. With the
addition of such “panel bending” crease lines, we describe the origami sheet
as a triangulated mesh. This enables the use of “bar-hinge” type models, in
which the complex deformation of an origami sheet to be captured in a reduced
description where the degrees of freedom are the positions of the vertices.

This “bar-hinge” approximation has been very helpful in identifying the
softest available motions for the large-scale origami sheet. Just as the very
small comparative energy penalty of the folding at the major creases allowed
us to identify the rigid folding motion, we may further choose to neglect the
small energetic penalties that may arise from the panel bending. In this
idealized case, the origami sheet resembles a truss, or a ball-spring network,
with vertices connected by springs as in Fig. 6c, and without any torsional
energy penalties coming from triangular panel orientations. In this limit, one
may use constraint counting and symmetry arguments to show that three
infinitesimal zero energy modes of deformation must exist in the infinite
system McInerney  _et al._ (2020). One of these modes is the quasi-planar
rigid folding motion established above, which extends nonlinearly. The other
two only achieve linearly small values, and are known colloquially as the
“twist” and “bend” modes, according to their qualitative effect on the sheet
midplane. We here review and establish the nature of these three modes via
their effect on midplane shape, as captured by the fundamental forms.

The first mode, referred to here as the “planar” mode, is parameterized by the
folding angle $\theta$, which may extend to nonlinear values on the order of
$\pi$. At the coarse scale, the effect of this folding motion is captured in
the first fundamental form via

$g_{ij}\rightarrow
g^{*}_{ij}(\theta)\equiv\frac{\mathbf{l}^{(i)}(\theta)\cdot\mathbf{l}^{(j)}(\theta)}{|\mathbf{l}^{(i)}(\theta_{0})||\mathbf{l}^{(j)}(\theta_{0})|}\,,$
(27)

where $\theta_{0}$ is the folding angle in the chosen reference state. This
equation is simply the special case of Eq. 25 in which the lattice vectors,
and thereby the metric, are now determined by this folding angle. In the case
of the Morph, the lattice vectors are orthogonal for all values of the
folding, and this general formula therefore simplifies to

$g^{*}_{ij}(\theta)\rightarrow\begin{bmatrix}\lambda_{1}(\theta)^{2}&0\\\
0&\lambda_{2}(\theta)^{2}\end{bmatrix}\,,$ (28)

as shown in main text Eq. 1, where

$\displaystyle\lambda_{1}(\theta)$
$\displaystyle=\frac{|\mathbf{l}^{(1)}(\theta)|}{|\mathbf{l}^{(1)}(\theta_{0})|}$
(29) $\displaystyle\lambda_{2}(\theta)$
$\displaystyle=\frac{|\mathbf{l}^{(2)}(\theta)|}{|\mathbf{l}^{(2)}(\theta_{0})|}$
(30)

are the principal stretches induced by the planar mode and, for the Morph,
$\mathbf{l}^{(1)},\mathbf{l}^{(2)}$ are the functions defined above in Eqs.
24&23.

To understand the second and third modes, referred to here (as well as in
other works) as the “twist” and “bend” modes, we must consider the second
fundamental form $h_{ij}$. As noted above, this tensor is symmetric, and
therefore determined by three scalars, $h_{11}$, $h_{12}$ and $h_{22}$. As
noted in previous work, the twist and bend modes may be fruitfully described
in terms of dihedral angles of the triangulated sheet, or, in terms of
“vertex” and “face” amplitudes McInerney  _et al._ (2022). However, in either
description each mode may only achieve linearly small values before they begin
to cause sheet stretching. Further, as the only nonlinear mode available to
the origami sheet $\theta$, does not contribute to this curvature tensor, the
curvature must also only extend to linear order in the detailed degrees of
freedom. Therefore, the components of this tensor may be written in in a basis
of three amplitudes

$h_{ij}\rightarrow
h^{*}_{ij}(\theta,t,b)=t\hat{t}_{ij}(\theta)+b\hat{b}_{ij}(\theta)+\kappa_{3}\hat{h}^{(3)}_{ij}(\theta)$
(31)

chosen so that the first two terms capture the twist and bend fields, while
the third term captures the component of curvature that does not arise from
bend or twist and therefore we will have $\kappa_{3}=0$ in our hyperribbon
theory. Our task is therefore to determine the nature of the curvature
subspace we are restricted to, defined by the basis $\hat{t}(\theta)$ and
$\hat{b}(\theta)$.

To determine this subspace, we first allocate the off diagonal component of
the curvature to the twist mode by setting

$\hat{t}(\theta)=\begin{bmatrix}0&1\\\ 1&0\end{bmatrix}\,.$ (32)

This indicates that the twist mode will generate no curvature along the
lattice vector directions, in accordance with prior twist mode descriptions
Schenk and Guest (2013); McInerney  _et al._ (2022). Note that this definition
removes any dependence of $\hat{t}$ on the folding angle $\theta$.

Within this definition of the twist, the generic form for the bend mode
curvature must take the form

$\hat{b}=\begin{bmatrix}b_{1}(\theta)&0\\\ 0&b_{2}(\theta)\end{bmatrix}\,.$
(33)

Rather than requiring additional computation of the unit cell structure, the
functions $b_{1}$ and $b_{2}$ may be related to the principal stretches
$\lambda_{1}$ and $\lambda_{2}$. Incrementation of the bend mode generates
(and defines) the out-of-plane Poisson’s ratio

$\nu_{\text{out}}\equiv-\frac{\delta h_{22}}{\delta
h_{11}}=-\frac{b_{2}(\theta)}{b_{1}(\theta)}$ (34)

while incrementation of the planar mode defines the in-plane Poisson’s ratio

$\nu_{\text{in}}\equiv-\frac{\delta g_{22}}{\delta g_{11}}\,.$ (35)

Knowing the form of the metric, and its increment, we may write this in terms
of the principal stretches

$\nu_{\text{in}}(\theta)=-\frac{\lambda_{2}(\theta)\lambda_{2}^{\prime}(\theta)}{\lambda_{1}(\theta)\lambda_{1}^{\prime}(\theta)}$
(36)

where prime denotes derivative with respect to the argument $\theta$. It has
been well-established that folding patterns in this category have
$\nu_{\text{in}}=-\nu_{\text{out}}$ Schenk and Guest (2013); Wei  _et al._
(2013); Pratapa  _et al._ (2019); McInerney  _et al._ (2022); Nassar  _et al._
(2022) and we may use this relation to write

$\displaystyle b_{1}(\theta)=$
$\displaystyle\lambda_{1}(\theta)\lambda_{1}^{\prime}(\theta)$ (37)
$\displaystyle b_{2}(\theta)=$
$\displaystyle-\lambda_{2}(\theta)\lambda_{2}^{\prime}(\theta)$ (38)

without any loss of generality. We are free to choose to put the negative sign
on the $b_{2}$ term as this choice may be absorbed into our definition of the
bend mode amplitude itself, along with any other fold-dependent multiplicative
factor applied to both $b_{1}$ and $b_{2}$ equally. Again, due to the symmetry
of the origami pattern, this bend mode does not contribute to the first
fundamental form to first order. Despite matching qualitatively with existing
literature, these particular amplitude definitions for the twist and bend
modes are unique to this work, having been chosen for analytic convenience in
our coarse description.

###  B.2 Spatially varying planar, twist, and bend modes

Following the central statement of the main text regarding spatially varying
generalizations of these local low-energy modes, we promote these scalar mode
parameters to fields
$\\{\theta,b,t\\}\rightarrow\\{\theta(X_{1},X_{2}),b(X_{1},X_{2}),t(X_{1},X_{2})\\}$
on the reference space (i.e. surface parameterization). One may think of these
fields, and their gradients, as constituting a recipe which determines the
change in crease folding and panel bending in each unit cell with respect to
the reference state. It therefore follows that the field modes constitute a
recipe to determine the fundamental forms across the sheet. In principal, all
of the field values as well as their gradients may contribute to determine the
fundamental forms. We are interested in deformations where the gradients are
small, as are the dimensionless twist and bend modes themselves
$a\,b(X_{1},X_{2}),a\,t(X_{1},X_{2})$. To lowest order approximation in these
small quantities, we are therefore concerned with the influence of the
nonlinear planar mode $\theta(X_{1},X_{2})$ as well as the influence of the
linear effects of the planar mode gradients $\nabla_{X_{1},X_{2}}\theta$ and
the twist and bend. However, to the symmetry of the Morph folding pattern upon
mirroring in the $X_{1}$ and $X_{2}$ directions, the gradients of the planar
mode may not contribute to the first fundamental form. We are therefore left
with the simple generalization of the above relations for the effect of the
uniform modes on the fundamental forms

$\displaystyle g_{ij}(X_{1},X_{2})$ $\displaystyle\rightarrow
g^{*}_{ij}(\theta(X_{1},X_{2}))$ (39) $\displaystyle h_{ij}(X_{1},X_{2})$
$\displaystyle\rightarrow
h^{*}_{ij}(\theta(X_{1},X_{2}),t(X_{1},X_{2}),b(X_{1},X_{2}))\,.$ (40)

Corrections to these expressions may arise at higher order, proportional to
second gradients of the folding angle $\theta$, as well as first gradients of
the twist and bend.

With Eqs 39&40, we have determined that the possible 6 local shape degrees of
freedom $g_{11},g_{12},g_{22},h_{11},h_{12},h_{22}$ of the origami midplane
will be determined entirely by three parameters (the local soft modes), and
therefore must be restricted to a 3 dimensional subspace of this possible 6
dimensional space. We may characterize this 3 dimensional subspace by first
considering the nonlinear mode $\theta(X_{1},X_{2})$. In the absence of the
other two mode-fields, the origami sheet must lie on a 1 dimensional nonlinear
manifold contained within the 3 dimensional shape space. Around each of these
points, the possible inclusion of the twist and bend fleshes out this one
dimensional manifold to become locally 3 dimensional. This structure in the 6
dimensional shape space which is locally long in one direction, short in
another (2) directions, and flat in the additional directions is reminiscent
of a three dimensional ribbon structure, and we therefore refer to this as the
hyperribbon.

###  B.3 Compatibility and the Gauss-Codazzi equations

As identified above, the spatially varying planar, bend, and twist mode fields
are expected to capture the response of the origami sheet. One might imagine
that any spatially varying pattern of these modes will correspond to a
candidate soft response. However, as also mentioned above, we may think of
these modes as a recipe determining the changes in the shape of the sheet
midplane. It is well established that in the continuum picture of such a
sheet, the strains and curvatures of the midplane must obey compatibility
relations, or else when we attempt to follow the recipe from the modes to
determine shape of the deformed sheet we will encounter geometric
inconsistencies. The problem of closure must be addressed to guarantee that a
particular mode pattern is geometrically valid, and therefore corresponds to a
possible continuous sheet configuration in three dimensions.

While approaches to the compatibility of the sheet may be analyzed using loop
closure relations, it is known that this question is readily addressed by the
Gauss-Codazzi-Mainardi-Peterson equations. That the first and second
fundamental forms must satisfy these equations is a requirement for
compatibility. Compatibility is required to ensure that the sheet may be
embedded in 3-dimensional space (i.e. that such a surface exists having the
corresponding pattern of fundamental forms). Further, according to the theorem
of Bonnet Toponogov and Rovenski (2006), this constitutes a sufficient set of
conditions to guarantee compatibility and therefore also guarantees a
corresponding low-energy configuration of the origami sheet.

The Gauss-Codazzi equations are captured in compact form in two parts. First
by the Gauss equation, which is a single scalar equation:

$R=2\kappa$ (41)

where

$\kappa\equiv\text{det}[h^{i}_{j}]=2\text{det}[g^{ik}h_{kj}]$ (42)

is the Gaussian curvature,

$R\equiv g^{ik}R^{j}_{ijk}$ (43)

is the (Ricci) scalar curvature defined in terms of the the inverse metric
$g^{ij}$ (satisfying $g^{ij}g_{jk}=\delta^{i}_{k}$) and the Riemann curvature
tensor capturing the intrinsic curvature of the sheet

$R^{l}_{ijk}\equiv\partial_{j}\Gamma^{l}_{ik}-\partial_{k}\Gamma^{l}_{ij}+\Gamma^{l}_{js}\Gamma^{s}_{ik}-\Gamma^{l}_{ks}\Gamma^{s}_{ij}$
(44)

which is in turn defined using the Christoffel symbol

$\Gamma^{l}_{ij}\equiv\frac{1}{2}g^{lk}\left(\partial_{j}g_{ki}+\partial_{i}g_{kj}-\partial_{k}g_{ij}\right)\,.$
(45)

Complementary to the Gauss equation are the Codazzi equations, which are
captured by the set of first order PDEs

$\displaystyle\epsilon^{ij}D_{i}h_{jk}=0$ (46)

which must hold for $k=\\{1,2\\}$. Here $D_{i}$ is the covariant derivative
defined by

$D_{i}h_{jk}\equiv\partial_{i}h_{jk}-\Gamma^{l}_{ij}h_{lk}-\Gamma^{l}_{ik}h_{lj}\,,$
(47)

and $\epsilon^{ij}$ is the totally antisymmetric tensor defined by component
as $\epsilon^{11}=\epsilon^{22}=0$ and $\epsilon^{12}=-\epsilon^{21}=1$.

These equations, for a general surface, impose three PDE constraints on the
six otherwise independent fields
$g_{11}(X_{1},X_{2}),g_{12}(X_{1},X_{2}),g_{22}(X_{1},X_{2}),h_{11}(X_{1},X_{2}),h_{12}(X_{1},X_{2}),h_{22}(X_{1},X_{2})$.
However, for the soft origami deformations, the six fundamental form fields
are determined by the three mode fields
$\\{\theta,b,t\\}\rightarrow\\{\theta(X_{1},X_{2}),b(X_{1},X_{2}),t(X_{1},X_{2})\\}$.
We therefore rewrite and simplify the GC equations by inserting the
fundamental forms from Eqs. 28&31. This process is tedious but
straightforward, the results detailed below.

For the Eq. 41 above, we first compute the Gaussian curvature
$\kappa=\text{det}[g^{ik}h_{kj}]$ which is the product of the sheet principal
curvatures. This reduces to

$\kappa=\frac{1}{[\lambda_{1}(\theta)\lambda_{2}(\theta)]^{2}}\left(b^{2}b_{1}(\theta)b_{2}(\theta)-t^{2}\right)\,.$
(48)

Next, the scalar curvature, which is the most tedious, may be expanded and
simplified before inserting our particular forms

$\displaystyle R=$ $\displaystyle
g^{ik}\left[\partial_{j}\Gamma^{j}_{ik}-\partial_{k}\Gamma^{j}_{ij}+\Gamma^{j}_{js}\Gamma^{s}_{ik}-\Gamma^{j}_{ks}\Gamma^{s}_{ij}\right]$
(49) $\displaystyle=$
$\displaystyle\frac{1}{2}g^{ik}\Big{\\{}(2\partial_{i}g_{kl}-\partial_{l}g_{ik})\partial_{j}g^{jl}-\partial_{i}g_{jl}\partial_{k}g^{jl}$
$\displaystyle\qquad+2g^{jl}\left(\partial_{i}\partial_{j}g_{kl}-\partial_{i}\partial_{k}g_{jl}\right)$
$\displaystyle\qquad+\frac{g^{jm}g^{ln}}{2}\Big{[}2\partial_{l}g_{mj}\partial_{i}g_{nk}-\partial_{l}g_{mj}\partial_{n}g_{ik}$
$\displaystyle\qquad-2\partial_{l}g_{mk}\partial_{i}g_{nj}+\partial_{l}g_{mk}\partial_{n}g_{ij}\Big{]}\Big{\\}}\,.$

After plugging in, this reduces to

$\displaystyle R$
$\displaystyle=\frac{2}{\lambda_{1}\lambda_{2}}\Bigg{[}\frac{\lambda_{1}^{\prime}\lambda_{2}^{\prime}(\partial_{2}\theta)^{2}}{\lambda_{2}^{2}}+\frac{\lambda_{1}^{\prime}\lambda_{2}^{\prime}(\partial_{1}\theta)^{2}}{\lambda_{1}^{2}}$
$\displaystyle-\frac{\lambda_{2}^{\prime\prime}(\partial_{1}\theta)^{2}}{\lambda_{1}}-\frac{\lambda_{1}^{\prime\prime}(\partial_{2}\theta)^{2}}{\lambda_{2}}-\frac{\lambda_{2}^{\prime}\partial_{1}\theta}{\lambda_{1}}-\frac{\lambda_{1}^{\prime}\partial_{2}\theta}{\lambda_{2}}\Bigg{]}$
(50)

where, as a reminder, $\lambda_{1}$ and $\lambda_{2}$ and their derivatives
(primes) are both functions of the mechanism parameter field $\theta$ even
though this dependence is suppressed here for ease of presentation. The Gauss
equation then may be assembled as a condition on the fields $\theta,b,t$ by
placing these relations Eq. 50&48 into Eq. 41, leading to main text Eq. 4.

Similar to the Gauss equation, the Codazzi equations reduce to

$\displaystyle 0=$
$\displaystyle\left(\partial_{1}t-\lambda_{1}\lambda_{1}^{\prime}\partial_{2}b-b(\lambda_{1}\lambda_{1}^{\prime\prime}+\lambda_{1}^{\prime}\lambda_{1}^{\prime})\partial_{2}\theta\right)$
$\displaystyle+\left(\frac{\lambda_{2}^{\prime}}{\lambda_{2}}-\frac{\lambda_{1}^{\prime}}{\lambda_{1}}\right)t\partial_{1}\theta+\left(\frac{\lambda_{1}^{\prime}}{\lambda_{1}}+\frac{\lambda_{2}^{\prime}}{\lambda_{2}}\right)b\lambda_{1}\lambda_{1}^{\prime}\partial_{2}\theta$
(51)

and

$\displaystyle 0=$
$\displaystyle-\left(\partial_{2}t+\lambda_{2}\lambda_{2}^{\prime}\partial_{1}b+b(\lambda_{2}\lambda_{2}^{\prime\prime}+\lambda_{2}^{\prime}\lambda_{2}^{\prime})\partial_{1}\theta\right)$
$\displaystyle+\left(\frac{\lambda_{2}^{\prime}}{\lambda_{2}}-\frac{\lambda_{1}^{\prime}}{\lambda_{1}}\right)t\partial_{2}\theta+\left(\frac{\lambda_{1}^{\prime}}{\lambda_{1}}-\frac{\lambda_{2}^{\prime}}{\lambda_{2}}\right)b\lambda_{2}\lambda_{2}^{\prime}\partial_{1}\theta\,.$
(52)

Combined with the Gauss equation, this imposes three PDE relations on the
three soft mode fields. This balance of equations and fields indicates that
the solution in the interior will be determined by the boundary conditions.
Specifically, the geometry of the soft mode will be determined by boundary
conditions on the fields $\theta,b,t$ Nonlinearly these may be restricted by
the lack of known existence and uniqueness theorems, yet the number of
possible solutions in the interior are expected to be proportional to the
boundary. The soft modes are therefore infinite in dimension, yet
subextensive.

###  B.4 Linear Deformations

The equations of compatibility determined in the previous section govern the
soft response of the origami metamaterial: each solution to these equations
guarantees a corresponding soft pattern of response which may arise under
loading. Therefore our task is to locate the space of solutions to these
reduced Gauss and Codazzi equations. For nonlinear strains, analytic solutions
are not obvious. However, linearizing these equations in the perturbative
modes $t(X_{1},X_{2}),b(X_{1},X_{2}),\delta\theta(X_{1},X_{2})$ allows closed-
form solutions to be obtained. Here,
$\delta\theta(X_{1},X_{2})\equiv\theta(X_{1},X_{2})-\theta_{0}$ implicitly
defines the quasi-planar state of folding (i.e. mechanism reference state) we
are deforming from. Such analysis makes clear the spatial character of these
soft modes, as well as the dualilty indicated in the main text.

Using $\lambda_{1}(\theta_{0})=\lambda_{2}(\theta_{0})=1$, the linearized
equations appear as

$\displaystyle 0$
$\displaystyle=(\partial^{2}_{1}-\frac{1}{\nu_{\text{in}}}\partial^{2}_{2})\delta\theta$
(53) $\displaystyle 0$
$\displaystyle=\partial_{1}t-\lambda_{1}^{\prime}\partial_{2}b$ (54)
$\displaystyle 0$
$\displaystyle=\partial_{2}t+\lambda_{2}^{\prime}\partial_{1}b$ (55)

where

$\nu_{\text{in}}=\frac{-\lambda_{2}^{\prime}}{\lambda_{1}^{\prime}}$ (56)

is the in-plane linear Poisson’s ratio of the planar folding mode around this
folding state $\theta_{0}$. These equations may be rewritten in a suggestive
form

$\displaystyle 0$ $\displaystyle=\partial_{w}\partial_{\bar{w}}\delta\theta$
(57) $\displaystyle 0$ $\displaystyle=\partial_{\bar{w}}H$ (58) $\displaystyle
0$ $\displaystyle=\partial_{w}\bar{H}$ (59)

where we have introduced new variables for the surface parametrization

$\displaystyle w$ $\displaystyle=X_{1}+\frac{1}{\gamma}X_{2}$ (60)
$\displaystyle\bar{w}$ $\displaystyle=X_{1}-\frac{1}{\gamma}X_{2}$ (61)

and another set of new variables capturing the soft curvature fields

$\displaystyle H$ $\displaystyle=t-\gamma\lambda_{2}^{\prime}b$ (62)
$\displaystyle\bar{H}$ $\displaystyle=t+\gamma\lambda_{2}^{\prime}b\,.$ (63)

Here, $\gamma=\sqrt{\frac{1}{\nu_{\text{in}}}}$ is defined for convenience and
the derivatives

$\displaystyle\partial_{w}$
$\displaystyle=\frac{1}{2}\left(\partial_{1}+\gamma\partial_{2}\right)$ (64)
$\displaystyle\partial_{\bar{w}}$
$\displaystyle=\frac{1}{2}\left(\partial_{1}-\gamma\partial_{2}\right)$ (65)

are implied by the requirements that
$\partial_{w}w=\partial_{\bar{w}}\bar{w}=1$ and
$\partial_{w}\bar{w}=\partial_{\bar{w}}w=0$. This compatibility theory has
notably decoupled the fields that generate curvature $t,b$ from the planar
mode $\theta$, into two compatibility problems that may be solved
independently. In the linear limit, these therefore divide into the “planar”
soft modes (those involving $\theta$) which do not generate curvature and the
“non-planar” modes involving $t$ and $b$ which do not influence the strain on
the surface.

The choices of transformed variables in Eqs. 60,61,62&63 are based largely on
that introduced in Czajkowski and Rocklin (2022), in which a related formalism
enabled general closed-form solutions to be written for the planar unimode
problem. Similarly, to solve for the non-planar soft modes, we note that all
solutions to Eqs. 58&59 must be written as $H(w,\bar{w})\rightarrow H(w)$ and
$\bar{H}(w,\bar{w})\rightarrow\bar{H}(\bar{w})$. However, these must also be
required to generate real-valued fields for $t=\frac{1}{2}(\bar{H}+H)$ and
$b=\frac{1}{2\lambda_{2}^{\prime}\gamma}(\bar{H}-H)$. Tuning the origami
metamaterial from anauxetic $\nu_{\text{in}}>0$ to auxetic $\nu_{\text{in}}<0$
changes the spatial coordinates $w,\bar{w}$ from pure real to complex. For the
auxetic origami sheet, we require that $H(w)^{*}=\bar{H}(\bar{w})$, while for
the anauxetic examples, these functions are simply required to be real-valued
themselves. Here, ∗ is the conventional complex conjugation.

The utility of this analytic result is displayed in the case of origami sheets
with the topology of the disc. In the auxetic case all solutions may be
assembled from analytic expansion $H(w)=\sum_{n}C_{n}w^{n}$ with arbitrarily
chosen complex coefficients, which then implies the expansion
$\bar{H}(w)=\sum_{n}D_{n}\bar{w}^{n}$, with $D_{n}=C_{n}^{*}$. Similarly, the
solutions for the anauxetic case are captured by the same sums, but with the
coefficients required to be real-valued.

Having solved the compatibility of the twist and bend modes, we turn our
attention to the planar mode. This particular question of planar motion
compatibility is again a previously solved problem. As shown in Czajkowski and
Rocklin (2022) the compatibility may also be usefully written

$\displaystyle\partial_{\bar{w}}\varepsilon$ $\displaystyle=0$ (66)
$\displaystyle\partial_{w}\bar{\varepsilon}$ $\displaystyle=0$ (67)

in terms of the $w,\bar{w}$ coordinates. Here, we have defined variables
capturing linear strain

$\displaystyle\varepsilon$ $\displaystyle=A\varepsilon_{\theta}-\gamma\phi$
(68) $\displaystyle\bar{\varepsilon}$
$\displaystyle=A\varepsilon_{\theta}+\gamma\phi$ (69)

where $\varepsilon_{\theta}$ is the (field) amplitude of strain arising from
the linear folding mechanism $\delta\theta$ and
$A=\sqrt{2/(1+\nu_{\text{in}}^{2})}$ is a mechanism strain normalization, and
$\phi$ is the coarse rotation of each continuum sheet element in the plane.
More explicitly, the mechanism strain field is defined so that the observed
unsymmetrized coarse strain $S_{ij}=\partial_{j}u_{i}$ from planar coarse
displacements $u_{1},u_{2}$ will appear as

$\mathrm{\varepsilon}=\varepsilon_{\theta}A\begin{bmatrix}1&0\\\
0&-\nu_{\text{in}}\end{bmatrix}+\phi\begin{bmatrix}0&-1\\\
1&0\end{bmatrix}\,.$ (70)

Finally, $\varepsilon_{\theta}=\delta\theta\frac{\lambda_{1}^{\prime}}{A}$
relates these scalar variables of planar strain.

The upshot of this additional formalism for linear deformations is that the
duality between the planar deformations and non-planar deformations is now
made explicit. Solutions for the curvature-inducing fields $H,\bar{H}$ may
readily be converted into solutions for the planar deformation fields
$\varepsilon,\bar{\varepsilon}$ and vice versa. This is accomplished with an
arbitrary multiplicative factor $\zeta$ with units of inverse length so that
the nonplanar mode corresponds to the planar mode via
$(H,\bar{H})\rightarrow\zeta*(\varepsilon,\bar{\varepsilon})$. Of course, this
may also be accomplished using the soft mode fields
$\varepsilon_{\theta},\phi,t,b$ directly. Furthermore, elements of each of the
dual soft mode spaces may be superimposed to span the soft deformations
available to the origami sheet.

It is straightforward to show that each field $\varepsilon_{\theta},\phi,t,b$
each individually satisfies the same second order PDE:
$\partial_{\bar{w}}\partial_{w}\theta=\partial_{\bar{w}}\partial_{w}\phi=\partial_{\bar{w}}\partial_{w}b=\partial_{\bar{w}}\partial_{w}t=0$.
To make this a bit more tangible, note that the second derivative
$\partial_{\bar{w}}\partial_{w}$ reduces to
$\frac{1}{4}(\partial^{2}_{1}-\frac{1}{\nu_{\text{in}}}\partial^{2}_{2})$,
which is a distorted Laplacian. Despite neglecting the connection between
$\phi$ and $\varepsilon_{\theta}$ and that between $t$ and $b$, these
relations capture the (identical) spatial patterns available to each mode.
Each mode may be written as the sum of two analytic functions i.e.
$\varepsilon_{\theta}=\frac{1}{2A}(\varepsilon(w)+\bar{\varepsilon}(\bar{w}))$.
Again, these functions may be independent real analytic functions in the
anauxetic case, and must be complex conjugates in the auxetic case. As
explored in more depth in Czajkowski and Rocklin (2022), this means that the
auxetic patterns will appear as boundary modes, which decay exponentially into
the bulk, while the anauxetic modes will appear as periodic functions composed
of plane waves. These characteristic behaviors are separated by an exceptional
point in the Poisson’s ratio at $\nu_{\text{in}}=0$, where the coordinate
transformations $w$ and $\bar{w}$ are identical and the mode becomes uniaxial.

##  Appendix C Simulation of Origami Sheets

Here we describe the process used to generate numerically force-balanced
states of origami sheets under point loading conditions.

###  C.1 Numerical investigations in a reduced model using the MERLIN software

To investigate the bulk force-balanced states available to the periodic fold
patterns in the Morph family, we employ the MERLIN2 software Liu and Paulino
(2017, 2018) developed by the Paulino research group at Georgia Tech. This
software, which runs in MATLAB, treats the origami sheet using a bar-hinge
model with the vertex positions (orange points in Fig. 6c) as the degrees of
freedom. To generate data for this manuscript, we employ the N4B5
triangulation scheme, in which a crease line is added to each panel along the
shorter of the diagonals as described below. This software has been highly
accurate in describing the deformation of real experimental origami sheets Liu
and Paulino (2017). As described in the main text, deformations of these
vertices incur energy penalties of three categories: bar stretching, crease
folding, and panel bending, captured by the terms

$E=E_{s}+E_{f}+E_{p}$ (71)

respectively.

Bar stretching penalizes changes in distance between any pair of nodes
connected by an edge such as the green and blue lines in Fig. 6c. An example
of such a bond connects node $u$ to $v$. The vector $\mathbf{r}_{uv}$ connects
these nodes in space and based on this the stretching energy penalty is

$E_{s}=\sum_{\text{bonds}\langle
uv\rangle}\frac{1}{2}a_{\text{bar}}d_{uv}W(\frac{(|\mathbf{r}_{uv}|-d_{uv})}{d_{uv}})$
(72)

where $a_{\text{bar}}$ is a cross-sectional area of the bar element, and
$d_{uv}$ is the distance between the nodes in the undeformed state (i.e.
length of the bar element). The bar energy density function

$W(x)=k_{s}\left(\sqrt{2x+1}-\frac{(2x+1)^{5/2}}{5}\right)$ (73)

captures an Ogden nonlinear constitutive model Liu and Paulino (2017). This
nonetheless resembles a system of conventional linear Hookean springs

$E_{s}\sim\sum_{\text{bonds}\langle
uv\rangle}\frac{a_{\text{bar}}k_{s}}{2d_{uv}}(|\mathbf{r}_{uv}|-d_{uv})^{2}$
(74)

for small strains of the bar elements. Even for the nonlinear macroscopic
deformations explored in main text Fig. 1 and 2, the bar elements are only
strained by linearly small amounts. Note that $k_{s}$ has the units of a
compression modulus or Young’s modulus ($J/m^{3}$).

The panel bending is intended to capture the effect indicated in the name.
Each parallelogram panel is itself a thin plate, which is therefore
susceptible to isometric deformations. These deformations are approximated by
the addition of a crease along the shorter of the two diagonals of the
parallelogram. The addition of this crease defines two triangular panels of
the parallelogram, and the panel bending energy penalizes the difference
between orientation of the normal direction of each panel, as indicated by the
dihedral angle $\rho_{p}$ in Fig. 6c. The energy penalty is computed from
these dihedral angles of panel bending via

$E_{p}=\sum_{\text{panels}-p}d_{p}k_{p}(\rho_{p}-\pi)^{2}\,,$ (75)

where the sum is taken over panels $p$, $d_{p}$ is the rest length of the
panel diagonal crease line and $k_{p}$ is the modulus of panel bending, which
has the units of line tension.

This crease line across the panel diagonal is also included as a panel
stretching spring in Eq. 74. This is essential to prevent the simulations from
accessing unrealistic zero energy deformations which would strain the
otherwise stiff parallelogram panels they aim to capture.

The energy of crease folding is captured in a similar manner to the panel
bending, but with a distinct rotational stiffness and the dihedral angles
having a rest position distinct from $\pi$. To calculate this energy requires
knowledge of each of the rest state dihedral angles $\rho_{q}^{0}$ at each
crease $q$ (i.e. given the initial quasi-planar state of the morph, what is
the angle of the crease $q$?). From this the energy is calculated via

$E_{f}=\sum_{\text{creases}q}d_{q}k_{f}(\rho_{q}-\rho_{q}^{0})^{2}\,,$ (76)

where the sum is taken over all the (bulk) crease segments of the folding
pattern, and $k_{f}$ is the modulus of crease folding.

###  C.2 Applied loading and constitutive parameters

The main text data is categorized into two types of constitutive models,
designated the “rigid” versus the “realistic” origami sheets. These
designations determine the combinations of the sheet constitutive parameters
$k_{s},k_{p},k_{f},A$ used in the simulations (as defined in Sec. C.1). The
primary difference between these is the presence of a finite stiffness against
folding at the creases in the “realistic” version, which is absent in the
“rigid” origami. They are therefore named according to their intended
interpretation, where “realistic” emulates a folded sheet of material, while
“rigid” emulates origami assembled from rigid panels connected by
comparatively energy-free hinges. For all the simulation data herein, lengths
are measured against the uniform crease length $a_{1}=a_{2}=c=1$.

To obtain the so-named “realistic” origami sheet configurations displayed in
Figs. 1&2, we employ the “auto” mode within the MERLIN2 code. Here, the
constitutive parameters for the bar-hinge model are specified indirectly, via
the specification of the Poisson’s ratio $\nu_{\text{sheet}}$ and Youngs
modulus $Y$ of the thin sheet, as well as a quantity which controls the
folding stiffness of the panel diagonals relative to the folded crease lines
$L^{*}=\frac{K_{b}}{k_{f}}$. The factor

$K_{b}=\frac{Yl_{t}^{3}}{12(1-\nu_{\text{sheet}}^{2})}$ (77)

is the bending modulus of the constitutive thin sheet, which is determined by
$Y$ and $\nu_{\text{sheet}}$, and in turn determines the stiffness of the
panel bending and crease line bending. We assign these parameters to values
which are appropriate for conventional thin sheets, choosing
$\nu_{\text{sheet}}=0.3$, $l_{t}=0.01$, $L^{*}=20$. The two origami patterns
simulated in main text Fig. 1 have crease geometry set by $\beta=1$,
$\alpha=\pi-\beta$. The folding is set to $\gamma=0.01$ for the top imitation
of an “uncreased sheet” in Fig. 1, which appears to be a sufficiently small
folding to remove the exotic effects of the crease pattern, while the lower
folded portion is set to a folding angle $\gamma=0.5$. Recall the precise
definition of $\gamma$ presented in App. A. The loading here is chosen to be
incompatible with a cylindrical configuration, expected for the isometries of
least bending and therefore least (small) energy. Then for main text Fig. 2,
the identical origami sheet is simulated as the folded example from Fig. 1,
but with different (arbitrary) displacement loading conditions applied to the
four corners of the sheet. For all of the simulations herein, the energy scale
is determined by choosing the Young’s modulus $Y=1$.

For the “rigid” data, we utilize the “manual” setting for the MERLIN2
simulation framework, allowing us to more directly set the parameters
$k_{f},k_{p},k_{s}$ and the cross-sectional area for the bar elements
$a_{\text{bar}}$. In this “rigid” data, the defining feature is that the
folding modulus of the major creases $k_{f}$ is set to zero, thereby realizing
the ideal mechanism limit, in which the planar mode costs zero energy. The
scale of energy is determined by the chosen panel bending modulus $k_{p}=1$.
Aside from these choices, and aside from the parameters we explicitly vary as
indicated, we use representative values for a sheet of paper, or similar
material, with thickness $l_{t}$. For a common piece of origami folded from
copy paper, the thickness (relative to the crease length) may be $\sim 0.01$.

The remaining undetermined parameters required for these simulations in the
manual setting are $k_{s},k_{p}$ and $a_{\text{bar}}$. To provide context for
these and estimate reasonable values, we first consider the bar elements,
which capture the stretching of the sheet in its own plane. In this case the
cross sectional area of a bar-element approximation of such phenomena will be
approximately $a_{\text{bar}}\sim l_{t}l_{\nu}$, and each “rigid” simulation
is hence run with $a_{\text{bar}}=0.01$. We may roughly approximate $k_{s}$
with the Young’s modulus $Y$. In contrast, the stiffness against bending of a
panel will scale as $k_{p}\sim\frac{l_{t}^{3}Y}{l_{\nu}}$, where the
denominator contains the linear dimension of the curved portion of the panel,
which we presume to be approximated by the edge length. Then the ratio between
the stretching and panel bending moduli will be of the order
$\frac{k_{s}}{k_{p}}\sim\frac{t^{3}}{l_{\nu}}\sim 10^{6}$ in the natural
length units of the system (chosen to be $l_{\nu}$). Or, connecting to the
dimensionless parameter in Fig. 5, $\frac{k_{s}}{ak_{p}}\sim 10^{6}$.
Therefore, with the the length $a_{1}$ set to unity, the default value of
$k_{s}$ in simulations will be $10^{6}$. This dramatic discrepancy between the
energy cost of deformations which stretch the constitutive material versus
those which simply bend (i.e. fold) reflects the isometric limit.

The dataset displayed in main text Fig. 3 is obtained using simulations of
origami sheets which extend by 20 unit cells in each lattice direction. As
suggested by the above analysis, the stretching modulus is set to
$k_{s}=10^{6}$ for these “rigid” simulations. These are run with reference
configuration geometry defined by the parameters in Table. 1. Each is run with
an identical displacement loading applied to the vertices at the four corners
of the quasi rectangular origami sheet. In this loading the vertices are
displaced according to the vectors
$\\{\mathbf{d}_{i}\\}=\\{\\{0.8,2.4,-2.0\\},\\{1.4,-0.8,2.4\\},\\{-2.8,2.4,-3.4\\},\\{0.6,-4.0,3.0\\}\\}$
(listed counter-clockwise from the bottom-left in the coordinate system
defined by the orthonormal vectors of the initial configuration
$\hat{l}^{(1)},\hat{l}^{(2)},\hat{\mathbf{N}}$ respectively). The simulation
will attempt to gradually apply these displacements in $70$ equal increments,
occasionally breaking some increments into finer steps in order to ensure
convergence of the force-balancing algorithm.

# | $\alpha$ | $\beta$ | $\theta_{0}$
---|---|---|---
1 | 1.4 | 1.0 | 0.5
2 | 1.6 | 1.0 | 0.5
3 | 1.7 | 1.0 | 0.5
4 | 2.01 | 1.0 | 0.5
5 | 2.1 | 1.0 | 0.5
6 | 1.4 | 1.0 | 0.7
7 | 1.5 | 1.0 | 0.7
8 | 2.01 | 1.0 | 0.7
Table 1: Origami reference state geometries for the nonlinear data displayed
in main text Fig. 3

To generate the data in main text Fig. 4, we replicate the approach to the
data in Fig. 3 for a modified set of sheet configurations. Again the modulus
of folding $k_{f}$ is zero while the modulus of panel bending sets the energy
scale $k_{p}=1$ and the modulus of panel stretching remains $k_{s}=10^{6}$ at
the previous value. Here, numerical data is collected for systems with unit
cell geometry defined by the angles $\theta_{0}=0.5$, $\beta=1.0$, and a
series of different $\alpha=\\{1.4,1.45,1.5,1.6,1.7,1.8,2.01,2.1,2.2\\}$,
generating the variety of lines, colored according to the reference state
Poisson’s ratio this implies. For each of these sheet configurations, data is
run at the system sizes $N=\\{10,20,30,40\\}$ and each of these is run at each
of the loading conditions below

| center-bottom | center-left | center-top | center-right
---|---|---|---|---
load 1 $\frac{1}{N}\mathbf{u}$ | $\\{0.00015,-0.001,0.00075\\}$ | $\\{0.0002,0.0006,-0.0005\\}$ | $\\{0.00035,-0.0002,0.0006\\}$ | $\\{-0.0007,0.0006,-0.00085\\}$
load 2 $\frac{1}{N}\mathbf{u}$ | $\\{0.0002,0.0006,-0.0005\\}$ | $\\{0.00015,-0.001,0.00075\\}$ | $\\{0.00035,-0.0002,0.0006\\}$ | $\\{-0.0007,0.0006,-0.00085\\}$
load 3 $\frac{1}{N}\mathbf{u}$ | $\\{0.0002,0.0006,-0.0005\\}$ | $\\{0.00035,-0.0002,0.0006\\}$ | $\\{0.00015,-0.001,0.00075\\}$ | $\\{-0.0007,0.0006,-0.00085\\}$
load 4 $\frac{1}{N}\mathbf{u}$ | $\\{0.0002,0.0006,-0.0005\\}$ | $\\{0.00035,-0.0002,0.0006\\}$ | $\\{-0.0007,0.0006,-0.00085\\}$ | $\\{0.00015,-0.001,0.00075\\}$
(a) Loads applied to vertices located at the center of each system edge

| bottom-left | bottom-right | top-left | top-right
---|---|---|---|---
load 5 $\frac{1}{N}\mathbf{u}$ | $\\{0.00015,-0.001,0.00075\\}$ | $\\{0.0002,0.0006,-0.0005\\}$ | $\\{0.00035,-0.0002,0.0006\\}$ | $\\{-0.0007,0.0006,-0.00085\\}$
load 6 $\frac{1}{N}\mathbf{u}$ | $\\{0.0002,0.0006,-0.0005\\}$ | $\\{0.00015,-0.001,0.00075\\}$ | $\\{0.00035,-0.0002,0.0006\\}$ | $\\{-0.0007,0.0006,-0.00085\\}$
load 7 $\frac{1}{N}\mathbf{u}$ | $\\{0.0002,0.0006,-0.0005\\}$ | $\\{0.00035,-0.0002,0.0006\\}$ | $\\{0.00015,-0.001,0.00075\\}$ | $\\{-0.0007,0.0006,-0.00085\\}$
load 8 $\frac{1}{N}\mathbf{u}$ | $\\{0.0002,0.0006,-0.0005\\}$ | $\\{0.00035,-0.0002,0.0006\\}$ | $\\{-0.0007,0.0006,-0.00085\\}$ | $\\{0.00015,-0.001,0.00075\\}$
(b) Loads applied to vertices located at the system corners

Table 2: Displacement loading conditions used to obtain numerical data
displayed in main text Figs. 4&5

which are sufficiently small loadings that this data corresponds with the
mechanics of linear deformations. Note that the values listed represent the
imposed displacement divided by the number of unit cells across $(N)$, thereby
indicating the manner in which the loads are scaled with system size. However,
while this system size feature in our code is important for the convergence of
nonlinear loading, this has negligible effect for the linear data of main text
Figs. 4&5. Data in these figures is generated using the analysis methods in
Sec. D, and averaging these results over the different loads in Table 2.

To generate the configuration data analyzed in main text Fig. 5 follows a very
similar procedure as that for Fig. 4, yet for a different series of
parameters. In this case, only the system size $N=20$ unit cells across is
considered and we probe instead the dependence on the panel stretching
stiffness. The data points in main text Fig. 5 are then generated at each of
the values
$k_{s}=\\{10^{1},10^{2},10^{3},10^{4},10^{5},10^{6},10^{7},10^{8}\\}$ and at
each of the loads and folding geometries from the above description of main
text Fig. 4. Again, averaging the fitting error data over the different loads
generates the displayed values. The method of quantifying the error in the
fitting is explained in Section. D.2 below.

##  Appendix D Analysis of Origami Sheets

Here we detail the methods to analyze the simulation data for a force-balanced
origami sheet in both the nonlinear and linear regimes. This will include the
hyperribbon error and the fitting to soft mode error.

###  D.1 Adherence to shape hyperribbon

The main text includes quantification of the observed data obtained in
Section. C via the quality of adherence of the sheet midplane shape to the
predicted hyperribbon described in Section. B.2. To quantify such adherence,
we must first evaluate the fundamental forms of the sheet midplane which
quantify the coarse shape locally. To do so, we follow previously employed
methods for the evaluation of nonlinear deformation in lattice metamaterials
Czajkowski  _et al._ (2022); Czajkowski and Rocklin (2022). We begin with the
evaluation of the metric. As the lattice system is spatially discrete, we do
not directly measure a continuous tensor field, but rather estimate local
tensor values at each unit cell grid index $(n_{1},n_{2})$. As described in
the beginning of Sec. B, the metric of deformation for an initially uniform
lattice may be estimated using the reference (initial)
$\mathbf{l}^{(1)}_{0},\mathbf{l}^{(2)}_{0}$ and deformed (final)
$\mathbf{l}^{(1)},\mathbf{l}^{(2)}$ lattice vectors via Eq. 25. To extract the
lattice vectors $\mathbf{l}^{(1)}(n_{1},n_{2}),\mathbf{l}^{(2)}(n_{1},n_{2})$
at the arbitrary unit cell index $(n_{1},n_{2})$, a common and simple choice
is to use the vector which connects the position of the $\mu$-th node
$\mathbf{R}_{\mu}(n_{1},n_{2})$ in this unit cell to the position of its
counterpart in the next unit cell $\mathbf{R}_{\mu}(n_{1}+1,n_{2})$. However,
the asymmetry in such a choice yields error in the estimation of the metric at
the order of the lattice lengths themselves and the gradients of the metric.
Conveniently, a symmetric choice to define these lattice vectors, averaging
“backward” and “forward” estimates of such together, eliminates this error in
favor of smaller, higher order errors. In addition, as noted in Ref. McInerney
_et al._ (2022), these lattice vectors must be averaged over the nodes of the
unit cell to avoid spurious symmetry-violating extra terms in the estimation
of these fundamental forms. With this, our formulas for extraction of the
lattice vectors from the numerical data are

$\displaystyle\mathbf{l}^{(1)}(n_{1},n_{2})$
$\displaystyle\sim\frac{1}{4}\sum_{\mu\in\text{unit
cell}}\frac{1}{2}\left(\mathbf{R}_{\mu}(n_{1}+1,n_{2})-\mathbf{R}_{\mu}(n_{1}-1,n_{2})\right)$
(78) $\displaystyle\mathbf{l}^{(2)}(n_{1},n_{2})$
$\displaystyle\sim\frac{1}{4}\sum_{\mu\in\text{unit
cell}}\frac{1}{2}\left(\mathbf{R}_{\mu}(n_{1},n_{2}+1)-\mathbf{R}_{\mu}(n_{1},n_{2}-1)\right)\,.$
(79)

With this, relying on the definition established in Eq. 25, the metric

${g}^{(\text{obs})}_{ij}(n_{1},n_{2})=\frac{\mathbf{l}^{(i)}(n_{1},n_{2})\cdot\mathbf{l}^{(j)}(n_{1},n_{2})}{|\mathbf{l}^{(i)}_{0}||\mathbf{l}^{(j)}_{0}|}$
(80)

is readily extracted for each unit cell in the simulation data; and similarly
for the strain using Eq. 26. Here, the superscript indicates that this is the
observed value of the metric, so estimated from the data.

The lattice vectors extracted using the recipes in Eqs. 78&79 are also used to
estimate the local value of the curvature tensor, and again the symmetric
definition averaged over the unit cell eliminates significant leading error.
First the lattice vectors are used to first define the local normal of the
smoothed sheet

$\hat{\mathbf{N}}(n_{1},n_{2})=\frac{\mathbf{l}^{(1)}(n_{1},n_{2})\times\mathbf{l}^{(2)}(n_{1},n_{2})}{|\mathbf{l}^{(1)}(n_{1},n_{2})\times\mathbf{l}^{(2)}(n_{1},n_{2})|}$
(81)

which is useful in defining the curvature. Further, the gradient of this
normal with respect to the sheet parametrization coordinates must be known.
This is achieved by again using symmetric derivatives so that we write

$\displaystyle\partial_{1}\hat{\mathbf{N}}(n_{1},n_{2})$
$\displaystyle=\frac{1}{2}(\hat{\mathbf{N}}(n_{1}+1,n_{2})-\hat{\mathbf{N}}(n_{1}-1,n_{2}))\partial_{2}\hat{\mathbf{N}}(n_{1},n_{2})$
$\displaystyle=\frac{1}{2}(\hat{\mathbf{N}}(n_{1},n_{2}+1)-\hat{\mathbf{N}}(n_{1},n_{2}-1))\,.$
(82)

With this, the recipe for the curvature tensor components

$h_{ij}=\begin{bmatrix}-\mathbf{l}^{(1)}\cdot\partial_{1}\hat{\mathbf{N}}&-\mathbf{l}^{(2)}\cdot\partial_{1}\hat{\mathbf{N}}\\\
-\mathbf{l}^{(1)}\cdot\partial_{2}\hat{\mathbf{N}}&-\mathbf{l}^{(2)}\cdot\partial_{2}\hat{\mathbf{N}}\end{bmatrix}$
(83)

where the dependence on the unit cell index location is now suppressed for
convenience, and index (lhs) versus matrix (rhs) representation are mixed for
simplicity of notation.

Once the fundamental forms are determined across the origami sheet, we must
determine the quality of adherence of these tensors to the hyperribbon. This
is achieved by first obtaining the closest value of the folding parameter
$\theta$ and then evaluating the twist and bend components of the curvature
tensor given this value of $\theta$. The folding angle is determined most
readily using the local value of $\mathbf{l}^{(2)}$. In the case of a
uniformly applied folding motion, this lattice vector length is determined by
the formula $|\mathbf{l}^{(2)}|=2a\cos(\theta)$ as in Section. A Eq. 23.
Inverting this equation for the folding angle gives a sufficiently accurate
estimate for the folding angle.

Having obtained the folding angle, it is possible to then also determine the
twist and bend mode amplitudes locally. The local twist mode is readily
obtained as the only off-diagonal component of the curvature tensor, so that
$t(n_{1},n_{2})=h_{12}(n_{1},n_{2})$. Note that using the other off-diagonal
component is equally valid $t(n_{1},n_{2})=h_{21}(n_{1},n_{2})$ and in
practice, the average pf these two estimates is employed to remove bias.
Extracting the bend component is similarly achievable using either of the
diagonal components of curvature $h_{11}$ or $h_{22}$. Explicitly, either of
the formulas

$\displaystyle b(n_{1},n_{2})$
$\displaystyle=\frac{h_{11}(n_{1},n_{2})}{\lambda_{1}(\theta(n_{1},n_{2}))\lambda_{1}^{\prime}(\theta(n_{1},n_{2}))}$
$\displaystyle b(n_{1},n_{2})$
$\displaystyle=-\frac{h_{22}(n_{1},n_{2})}{\lambda_{2}(\theta(n_{1},n_{2}))\lambda_{2}^{\prime}(\theta(n_{1},n_{2}))}$

where $\lambda_{1}$ and $\lambda_{2}$ are the known principal stretch
functions introduced in Sec. B.1. Again, in practice each choice of estimation
method is roughly equivalent and the two estimates are averaged to remove
bias.

Finally, having obtained local estimates of the soft mode fields, we can
quantify the accuracy of this hyperribbon description. To do this, we then
produce local estimates of the projection of the fundamental forms onto the
hyperribbon ${g}^{(\text{proj})}_{ij}(n_{1},n_{2})$ and
${h}^{(\text{proj})}_{ij}(n_{1},n_{2})$ at each unit cell. This is achieved by
reconstructing the fundamental forms from the soft mode field data using the
discrete version of Eqs. 39&40

$\displaystyle{g}^{(\text{proj})}_{ij}(n_{1},n_{2})$
$\displaystyle=g^{*}_{ij}(\theta(n_{1},n_{2}))$ (84)
$\displaystyle{h}^{(\text{proj})}_{ij}(n_{1},n_{2})$
$\displaystyle=h^{*}_{ij}(\theta(n_{1},n_{2}),t(n_{1},n_{2}),b(n_{1},n_{2}))\,,$
(85)

where the functions $g^{*}_{ij}$ and $h^{*}_{ij}$ are defined in Sec. B.1 to
capture the hyperribbon. With this, the normalized tensor error in either the
metric or the curvature may (first) be readily constructed from these two
tensor fields via the tensor error function

$\Delta[t^{(1)}_{ij},t^{(2)}_{ij}]=\sqrt{\frac{\sum_{n_{1},n_{2}}\sum_{ij}(t^{(1)}_{ij}(n_{1},n_{2})-t^{(2)}_{ij}(n_{1},n_{2}))^{2}}{\sum_{n_{1},n_{2}}\sum_{ij}t^{(1)}_{ij}(n_{1},n_{2})t^{(1)}_{ij}(n_{1},n_{2})}}\,.$
(86)

which appropriately yields the fractional difference between two discrete
tensor grids $t^{(2)}_{ij}(n_{1},n_{2})$ and $t^{(2)}_{ij}(n_{1},n_{2})$ which
are close in value. It is important to notice that this error evaluation
method is most appropriate for tensors which are zero in the undeformed state.
Otherwise, the quality of error is artificially inflated for small
deformations. The curvature tensor satisfies this condition, and therefore the
appropriate dimensionless measure of the curvature portion of the hyperribbon
error will be
$\Delta_{h}=\Delta[{h}^{(\text{obs})}_{ij},{h}^{(\text{proj})}_{ij}]$.
However, to quantify the fractional non-hyperribbon component of the metric we
instead use the deviational part of this tensor, which is the strain. The
strain is constructed by
$\varepsilon_{ij}=\frac{1}{2}\left(g_{ij}-g_{ij}^{(0)}\right)$, where the
initial metric $g_{ij}^{(0)}$ in the case here is the identity $\delta_{ij}$.
Then this portion of the hyperribbon error corresponding to the metric is
$\Delta_{\varepsilon}=\Delta\left[\frac{1}{2}\left({g}^{(\text{obs})}_{ij}-\delta_{ij}\right),\frac{1}{2}\left({g}^{(\text{proj})}_{ij}-\delta_{ij}\right)\right]$.

###  D.2 Analytic soft mode fitting

The theory presented in the main text includes not only the assertion the
local fundamental forms be characterized by the three mode amplitudes
(hyperribbon), but also prescribes the manner in which these modes will be
distributed across the sheet spatially. In the limit of linear deformations,
we indicate that these modes must be controlled by sheared analytic functions,
detailed in Sec. B.4. In order to evaluate this hypothesis, in this section we
present an analysis procedure used to take a sheet geometric configuration and
extract a single dimensionless number describing how close the state is to
such a sheared analytic mode. This method is used to produce the data
displayed in main text Fig. 5. As with the above hyperribbon error analysis,
this is performed separately for the strain and curvature portions of the
fundamental form variations. Note that this analysis closely follows that
established in Ref. Czajkowski and Rocklin (2022).

This fitting procedure first follows the analysis in Sec. D.1 to obtain a
numeric estimate of the local mode values
$\theta(n_{1},n_{2}),t(n_{1},n_{2}),b(n_{1},n_{2})$ for each unit cell grid
index $(n_{1},n_{2})$ in the sheet. Each unit cell is affiliated with a
location

$\mathbf{R}(n_{1},n_{2})\equiv\frac{1}{4}\sum_{\mu\in\begin{smallmatrix}\text{unit}\\\
\text{cell}\end{smallmatrix}}\mathbf{R}_{\mu}(n_{1},n_{2})$ (87)

by summing over the four vertex locations within the unit cell, again similar
to that described in the previous section for the lattice vectors. Here,
$\mathbf{R}_{\mu}(n_{1},n_{2})$ is the location of the $\mu$-th of the 4
vertices of the $(n_{1},n_{2})$ unit cell. Equipped with a sufficient set of
field and position data, we set out to fit each field individually to an
appropriate analytic form. As detailed in Section B.4, linear combinations of
the curvature fields,

$\displaystyle H$ $\displaystyle=t-\gamma\lambda_{2}^{\prime}b=p_{1}(w)$ (88)
$\displaystyle\bar{H}$
$\displaystyle=t+\gamma\lambda_{2}^{\prime}b=p_{2}(\bar{w})\,,$ (89)

will be captured by analytic functions of the sheared reference-frame spatial
variables $w,\bar{w}$ (termed “sheared analytic functions”) for the lowest-
energy deformations. Linear combinations of the strain variables are captured
similarly

$\displaystyle\varepsilon$
$\displaystyle=A\varepsilon_{\theta}-\gamma\phi=q_{1}(w)$ (90)
$\displaystyle\bar{\varepsilon}$
$\displaystyle=A\varepsilon_{\theta}+\gamma\phi=q_{2}(\bar{w})\,.$ (91)

Values for these spatial variables $\\{w(n_{1},n_{2}),\bar{w}(n_{1},n_{2})\\}$
are easily obtained from the initial unit cell position data
$\\{\mathbf{R}_{0}(n_{1},n_{2})\\}$ via the recipes
$w(\mathbf{R}_{0})=\mathbf{R}_{0}\cdot\mathbf{l}^{(1)}_{0}+\sqrt{\nu_{\text{in}}}\,\mathbf{R}_{0}\cdot\mathbf{l}^{(2)}_{0}$
and
$\bar{w}(\mathbf{R}_{0})=\mathbf{R}_{0}\cdot\mathbf{l}^{(1)}_{0}-\sqrt{\nu_{\text{in}}}\,\mathbf{R}_{0}\cdot\mathbf{l}^{(2)}_{0}$.
Here, we borrow from the local basis notation constructed using the lattice
vector directions in the previous section (D.1). Choosing the coordinate
system so that the origami sheet midplane initially lies in the x-y plane,
these may be pictured as the x and y spatial components. Note that because
these are reference frame coordinates they do not change with the sheet
deformation, and are always defined using $\\{\mathbf{R}_{0}(n_{1},n_{2})\\}$.

Solving directly for the mechanism strain amplitude $\varepsilon_{\theta}$ in
Eqs. 90&91 we see that this will be captured by a particular analytic form

$\tilde{\varepsilon_{\theta}}(w,\bar{w})=\frac{1}{2A}(q_{1}(w)+q_{2}(\bar{w})).$
(92)

Ignoring the rotation field $\phi$, this captures the essential restriction on
the possible spatial forms that this mechanism strain may take through space.
The metric of deformation is independent of this rotation field and evaluating
the adherence of the first fundamental form variations to our predictions is
accomplished by fitting the spatial data for this mechanism strain to the form
above. This procedure appears different on each side of the exceptional point
at $\nu_{\text{in}}=0$ and it is therefore useful to include these details
more carefully. We consider the case of the mechanism strain pattern first,
with the procedure for the twist and bend fields following a nearly identical
procedure as discussed near the end of this section. In what follows we will
absorb the factor of $\frac{1}{2A}$ into the definition of the strain modes,
so that the procedure applies equally to each of the modes we will fit
$(\varepsilon_{\theta},t,b)$ which are each fit to the form
$(q_{1}(w)+q_{2}(\bar{w}))$.

####  D.2.1 Auxetic Case

For origami sheets with auxetic mechanism $\nu_{\text{in}}<0$, the quantity
$\gamma=\sqrt{1/\nu_{\text{in}}}$ is complex-valued, and so are the position
data $\\{w(n_{1},n_{2}),\bar{w(n_{1},n_{2})}\\}$. Therefore, analytic
functions $q_{1},q_{2}$ will necessarily also be complex valued. The
combination of these two functions to produce the local mechanism strain
amplitude, however, is required to be real-valued. Therefore, in the auxetic
case it will be required that these functions are complex conjugates of one
another. As these functions admit expansions,

$\displaystyle q_{1}(w)$ $\displaystyle=\sum_{n=0}^{N_{c}}C_{n}w^{n}$ (93)
$\displaystyle q_{2}(\bar{w})$
$\displaystyle=\sum_{n=0}^{N_{c}}D_{n}\bar{w}^{n}$ (94)

and since the bar operation on $w$ is here equivalent to complex conjugation,
the requirement that $\varepsilon_{\theta}$ be real-valued is captured by
$C_{n}=D^{*}_{n}$. Therefore, finding the closest fit for the mechanism strain
data is a process of varying the $C_{n}$ complex coefficients to minimize the
fitting error function

$\delta_{\text{fit}}[\\{\varepsilon_{\theta}^{(k)}\\},\tilde{\varepsilon_{\theta}}()]=\sum_{k}\left|\varepsilon_{\theta}^{(k)}-\tilde{\varepsilon_{\theta}}(w_{k},\bar{w}_{k})\right|^{2}$
(95)

where all interior unit cells are now indexed by a single number $k$ so that
$\mathbf{R}_{k}=\mathbf{R}(n_{1}(k),n_{2}(k))$, and again a boundary layer of
unit cells which are two or less away from the system edge are excluded from
the analysis to minimize such boundary effects. For the data explored here, we
choose to have $N_{c}=20$.

To solve for the coefficients which minimize the fitting error function, we
first break down into real and imaginary parts $C_{n}=A_{n}+iB_{n}$. In this
case the function becomes

$\delta_{\text{fit}}=\sum_{k}\left[\varepsilon_{\theta}^{(k)}-\sum_{n}(2A_{n}\mathrm{Re}[w_{k}^{n}]-2B_{n}\mathrm{Im}[w_{k}^{n}])\right]^{2}$
(96)

and the minimizing coefficients are then the solution of the set of equations

$\displaystyle\frac{\partial\delta_{\text{fit}}}{\partial A_{m}}=$
$\displaystyle
0=-2\sum_{k}\left[\varepsilon_{\theta}^{(k)}-\sum_{n}(2A_{n}\mathrm{Re}[w_{k}^{n}]-2B_{n}\mathrm{Im}[w_{k}^{n}])\right]2\mathrm{Re}[w_{k}^{m}]$
(97) $\displaystyle\frac{\partial\delta_{\text{fit}}}{\partial B_{m}}=$
$\displaystyle
0=2\sum_{k}\left[\varepsilon_{\theta}^{(k)}-\sum_{n}(2A_{n}\mathrm{Re}[w_{k}^{n}]-2B_{n}\mathrm{Im}[w_{k}^{n}])\right]2\mathrm{Im}[w_{k}^{m}]$
(98)

which must be satisfied simultaneously for every $m$. This is a linear algebra
problem which may be put conveniently into a matrix form

$\begin{bmatrix}\mathbf{f}_{1}\\\
\mathbf{f}_{2}\end{bmatrix}=\begin{bmatrix}\mathrm{M}_{1A}&\mathrm{M}_{1B}\\\
\mathrm{M}_{2A}&\mathrm{M}_{2B}\end{bmatrix}\cdot\begin{bmatrix}\mathbf{A}\\\
\mathbf{B}\end{bmatrix}.$ (99)

Here, $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$ are vectors having components
$(f_{1})_{m}=\sum_{k}\varepsilon_{\theta}^{(k)}\mathrm{Re}[w_{k}^{m}]$ for
$m={0,1,2,...N_{c}-1}$ and
$(f_{2})_{m}=-\sum_{k}\varepsilon_{\theta}^{(k)}\mathrm{Im}[w_{k}^{m}]$ for
$m={1,2,...N_{c}-1}$, while
$\mathrm{M}_{1A},\mathrm{M}_{2A},\mathrm{M}_{1B},\mathrm{M}_{2B}$ are matrices
defined by

$\displaystyle(\mathrm{M}_{1A})_{mn}$
$\displaystyle=2\sum_{k}\mathrm{Re}[w_{(k)}^{m}]\mathrm{Re}[w_{(k)}^{n}]$
(100) $\displaystyle m=\\{0,1,...N_{c}-1\\}\quad n=\\{0,1,...N_{c}-1\\}$ (101)
$\displaystyle(\mathrm{M}_{2A})_{mn}$
$\displaystyle=2\sum_{k}\mathrm{Im}[w_{(k)}^{m}]\mathrm{Re}[w_{(k)}^{n}]$
(102) $\displaystyle m=\\{1,2,...N_{c}-1\\}\quad n=\\{0,1,...N_{c}-1\\}$ (103)
$\displaystyle(\mathrm{M}_{1B})_{mn}$
$\displaystyle=2\sum_{k}\mathrm{Re}[w_{(k)}^{m}]\mathrm{Im}[w_{(k)}^{n}]$
(104) $\displaystyle m=\\{0,1,...N_{c}-1\\}\quad n=\\{1,2,...N_{c}-1\\}$ (105)
$\displaystyle(\mathrm{M}_{2B})_{mn}$
$\displaystyle=2\sum_{k}\mathrm{Im}[w_{(k)}^{m}]\mathrm{Im}[w_{(k)}^{n}]$
(106) $\displaystyle m=\\{1,2,...N_{c}-1\\}\quad n=\\{1,2,...N_{c}-1\\}$ (107)

and the vectors of coefficients are

$\displaystyle\mathbf{A}^{T}$ $\displaystyle=\\{A_{0},A_{1},...A_{N_{c}-1}\\}$
(108) $\displaystyle\mathbf{B}^{T}$
$\displaystyle=\\{B_{1},B_{2},...B_{N_{c}-1}\\}.$ (109)

Note that the element $B_{0}$ has no impact on the values of the function
$\tilde{\varepsilon_{\theta}}(w,\bar{w})$, and we have already omitted the
corresponding row and column from the elements in Eq. 99 leaving $B_{0}$
undetermined, and avoiding numerical issues with solving the linear algebra
problem. This linear system of equations (Eq. 99) is solved in Mathematica,
yielding the coefficients $C_{n}$ that define the closest fit function.

The above procedure obtains the closest fit sheared analytic mode matching to
the mechanism strain field $\tilde{\varepsilon_{\theta}}()$. The procedure is
carried out identically to procure predictive fittings to the twist and the
bend fields each, respectively yielding the closest fit functions
$\tilde{t}()$ and $\tilde{b}()$. It is worth noting that this procedure does
not as carefully probe that the combination field $H\equiv
t-\sqrt{\frac{1}{\nu_{\text{in}}}}\lambda_{2}^{\prime}b$ precisely adheres to
a single sheared analytic mode. Rather, this process confirms that the twist
and bend each individually obey the conditions entailed by the sheared
analytic solutions. This is like checking a complex function for analyticity
by instead checking that the real and imaginary components each obeyed the
laplace equation, without checking the harmonic conjugacy directly. Further
investigations may improve upon this method. However, the quality of the data
and other informal measures indicate that these fields indeed obey the sheared
analytic generalization of harmonic conjugacy.

Given these closest fit functions
$\tilde{\varepsilon_{\theta}}(),\tilde{b}(),\tilde{t}()$, we would like to now
evaluate how close the estimates are to the target data they are intended to
approximate. This is accomplished similar to in the previous section, by first
extracting the corresponding strain

$\varepsilon^{(\text{pred})}_{ij}(n_{1},n_{2})=\tilde{\varepsilon_{\theta}}\left(w(n_{1},n_{2}),\bar{w}(n_{1},n_{2})\right)A\begin{bmatrix}1&0\\\
0&-\nu_{\text{in}}\end{bmatrix}$ (110)

and the corresponding curvature

$h^{(\text{pred})}_{ij}(n_{1},n_{2})=\tilde{t}\left(w(n_{1},n_{2}),\bar{w}(n_{1},n_{2})\right)\begin{bmatrix}0&1\\\
1&0\end{bmatrix}+\tilde{b}\left(w(n_{1},n_{2}),\bar{w}(n_{1},n_{2})\right)\begin{bmatrix}b_{1}(\theta_{0})&0\\\
0&b_{2}(\theta_{0})\end{bmatrix}.$ (111)

Then the error in the fitting for the strain is accomplished by reusing the
error function (Eq. 86) from Section D.1, this time using the analytic
predictions of the tensors rather than the projections into the hyperribbon.
For the case of the curvature, then, the fractional error in the analytic
fitting is
$\Delta^{(fit)}_{h}=\Delta[h^{(\text{pred})}_{ij},{h}^{(\text{obs})}_{ij}]$,
which takes two grids of tensors and reduces them to their root-mean-square
difference normalized by their root-mean-square projection onto one another.
Similarly, this fractional error in fitting the strain is
$\Delta^{(fit)}_{\varepsilon}=\Delta[\varepsilon^{(\text{pred})}_{ij},{\varepsilon}^{(\text{obs})}_{ij}]$,
where
${\varepsilon}^{(\text{obs})}_{ij}=\frac{1}{2}\left({g}^{(\text{obs})}_{ij}-\delta_{ij}\right)$.

####  D.2.2 Anauxetic Case

The procedure for closest fitting to the soft modes in the case of the
anauxetic mechanism is distinct yet quite similar. The important distinction
comes from $\gamma=\sqrt{\frac{1}{\nu_{\text{in}}}}$ which is now real-valued;
therefore the coordinates $w$ and $\bar{w}$ are also real-valued everywhere.
Therefore instead of the inter-relation between complex coefficients of the
analytic expansions

$\displaystyle q_{1}(w)$ $\displaystyle=\sum_{n=0}^{N_{c}}C_{n}w^{n}$ (112)
$\displaystyle q_{2}(\bar{w})$
$\displaystyle=\sum_{n=0}^{N_{c}}D_{n}\bar{w}^{n}$ (113)

we instead have real and independent coefficients $\\{C_{n}\\}$ and
$\\{D_{n}\\}$. Similar to the auxetic case, we again minimize the fitting
error function, which now takes the form

$\delta_{\text{fit}}=\sum_{k}\left[\varepsilon_{\theta}^{(k)}-\sum_{n}C_{n}w_{k}^{n}-\sum_{n}D_{n}\bar{w}_{k}^{n}\right]^{2}.$
(114)

We then take partial derivatives to minimize, again leading to a set of linear
equations, captured in matrix form

$\begin{bmatrix}\mathbf{f}_{1}\\\
\mathbf{f}_{2}\end{bmatrix}=\begin{bmatrix}\mathrm{M}_{1C}&\mathrm{M}_{1D}\\\
\mathrm{M}_{2C}&\mathrm{M}_{2D}\end{bmatrix}\cdot\begin{bmatrix}\mathbf{C}\\\
\mathbf{D}\end{bmatrix}.$ (115)

Here, $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$ are vectors having components
$(f_{1})_{m}=\sum_{k}\varepsilon_{\theta}^{(k)}\mathrm{Re}[w_{k}^{m}]$ for
$m={0,1,2,...N_{c}-1}$ and
$(f_{2})_{m}=-\sum_{k}\varepsilon_{\theta}^{(k)}\mathrm{Im}[w_{k}^{m}]$ for
$m={1,2,...N_{c}-1}$, while
$\mathrm{M}_{1A},\mathrm{M}_{2A},\mathrm{M}_{1B},\mathrm{M}_{2B}$ are matrices
defined by

$\displaystyle(\mathrm{M}_{1C})_{mn}$
$\displaystyle=\sum_{k}w_{(k)}^{m}w_{(k)}^{n}$ (116) $\displaystyle
m=\\{0,1,...N_{c}-1\\}\quad n=\\{0,1,...N_{c}-1\\}$ (117)
$\displaystyle(\mathrm{M}_{2C})_{mn}$
$\displaystyle=\sum_{k}\bar{w}_{(k)}^{m}w_{(k)}^{n}$ (118) $\displaystyle
m=\\{1,2,...N_{c}-1\\}\quad n=\\{0,1,...N_{c}-1\\}$ (119)
$\displaystyle(\mathrm{M}_{1D})_{mn}$
$\displaystyle=\sum_{k}w_{(k)}^{m}\bar{w}_{(k)}^{n}$ (120) $\displaystyle
m=\\{0,1,...N_{c}-1\\}\quad n=\\{1,2,...N_{c}-1\\}$ (121)
$\displaystyle(\mathrm{M}_{2D})_{mn}$
$\displaystyle=\sum_{k}\bar{w}_{(k)}^{m}\bar{w}_{(k)}^{n}$ (122)
$\displaystyle m=\\{1,2,...N_{c}-1\\}\quad n=\\{1,2,...N_{c}-1\\}$ (123)

and the vectors of coefficients are

$\displaystyle\mathbf{C}^{T}$ $\displaystyle=\\{C_{0},C_{1},...C_{N_{c}-1}\\}$
(124) $\displaystyle\mathbf{D}^{T}$
$\displaystyle=\\{D_{1},D_{2},...D_{N_{c}-1}\\}.$ (125)

Again, one element and row have been omitted to account for the single
indeterminacy. In this case it is because $C_{0}$ and $D_{0}$ are redundant in
determining the uniform part of the resulting soft mode. In this case it is a
choice to omit $D_{0}$.

Solving these equations numerically in Mathematica yields the prediction for
the modes, and analyzing the fractional error in this fitting follows the same
procedure as the auxetic case, leading to the data in main text Fig. 5.

## References

  * Jellett (1849) J. H. Jellett, The Transactions of the Royal Irish Academy 22, 343 (1849).
  * Cerda and Mahadevan (1998) E. Cerda and L. Mahadevan, Physical Review Letters 80, 2358 (1998).
  * Witten (2007) T. A. Witten, Reviews of Modern Physics 79, 643 (2007).
  * MIURA and Koryo (1985) MIURA and Koryo, The Institute of Space and Astronautical Science report 618, 1 (1985).
  * Tachi (2009) T. Tachi, Journal of the International Association for Shell and Spatial Structures 50, 2287 (2009).
  * Schenk and Guest (2013) M. Schenk and S. D. Guest, Proceedings of the National Academy of Sciences of the United States of America 110, 3276 (2013).
  * Pratapa  _et al._ (2019) P. P. Pratapa, K. Liu, and G. H. Paulino, Physical Review Letters 122, 155501 (2019).
  * Dieleman  _et al._ (2019) P. Dieleman, N. Vasmel, S. Waitukaitis, and M. van Hecke, Nature Physics , 1 (2019).
  * Feng  _et al._ (2020) F. Feng, X. Dang, R. D. James, and P. Plucinsky, Journal of the Mechanics and Physics of Solids 142, 104018 (2020).
  * McInerney  _et al._ (2022) J. McInerney, G. H. Paulino, and D. Z. Rocklin, Proceedings of the National Academy of Sciences 119, e2202777119 (2022).
  * Mahadevan and Rica (2005) L. Mahadevan and S. Rica, Science 307, 1740 (2005).
  * Audoly and Boudaoud (2008) B. Audoly and A. Boudaoud, undefined 56, 2444 (2008).
  * Christensen  _et al._ (2015) J. Christensen, M. Kadic, O. Kraft, and M. Wegener, MRS Communications 5, 453 (2015).
  * Bertoldi  _et al._ (2017) K. Bertoldi, V. Vitelli, J. Christensen, and M. van Hecke, Nature Reviews Materials 2, 17066 (2017).
  * Brittain  _et al._ (2001) S. T. Brittain, O. J. Schneller, H. Wu, S. Whitesides, and G. M. Whitesides, Journal of Physical Chemistry B 105, 347 (2001).
  * Temmerman  _et al._ (2007) N. D. Temmerman, M. Mollaert, T. V. Mele, and L. D. Laet, http://dx.doi.org/10.1260/026635107782218868 22, 161 (2007).
  * Randhawa  _et al._ (2011) J. S. Randhawa, K. E. Laflin, N. Seelam, and D. H. Gracias, Advanced Functional Materials 21, 2395 (2011).
  * Rogers  _et al._ (2016) J. Rogers, Y. Huang, O. G. Schmidt, and D. H. Gracias, MRS Bulletin 41, 123 (2016).
  * Li  _et al._ (2019) S. Li, H. Fang, S. Sadeghi, P. Bhovad, and K. W. Wang, Advanced Materials 31, 1805282 (2019).
  * Okuzaki  _et al._ (2008) H. Okuzaki, T. Saido, H. Suzuki, Y. Hara, and H. Yan, Journal of Physics: Conference Series 127, 012001 (2008).
  * Hawkes  _et al._ (2010) E. Hawkes, B. An, N. M. Benbernou, H. Tanaka, S. Kim, E. D. Demaine, D. Rus, and R. J. Wood, Proceedings of the National Academy of Sciences of the United States of America 107, 12441 (2010).
  * Jeong and Lee (2018) D. Jeong and K. Lee, Robotica 36, 261 (2018).
  * Lazarus  _et al._ (2019) N. Lazarus, G. L. Smith, and M. D. Dickey, Advanced Intelligent Systems 1, 1900059 (2019).
  * Son  _et al._ (2022) H. Son, Y. Park, Y. Na, and C. K. Yoon, Polymers 2022, Vol. 14, Page 4235 14, 4235 (2022).
  * An and Rus (2012) B. An and D. Rus, Proceedings - IEEE International Conference on Robotics and Automation , 3299 (2012).
  * Felton  _et al._ (2014) S. Felton, M. Tolley, E. Demaine, D. Rus, and R. Wood, Science 345, 644 (2014).
  * Miyashita  _et al._ (2016) S. Miyashita, S. Guitron, K. Yoshida, S. Li, D. D. Damian, and D. Rus, in  _Proceedings - IEEE International Conference on Robotics and Automation_, Vol. 2016-June (Institute of Electrical and Electronics Engineers Inc., 2016) pp. 909–916.
  * Misseroni  _et al._ (2022) D. Misseroni, P. P. Pratapa, K. Liu, and G. H. Paulino, Extreme Mechanics Letters 53, 101685 (2022).
  * Waitukaitis  _et al._ (2015) S. Waitukaitis, R. Menaut, B. G. G. Chen, and M. V. Hecke, Physical Review Letters 114, 055503 (2015).
  * Stern  _et al._ (2017) M. Stern, M. B. Pinson, and A. Murugan, Physical Review X 7, 041070 (2017).
  * Chen and Santangelo (2018) B. G. G. Chen and C. D. Santangelo, Physical Review X 8, 011034 (2018).
  * Lee-Trimble  _et al._ (2022) M. E. Lee-Trimble, J. H. Kang, R. C. Hayward, and C. D. Santangelo, Soft Matter 18, 6384 (2022).
  * Czajkowski and Rocklin (2021) M. Czajkowski and Z. Rocklin, Bulletin of the American Physical Society , B17 (2021).
  * Czajkowski  _et al._ (2022) M. Czajkowski, C. Coulais, M. van Hecke, and D. Z. Rocklin, Nature Communications 2022 13:1 13, 1 (2022).
  * Czajkowski and Rocklin (2022) M. Czajkowski and D. Z. Rocklin, “Duality and sheared analytic response in mechanism-based metamaterials,” (2022), arXiv:2205.10751 [cond-mat.soft] .
  * Zheng  _et al._ (2022) Y. Zheng, I. Niloy, P. Celli, I. Tobasco, and P. Plucinsky, Physical Review Letters 128, 208003 (2022).
  * Nassar  _et al._ (2017) H. Nassar, A. Lebée, and L. Monasse, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473, 20160705 (2017).
  * Nassar  _et al._ (2022) H. Nassar, A. Lebée, and E. Werner, Extreme Mechanics Letters 53, 101722 (2022).
  * Xu  _et al._ (2023) H. Xu, I. Tobasco, and P. Plucinsky, “Derivation of an effective plate theory for parallelogram origami from bar and hinge elasticity,” (2023), arXiv:2311.10870 [cond-mat.soft] .
  * Wei  _et al._ (2013) Z. Y. Wei, Z. V. Guo, L. Dudte, H. Y. Liang, and L. Mahadevan, Physical Review Letters 110, 215501 (2013).
  * Toponogov and Rovenski (2006) V. Toponogov and V. Rovenski,  _Differential Geometry of Curves and Surfaces: A Concise Guide_ (Birkhäuser Boston, 2006).
  * Filipov  _et al._ (2017) E. T. Filipov, K. Liu, T. Tachi, M. Schenk, and G. H. Paulino, International Journal of Solids and Structures 124, 26 (2017).
  * Liu and Paulino (2017) K. Liu and G. H. Paulino, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473 (2017), 10.1098/RSPA.2017.0348.
  * Pretko and Radzihovsky (2018) M. Pretko and L. Radzihovsky, Physical Review Letters 120, 195301 (2018).
  * Fruchart  _et al._ (2020) M. Fruchart, Y. Zhou, and V. Vitelli, Nature 577, 636 (2020).
  * McInerney  _et al._ (2020) J. McInerney, B. G. G. Chen, L. Theran, C. D. Santangelo, and D. Z. Rocklin, Proceedings of the National Academy of Sciences of the United States of America 117, 30252 (2020).
  * Belcastro and Hull (2002) S. M. Belcastro and T. C. Hull, Linear Algebra and Its Applications 348, 273 (2002).
  * Liu and Paulino (2018) K. Liu and G. Paulino, Origami 7, 1167 (2018).