Branches of triangulated origami near the unfolded state

TL;DR

This paper develops a second-order rigidity theory for triangulated origami near the unfolded state, revealing multiple folding pathways and the likelihood of misfolding, which impacts self-folding design strategies.

Contribution

It introduces a second-order infinitesimal rigidity framework for unfolded triangulated origami and analyzes the complex branching of nearly-unfolded configurations.

Findings

Number of folding branches grows exponentially with vertices.

Branches are characterized by internal vertices popping up or down.

High pathway multiplicity suggests difficulty in avoiding misfolded states.

Abstract

Origami structures are characterized by a network of folds and vertices joining unbendable plates. For applications to mechanical design and self-folding structures, it is essential to understand the interplay between the set of folds in the unfolded origami and the possible 3D folded configurations. When deforming a structure that has been folded, one can often linearize the geometric constraints, but the degeneracy of the unfolded state makes a linear approach impossible there. We derive a theory for the second-order infinitesimal rigidity of an initially unfolded triangulated origami structure and use it to study the set of nearly-unfolded configurations of origami with four boundary vertices. We find that locally, this set consists of a number of distinct "branches" which intersect at the unfolded state, and that the number of these branches is exponential in the number of vertices. We find numerical and analytical evidence that suggests that the branches are characterized by choosing each internal vertex to either "pop up" or "pop down". The large number of pathways along which one can fold an initially-unfolded origami structure strongly indicate that a generic structure is likely to become trapped in a "misfolded" state. Thus, new techniques for creating self-folding origami are likely necessary; controlling the popping state of the vertices may be one possibility.

Tables

Table 1. Table 1: Summary of random triangulation computations: V i subscript 𝑉 𝑖 V_{i} is the number of internal vertices of the triangulation, the column “systems with V t > 0 subscript 𝑉 𝑡 0 V_{t}>0 ” gives the number of triangulations found that include triangulated triangles (in all cases these were flat octahedra), the column “MV duplicates” gives the number of triangulations that included at least one pair of branches with the same mountain and valley fold assignments. In all cases, the number of branches with multiplicity was 2 V i superscript 2 subscript 𝑉 𝑖 2^{V_{i}} , the number of distinct branches was 2 V i − V t superscript 2 subscript 𝑉 𝑖 subscript 𝑉 𝑡 2^{V_{i}-V_{t}} and the branches were in 2 to 1 correspondence with their vertex sign patterns (defined in Section III.2 ).

$V_i$	triangulations generated	precision used	systems with $V_t>0$	MV duplicates
2	100	500	0	0
3	5000	690	0	0
4	1000	690	2	0
5	1000	690	1	17
6	1000	690	1	28
7	300	690	0	17
8	50	690	0	8

Equations

$$|{\mathbf{X}}_n-{\mathbf{X}}_m{|}^2-L_{(n,m)}^2=0,$$

$$2\frac{{\mathbf{U}}_n-{\mathbf{U}}_m}{\left|{\mathbf{U}}_n-{\mathbf{U}}_m\right|}\cdot \left({\mathbf{u}}_n-{\mathbf{u}}_m\right)+\frac{{\left(h_n-h_m\right)}^2}{\left|{\mathbf{U}}_n-{\mathbf{U}}_m\right|}\approx 0.$$

$$\left[\mathbf{C}\left(\begin{array}[]{c}\mathbf{u}_{1}\\ \vdots\end{array}\right)\right]_{(n,m)}=\frac{\mathbf{U}_{n}-\mathbf{U}_{m}}{\left|\mathbf{U}_{n}-\mathbf{U}_{m}\right|}\cdot\left(\mathbf{u}_{n}-\mathbf{u}_{m}\right).$$

$$2[\mathbf{C}\cdot \mathbf{u}{]}_{(n,m)}+\frac{{\left(h_n-h_m\right)}^2}{\left|{\mathbf{U}}_n-{\mathbf{U}}_m\right|}\approx 0.$$

$$\sum _{(n,m)}\frac{{\sigma }_{(n,m)}}{|{\mathbf{U}}_n-{\mathbf{U}}_m|}{\left(h_n-h_m\right)}^2=0,$$

$$0$$

$${\tilde{\sigma }}_{I,I+1}$$

$${\tilde{\sigma }}_I$$

$$\sum _{n,m}{Q}_{nm}{h}_n{h}_m=0.$$

$$D=V_e-3.$$

$$B_{distinct}=2^{V_i-V_t}.$$

$$E[x]=\frac{1}{2}\sum _{(n,m)}{k}_{(n,m)}({\theta }_{(n,m)}-{\bar{\theta }}_{(n,m)}{)}^2,$$

$$\sum _{i=1}^{N_F}\frac{\partial E}{\partial {\theta }_i}\cdot \Delta {\theta }_i\ge 0.$$

$${\psi }_i=\frac{\pi }{2}-\frac{h_i-h_0}{L_i}.$$

$$\cos {\alpha }_{i,i+1}=\cos {\psi }_i\cos {\psi }_{i+1}+\sin {\psi }_i\sin {\psi }_{i+1}\cos {\beta }_{i,i+1}.$$

$$\delta {\beta }_{i,i+1}=\frac{2\delta {\psi }_i\delta {\psi }_{i+1}-\delta {\psi }_i^2\cos {\alpha }_{i,i+1}-\delta {\psi }_{i+1}^2\cos {\alpha }_{i,i+1}}{2\sin {\alpha }_{i,i+1}}.$$

$$\sum _i\delta {\beta }_{i,i+1}=K.$$

$$\delta {\beta }_{i,i+1}$$

$$0$$

$${\sigma }_{i,i+1}=-\csc {\alpha }_{i,i+1}\frac{\Delta L_{i,i+1}}{L_i{L}_{i+1}}$$

$${\sigma }_i$$

$${\sigma }_{i,i+1}=-\csc {\eta }_{i,i+1}/L_i,$$

$${\sigma }_{i-1,i}=-\csc {\eta }_{i-1,i}/L_i,$$

$$F_{rim,\parallel }$$

$$=F_{rim,\parallel ,-}+F_{rim,\parallel ,+}$$

$$F_{rim,\parallel ,-}$$

$$F_{rim,\parallel ,+}$$

$$\csc {\alpha }_{i-1,i}/L_{i-1}-\cot {\alpha }_{i-1,i}/L_i.$$

$$\sum _{j=1}^3{\sigma }_j{\mathbf{u}}_j=0.$$

$${\sigma }_1/\sin {\alpha }_{2,3}={\sigma }_2/\sin {\alpha }_{3,1}={\sigma }_3/\sin {\alpha }_{1,2}$$

$$\frac{{\sigma }_1}{\sin {\alpha }_{2,3}}=$$

$$-\frac{\cot {\alpha }_{1,2}}{L_1}-\frac{\cot {\alpha }_{3,1}}{L_1})$$

$$=$$

Full text

I Introduction

The development of responsive materials has paved the way to the fabrication of self-folding structures, based on origami, in which flat sheets of a material can be folded along a discrete network of creases into a targeted three-dimensional configuration Silverberg et al. (2015); Waitukaitis et al. (2015); Na et al. (2015); Sussman et al. (2015); Ryu et al. (2012); Tolley et al. (2014); Liu et al. (2014); Tachi and Hull (2016). The creases and vertices formed at their junctions together form a kind of geometric “program” which determines the shape from the strong constraints on how a flat sheet can fold into space. This attractive design paradigm suggests the use of origami as the foundation for mechanical metamaterials Wei et al. (2013); Silverberg et al. (2014); Tachi (2010); Evans et al. (2015); Lv et al. (2014); Schenk and Guest (2013) and deployable structures Miura (1980); Pellegrino (2001). Yet, flexibility is both a blessing and a curse: a single origami crease pattern can admit many different folding pathways Hull (2002); Ballinger et al. (2015); Abel et al. (2016) and, indeed, manipulating a nearly-unfolded origami structure with one’s hands (e.g. the “map-folding problem”) illustrates the competition between pathways that can impede folding to a specific desired configuration Silverberg et al. (2014); Waitukaitis et al. (2015); Tachi and Hull (2016); Brunck et al. (2016). Furthermore, experiments on self-folding gel origami do not always fold into the expected, programmed shape Silverberg et al. (2015).

Let us now fix terminology so that we can discuss these issues more precisely. An origami structure refers to a system of rigid flat plates joined pairwise by ideal hinges, or creases. Origami structures considered here will always arise from a plane polygon decorated with a network formed by the creases and their junctions at vertices, which we call the crease pattern. Origami structures can take on a variety of configurations in 3D space, which are uniquely specified by the positions of their vertices.

To better understand the phenomena of multiple folding pathways and misfolding, it is useful to distinguish two notions of floppiness in an origami structure: (1) the number of degrees of freedom, $D$, which is the dimensionality of the space of motions and scales with the number of boundary sides of a generic origami crease pattern Tachi (2010); Chen et al. (2016); Santangelo (2017), and (2) the number of distinct branches $B$, or folding pathways. As we discuss later, the flat unfolded configuration of an origami structure is a singularity in the space of origami configurations where $B$ branches of dimension $D$ intersect. Consider the triangulated origami structure in Fig. 1 as an example; the figure shows some of the allowed configurations as a function of three of the fold angles. In this example, $D=1$ so the configuration space is locally curve-like almost everywhere; yet, $B=4$, which can be seen as four distinct one-dimensional branches intersecting at the unfolded configuration in the center.

This paper addresses the following questions: how many distinct branches does a generic, triangulated origami crease pattern have, and how can those branches be distinguished? We focus on triangulated crease patterns because triangulated structures are marginally rigid while being maximally flexible Chen et al. (2016), and because triangulated origami can encode the kinematics of origami with bendable faces more generally Silverberg et al. (2014).

To answer these questions, we first show that the small deformations around any initially unfolded configuration of a triangulated origami structure can be described by the simultaneous solutions of a system of $V_{i}$ quadratic equations (Section II) in the $V_{i}+V_{e}$ vertical displacements at the vertices, where $V_{i},V_{e}$ are the number of internal and external vertices, respectively. We will give two interpretations for each of these equations, one coming from statics, showing that we have one equation for each self stress in the system, and one coming from kinematics, showing that each equation enforces the vanishing of Gaussian curvature at an internal vertex. We use this formalism to review the geometry of nearly-unfolded $n$-fold single-vertex origami structures and give a new proof of the fact that their configuration spaces look like $(n-3)$-dimensional double cones Kapovich and Millson (1997), where the two nappes are distinguished by whether the vertex “pops” up or down Abel et al. (2016); Tachi and Hull (2016).

Moving on to the case of triangulated origami with multiple vertices, we restrict our attention in this paper to triangulated origami with four boundary vertices, where the number of degrees of freedom $D=1$. We provide in Section III numerical evidence from a model of random origami squares that the number of branches, $B$, is generically $2^{V_{i}}$, and they are with high probability all distinct. We find a small number of exceptions (appearing with frequency $\sim 1/1000$) that can all be identified from the crease pattern. The branches are not necessarily distinguished by the mountain and valley assignments of the folds, that is, which folds have dihedral angle larger or smaller than $\pi$ respectively. However, we find that pairs of branches appear to be in one-to-one correspondence with pairs of vertex sign patterns, which are assignments of $\pm 1$ to each internal vertex specifying their popping state.

In Section IV we show that a special class of triangulations (roughly, those formed from a sequence of adding degree-3 vertices to the boundary) do satisfy $B=2^{V_{i}}$. We conclude with a discussion in Section V with a discussion on the implications of our results. In particular, we note that the exponential number of branches of a generic origami crease pattern interferes with typical designs for self-folding origami, thus requiring a deeper understanding of how to engineer the origami configuration space topology. Specifically, our results suggest that methods for controlling the popping state of vertices should be investigated.

II Analytical Methods

II.1 Model and second-order deformations

Our kinematic model for origami consists of a triangulated network of springs joining vertices in two dimensions which can, upon deformation, come out of the plane (Fig. 2). We will only consider networks that are planar triangulations of disks (polygons). The edges that are in the interior of the disk will be called folds, since we think of the network as a representation of an origami crease pattern, and since they separate pairs of triangles whose relative orientations differ by some dihedral angle at that edge. We refer to configurations of the origami structure as folded if not all of the dihedral angles between adjacent faces are equal to $\pi$, and unfolded otherwise. Because the network is made from triangles, the angles at vertices of faces between adjacent edges will be preserved as long as lengths are preserved. Furthermore, the triangular faces cannot bend, making this a good model for rigid origami. We label the vertices of the origami structure by an integer, and label edges by a pair $(n,m)$ when the fold joins vertex $n$ to vertex $m$. Using this notation, the kinematic constraints are given by the equations,

[TABLE]

for each pair of vertices, $(n,m)$, joined by a spring of equilibrium length $L_{(n,m)}$. These equations define the configuration space $\mathcal{C}$ of the origami. Note that while self-intersections are allowed, we will work in neighborhoods of $\mathcal{C}$ consisting of configurations which do not self-intersect.

Fig. 1 shows numerical solutions to these length equations for a simple crease pattern. We have chosen to draw $\mathcal{C}$ in this figure using fold angles as coordinates here since those variables are more intrinsic (for instance, they do not see the overall position and orientation of the structure). Near the unfolded state, which we take as the origin, one can linearize the fold angles as functions of the displacements. Therefore the shape of $\mathcal{C}$ near the origin looks the same (up to this linear map) in either set of coordinates. This representation is also useful when thinking about the self-folding paradigm, which we will return to in the concluding section.

The two notions of floppiness described in the introduction have natural interpretations in terms of the geometry of $\mathcal{C}$. The number of degrees of freedom $D$ is the dimension of the configuration space. Note that the configuration space may have singularities (and in this work, this is the case of particular interest), so the notion of “dimension” becomes subtle (Harris (1992), Lecture 11). For our purposes, it will suffice to say that the dimension of the configuration space is the dimension at any nonsingular point.

The number of branches $B$ is a property of a singular point of $\mathcal{C}$. For instance, at a singular point consisting of the intersection of multiple distinct curves or surfaces, each one of those would consist of a branch. A general definition of a branch (as an irreducible component of the analytic germ at the singularity) would take us a bit too far afield into singularity theory Greuel et al. (2007); we give a more concrete definition for our case later in Section II.3.

We are interested in deformations of unfolded origami, where the vertices all lie in a single plane and the faces do not overlap. Without loss of generality, we will assume the initial unfolded configuration lies in the $xy-$plane. We write the position of vertex $n$ as $\mathbf{X}_{n}=\mathbf{U}_{n}+\mathbf{u}_{n}+h_{n}\mathbf{\hat{\mathbf{z}}}$, where $\mathbf{U}_{n}$ are the equilibrium positions of the vertex in the $xy-$plane, $\mathbf{u}_{n}$ is a vertex displacement in the $xy-$plane and $h_{n}$ is a vertical displacement out of the plane. Expanding Eq. (1) to lowest order in the displacements yields

[TABLE]

Because the vertical displacement decouples from the in-plane displacement and the linear terms in $h_{n}$ vanish, any displacement with $\mathbf{u}_{n}=0$ for all $n$ preserves lengths to first order (i.e. any displacement consisting only of height changes is a first-order flex or motion). But by stopping here we have not captured enough information to see the branches, as the lowest-order information lies in the quadratic terms of Eq. (2). Since the term quadratic in height leads to a change in bond lengths of the same order as the linear term in the in-plane displacements, we can safely neglect terms of order $\mathcal{O}(\mathbf{u}_{n}^{2})$.

The first term of Eq. (2) governs the infinitesimal displacements of the in-plane degrees of freedom. We can rewrite this expression by concatenating the in-plane displacements into a vector $(\mathbf{u}_{1},\mathbf{u}_{2},\cdots)$, and define an in-plane compatibility matrix such that row $(n,m)$ of $\mathbf{C}$ is defined by the equation

[TABLE]

The matrix $\mathbf{C}$ has pairs of columns indexed by vertices $n$ and rows indexed by the unique folds $(n,m)$. Formally, it maps vectors of in-plane deformations to vectors of in-plane spring displacements and governs the linear deformations of the unfolded configuration of the origami structure that keep it in the $xy$-plane. Using $\mathbf{C}$, Eq. (2) becomes

[TABLE]

Since the in-plane deformations and out-of-plane deformations are decoupled in Eq. (4), the in-plane motions are governed (at this order) by the kernel of $\mathbf{C}$. And because 2D triangulated networks are generically rigid in the plane, by a counting argument and Laman’s theorem Graver et al. (1993), $\textrm{ker}~{}\mathbf{C}$ is generically generated by translations and rotations only. As these are not of much interest to us here, from now on, we will only consider how Eq. (4) constrains out-of-plane deformations.

To extract constraints on the vertical displacements $h_{n}$ from Eq. (4), we make use of the self stresses of the network. These are row redundancies in $\mathbf{C}$, defined by $\bm{\sigma}^{T}\cdot\mathbf{C}=\mathbf{0}^{T}$, where $\bm{\sigma}$ is a stress vector with one component per fold and T is the transpose 111Note that unlike the dimension of the space of linear deformations, the space of self stresses does not depend on whether the unfolded configuration is considered to be embedded in an ambient 2D or 3D space.. Taking the dot product of $\bm{\sigma}$ with both sides of Eq. (2), we obtain

[TABLE]

where $\sigma_{(n,m)}$ is the component of $\bm{\sigma}$ along the fold $(n,m)$ and the sum is taken over all folds in the network. Thus, each self stress of the unfolded configuration of the origami gives an equation (5) that constrains the vertical displacements $h_{n}$ to second order. This is a special case of the more general formalism of Connelly and Whiteley Connelly and Whiteley (1996) who also show that all second-order constraints are generated by the self stresses (that is, these conditions which are necessary for deformations to preserve lengths to second-order are also sufficient).

In most of the rest of this paper, we will be considering solutions to the system of equations coming from applying Eq. (5) to an independent basis of self stresses of triangulated crease patterns. This system defines a certain subspace in the space of all vertical displacements, whose dimension we will discuss in Section II.3. By considering only vertical displacements as our variables we are implicitly removing the in-plane degrees of freedom, and in particular, the in-plane translations and rotation about the $z$-axis will not come into play.

The solutions $h_{n}$ to Eq. (5) only give an approximation to the configuration space at the unfolded state, and may not be a wholly accurate picture, even qualitatively (Harris (1992), Lecture 20). The issue is that while the non-existence of folded solutions at second-order proves that there can be no folded configurations (second-order rigidity implies rigidity) Connelly and Whiteley (1996); Connelly (1980), a second-order solution may not be a solution at all orders. In all cases that we checked, e.g. in making Fig. 1 and the Supplementary Movies, solutions to Eq. (5) did seem to correspond to true solutions of Eq. (1) (see Section III). And for single-vertex origami, the second-order deformations from the unfolded state can be shown to extend to actual rigid motions Kapovich and Millson (1997); Demaine et al. (2016), see also Appendix C). We also discuss in Section IV some special multiple-vertex crease patterns where second-order motions also extend to true motions. However, we know of no such general guarantee, and there are (non-triangulated) origami examples where second-order solutions do not correspond to points lying in the true configuration space Demaine et al. (2016). Nonetheless, the displacements allowed here can only change the stretching energy to at most sixth-order, so we will mostly ignore this issue in what follows and simply refer to our second-order approximations as configurations. Because our analysis is nonlinear, we will end up relying on several more such hypotheses which we have chosen to deal with post hoc by checking that they are satisfied in our numerics, rather than seeking a rigorous proof here.

II.2 Wheel stresses, Gaussian curvature and single-vertex origami

Our first observation is that an unfolded origami structure with $V_{i}$ internal vertices has at least $V_{i}$ self stresses. To construct them, we first isolate the faces around each internal vertex and consider the mechanics of the isolated vertex stars apart from the rest of the origami structure (Fig. 3a). Those faces and their edges make a spoked wheel of folds emerging from a single internal vertex and meeting the vertices of a polygon. If there are $N$ spokes, the $2N+2$ in-plane positions of the vertices are subject to $2N$ constraints. Since there are three planar Euclidean motions, there must be generically one self stress (Fig. 3b). This wheel stress is preserved if we embed it into the larger structure by setting the remaining components of $\sigma_{(n,m)}$ to zero.

Interestingly, the second-order constraints arising from using the wheel stresses in Eq. (5) also have a natural geometric interpretation in terms of the Gaussian curvature at each internal vertex, measured by the sum of the angles between adjacent folds around it. To see this, note that we can also generate constraints by enforcing the condition that the Gaussian curvature remains zero at each internal vertex after deformation. We label the folds around each internal vertex with an index $I$ which we take modulo the number of folds meeting at the vertex. Then let $\alpha_{I,I+1}$ be the planar angle between folds $I,I+1$ and let $\psi_{I}$ be the angle the $I^{th}$ fold makes with respect to the $xy-$plane. Spherical trigonometry yields the constraint

[TABLE]

valid up to quadratic order in the $\psi_{I}$.

To lowest order, we have $\psi_{(n,m)}=(h_{n}-h_{m})/|\mathbf{U}_{n}-\mathbf{U}_{m}|$ in terms of the height displacements. We thus have an equation for each internal vertex in the form of Eq. (5) from which we can read off a candidate self stress $\sigma_{(n,m)}$. If we denote the self stress on the outside edges (on the rim of the wheel) of each vertex by $\tilde{\sigma}_{I,I+1}$ and the self stress on the spokes $I$ by $\tilde{\sigma}_{I}$ (Fig. 3b), then we obtain

[TABLE]

where $L_{I}$ is the length of fold $I$ and $\Delta L_{I,I+1}^{2}=L_{I}^{2}+L_{I+1}^{2}-2L_{I}L_{I+1}\cos\alpha_{I,I+1}$ (Appendix A). In Appendix B, we prove that Eqs. (II.2) do give the coefficients of a self stress and thus that the wheel stress constraint is precisely Gaussian curvature preservation.

We now discuss the configuration space of single-vertex origami structures. Our analysis here will play a big role in our later treatment of multiple-vertex structures. Let $Q_{nm}$ be the symmetric matrix corresponding to $\bm{\sigma}$ so that the quadratic form in Eq. (5) is written in terms of the vector of vertical displacements at vertices $\mathbf{h}$,

[TABLE]

If the vertex associated with $Q_{nm}$ has $N$ folds, we find that the $(N+1)\times(N+1)$ matrix $Q_{nm}$ has $N-2$ nonzero eigenvalues, exactly one of which is negative (Ref. Kapovich and Millson, 1997 and Appendix C). This means that a single $N-$fold vertex has $N-2$ first-order motions once translations and rotations have been removed. The second-order motions are the solutions to Eq. (8), which defines a $(N-3)$-dimensional surface in the linear space of first-order motions, namely, the null-cone of this quadratic form. Since $Q_{nm}$ has one negative eigenvalue, its null-cone has two conical components (nappes) that meet at the unfolded state (Fig. 4a). Topologically, the nappes are cones over $(N-4)$-dimensional spheres.

Let $\mathbf{h}_{-}$ be the eigenvector corresponding to this negative eigenvalue. Points on the two nappes can be distinguished by the sign of their dot products with $\mathbf{h}_{-}$ since the plane normal to $\mathbf{h}_{-}$ separates the nappes. This eigenvector gives a set of displacements that maximizes the change in the Gaussian curvature. Indeed, it can be shown from the formulae given in Appendix C that the component of largest magnitude in this eigenvector is the displacement at the vertex itself and the neighboring vertices are moved by smaller amounts in the opposite direction. Such a displacement (Fig. 4b) leads to a conical deformation at the vertex. This suggests that the difference between rigid origami configurations in the two nappes is related to whether the vertex is buckled up or down (relative to the upwards normal of the origami sheet). We make this more precise in the rest of this subsection.

The trace of an origami vertex is defined to be the spherical polygon obtained by intersecting the origami with a small sphere centered at the vertex; it is non-self-intersecting for a vertex sufficiently close to being unfolded and thus cuts the sphere into two pieces, corresponding to the upper and lower sides of the origami sheet. Since the trace of a folded origami vertex lies completely in an open hemisphere Streinu and Whiteley (2005), one of those pieces will have area less than $2\pi$. If that piece corresponds to the upper side of the origami, the configuration is called popped down (as the vertex “points” towards the lower side of the sheet) and otherwise it is called popped up Abel et al. (2016) (Fig. 5). Since configurations of these two types meet only in the unfolded state, one of the nappes of the double cone configuration space consists of popped up configurations and the other consists of popped down configurations, so we can use dot products with $\mathbf{h}_{-}$ to distinguish them computationally. Note that Ref. Abel et al. (2016) give a simpler definition, where popped down (up) vertices are those whose edge vectors are all in the northern (southern) hemisphere. Our definition is a rotation-invariant generalization that is better suited for considering configurations at vertices that are part of larger multiple-vertex origami.

II.3 Consequences for multiple-vertex origami and the definition of branches

Eqs. (5) provide a way to count the number of infinitesimal degrees of freedom of an arbitrary triangulated origami structure (subject to the caveats described at the end of Section II.1), and provide information on the number of distinct ways of folding a given crease pattern from the unfolded state.

Suppose we have an unfolded origami structure with $V_{i}$ internal vertices and $V_{e}$ boundary vertices. There are $V_{i}+V_{e}$ linear degrees of freedom corresponding to vertical displacements, but $V_{i}$ quadratic equations constraining them. We are assuming here that all folds are incident to at least one internal vertex (i.e. we cannot disconnect the crease pattern by cutting along any one fold as in Fig. 6). Therefore, the number of nontrivial degrees of freedom for a generic unfolded triangulation, should be given by

[TABLE]

The term $3$ arises from removing the three remaining out-of-plane Euclidean motions (this can be done by e.g. pinning the vertices of an arbitrary triangle to the $xy-$plane). Eq. (9) recovers the count for the degrees of freedom for a generic (folded) triangulated origami derived from a linear analysis Chen et al. (2016). The unfolded configuration admits $V_{i}+V_{e}-3=D+V_{i}$ nontrivial linear motions, so this linear analysis fails, though we see that counting quadratic constraints as we do here leads to the expected number $D$. Geometrically, the unfolded configuration is a singular point of the configuration space, where the dimension of the tangent space exceeds the dimension at other nearby points.

More precisely, $D$ should be the local dimension of the configuration space at nonsingular points. However, since the constraint equations are nonlinear, our derivation of Eq. (9) is not rigorous. To give a proof, we would also need to show that the hypersurfaces arising from Eq. (5) for each of the $V_{i}$ internal vertices intersect transversely (Harris (1992), Example 11.8). This should be true generically, and in all of the numerical examples considered in this paper this is indeed the case.

Before turning to our discussion of distinct folding pathways and branches, we consider first the case $V_{e}=3$. Such crease patterns are triangulated triangles (see the upper left triangle and its red interior folds in Fig. 6). Triangulated triangles are equivalent to planar projections (Schlegel diagrams) of triangulations of spheres, e.g. a degree-3 vertex lies in the center of a projected tetrahedron. Our count $D=0$ suggests that these should have 0 degrees of freedom, i.e. that they should be rigid. Gluck proved that generic triangulated spheres are rigid Gluck (1975). In our case, the triangulated triangles are flat and thus nongeneric, but Connelly proved that these are rigid at second-order in 3D as well Connelly (1980). Thus, in this case, the configuration space is simply a point at the unfolded state (with multiplicity, as we will see later).

To better understand the neighborhood of the unfolded state and to define the notion of “branches” when $D\geq 1$, we consider the solutions of Eqs. (5), which are sets of vectors $\mathbf{h}$ in $\mathbb{R}^{V_{i}+D}$. Since the quadratic equations are all homogeneous in the vertex heights, any vector $\lambda\mathbf{h}$ solves Eqs. (5) if $\mathbf{h}$ does, for any real number $\lambda$.

Therefore, we are led to consider solutions of Eqs. (5) in projective space $\mathbb{RP}^{V_{i}+D-1}$, where a height vector $h$ is identified with $\lambda h$ for any nonzero real number $\lambda$. Let $B$ be the number of connected components of the solution set (counted with multiplicity) in $\mathbb{RP}^{V_{i}+D-1}$ and let us assume that $B>0$. Roughly speaking each of these connected components is generically a $(D-1)$-dimensional component of the intersection of a small sphere centered at the unfolded state with the origami configuration space, where components related by the $z$-reflection symmetry $h\mapsto-h$ are identified. In this paper we will consider mostly the case where $D=1$ where these components are simply points. Back in the space of first-order deformations $\mathbb{R}^{V_{i}+D}$, each of these components induces a double cone over some reflection-related pair of components on this sphere, and all of the origins of these cones (singular when $D>1$) intersect at the unfolded state. We will refer to these $B$ cones as branches. When $D=1$, these cones are simply lines.

Note that by our discussion in Section II.2, single-vertex origami structures have a single branch. For the example in Fig. 4a, where $V_{i}=1,D=2$, the branch is a double cone over a closed curve in $\mathbb{RP}^{2}$.

III Numerical results

III.1 Counting branches

We now discuss the case of $D=1$ in more concrete terms. In that case, each of the branches corresponds to a curve that passes through the unfolded state, as in Fig. 1. To understand the typical number of branches $B$, we consider a class of random, one-degree-of-freedom origami, generalizing the example of Fig. 1. We generated random origami structures by computing a Delaunay triangulation on the point set consisting of $V_{i}$ random points uniformly distributed within a square, together with the corners of the square. We omitted configurations having 3-fold vertices since such vertices are always rigid. We also found no triangulations with edges connecting opposite corners of the square (such as in Fig. 6). Such edges would break the square into two rigid triangulated triangles.

Since there is a vertex at each corner of the square, $V_{e}=4$ and so Eq. (9) yields $D=1$, no matter how the interior vertices and edges are arranged. We fix the overall position and orientation by setting the height of the vertices of one triangle to zero. To remove the scaling symmetry from the homogeneous system coming from Eqs. (5), we add a normalizing equation $\sum_{j}h_{j}^{2}=1$, resulting in $V_{i}+1$ quadratic polynomials that must be simultaneously solved in terms of $V_{i}+1$ vertex heights. We solved these systems in Mathematica 11, which uses a homotopy continuation algorithm for numerical root finding of polynomial systems Morgan and Sommese (1987). Note that the solutions of these polynomial systems come in pairs related by multiplication by -1 corresponding to $z$-reflection symmetry (as discussed above, each branch is a line through the origin and these intersect the normalizing unit sphere twice), so the number of branches $B$ is half the number of real solutions.

We generated several thousand random triangulated origami squares with $V_{i}=2$ to 8 and computed and analyzed their branches (Table 1). Accurate solutions seem to require very high precision arithmetic, especially as $V_{i}$ becomes larger; to ensure good results we used up to $690$ digits of precision and verified the resulting solutions, which allowed us to find solutions up to $V_{i}=8$.

Using these second-order branches, we numerically computed approximate configurations (solutions to the length equations in Eq. (1)) to create our figures and movies by using the solutions to perturb the unfolded configuration and minimizing a stretching energy (sum of squared differences of edge lengths to the lengths in the unfolded state) on the coordinates until it was zero to high accuracy.

We now give a little bit of background on systems of polynomial equations to give context for our main results. In general, there is little one can say about the simultaneous roots of a completely arbitrary system of polynomials, especially if one is interested in real roots. One result, known as Bézout’s theorem, states that if the solutions to a system of polynomial equations are isolated points, then the number of complex solutions, counted with multiplicity and including points “at infinity”, is equal to the product of the degrees of the equations (Harris (1992), Lecture 18). For our systems of $V_{i}+1$ quadratic polynomials in $V_{i}+1$ height variables, this yields $2^{V_{i}+1}$ roots. However, we are only interested in real and finite roots. Our systems have the property that all coefficients are real, but this merely guarantees that nonreal roots come in complex-conjugate pairs; there might still be no real roots.

Now we state the first of our main numerical findings: with only a few exceptions we will discuss shortly, all the solutions of our system are distinct and isolated and amazingly, all $2^{V_{i}+1}$ of them are real! Since the branches correspond to $\pm$ pairs of solutions, we therefore have $B=2^{V_{i}}$.

As mentioned, not all crease patterns lead to $2^{V_{i}}$ distinct branches. We found a very small number of systems ($2$ with $V_{i}=4$, $1$ with $V_{i}=5$ and $1$ with $V_{i}=6$, see Table 1) with fewer branches; in these cases all branches came with some multiplicity. In these cases, we were able to identify triangulated triangles within the crease pattern. Let $V_{t}$ be the number of vertices in the interior of all such triangulated triangles (if we had not excluded crease patterns with degree-3 vertices from our computations, we would count these in $V_{t}$). These vertices must remain unfolded in all branches of configurations, and it follows that their height variables do not contribute to distinct roots of the quadratic equations, but rather only give rise to multiplicity. We found that in such cases, the distinct roots of the quadratic equations had multiplicity $2^{V_{t}}$. Thus we may account for this effect with the (conjectural) formula

[TABLE]

Let us call the $V_{i}-V_{t}$ internal vertices which are not interior to triangulated triangles foldable vertices. Note that the heights of non-foldable vertices are determined by the foldable ones by the linear condition that each triangulated triangle is planar, so in the space of linear motions $\mathbb{R}^{V_{i}+1}$ the branches must lie in a lower-dimensional $\mathbb{R}^{V_{i}-V_{t}+1}$.

We now briefly discuss the pattern of mountain and valley folds in the branches. Following a short Euler characteristic argument, the number of folds (internal edges) for a triangulated square is $3(V_{i}-V_{t}+1)$, so the potential number of distinct mountain-valley (MV) assignments up to a global sign change is $2^{3(V_{i}-V_{t}+1)}/2$. Naïvely, the fraction of these MV assignments that we should expect to find among the branches is at most $(1/4)^{V_{i}-V_{t}-1}$, which approaches zero exponentially in the number of vertices. However, there are combinatorial consistency constraints on the MV patterns, along the lines of those derived in Ref. Abel et al. (2016) for single-vertex origami, so this is certainly an overestimate. We did not attempt to work out these consistency conditions, but as some evidence that they play a role, in our computer-generated examples with $V_{i}=5,6,7,8$ we found an increasing number of crease patterns where multiple branches have coinciding MV assignments (Table 1 and Section III.3). This phenomenon is well-known in the origami community Hull and an illuminating example is the 6-fold origami vertex with alternating mountain and valley folds (MVMVMV in cyclic order around the vertex) called the “waterbomb structure” Brunck et al. (2016); Tachi and Hull (2016). Since the distribution of mountain and valley folds do not distinguish different branches in configuration space, and further, since it is hard to guess which MV assignments are allowed, the question remains, is there anything that does distinguish those branches from each other?

III.2 Vertex sign patterns

Given a folded configuration of an origami square with multiple internal vertices, we can ask whether each vertex is popped up or down. This data is encoded as vertex sign patterns, assignments of $+1$ or $-1$, to the foldable vertices if they are popped up or down, respectively. (Since the $V_{t}$ vertices lying within triangulated triangles are always unfolded, we could extend the vertex sign pattern to these vertices by assigning them the value 0). As there are therefore $2^{V_{i}-V_{t}}$ choices of vertex sign patterns, it is natural to hope that there is a one-to-one correspondence with the branches. However, branches consist of pairs of solutions related by the $z$-reflection symmetry, so we must identify vertex sign patterns related by a global sign change and we are left with only $2^{V_{i}-V_{t}-1}$ equivalence classes. Note that in a single 4-fold vertex origami, we have only one vertex sign pattern up to sign and two branches, so instead the best we can hope for is a one-to-one correspondence between pairs of branches and sign-related pairs of vertex sign patterns.

We check this correspondence computationally as follows: we first determine the eigenvector with negative eigenvalue of $Q_{nm}$ for each internal vertex $n$, $\mathbf{e}_{n}$ (shown as an arrow in Fig. 4). Here we ensure that $+\mathbf{e}_{n}$ corresponds to the popping up deformation and we extend the vector with zero components so that its dimension is the same as that of $\mathbf{h}$. Then the vertex sign pattern is defined by $\sigma_{n}=\textrm{sgn}[\mathbf{e}_{n}\cdot\mathbf{h}]$. Both $\mathbf{h}$ and $-\mathbf{h}$ correspond to the same branch, so we associate to each branch a pair of vertex sign patterns related by a global sign change. Our second main numerical finding is the following: remarkably, in all of our computed examples (summarized in Table 1), there are exactly two branches with each such pair of sign patterns, in agreement with our guess above!

If we look at the branches as lines intersecting in the singular unfolded state in the second-order configuration space $\mathbb{R}^{V_{i}-V_{t}+1}$, each of the $\mathbf{e}_{n}$ defines a hyperplane separating the popped-up configurations at $n$ from those popped-down there (Fig. 7). The set of all such planes arising from the $V_{i}-V_{t}$ foldable vertices divides $\mathbb{R}^{V_{i}-V_{t}+1}$ into $2^{V_{i}-V_{t}}$ chambers, each labeled by a different vertex sign pattern. Each of these chambers is topologically the product of an orthant of $\mathbb{R}^{V_{i}-V_{t}}$ with a real line. For instance in Fig. 7, each chamber is topologically the product of a quarter-plane with a line, a “wedge-shaped” region of 3D space (the projection has been specially chosen so that these chambers are seen “edge-on”, otherwise the dividing planes obstruct the view). In these terms, our observation is that the $2^{V_{i}-V_{t}+1}$ distinct rays of the branches always seem to be distributed so that two rays lie in each of these chambers, implying that the solutions to the coupled system of nonlinear equations Eq. (5) are controlled to some extent by what happens at each vertex. Again, from the perspective of random solutions to real quadratic equations, one might have expected some of these branches to be complex conjugate pairs and that the real branches would be distributed much more unevenly in the chambers. Simple tests with random perturbations of the coefficients of our equations (so that they no longer come from realizable crease patterns) confirm this expectation – after perturbing, the solutions of most systems have many complex-conjugate pairs and real solutions are not equidistributed.

III.3 Branches with non-unique Mountain-Valley assignments

We now have two pieces of combinatorial data associated to each branch: first, the well-known assignment of which folds are mountain and which are valleys (MV assignment) and second, the vertex sign pattern (modulo sign). For most of the origami crease patterns we computed (Table 1), the branches all have different MV assignments. However, this is not always the case, as mentioned at the end of Section III.1. When pairs of branches have identical MV assignments, they can usually be distinguished by the vertex sign patterns; in particular there is usually exactly one vertex sign that differs. However, interestingly, we did find a number of examples where the two branches with the MV assignment also have the same vertex sign pattern. In the rest of this section we will show examples of these occurrences.

In Fig. 8 we show a typical origami crease pattern from our data exhibiting two noncongruent branches with coincident MV assignments. The branches in this example (and most of the other “MV-coincident pairs” we found) can be distinguished by the popping state of exactly one vertex. It follows from Corollary 1 of Ref. Abel et al. (2016), that a vertex which can both pop up and pop down must have degree at least 6, and indeed, must contain both a mountain “bird’s foot” and a valley bird’s foot as subsets of the folds around the vertex. Here, a bird’s foot is a sequence of four not-necessarily adjacent folds $c_{1},c_{2},c_{3},c_{4}$ in counter-clockwise order around the vertex such that the angles between $c_{1},c_{2},c_{3}$ are between 0 and $\pi$ and $c_{1},c_{2},c_{3}$ have the same sign (all mountains or all valleys) and $c_{4}$ has the opposite sign. Note that the waterbomb vertex contains both a mountain bird’s foot and a valley bird’s foot.

More interestingly, we found a few examples (1 configuration with $V_{i}=5$, 2 with $V_{i}=6$ and 1 with $V_{i}=7$, out of the configurations computed for Table 1) where MV-coincident branches also had the same vertex sign patterns. The example with $V_{i}=5$ is shown in Fig. 9.

In all the examples we computed where branches had coincident MV assignments, they either had exactly one high-degree vertex whose popping state distinguished the branches, or none. We do not dare venture to guess about the relative frequency of such examples as $V_{i}$ gets large.

IV H1 triangulations and butterflies

In this section we describe a class of $D=1$ triangulations whose configuration spaces are particularly easy to analyze. We will show that near the unfolded state, their configuration spaces consist of $2^{V_{i}}$ intersecting 1D branches. In contrast to most of the rest of our results, we will be discussing actual configurations here, not just second-order approximations.

Henneberg moves Graver et al. (1993); Tay and Whiteley (1985); Nixon and Ross (2014) are rigidity-preserving transformations of graphs that are useful in rigidity theory. We will only consider the Henneberg-1 move, which in three dimensions is just the addition of a new vertex and attaching it to the original graph with three new edges. In particular, we are interested in triangulations that can be built by repeatedly adding such degree-3 vertices to the boundary. To be more precise, suppose we have some triangulation $T$. First choose three vertices $v_{1},v_{2},v_{3}$ on the boundary of $T$ that are adjacent in the cyclic ordering. Then add a new vertex $w$ and join it with new edges to $v_{1},v_{2},v_{3}$. We will say that $T$ and the resulting triangulation $T^{\prime}$ are related by an H1 move and will not refer to general Henneberg-1 moves any more. See Fig. 10(b), where the vertex labeled (1) has been added by an H1 move to the rest of the triangulation; note how (1) and the three edges incident to it form a pair of adjacent red triangles on the boundary of the triangulation. We will require that $v_{2}$ has degree at least 3, as otherwise, it would become a degree-3 vertex in the interior of $T^{\prime}$, and thus be non-foldable.

The key property of H1 moves is that given a configuration of $T$ where the faces around $v_{2}$ are in a folded configuration, there are always two distinct folded configurations of $T^{\prime}$. To see this, we consider the motion of butterflies. A butterfly is the 1-degree-of-freedom rigid origami consisting of two triangular faces joined by a shared fold (Fig. 10a). Its configuration space is topologically a circle, parametrized by the dihedral angle at the shared fold. Consider the distance $d$ between the two non-shared vertices of the butterfly. As the dihedral angle varies from 0 to $2\pi$, the distance $d$ increases monotonically from its minimum value $d_{min}$ in the flat folded state, to its maximum value $d_{max}$ when the butterfly is flat and unfolded, and then decreases monotonically again to $d_{min}$. Thus for each value $d_{min}<d<d_{max}$, there are two distinct configurations of the butterfly that are related by reflection. Now the claim follows since an H1 move can be viewed as gluing a butterfly at $v_{1},v_{2},v_{3}$. This is essentially the 3D version of a construction for graphs embedded in the plane described in Ref. Borcea and Streinu (2004). Note that when the faces around $v_{2}$ are unfolded, then $d$ is maximized and we can only attach the butterfly in its unfolded state.

The basic idea in the rest of this section will be that triangulations that are reducible by reverse H1 moves to a seed triangulation with a simple configuration space can also be understood easily. Unfortunately, it seems that very few triangulations are reducible at all by a reverse H1 move. We generated 20000 triangulations with $V_{i}=3$ to 8 by the method described in Section III and found that an increasing fraction of triangulations had no degree-3 vertices at the boundary (starting from 0% at $V_{i}=3$ to 88.5% at $V_{i}=8$). Furthermore, the number of triangulations in this data set that can be reduced further with more reverse H1 moves appears to decay exponentially.

Nonetheless, we now narrow our focus to triangulations that can be constructed from a sequence of H1 moves from a butterfly (the seed). We will call these H1 triangulations (Fig. 10b). In line with our observations in the last paragraph, we found that the fraction of H1 triangulations decreased roughly exponentially from 1 at $V_{i}=3$ to 0.0077 at $V_{i}=8$. (It’s not hard to check that under our restriction of no interior degree-3 vertices, all $V_{i}=1,2,3$ triangulations are H1).

Note that each reverse H1 move results in the deletion of one boundary vertex and the conversion of one internal vertex to a boundary vertex. This means that a H1 triangulation with $V_{i}$ internal vertices can be decomposed into a sequence of $V_{i}$ butterflies. Also, since the seed butterflies always have four boundary vertices, every H1 triangulation also has four boundary vertices. We will assume in this section that the triangulations have no interior degree-3 vertices and also that there are no collinear folds meeting at vertices (there are no “crosses” in the terminology of Ref. Abel et al. (2016)), as these nongeneric collinearities can rigidify folds in a flat state.

Let us now discuss the configuration space of H1 triangulations near the unfolded state. As discussed above, the seed butterfly has a configuration space that is a circle. Now consider the butterfly that is attached by the first H1 move. We showed that there are two distinct ways of attaching it when the seed is nonflat, and one way of attaching it when the seed is unfolded. We can always choose a small enough interval around the unfolded state in the configuration space of the seed where the dihedral angle of the new butterfly never reaches $\pi$. This implies that a neighborhood of the unfolded state of the configuration space of the two butterflies together consists of two intersecting lines, topologically a letter X, since we have two choices over every nonzero initial dihedral angle, glued together at the unfolded state. Continuing, one sees that for each H1 move, the number of 1D branches doubles, and they all still intersect at the unfolded state. This results in $2^{V_{i}}$ branches, as desired.

We have not been able to show that all vertex sign patterns are realized twice, as this seems to require some careful analysis of the configurations of the boundary vertices and when they result in popped up / down states of interior vertices after an H1 move. However, one can generalize the arguments in this section to H1-like triangulations where the seed is not just a butterfly but some other simple triangulation, e.g. a single-vertex origami. We hope to elaborate on this elsewhere.

V Discussion

Our results have broad importance in using origami techniques to manufacture shapes from flat substrates. In the standard paradigm for self-folding origami, an initially unfolded sheet is “programmed” by setting the equilibrium dihedral angle of each fold to a nonzero value Silverberg et al. (2014); Waitukaitis et al. (2015); Tachi and Hull (2016); Brunck et al. (2016). One illustrative self-folding energy functional takes the form

[TABLE]

where $\theta_{(n,m)}$ is the folding angle of fold $(n,m)$, $\bar{\theta}_{(n,m)}$ is the programmed equilibrium value, $k_{(n,m)}$ is an angular spring constant for the fold, and the sum is taken over all folds $(n,m)$. One can visualize geometrically this functional in Fig. 1 as a generalized squared-distance function from the point $\bar{\theta}$ corresponding to the programmed folding angles. Based on the branched structure of configuration space we have found (as exemplified by Fig. 1), once an origami structure begins folding along the wrong branch, it is potentially very difficult to return it to the desired branch.

This possibility motivated Tachi and Hull, in Ref. Tachi and Hull (2016) to introduce the notion that the bending moments driving a self-folding origami structure should drive the pattern in a direction perpendicular to any undesirable branches in configuration space; i.e. the gradient of the energy functional at the unfolded state should project onto only one branch. Their condition can be justified from the point of view of an analytic energy landscape. In order to avoid mis-folding, the energy should not decrease along any undesired branches. If we assume that the energy of a structure depends only on the angles of the folds, $E(\theta_{1},\theta_{2},\cdots,\theta_{N_{F}})$, then naturally we require that any infinitesimal change in fold angles $\Delta\theta_{i}$ along an undesired configuration branch satisfy

[TABLE]

Since branches are symmetric under $\theta\rightarrow-\theta$, $\partial E/\partial\theta_{i}$ must then be perpendicular to all undesired branches in the $N_{F}$-dimensional space of fold angles.

In order for this to be satisfied easily, one could hope that under some circumstances, an origami structure would have very few branches; however our results suggest that in generic triangulated patterns with $D=1$, the number of branches depends only on $V_{i}$ and not the geometry. If we combine the fact that in a triangulation, the number of internal folds is $N_{F}=3V_{i}+const$ with our finding that the number of branches is $2^{V_{i}}$, we have thus shown that it is impossible to satisfy the above orthogonality constraint for even modest numbers of internal vertices! Nevertheless, it is an interesting open question whether this condition could be satisifed approximately, i.e. whether, for sufficiently large origami crease patterns, we might have most branches within 99% of being orthogonal, as this happens for random vectors in high dimensions Cai et al. (2013). Some preliminary results have been collected in Appendix D which show that in fact branches tend to be less orthogonal in fold-angle space than a uniformly random set of vectors. It remains to be seen how this impacts the programmability of self-folding origami structures.

We also note that our methods can be applied even to non-triangulated origami. In that case, we would first triangulate the origami, then add additional quadratic constraints to rigidify certain folds. The dimension of the configuration space becomes $D=V_{e}-4-E$, where $E$ is the number of rigidified folds in the triangulation that are not in the original pattern. When this formula becomes negative, the constraints in the system become redundant and we instead take $D=0$. We would also like to understand what happens when $D>1$; it may be that when the dimension of the branches increases, that the number of them no longer scales exponentially. In a recent paper, Stern et al Stern et al. (2017) studied examples with quadrilateral faces where $V_{e}$ was large and $D=0$. They studied branches defined to be minimal-energy directions (as opposed to zero-energy, as in this work) away from the unfolded state. Intriguingly, the scaling of their $B$ is exponential in $\sqrt{V_{i}}$. It would be very interesting to understand the crossover between the regime studied in that work and this one.

Proving that each sign-related pair of vertex sign patterns corresponds to a unique pair of folding branches remains elusive. Our problem of finding all branches fits into the broader context of enumeration problems in geometric rigidity such as finding all realizations of isostatic graphs in the plane Borcea and Streinu (2004); Capco et al. (2017) and enumerating all rigid clusters of sticky hard spheres Holmes-Cerfon (2016). The vertex sign patterns seem to impose some structure on the branches in configuration space that leads to the much more tractable formula $B=2^{V_{i}}$ than e.g. the recent recursion formula for the number of complex realizations of isostatic graphs in the plane Capco et al. (2017).

Understanding the relation between vertex sign patterns and the geometry of the branches better may be useful for developing robust self-folding structures. In particular it would be interesting to see whether a self-folding paradigm based not only on preferred dihedral angles but also on vertex popping states could be easier to control experimentally, e.g. with conical indenters above and below the sheet, or perhaps with some kind of actuation or swelling that breaks the up-down symmetry of the sheet near each vertex.

Finally, since MV assignments seem to be more frequently ambiguous, particularly for origami with many internal vertices (see Table 1 where the fraction of crease patterns with branches having duplicate MV assignments seems to be increasing), we suggest that perhaps origami crease patterns for folding paper origami should be given with vertex popping states as well.

Acknowledgements.

We thank the kind hospitality of the Kavli Instutute of Theoretical Physics in Santa Barbara, CA and thank Tom Hull, Louis Theran, Menachem Stern and Arvind Murugan for helpful discussions. BGC thanks ICERM in Providence, RI and the participants and staff during the 2016 “Topology in Motion” semester program there for support and encouragement. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1125915, and also Grant No. EFRI ODISSEI-1240441.

Appendix A Deriving a quadratic constraint for each internal vertex from the vanishing of the Gaussian curvature

One way to express the Gaussian curvature constraint around a single vertex is in terms of the interior angles of the spherical triangle made by a pair of adjacent folds and the $\hat{\mathbf{z}}$ axis (Fig. 11). We denote $\alpha_{i,i+1}$ as the planar angle between adjacent folds, which becomes the length of one side of the triangle. Similarly, we denote $\psi_{i}$ as the angle between $\hat{\mathbf{z}}$ and the $i^{th}$ fold.

Let $h_{0}$ be the height of the central vertex, and $h_{i}$ be the heights of the vertices around the interior vertex at $h_{0}$ and $L_{i}$ are the lengths of the folds from $h_{0}$ to $h_{i}$. We assume that all quantities are periodic in the index and numbered in counterclockwise order. Then we have angles,

[TABLE]

Finally, define the interior angle of the spherical triangle between $\hat{\mathbf{z}}$ and the two folds as $\beta_{i,i+1}$ (Fig. 11). The relationship between these angles is given by the spherical law of cosines,

[TABLE]

Since the vertices we consider are nearly unfolded, we will expand around an unfolded configuration in which $\psi_{i}=\pi/2$ and $\beta_{i,i+1}=\alpha_{i,i+1}$. Therefore, if $\psi_{i}=\pi/2+\delta\psi_{i}$, $\beta_{i,i+1}=\alpha_{i,i+1}+\delta\beta_{i,i+1}$, expanding to quadratic order yields

[TABLE]

The sum of the interior angles $\sum_{i}\alpha_{i,i+1}=2\pi-K$, where $K$ is the deficit angle, equivalent to the discrete Gaussian curvature of the vertex; however, $\sum_{i}\beta_{i,i+1}=2\pi$ no matter $K$. Therefore, for small $K$, we require that

[TABLE]

When the origami is unfolded, $K=0$ at each vertex.

Now we can rewrite the angle $\delta\beta_{i,i+1}$ in terms of the heights as

[TABLE]

where $\Delta L_{i,i+1}=\sqrt{L_{i}^{2}+L_{i+1}^{2}-2L_{i}L_{i+1}\cos\alpha_{i,i+1}}$ is the distance between the vertices at $h_{i}$ and $h_{i+1}$.

Finally, we have

[TABLE]

This is precisely the form of the quadratic constraint that we get from the self stresses in Eq. (5). In the next appendix, we verify that the coefficients here do correspond to a self stress.

Appendix B Verification that the Gaussian curvature constraint yields a self stress

Using Eq. (A), we can immediately read off the components of the stress that would give the quadratic constraint in Eq. (A). In particular, the stress in the edges along the “rim” of the wheel should be

[TABLE]

and on the spoke edges

[TABLE]

What remains to to verify that these do in fact satisfy the equation $\mathbf{\sigma}^{T}\cdot\mathbf{C}=\mathbf{0}^{T}$ defining self stresses. Recall that self stresses are assignments of tensions and compressions to edges that preserve force balance at each vertex. Therefore we must check force balance at the spoke vertices (labeled $i=1$ to $N$) and the hub vertex 0.

B.1 Force balance at the spokes

We first show force balance at spoke vertex $i$. Recall that the position of vertex $j$ is the vector $\mathbf{U}_{j}$ so that $L_{i}=|\mathbf{U}_{i}-\mathbf{U}_{0}|$ and $\Delta L_{i,i+1}=|\mathbf{U}_{i+1}-\mathbf{U}_{i}|$ (Fig. 12). Let’s first check the forces in the direction perpendicular to the spoke edge vector ($\mathbf{U}_{i}-\mathbf{U}_{0}$). We only need to use the stresses along the outer edges (Eq. (19)) for this.

First we use the law of sines on the triangle with sides $L_{i},L_{i+1},\Delta L_{i,i+1}$ to transform Eq. (19) to

[TABLE]

where $\eta$ is the angle opposite $L_{i+1}$. A similar argument shows that

[TABLE]

where $\eta_{i-1,i}$ is opposite $L_{i-1}$.

The magnitudes of the forces perpendicular to the spoke vector are given by $\sigma_{i,i+1}\sin\eta_{i,i+1}=-\sigma_{i-1,i}\sin\eta{i-1,i}$ so there is force balance along this direction at vertex $i$.

Next we consider the force at vertex $i$ in the direction parallel to the spoke vector. Using Eqs. (21) and (22), we see that the contribution from the two rim edges to the force parallel to the spoke vectors may be written

[TABLE]

Consider the term $F_{rim,\parallel,-}$. This is minus the reciprocal of $L_{i}\tan\eta_{i-1,i}$ which is the side length $A$ of a certain right triangle with $L_{i}$ as one of its legs (Fig. 12). There is a similar interpretation for the term $F_{rim,\parallel,+}$.

Now consider the contribution from the spoke edge given by Eq. (B). We will split this expression into two pieces, which will each cancel one of the two terms $F_{rim,\parallel,\mp}$. The first piece is

[TABLE]

The first term of Eq. (26) is the reciprocal of $L_{i-1}\sin\alpha_{i-1,i}$, which is the side length $h$ of a right triangle with $L_{i-1}$ as the hypotenuse. The second term is the reciprocal of $L_{i}\tan\alpha_{i-1,i}$ which is the side length $B$ of a (different) right triangle with $L_{i}$ as one of its legs. See Fig. 12 where these sides are depicted as dashed lines.

The three side lengths $A,B,h$ are related as in the “Crossed Ladders problem” Wikipedia (2016): for parallel line segments $A,B,h$ in the configuration shown in Fig. 12, we have the identity $1/A+1/B=1/h$. Therefore the first part of the force along the spoke edge (Eq. (26)) cancels $F_{rim,\parallel,-}$ (Eq. (24)).

It should be clear that the argument of the preceding three paragraphs can also be carried out for the second piece of Eq. (B), to show that it cancels $F_{rim,\parallel,+}$.

B.2 Force balance at the hub

For force balance at the hub vertex, we proceed with an induction on the number of spokes. For the base case with 3 spokes, we must show that the vectors with lengths $\sigma_{1},\sigma_{2},\sigma_{3}$ directed along the spoke directions sum to zero. Let $\mathbf{u}_{j}=(\mathbf{U}_{j}-\mathbf{U}_{0})/|\mathbf{U}_{j}-\mathbf{U}_{0}|$ be the unit vector along the $j$th spoke edge. Then this is:

[TABLE]

We interpret this equation as the condition for the closure of the polygonal path with sides $\sigma_{j}\mathbf{u}_{j}$. Generally, when a vertex is in equilibrium, the forces acting on it can be seen as the sides of a closed force polygon. Note that in this force triangle the angles $\alpha$ between the spoke edges become the “turning angles”, so that e.g. the angle opposite side $\sigma_{1}$ in the force triangle is $\pi-\alpha_{2,3}$. We can show that this triangle is closed using the (converse of the) law of sines, which amounts to proving that the following are all equal:

[TABLE]

We will just show the first equality, as the proof of the others is exactly the same.

[TABLE]

We will need the following trigonometric identity, equivalent to the sine addition rule:

[TABLE]

Note that $-\csc\alpha_{3,1}=\csc(\alpha_{1,2}+\alpha_{2,3})$ since $\alpha_{1,2}+\alpha_{2,3}+\alpha_{3,1}=2\pi$. Therefore by applying Eq. (31) twice, we have,

[TABLE]

For the induction step, assume that we have shown that we have force balance at the hub vertex in any wheel graph with $n$ spokes with stresses given by Eq. (B). Consider the star subgraph $G_{n+1}$ of the wheel graph, consisting of $n+1$ spoke vertices connected to the hub; we will reduce the force balance at its hub to force balance in the star graph $G_{n,i}$ formed by removing (an arbitrary) spoke vertex $i$ and the edge joining it to the hub. As in the argument for the base case, we will work with force polygons and prove that the closure of the force polygon $P_{n,i}$ of $G_{n,i}$ implies the closure of the force polygon $P_{n+1}$ of $G_{n+1}$.

Note that the expression in Eq. (B) for spoke $j$ depends only on the lengths $L_{j-1},L_{j},L_{j+1}$ and the angles $\alpha_{j-1,j}$ and $\alpha_{j,j+1}$. This means that the self stress evaluated on corresponding spokes of $G_{n+1}$ and $G_{n,i}$ are identical except at the edges $i-1,i,i+1$. Thus $P_{n,i}$ and $P_{n+1}$ are identical except at those edges too. We will prove that the situation is as depicted in Fig. 13. There, the edges of $P_{n+1}$ are depicted as thick lines, with the dashed edges $i-1,i+1$ of $P_{n,i}$ overlaid.

To get started, we observe that the closure of $P_{n,i}$ implies there is a vertex where the edges corresponding to spokes $i-1$ and $i+1$ (dashed in Fig. 13) intersect. We will denote lengths of edge $i-1$ and $i+1$ in $P_{n,i}$ are $\sigma^{0}_{i-1}$ and $\sigma^{0}_{i+1}$, respectively. In $P_{n+1}$, edges $i-1$ and $i+1$ can be taken to lie on the corresponding edges of $P_{n,i}$ but they will in general have different lengths, which are simply $\sigma_{i-1}$ and $\sigma_{i+1}$. For $P_{n+1}$ to be closed, edge $i$ must begin at the end point of $i-1$ and end at the starting point of $i+1$. This means that there is a closed triangle with side lengths $\Delta\sigma_{i-1}=\sigma_{i-1}-\sigma^{0}_{i-1}$, $\Delta\sigma_{i+1}=\sigma_{i+1}-\sigma^{0}_{i+1}$, and $\sigma_{i}$ and whose angles are determined by the angles $\alpha_{i-1,i}$ and $\alpha_{i,i+1}$ between the spoke edges $i-1,i$ and $i+1$. Note that the angle $\beta_{i}=\pi-\alpha_{i-1,i}-\alpha_{i,i+1}$.

To prove that this triangle is closed, we will again use the law of sines. For convenience, here are the formulas for the lengths:

[TABLE]

We will only show $\sigma_{i}/\sin\beta_{i}=\Delta\sigma_{i-1}/\sin\alpha_{i,i+1}$, as the proof of the other identity is the same.

[TABLE]

where in the second and third equalities we have applied Eq. (31). This completes the proof that Eqs. (19) and (B) define a self stress on the wheel graph.

Appendix C Vertex Quadratic Forms

Let us consider the quadratic form in Eq. (A). This form has the interpretation as the energy if the wheel self stress is applied as a pre-stress Connelly and Whiteley (1996). Note that with our sign convention, this pre-stress places the spokes under compression and the outer edges under tension.

Using the same notation as in Section II, where $h_{0}$ is the height of the vertex and $h_{1}$ through $h_{N}$ are the heights of the adjacent vertices, Eq. (A) then reads

[TABLE]

where

[TABLE]

The corresponding matrix has three null directions, corresponding to the vertical translation and rotations about the $x-$ and $y-$axes of the entire origami structure. We can remove these by setting $h_{0}=h_{1}=h_{N}=0$ explicitly. Then we have

[TABLE]

This gives a tridiagonal matrix whose determinant, $\chi(N)$, can be computed to be

[TABLE]

This can be proven by induction using the continuant,

[TABLE]

When all the angles between vertices are identical, $\alpha_{n,n+1}\equiv\alpha$ and the lengths of the folds are also the same, $L_{n}=L$, the quadratic form of Eq. (45) is also Toeplitz, and its eigenvalues can be determined explicitly by

[TABLE]

for $m$ ranging from $1$ to $N-1$. For $N$ folds, $\alpha=2\pi/N$. Thus,

[TABLE]

Consequently, $\lambda_{1}<0$ and the other eigenvalues are positive. Now consider changing the angles smoothly. Since no eigenvalue can change sign without the determinant $\chi=0$, we see that Eq. (46) implies that there is always one negative eigenvalue so long as the angles between any pair of adjacent folds remain between 0 and $\pi$. When $N\geq 4$ there are thus eigenvalues of both signs, and there is always a solution to Eq. (45) for a single vertex; the vertex is infinitesimally rigid when all eigenvalues have the same sign (as happens when $N=3$, for instance).

We can get some physical intuition for the distribution of eigenvalues as follows. Under the wheel pre-stress corresponding to this quadratic energy, the eigenvector with negative eigenvalue should correspond to an out-of-plane displacement that is maximally unstable. We can imagine that there is one that increases the lengths of the compressed spokes while not increasing the lengths of the outer edges too much, as such a motion will be destabilized by the stress. The other $N-3$ motions correspond to out-of-plane displacements that are stabilized by the pre-stress, i.e. those that primarily increase the lengths of the outer edges.

The discussion above recovers a special case of a result of Kapovich and Millson Kapovich and Millson (1997) who studied the configuration space of origami vertices using techniques from the deformation theory of representations of $SO(3)$. In fact, they considered unfolded origami vertices that may fold back on themselves (allowing $\alpha_{j,j+1}<0$ for some $j$) or have a different winding around the vertex ($\sum_{j=1}^{N}\alpha_{j,j+1}=2\pi w$, for some integer $w$ not necessarily equal to 1). The relevant result of their paper is the following

Theorem 1** (Kapovich and Millson, 1997)**

Theorem 1.1(iii). Assume all planar angles satisfy $0<|\alpha_{j,j+1}|<\pi$. The germ of the configuration space of an origami vertex with $f$ forward-tracks, $b$ back-tracks, and winding $w$ is isomorphic to the germ of the null-cone of a quadratic form with nullity 3, and signature $(f-2w-1,b+2w-1)$.

Here $f$ counts the number of forward-tracks, defined to be the planar angles with $\alpha>0$ and $b$ counts the number of back-tracks, defined to be those angles with $\alpha<0$.

We showed above the case of this theorem when $f=N$, $b=0$ and $w=1$, and in fact our sign convention for the sign of the quadratic form also agrees with theirs. (A change of sign swaps negative and positive eigenvalues, but of course leads to the same null-cone). Furthermore, it is easy to see that our observation that Eq. (46) only vanishes when one of the $\alpha$ is an integer multiple of $\pi$ is consistent with their statement that the signature changes when any of $f,b,w$ change. Indeed, Eq. (45) for the quadratic form applies in full generality, and so we can actually use it to give a complete proof of Theorem 1.

Proof of Theorem 1: The first claim about the germ being isomorphic to the germ of the null-cone essentially states that the quadratic form we have derived in Eq. (45) as the lowest order constraint on the configuration space does give an accurate picture of a neighborhood of the singular point (the unfolded state), i.e. that the second-order motions satisfying Eq. (45) extend to true motions. For this we refer to Theorem 4 of Demaine et al. (2016) which gives an elementary proof of this fact. (In fact, Kapovich and Millson (1997) prove the stronger result that there is an analytic isomorphism between neighborhood germs.)

The rest of the theorem addresses the signature of the quadratic form. The nullity of 3 corresponds to the global isometries mentioned earlier. So we just need to derive the expression $(f-2w-1,b+2w-1)$ for the number of (positive, negative) eigenvalues. Since the defining symmetric matrix is $(n+1)\times(n+1)$ and $f+b=n$ it is enough to show that the number of negative eigenvalues $N_{-}=b+2w-1$.

Suppose we have a real symmetric $k\times k$ matrix $M$, and let $\Delta_{j}$ is the determinant of the upper-left $j\times j$ submatrix of $M$ (the $j$th principal minor of $M$). The key observation is that, provided none of the $\Delta_{j}$ vanish, $N_{-}$ is equal to the number of sign changes of the sequence $\Delta_{0}=1,\Delta_{1},\dots,\Delta_{k}$ Ghys and Ranicki (2016).

Due to the tridiagonal nature of the stress matrix, $\chi(j)$ is also an expression for $\Delta_{j}$ of the stress matrix. So we just have to show that we get exactly $b+2w-1$ sign changes. Let us consider the ratio:

[TABLE]

This consists of two factors. The first factor $\csc\alpha_{j-1,j}$ is negative if and only if $\alpha_{j-1,j}<0$, i.e. when the $j$th planar angle is a back-track. The second factor is negative if and only if the two quantities $\sum_{n=1}^{j-1}\alpha_{n,n+1}$ and $\sum_{n=1}^{j}\alpha_{n,n+1}$ sandwich an integer multiple of $\pi$. Let us first assume that the second factor does not vanish (the first never will, by our assumptions on the $\alpha$’s).

We get a net negative sign in Eq. (50) if and only if one of the following two possibilities occurs. Possibility A, the $j$th planar angle is a back-track and the partial sum $\sum_{n=1}^{j}\alpha_{n,n+1}$ does not pass an integer multiple of $\pi$. If $C_{b}$ is the number of back-track crossings over integer multiples of $\pi$, this occurs $b-C_{b}$ times. Possibility B, the $j$th planar angle is a forward-track and the partial sum does cross an integer multiple of $\pi$. This occurs $C_{f}$ times, where $C_{f}$ is the number of forward-track crossings of integer multiples of $\pi$.

It follows that $2w-1$ is the net number of times $C_{f}-C_{b}$ that integer multiples of $\pi$ are crossed, since $w$ is the total winding number. Therefore we have $b-C_{b}+C_{f}=b+2w-1$ sign changes in total, as desired.

Finally, we treat the case where $\Delta_{j}$ vanishes due to the partial sum $\sum_{n=1}^{j}\alpha_{n,n+1}$ being equal to an integer multiple of $\pi$. However, as argued after Eq. (49), we can perturb the angles $\alpha$ to avoid these cases without changing the sign of the overall determinant (note that the angle sum in $\chi(N)$ is equal to $2\pi-\alpha_{N,1}$), and hence without changing the signature. Such a perturbation can also be chosen small enough so that $f,b,w$ do not change, so the same formula applies.

This concludes the proof of Theorem 1.

Appendix D High-dimensional geometry and random origami

To better understand how the branches of a random triangulated $D=1$ origami are distributed as directions in configuration space, we computed pair distribution functions between the branches, where the branches are intepreted as lines in the $3V_{i}+1$-dimensional space of infinitesimal changes in dihedral angle. (Note that the branches lie in a $V_{i}+1$-dimensional linear subspace, as not all sets of dihedral angles are induced by height displacements). Computationally, given a particular configuration, we take the numerically computed branches given as vectors of vertex heights and transform them into vectors of dihedral angles. Here we choose to work with both $\pm$ ends of the branch line, so that we have $2^{V_{i}+1}$ vectors. We then compute the $(2^{V_{i}+1})^{2}$ angles coming from dot products between all ordered pairs of branches. These angles are the distances between branches as points on the $V_{i}$-dimensional sphere. For each $V_{i}$ we compute these angles for all of the configurations in Table 1 to create histograms of the angles, or pair distribution functions.

In Fig. 14, we show the histogram of the angle between pairs of branches for $V=6$ using bins of size $\pi/128$ (black line). We compare the results to random points on a $V_{i}$-dimensional sphere (gray line). The data shows that random origami has a slight enhancement in the number of orthogonal branches but a more prominent enhancement at the two tails of the distribution. In particular, the results of Ref. Cai et al. (2013) imply that for random points on a sphere, the angle distribution approaches a Gaussian with variance going to zero. However, the tails of the distributions of angles between branches in random origami appear to decay exponentially, rather than as a Gaussian.

To indicate what happens as $V_{i}$ changes, Fig. 15 depicts the distributions for random origami with vertices from $V_{i}=2$ – $8$, again binning the results into bins of size $\pi/128$. Since the branches have no natural orientation, for each angle $\theta$, we also include the angle $\pi-\theta$; as a consequence, for small $V_{i}$, there is an enhancement for branches that are almost $\pi$ apart in angle, since each branch has an angle [math] and $\pi$ with itself. It is interesting to note that while the distributions for $V_{i}=6,7$ very nearly coincide, they differ significantly from that of $V_{i}=8$. We do not believe that this is due to a lack of data. Even though we computed fewer configurations at $V_{i}=8$, each configuration has twice the number of branches and hence four times as many angle pairs, so the number of data points going into the histograms is comparable.