\section{Introduction}
Origami has long surpassed its artistic vocation. Its early uses in architectural geometry date back to the 1920s \cite{wingler1978bauhaus, Lebee2015}. By now, origami techniques have been applied to the design of deployable structures in astronautics, robotics and biomedical engineering \cite{Kuribayashi2006, Schenk2011a, Zirbel2013, Callens2017}. Origami tessellations in particular make up versatile morphing 2D metamaterials that can access various low-energy high-curvature configurations \cite{Seffen2012,norman2009phd}.
Although inextensible in principle, the fine crease patterns of origami metamaterials enable them to effectively deform non-isometrically. In this letter, we determine the strains that are compatible with coarse-grained origami kinematics. Strain compatibility permits to fully characterize the 3D-embedded geometry when suitable boundary conditions are prescribed. In most cases however, the geometry further solves an elastostatic equilibrium for which we calculate the appropriate elasticity functional.

Morphing in origami metamaterials can be triggered in response to different loading conditions such as bending, twisting and stretching all of which couple in a manner unparalleled in other plates and shells; see Fig.~\ref{fig:models}. For instance, standard pure bending kinematics of plates imply that the curvatures in the bending plane and in the plane orthogonal to it, $\kappa_1$ and $\kappa_2$, are in a proportion equal and opposite to the in-plane Poisson's ratio $\nu$, i.e., $\nu=-\kappa_2/\kappa_1$. Thus auxetic plates tend to bend synclastically while anauxetic plates tend to bend anticlastically \cite{Naboni2016}. Origami metamaterials do not \cite{Schenk2011}. In fact, in what follows, we demonstrate that they typically exhibit the exact opposite behavior, i.e., $\nu=+\kappa_2/\kappa_1$. Other loads that induce morphing, such as ``pinching'', are unique to origami \cite{Grey2019}. Indeed, in plates and membranes, a pinch usually produces highly localized and singular deformations. By contrast, pinching an origami metamaterial brings on global and nearly-uniform deformations \cite{Seppecher2011, Filipov2015}. The consequence, as we shall prove, is that membrane energy penalizes, not strains, but their gradients, i.e., the Christoffel symbols~\cite{Mindlin1964, Gray2006, ciarlet2006}.

\begin{figure}[ht!]
    \centering
    \includegraphics[width=\linewidth]{figures/figure0.png}
    \caption{Two views of an origami-folded hyperboloid: A periodic ``Mars'' pattern is laser-etched into a PET sheet then hand-folded, circularly bent and fixated with adhesive tape. The pattern is auxetic but bends into a saddle. The letter explains the observed self-equilibrium geometry.}
    \label{fig:models}
\end{figure}

\begin{figure*}[ht!]
    \centering
    \includegraphics[width=\linewidth]{figures/figure1.pdf}
    \caption{Zigzag sums. (a) A zigzag sum in a periodic reference configuration: translating one zigzag along another sweeps a polyhedral surface corresponding to a crease pattern complete with a natural assignment of mountain and valley folds. (b) A super-cell in an abstract reference state: indices $(ij)$ are discrete coordinates for the vertices; the central node is at $\bm\xi\equiv\bm\xi^{(00)}$. (c) A super-cell in a physical folded state: vectors $(\t u,\t v, \t s,\t t)$ are $\bm\xi$-dependent non-dimensional vectors aligned with the creases; they fully characterize fold-only configurations as well as the tangent vectors $\t x_1$ and $\t x_2$. (d) A super-cell following an infinitesimal bending $\delta\t x$ of the panels: note how the parallelograms are slightly distorted now. (e) A global view of a finitely deformed zigzag sum $(\Sigma)$: each super cell is designated by its central node's coordinates $\bm\xi=(\xi_1,\xi_2)$ (red dot) and is obtained by the process illustrated in panels (b-d). (f) An orthogonal zigzag sum: initially orthogonal zigzag directions remain orthogonal during folding ($g_{12}\equiv 0$). (g) An oblique zigzag sum: initially non-orthogonal zigzag directions are sheared during folding; even then, $g_{12}$ remains constant.}
    \label{fig:zigzags}
\end{figure*}

Previous contributions have succeeded in describing several aspects of origami mechanics, including calculations of Poisson's ratios and bending moduli \cite{Wei2013, Schenk2013, Pratapa2019}. Hereafter, we derive an effective medium theory complete with its field equations of compatibility and equilibrium. The theory naturally accounts for the observations described above and shows how the different deformation modes couple kinematically and energetically. The subject of the theory is a class of origami patterns that we call ``zigzag sums''. These are generated by translating one zigzag line along another so as to sweep a periodic polyhedral surface $\Sigma$ whose edges and panels we interpret as an origami crease pattern; see Fig.
~\ref{fig:zigzags}a. Equivalently, $\Sigma$ is a periodic crease pattern with a unit cell composed of four parallelograms. Four quadrilaterals per unit cell is as simple as origami tessellations can be without being trivial; moreover, here, the quads are parallelograms. Nonetheless, zigzag sums span various well-known patterns (e.g., the Miura-ori, the eggbox, ``Sakoda's staircase'' and ``Barreto's Mars'' patterns \cite{Lang2018}) and encompass their oblique, non-developable and non-flat-foldable variants by maintaining two basic requirements: periodicity and parallelism. Meanwhile, the zigzag generators provide a handful of geometric seeds that can readily tune their deformation paths as well as the underlying elastic behavior.

\begin{figure*}[ht!]
    \centering
    \includegraphics[width=\linewidth]{figures/figure2.pdf}
    \caption{Convergence analysis. (a) Convergence speed: the maximum error in equation~\eqref{eq:Poisson} decreases like $r$ or, equivalently, like the reciprocal of the number of unit cells per row. The plot also shows that the maximum angular deflection $r\delta^{(ij)}$ decreases in the same fashion confirming the fold-to-bend orders of magnitude adopted in~\eqref{eq:scaling} and that, in the limit $r\to 0$, all panels are flat. The inset shows the unit cell adopted for this example. (b,c) Two iterates: snapshots of the way in which the zigzag sum converges to a smooth surface taken for $64$ and $128$ unit cells per row. Note that the zigzag sum is finitely and non-uniformly stretched and bent. Illustrated surfaces are constructed by ``propagating'' inextensibility constraints away from a given boundary; see~\cite{Nassar2018b} for details.}
    \label{fig:convergence}
\end{figure*}

\section{Geometric mechanics of zigzag sums}
The deformations of a zigzag sum $\Sigma$ are the outcome of two competing mechanisms \cite{Lechenault2014}. The first is the folding and unfolding of $\Sigma$ along its edges. The second is the isometric bending of its panels. Signature origami kinematics emerges when the bending stiffness of the panels is significantly larger than the folding stiffness of the edges.
Accordingly, we can assume that individual panels bend infinitesimally; $\Sigma$, on the other hand, can still embrace finitely curved surfaces by accumulating small deflections across many unit cells. Formally, we characterize the deformations of $\Sigma$ in terms of a continuous parametrization $\t x=\t x(\bm\xi)$ of its mid-plane, where $\t x$ is the current position of vertex $\bm\xi$, \highlight{with $\bm\xi=(\xi_1,\xi_2)$ being the effective curvilinear coordinate along the zigzags directions}; see Fig.~\ref{fig:zigzags}b-e. Consider then a non-elementary unit cell, i.e., a supercell, centered about vertex $\bm\xi$ and containing the neighboring vertices $\bm\xi^{(ij)}$ placed at $\t x^{(ij)}$. Then, the adopted origami kinematics imply that position
\[\label{eq:scaling}
    \t x^{(ij)} = \t X^{(ij)} + r^2\delta\t x^{(ij)}
\]
is attained by a fold-only motion $\bm\xi^{(ij)}\mapsto\t X^{(ij)}$ (Fig.~\ref{fig:zigzags}b-c) followed by a bending-induced perturbation $r^2\delta\t x^{(ij)}$ that is small compared to the characteristic size of the panels denoted $r$ (Fig.~\ref{fig:zigzags}c-d). In the limit of a tessellation where $\Sigma$ contains a large number of unit cells, i.e., for $r\to 0$, the plane locally tangent to $\Sigma$ is spanned by two vectors
\[
\begin{split}
\t x_1\equiv\frac{\partial\t x}{\partial\xi_1} &= \lim_{r\to 0}\frac{\t x^{(10)}-\t x^{(\bar10)}}{r}=\t u-\t v,\\
\t x_2\equiv \frac{\partial\t x}{\partial\xi_2} &= \lim_{r\to 0}\frac{\t x^{(01)}-\t x^{(0\bar1)}}{r}=\t s-\t t,
\end{split}
\]
where $\t u$, $\t v$, $\t s$ and $\t t$ are the crease vectors meeting at vertex $\bm\xi$ and normalized with respect to $r$. The apparent in-plane deformations of $\Sigma$ can therefore be quantified in terms of an effective metric tensor $\t g$ of components
\[\label{eq:constr1}
g_{\mu\nu}=\scalar{\t x_\mu,\t x_\nu} = \begin{bmatrix}
\norm{\t u-\t v}^2 & \scalar{\t u-\t v,\t s-\t t} \\
\scalar{\t u-\t v,\t s-\t t} & \norm{\t s-\t t}^2
\end{bmatrix}.
\]
Note that bending-related contributions disappear in the limit $r\to 0$ and can be disregarded for now.

It is noteworthy that $g_{12}=\scalar{\t u,\t s} + \scalar{\t v,\t t} - \scalar{\t u,\t t}-\scalar{\t v,\t s}$ is a motion constant (i.e., $\mathrm{d} g_{12}\equiv 0$) as it is equal to a combination of constant lengths and angles. Thus, if $g_{12}=0$ so that the zigzags $\xi_1$- and $\xi_2$-contours are initially orthogonal, then orthogonality is maintained in any subsequent motion (Fig.~\ref{fig:zigzags}f). However, in oblique tessellations, with $g_{12} = \norm{\t x_1}\norm{\t x_2}\cos\theta \neq 0$, stretching couples to shearing according to
\[
    \tan\theta\,\mathrm{d}\theta =
    \frac{\mathrm{d}\!\norm{\t x_1}}{\norm{\t x_1}}
    +
    \frac{\mathrm{d}\!\norm{\t x_2}}{\norm{\t x_2}}
    =\frac{\mathrm{d} g}{2g_{11}g_{22}}
\]
with $g=\det\t g$. Rather intuitively then, shearing, understood as aligning the zigzags ($\theta\to 0$ or $\pi$), always goes to reduce the pattern's effective area (Fig~\ref{fig:zigzags}g). As for the stretch ratios in directions $1$ and $2$, they are related to one another through a geometric Poisson's ratio $\nu\equiv
    -g_{11}/g_{22}\times\mathrm{d} g_{22}/\mathrm{d} g_{11}$.
A general expression for $\nu$ is pursued in the Supplemental Material \cite[App. A]{note}.

Most important is the observation that the fold-only motion $\bm\xi^{(ij)}\mapsto\t X^{(ij)}$ maintains crease parallelism and periodicity: this is truly where the fact that $\Sigma$ is made out of parallelograms is crucial. Hence, folding alone produces no curvature and, in that regard, it is bending that takes over now. This is most apparent in the expressions of the second-order derivatives which, in the limit $r\to 0$, read
\[
\begin{split}
\t x_{11}&=\delta\t x^{(20)}-2\delta\t x+\delta\t x^{(\bar{2}0)},\\
\t x_{22}&=\delta\t x^{(02)}-2\delta\t x+\delta\t x^{(0\bar{2})},\\
\t x_{12}&=\delta\t x^{(11)}-\delta\t x^{(\bar{1}1)}-\delta\t x^{(1\bar{1})}+\delta\t x^{(\bar{1}\bar{1})},
\end{split}
\]
where $\t x_{\mu\nu}\equiv\partial^2\t x/\partial\xi_\mu\partial\xi_\nu$. Therefore, the $\t x_{\mu\nu}$ at $\bm\xi$ are linear functions of the bending DOFs of the super cell centered at $\bm\xi$. \highlight{That being said, not all $\delta\t x^{(ij)}$ are compatible with the inextensibility constraints. Ultimately, it is possible to show that $(\t x_{11}, \t x_{22}, \t x_{12})$ belong to a linear subspace of $\mathbb{R}^3\times\mathbb{R}^3\times\mathbb{R}^3$ spanned by 4 DOFs, each attributed to a planarity defect $(\delta^{(ij)})$ of one of the four panels of the central unit cell (Fig~\ref{fig:convergence}a) \cite[App. B]{note}.} By the rank-nullity theorem, the $\t x_{\mu\nu}$ satisfy $9-4=5$ linear compatibility equations. Straightforward calculations show that the first two of these equations are $g_{12,1}=g_{12,2}=0$, i.e., they re-produce $\mathrm{d} g_{12}=0$. The other three relations are far more interesting; they combine into the vector identity \cite[App. B]{note}
\[\label{eq:identity}
\t x_{22} = \frac{\scalar{\t u,\t t\wedge\t s}\scalar{\t v,\t t\wedge\t s}}
{\scalar{\t t,\t u\wedge\t v}\scalar{\t s,\t u\wedge\t v}}\t x_{11}.
\]
The tangent components of $\t x_{22}$ and $\t x_{11}$ describe the strain gradients in the mid-plane of $\Sigma$ in direction $\t x_1$ and $\t x_2$ whereas the normal components describe the normal curvatures $\kappa_1$ and $\kappa_2$ of the mid-plane in the same directions. The above identity states that these respective quantities necessarily occur in equal proportions. Specifically, the in-plane and out-of-plane Poisson's ratios, $\nu$ and $-\kappa_2/\kappa_1$, of a zigzag sum $\Sigma$ are equal and opposite and are given by the current geometry of the folds as per
\[\label{eq:Poisson}
\nu = \frac{\kappa_2}{\kappa_1} = \frac{g_{11}}{g_{22}}\frac{\scalar{\t u,\t t\wedge\t s}\scalar{\t v,\t t\wedge\t s}}
{\scalar{\t t,\t u\wedge\t v}\scalar{\t s,\t u\wedge\t v}},
\]
with $\kappa_\mu=\scalar{\t x_{\mu\mu},\t n}/g_{\mu\mu}$ (no sum), $\t n = \t x_1\wedge\t x_2/\norm{\t x_1\wedge\t x_2}$ being the unit normal \cite[App. C]{note}.

Identity~\eqref{eq:Poisson} is of an asymptotic nature: in the limit $r\to 0$ of infinitesimal unit cells, the Poisson's ratios of a zigzag sum $\Sigma$ converge to one another. In other, more practical, terms, the equality is accurate within an error proportional to $r$. See Fig~\ref{fig:convergence} for a numerical illustration of the convergence process and the relevant error bounds. Identity~\eqref{eq:Poisson} is also universal: it holds whether $\Sigma$ is developable, flat-foldable, symmetric, rectangular, or not. Therein, the mixed products (e.g., $\scalar{\t u,\t t\wedge\t s}$) have their signs determined by the mountain-valley assignment of the concerned folds. This allows us to deduce the rather nice corollary: if a zigzag sum $\Sigma$ has an odd number of mountain folds meeting at a vertex then: $(i)$ it is necessarily auxetic; and $(ii)$ it necessarily bends into a saddle. In particular, all developable zigzag sums are auxetic and bend into saddles. Note that these results are only valid for nondegenerate configurations, i.e., where all folds are partially folded. Indeed, when a fold is flat, be it open or closed, mountain and valley assignments coalesce and the Poisson's ratio could become null, infinite or multivalued \cite{Pratapa2019}. Note also that the proven identity features the elongations and curvatures in two specific material directions $\xi_1$- and $\xi_2$-contours, be them orthogonal or not; in particular, the identity bears no immediate consequences on the Poisson's coefficient defined for two general orthogonal directions.

\begin{figure*}[ht!]
    \centering
    \includegraphics[width=\linewidth]{figures/figure3.pdf}
    \caption{``Mars''-folded origami pillars are one-sheeted hyperboloids. (a) A unit cell of the ``Mars'' pattern: the supplementary opposite angles ensure developability and flat-foldability; an equilateral pattern is further rectangular. (b) Action of the first DOF on the shape of the pillar: given a small enough aspect ratio $L_2/L_1$, the pillar can be inverted inside-out; here, the pattern has $N=100$ unit cells per meridian, and $8$ cells per parallel, with a flat-unfolded aspect ratio $L_2/L_1= 8\%$. (c) Action of the second DOF: the pattern flattens as the opening angle $\angle(\t u,\t v)$ at a given parallel is increased; here, the pattern has $N=100$ unit cells per meridian and $40$ unit cells per parallel with $L_2/L_1= 40\%$. (d) Bending energy: plot of normalized total bending energy as a function of the opening angle $\angle(\t u,\t v)$ at the equator for origami pillars centered about the equator with $L_2/L_1=40\%$; the discretely summed energy for $N=100$ is plotted against the continuum limit $N\to\infty$; data points $(i-iv)$ are the ones illustrated on panel (c); greyed out domain cannot be accessed with the current aspect ratio; the continuum theory satisfactorily predicts both energy and kinematical admissibility. Adopted values: $\alpha=\pi/3$, $\beta=\pi/2$.}
    \label{fig:hyperboloids}
\end{figure*}

\section{Self-equilibrium of origami pillars}
The equality of the Poisson's ratios greatly constrains the 3D geometries accessible to a zigzag sum $\Sigma$ \cite{Nassar2017a, Nassar2017e, Nassar2018b, zheng2021continuum}. Consider, for tractability, the ``Mars'' tessellation of Fig.~\ref{fig:hyperboloids}a. The tessellation lacks mirror symmetry but is otherwise developable, flat-foldable and equilateral. We find that compatible metrics and curvatures satisfy \cite[App. D]{note}
\[\label{eq:metricEx}
\begin{gathered}
g_{12}=0,\quad
g_{11}(4-g_{22}) =
4(\cos\alpha+\cos\beta)^2,\\
\frac{\kappa_2}{\kappa_1} = \nu =
-4\frac{(\cos\alpha+\cos\beta)^2}{g_{11}g_{22}}.
\end{gathered}
\]
This is in fact a system of non-linear PDEs weighing on the admissible configurations $\t x$. For instance, we explore the configurations of a rectangular domain of flat dimensions $L_1\times L_2$ that are folded and wrapped around an axis of symmetry into an origami pillar \cite{Lang2018}. These pillars are surfaces of revolution with a Cartesian parametrization of the form
\[
\t x(\bm\xi)=
(\rho(\xi_2)\cos(\omega\xi_1),\rho(\xi_2)\sin(\omega\xi_1),z(\xi_2)).
\]
Here, the axis of symmetry is the $z$-axis, $\xi_1$-contours are meridian lines and $\xi_2$-contours are parallel lines. Equations~\eqref{eq:metricEx} then simplify into a system of ODEs governing $z$ and $\rho$. Its solutions constitute a 2-DOF family of one-sheeted hyperboloids (Fig.~\ref{fig:hyperboloids}b,c) \cite[App. E]{note}. The range of motion of each pillar is bounded by the maximally folded and unfolded states of the pattern. Thus, for increasing aspect ratio $L_2/L_1$, the pillar will have a smaller range of motion, until it jams, or even becomes impossible to form without tearing.

Among all accessible pillars, the ones that are in self-equilibrium exhibit minimum levels of strain energy $\psi$. Strain energy is composed of two contributions: crease folding and panel bending. Crease folding energy takes the form of a classical, however nonlinear, membrane energy. In what follows, we neglect this contribution and focus on the less-explored influence of panel bending \cite{Grey2020}. For ``Mars''-folded origami pillars, centered about the equator, strain energy is exemplified on Fig.~\ref{fig:hyperboloids}d. It appears that unfolding the pattern, by increasing the angle $\angle(\t u,\t v)$, reduces bending energy. The unfolding flattens the pillar and reduces both curvatures $\kappa_1$ and $\kappa_2$. That trend continues up to a point where energy is minimum (near state $(iv)$) and beyond which any further unfolding, or flattening, of the pattern actually causes a dramatic increase in bending energy. This suggests that bending energy $\psi$ depends on some other deformation measures besides curvatures. As a matter of fact, the bending energy density $b$, where $\psi=\int_\Omega b\mathrm{d}\bm\xi$, $\Omega$ being the reference domain of the pattern, is a quadratic form of the planarity defects $\delta^{(ij)}$ of the panels. For their part, the defects are linear forms of the full, both out-of-plane and in-plane, components of the parametrization's second derivatives $\t x_{\mu\nu}$. The out-of-plane components $\Gamma_{3\mu\nu}\equiv\scalar{\t n,\t x_{\mu\nu}}$ are the coefficients of the second fundamental form; they quantify the curvatures and torsion of the embraced surface. The in-plane components $\Gamma_{\sigma\mu\nu}\equiv\scalar{\t x_\sigma,\t x_{\mu\nu}}$ are the Christoffel symbols; they quantify, not the strains, but their gradient, i.e., terms of the form $g_{\mu\nu,\sigma}$. Accordingly, $b=b(\bm\Gamma;\t g)$ is a metric-dependent quadratic form of the curvatures, torsion, and strain gradient. Detailed derivations of the expression of $b$ leading to Fig.~\ref{fig:hyperboloids}d can be found in the Supplemental Material \cite[App. F, G]{note}.

Back to the pillars, unfolding the pattern reduces both curvatures and strain-gradients. However, as the geometric rigidities associated with the curvatures (i.e., $\partial^2b/\partial\Gamma_{3\mu\nu}$) remain bounded, the ones associated with the strain gradient (i.e., $\partial^2b/\partial\Gamma_{\sigma\mu\nu}$) diverge for states that are close to being flat. Such states occur near the outer rims of sufficiently unfolded pillars; their presence further signals that the pillar has reached the boundary of the kinematically admissible domain, hence the energy blow-up observed on Fig.~\ref{fig:hyperboloids}d. Strain-gradient energy further dominates the response of any plane, non-uniform, state such as the ``ring'' observed midway through Fig.~\ref{fig:hyperboloids}b. Indeed, uniformly folded plane states do not engage panel bending and therefore have zero energy. By contrast, non-uniformly folded states, with gradients of folding angles, cannot be achieved without panel bending. When the state is plane, it has zero curvatures and torsion, and energy becomes function of the strain gradient alone \cite{Seppecher2011}. More generally, it is noteworthy that $\psi$ does not penalize strains or fold angles, however large, so long as they are uniform. In fact, $\psi$ does not even refer to a specific natural state or any specific natural fold angles in reference to which strains should be measured \cite{ciarlet2006}. Instead, it refers to a higher-order strain measure, namely $\bm\Gamma$, which quantifies not how much the pattern is folded, but rather how far it is from being uniformly folded.

\section{Conclusion}
In conclusion, we proposed a continuum theory of origami-folded shells where the crease pattern is a zigzag sum. The theory accounts for geometric nonlinearities and accurately predicts the states that are accessible under isometric folding thanks to a universal identity relating in-plane deformations and out-of-plane curvatures. The study further demonstrates that the energy density of origami-folded shells can heavily depend on the in-plane strain-gradient in a way that is unparalleled in classical shell theories. \highlight{It should be possible to extend the proposed theory, in one form or another, to more general tessellations that possess periodic folding motions. By contrast, alternative approaches are likely to be needed to tackle tessellations that do not fold periodically (e.g., Huffman grids and Yoshimura-like patterns  \cite{Evans2015a})}. Finally, it would be of interest to investigate how the proposed theory can inform shape programming and morphing planning of origami and kirigami structures \cite{tachi2013, Dudte2016a, Guseinov2020, Siefert2020, Jin2020, Feng2020, Li2021,  Yang2021, Melancon2021a} for robotics and 4D printing applications.

\section*{Acknowledgements}
HN acknowledges support by the NSF under CAREER award No. CMMI-2045881. AL acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant No. ANR-17-CE08-0039 (project ArchiMatHOS). EW acknowledges support by the NASA Missouri Space Grant Consortium.