Rigidly flat-foldable class of lockable origami-inspired metamaterials with topological stiff states

Origami crease patterns have inspired the design of reconfigurable materials that can transform their shape and properties through folding. Unfortunately, most designs cannot provide load-bearing capacity, and those that can, do so in certain directions but collapse along the direction of deployment, limiting their use as structural materials. Here, we merge notions of kirigami and origami to introduce a rigidly foldable class of cellular metamaterials that can flat-fold and lock into several states that are stiff across multiple directions, including the deployment direction. Our metamaterials rigidly fold with one degree of freedom and can reconfigure into several flat-foldable and spatially-lockable folding paths due to face contact. Locking under compression yields topology and symmetry changes that impart multidirectional stiffness. Additionally, folding paths and mixed-mode configurations can be activated in situ to modulate their properties. Their load-bearing capacity, flat-foldability, and reprogrammability can be harnessed for deployable structures, reconfigurable robots, and low-volume packaging.


S1. Kinematic analysis
To study the mobility of our unit kinematic chain ( Supplementary Fig. 1a), we use a surrogate bar network where the planar faces are replaced with pin-jointed bars ( Supplementary Fig. 1b, c). This assembly consists of n b inextensible bars, n j pin joints and rigid triangular panels; quad panels are replaced by a set of two connected rigid triangles, thus eliminating in-plane shear deformation. The rigidity of each quad, or triangular panel, is ensured by imposing two constraints: (i) bar inextensibility and (ii) panel coplanarity for the triangular faces making up a quad panel ( Supplementary Fig. 1c, d).
Supplementary Fig. 1 Unit chain with its representative bar and hinge kinematic model. a Geometric parameters describing our reconfigurable unit with N = 4. Superscript L assigned to a given variable in the first locked configuration mode of the unit chain, when all dihedral angels i are equal. b Equivalent network of rigid panels connected through flexible hinges replaced by c a bar and pin-joint mechanism described by parameters used to impose planarity condition d.
To determine the first set of constraints for panel rigidity, we first denote the position of joint i by vector x i (boldface used here for vectors and tensors) and joint j by vector x j ; the bar connecting joint i to joint j is given by the vector r ij = (x i -x j ) with length where ‖ ‖ represents the vector norm. With this notation, the first-order inextensibility condition for bar r ij can be written by equating the first derivative of the bar length to zero, 0, as where dx i and dx j are the infinitesimal displacements of the joints i and j, respectively.
The second set of constraints imposes panel planarity for the quad faces, to impede the rotation of triangular faces that make the quad panel. It is formulated as follows. Assume a quad panel connecting joints i, j, k, and l, and comprising two triangular panels p (connecting joints i, j and k) and q (connecting joints i, j and l), with the normal vectors n p = r ij × r ik and n q = r ij × r il ( Supplementary Fig. 1d). The angle between two adjacent panels p and q can be expressed as cos • . Equating the first differential of the angle  to zero, 0, gives the coplanarity condition for the quad panel as 0. (S2) The rigidity constraints, Eqs. (S1) and (S2), can be expressed in matrix form for all panels as where R is the rigidity matrix consisting of the inextensibility constraints (compatibility) matrix B and the planarity constraints matrix P, and dx is the infinitesimal displacement vector of the joints.
The relations above can now be used to determine the number of degrees of freedom (DoF) m, of the pinjointed triangulated network [1] as where n K is the number of external kinematic constraints, such as those confining the rigid motions, and r is the rank of the rigidity matrix R in Eq. (S3). In Eq. (S4), d represents the dimensions of the problem; for spatial mechanisms d is 3. For an infinitesimally rigid structure m is zero, while for infinitesimal mechanisms m > 0.
With the rigidity matrix, we can now use the singular value decomposition (SVD) to determine the kinematic paths that our unit can travel upon bifurcation and track its configuration using a predictorcorrector type incremental method [2].
to the predicted configuration C i′ as When a new path has been found, we must control the independency of the kinematic paths. We do so by comparing the newly found path with other paths. This strategy is necessary since using the predictorcorrector algorithm with the first-order analysis does not guarantee that the obtained path converges to a distinct finite mechanism. The predictor-corrector algorithm is implemented in MATLAB, and all the simulations are based on an in-house developed code.
The algorithm can execute large displacement simulations, but cannot find all possible postbifurcation kinematic paths. We resort to principles of Pólya Enumeration Theory that enables us to find all independent post-bifurcation kinematic paths, as explained in Supplementary Discussion S3.

S2. Additional constraints (imposed by unit chain stacking)
A single unit kinematic chain has multiple DoFs that lead to certain configurations, and the fully developed (flat) configuration is one among others. To reduce its DoF we resort to layer stacking, which leads our unit to share certain panels with those of the units placed above and below. Supplementary Fig. 2a shows a physical prototype where only the triangular panels are shared between two units, whereas the others are not shared. Here, the non-zero thickness of adjacent hinged panels causes an offset between the hinges.
Supplementary Fig. 2b shows the only possible rigid body motion the planes can undergo. The quad planes cannot rotate in the same direction, shown in Supplementary Fig. 2c. This is due to the internal forces exerted by planes and hinges. The outcome is that both plane interference and hinge offset geometrically constraint the unit to respect the rigid body motion. For example, the motion in Supplementary Fig. 2c is only possible if the quad planes can undergo flexural deformation ( Supplementary Fig. 2d), which violates our rigid motion assumption.
Supplementary Fig. 2 Layer stacking as a pathway to reduced mobility. a Flat configuration for a two-layered unit kinematic chain made of non-negligible thickness panels. b A kinematic path of the multilayered units, where two adjacent quad panels rotate in opposite directions under in-plane load f. c Example of inadmissible kinematic path for two quad panels rotating in identical direction under in-plane load f; this kinematic path violates the rigidity assumption of the panels and is feasible only if the neighboring panels can deform d.

S3. Pólya enumeration theorem
The post-bifurcation kinematic analysis given in the previous section cannot guarantee that all postbifurcation paths are found [3]. We now use principles of Pólya Enumeration Theory [4] to identify and count all independent post-bifurcation kinematic paths and available modes. The Pólya Enumeration Theorem is instrumental to identify the reconfiguration modes of our unit kinematic chain.
As stated in the main text of the paper, the problem of finding all independent regular modes of N n is equivalent to the classical necklace problem in combinatorics. There, the goal is to find the number of necklace arrangements we can construct from knowledge of N/2 colored beads, each painted in two distinct colors, e.g., white or black. In our case, the dihedral angle of our unit chain represents one color, e.g., white, and the other, e.g., black. Due to the mountain and valley assignment, the transformations are defined on hypothetical N/2 sided polygons since congruent triangular faces are located on every other corner of the original chain. Thus, the necklace problem of an N even-sided polygonal unit chain, N N , with constraints of mountains and valleys degenerates to an N/2 sided necklace.
We start by stretching the necklace with beads into the shape of a uniform -sided regular polygon with one bead at each corner. We define the set of the vertices of this polygon by S = { 1, 2, 3, …, } and the set of binary colors by R = { , }, and introduce set X = {all possible coloring arrangement of S by R}. Next, we define rigid body transformations of the necklace. We assume that the beads can rotate around the necklace, and the necklace can be flipped over. These operations are considered as rigid motions applied to the necklace. To illustrate the above, we consider N6 as an example for the case of 0˚ rotation. We note the vertex assignment is identical to the initial assignment. The transformation corresponding to its rigid motion is the identity operation, which maps every vertex to itself regardless of the number of times the transformation is applied. Its vertex counterpart can be written as (1)(2)(3). By using the Pólya's enumeration theorem, this rigid motion is made of three cycles each of length one. On the other hand, a 120˚ rotation applied to N6 is described by (132), which consists of 1 cycle of length three.

Supplementary Fig. 3 Schematic representation of the unit kinematic chains with its equivalent necklace.
Schematic representation of unit kinematic chain with N = 4, 6, 8 and 10 (left) and its equivalent necklace (right) with all its possible rigid motions below the first row for each unit chain. Two valley "dihedral angle pairs" are shown as an example for N4 only. Symmetries of the necklaces are also represented as permutation of the vertices.
Below the first row are the possible rigid body motions each unit chain can undergo. The alternating light and dark colors on the triangular faces of the kinematic chain refer to mountain and valley assignments, respectively. The equivalent necklace has numbered vertices and describes a given rigid motion. Below each rigid motion is a sequence of parentheses describing the rigid body motion in cycle notation. A black dot represents zero-degree rotation. A combination of dashed lines and circular arcs represent the rotation operations. A dashed line is used for reflection operations and indicates the mirror line.

Assume
∈ G is a rigid motion operation that can be described by a combination of its cycles. . The mode is lockable when at least two successive acute dihedral angles exist in the pattern sequence. For example, in the case of N n , the lockable mechanisms are (i.e., the one with the sequence of ), and . Supplementary

S4. Tessellations: from individual stacks of unit chains to a periodic cellular material
Each multilayered unit chain can be joined with others to form a periodic material system that can reconfigure along either lockable or flat-foldable modes. To examine their in-plane tessellation patterns (see Supplementary Fig. 4b), we refer to the periodic array our N n forms in 2 dimensions, with centroids falling onto a simple Bravais lattice of the base-centered orthorhombic family ( Supplementary Fig. 4a).
Here, we study ways of connecting multilayered unit chains to tile the x-y plane while preserving the characteristic kinematics. We do so by joining in-plane the triangular rigid panels of our multilayered unit chains, without making use of flexible hinges at their interface. This strategy is realized through two types of connections: 1) the whole triangular panel is shared between two multilayered unit chains, i.e., the two patterns overlap (see for example the tessellation pattern for N = 6 in Supplementary Fig. 1b); 2) the interface of the base edge of their triangular panels acts as connection; this means that the two triangular panels form a single parallelogram panel that is rigid (see for example, the tessellation pattern for N = 4 in Supplementary Fig. 1b). This leads us to conceive two classes of tessellation containing similar unit chains: i) Most packed tessellation patterns. This class describes the densest lattices that can be made from our primitive units, ii) Less packed tessellation patterns. This class includes realizations of other possible low-density lattices; it is analyzed through the notion of macro-chains.
Other types of connection between dissimilar units are possible, but are not examined here.

S4.1. Most packed tessellation patterns
The most compact form of tessellation can be described by a periodic array of N n multilayered unit chains in two dimensions, whose centroids (centers of mass) fall onto a simple Bravais lattice of the base-centered orthorhombic family with basis vectors {e1, e2}.
Let us consider a rhombus constructed from vertices that are the centroids of the four adjacent unit chains in a tessellated pattern ( Supplementary Fig. 4a). The geometric constraints imply that its vertex angle ≡ cos • impedes the overlap of unit chains along the short diagonal of the rhombus ( Supplementary Fig. 4a). Additionally, the vertex angle must be made of segments of (where ), i.e., , and its supplementary angle ′ must be made of ′ segments of , i.e., . All admissible angles can be obtained by solving the linear Diophantine equation [5] /2 (S12) which is written assuming the condition (see Supplementary Fig. 4a).
Supplementary Table 2 Fig. 4). In certain cases, e.g., hexagonal family ( ) with an odd (see for example, the case N = 18 in Supplementary Table 3

S4.2. Less-packed tessellation patterns using macro-chains
An alternative way to tessellate our pattern is to start from a macro-chain containing primitive units. A macrochain is constructed by the in-plane connection of a certain number of similar N n unit chains, i.e., units with prescribed N, to form a regular polygon; this is done by replacing each edge with one primitive unit, e.g., Supplementary Fig. 4c, d. The shape of macro-chains is governed by the geometry of their constituents. Their existence can be formulated as follows. We denote the i th positive factor of N with ℳ . For any given ℳ 4 one ℳ -sided regular polygonal macro-chain (N ℳ n ) exists provided either of ℳ or ℳ ℳ 2 becomes an even number. N ℳ n is a macro-chain connecting the centroids of the unit chains N n .
In-plane tessellations of macro-chains can be formulated using the method described above if they fall into one of the Bravais lattice families 2 . Supplementary Fig. 4b illustrates the smallest macro-chain, i.e., N n , and its tessellation that can be built using the unit kinematic chain N n . Supplementary Fig. 4c, d, e shows examples of possible tessellations one can imagine using the macro-chain of N n ; several more exist that are beyond the scope of this work. In the figures, orange lines indicate the connectivity between units in a macro-chain, and red lines describe the connectivity of macro-chains in the tessellation.
In the tessellations of N ℳ n macro-chains, the connection of the macro-chains may be much more complex than that of the single unit chains, as inter-chain connections can also be established, e.g., red dash-lines in Supplementary Fig. 4c. Here, the main bonds between macro-chains occur either through the bases of the isosceles triangulated panels that intersect the direction vectors e 1 and e 2 (see Supplementary   Fig. 4c, d), or through the overlapped unit-chains N n as shown in Supplementary Fig. 4e.
The results of the analysis above exemplify the breadth of the design space that can be achieved through the selection of alternative tessellation patterns. Each tessellation pattern is characterized by its own set of rigid foldability attributes, physical properties, and load-bearing capacity.
2 ℳ 4 results in the Bravais lattice of the tetragonal family and ℳ 3 results in the Bravais lattice of the hexagonal family. Supplementary Fig. 4 Most-packed and less-packed tessellations. In-plane view of most compact tessellations for N n along with three other possible less-packed tessellations for N n .

S5. Rigidity
To study the mobility of our unit, we replace the actual panels with their bar and hinge counterparts. Upon locking, certain bars and hinges coincide ( Supplementary Fig. 5). This holds only if, at the lock state, we assume there is full contact between edges and faces, and thus relevant bars and hinges condense into a single bar and a single hinge. Such an assumption underpins our rigidity analysis and is considered valid if the direction of the force is unchanged [6].
To analyze the rigidity of our unit chains one panel is constrained to eliminate rigid body motions.
Supplementary Table 3

S6. Geometric mechanics of representative unit cell (RUC)
During reconfiguration, our unit chains pass a kinematic bifurcation instant that provides access to dissimilar kinematic paths, each having modes that can be either flat-foldable and lockable. Each path has its own physical characteristics, e.g., geometry, relative density, and Poisson's ratio. In this section, we first derive the closed-form expressions that describe changes in properties upon folding for N n and N n units. Next, we derive the representative energy landscape relationships for a given lockable and flat-foldable modes, and generate energy-phase diagram that can explain the underlying physics of activation through in-plane confinement.

S6.1 RUC model with boundary and loading conditions for unit cell
We examine a representative unit cell for N n that can describe the behavior of its entire periodic tessellation. Supplementary Fig. 6 shows the RUC with its own geometric parameters and loading conditions in a Cartesian coordinate system. In modes and , the RUC is a square (red) with side L 1 = 2S 1 , in which S 1 is the length of the line connecting the mid-points of the bases of two adjacent triangles separated by a parallelogram (right of Supplementary Fig. 6b) and can be calculated as cos sin cos . (S13) Upon kinematic bifurcation, the square shape of the RUC switches to a parallelogram. One of its sides can be computed from Eq. (S13), while the other side L 2 = 2S 2 can be calculated by replacing in Eq. (S13) as cos sin cos . (S14) To obtain the skewness angle of the parallelogram angle in Supplementary Fig. 6b, we first calculate the angles and from tan cos sin cos (S15) tan cos sin cos and then, replace them in the relation
(S16) Then, the cross-section area of the RUC in the x-y plane can be obtained from sin (S17)

S6.2 RUC model with boundary and loading conditions for unit
We now examine N n unit and derive geometric relations that describe the deformations of its kinematic paths. In this case, the smallest unit with parallel edges is a rectangle. The length of the line connecting the midpoints of the two adjacent triangles are If we connect the centroids of every other triangle around a unit chain, we obtain an additional triangle.
The sides of this triangle relates to S1 and S2, and s through The parallelogram of the RUC is defined by two triangles. One of the sides of the parallelogram is while the other side is given by Now, using the cosine rule, we obtain

S7. Relative density
The relative density of our multimodal rigid-foldable materials upon reconfiguration (and transition between modes) can be expressed as a function of the dihedral angle and other geometric parameters of the unit chain for a given tessellation type. In addition, upon kinematic bifurcation multiple paths emerge, and a distinct relative density can be associated to each given kinematic path. Supplementary Fig. 8 shows the relevant geometric parameters for N n , here taken as example for the analysis. Supplementary Fig. 8 General geometric parameters of our tessellated units. a Geometric parameters of a unit kinematic chain; b Connecting the mid-point of the bases of the triangles yields an N-sided polygon, shown in red; c Geometric parameters of the Bravais lattice tessellation. Regardless of N, four units can always be combined to form a rhombus defined by the angle , which is obtained from Eq. (S12). Supplementary Fig. 8b shows that by connecting the mid-points of the bases of the triangles, we obtain a regular polygon (red decagon here). R is the radius of the circumscribing circle of the polygon. The side of the polygon, L, is given by which relates to R through The volume of the rhombus shown in Supplementary Fig. 8 with height h (recall that it is now a 3D structure) for modes / and / can be related to R using in the x-y plane. We observe the relative density of flat-foldable modes approaches infinity as the lattice approaches its flat-foldable configuration. Also, each deformation mode attains a certain minimum relative density during the reconfiguration process. As we expect, lock mode / has the highest density of all locked configurations, while / is the least dense flat-foldable mode.
Supplementary Fig. 9 Relative density upon reconfiguration for different lock and flat-fold modes. Relative density computed for N n and N n lattices during reconfiguration along their two kinematic paths. t is the panel thickness, and a is the length of the sides of the primitive regular polygon.

S8. Poisson's ratio
During folding our metamaterials exhibit auxetic in-plane behavior. The following derives expressions of the Poisson's ratio that emerge at given modes.

S8.1. Poisson's ratio of / and / modes
To obtain the expressions of the Poisson's ratio for N n , we begin by tessellating the rhombus of Supplementary Fig. 8. By doing so, we can express the increments of the strains, i.e., , and , in terms of the derivatives of the parameters ℒ , ℒ and the height of the N n unit. We first express ℒ , ℒ as a function of the radius of the circumscribing circle, R, and the skewness angle . Interestingly, we observe a change in the sign of the out of plane Poisson's ratios. This entails that passing the kinematic bifurcation, our folding metamaterials in / mode have omnidirectional auxetic behavior.

S8.2. Poisson's ratio of other modes
The Poisson's ratio in other modes depends on the lattice (i.e., N n ) and the kinematic path. We examine here N n and N n and present the results of the Poisson's ratio for other mode shapes.
To simplify the calculation, we refer to the RUC in Supplementary Figs. 6 and 7. In this case, the in-plane Poisson's ratio expression can be written as where the parameters , and are presented in Eq. (S13), Eq. (S14) and Eq. (S16) for N4 pattern. The expression for the Poisson's ratio of N6 pattern can be written as The parameter , and can be obtained from Eqs. (S23), (S24) and (S26). Supplementary Fig. 10 shows the in-plane Poisson's ratio calculated for N n and N n evolving along their kinematic paths. We observe that the in-plane Poisson's ratio of each lattice highly depends on the tessellation and kinematic path as opposed to that of the modes / and / .
dihedral angle, and are the in-plane forces. ⟘ and ⟘ are the distances of the boundaries of RUC where and are applied ( Supplementary Fig. 6), and can be defined as Moreover, the initial distances of the RUC boundaries are ⟘ and ⟘ . Eq. (S40) is valid for all in-plane loading conditions of N n .

S9.2. Energy analysis of under in-plane biaxial loading
Starting from a flat configuration (i.e., 0, 0) and applying biaxial loads, we can find a closed form expression of the critical bucking load that is necessary to move from a flat configuration to mode. By taking the first variation of the energy expression with 0 in the first lock mode and setting 180°, we obtain the critical load as This boundary is shown by line in Supplementary Fig. 11. Passing the bifurcation, N n buckles into its first mode. A load increase makes the metamaterial to continue its deformation in the mode until reaches the following limit, which corresponds to the kinematic bifurcation at 90° (line in Supplementary Fig. 11) cos . (S43) To determine the post-bifurcation behavior, we compare the slope of the energy for dissimilar modes at 90°. The path taken by our metamaterial is the one with the lowest gradient.
To obtain the boundaries of and modes, i.e., lines and in Supplementary Fig. 11, we note that the mode is activated when the slope of the energy at 90° is below the slope of the energy in mode. This condition can be written as (S45) The equality signs in Eq. (S45) define the boundary of ( ) with (lines and in Supplementary   Fig. 11). If the compressive forces are equal while passing the bifurcation (line in Supplementary Fig. 11), the metamaterial shifts to mode . The reconfiguration in mode continues until the loads reach a limit defined by the lock angle. This limit is obtained by taking the variation of Eq. (S40) and solving for as 8 sin sin cos sin cos sin cos .

(S46)
This limit is shown by line in Supplementary Fig. 11.

S9.3. Energy analysis of under in-plane biaxial loading
To obtain the reconfiguration behavior of N n we use Eq. (S40). The length of the sides of RUC in Supplementary Fig. 7 can be obtained by using Supplementary Fig. 12 shows the phase diagram of N6. Initially, the metamaterial is assumed to be flat. By compressing it, a certain limit is reached, and the metamaterial buckles to (line in Supplementary Fig. 12). This limit can be obtained by taking the first variation of Eq. (S40) at 180° as
This deformation mode is maintained until the metamaterial reaches 90° (line in Supplementary Fig. 12). The equation of line can be obtained from the first variation of Eq. (S40) by setting 90°. Accordingly, this limit can be expressed as

√3
cos . (S49) This also applies to the case of N4. At the bifurcation point, the ratio of the forces plays a major role in governing the resulting deformation path. The boundaries of , , and modes can be obtained by comparing the slope of the corresponding energy expressions. The smallest slope defines the deformation after the kinematic bifurcation. Therefore, we can write which can be simplified as (S57) S10. Experiments

S10.1. Base material characterization
To characterize the Young's modulus of the base paperboard material we fabricated rectangular samples of dimensions 240 mm × 5 mm (30 mm on each side for the grip) as suggested in literature [7] for a better outcome than that provided by standardized methods, i.e., ISO 1924-2 (see the left-hand-side of Supplementary Fig. 13). The stress-strain response of the base paperboard material of our proof-of-concept prototype was characterized using in-plane uniaxial tensile tests conducted on rectangular samples. Tests were carried out under in-plane uniaxial tension along two main directions, i.e., machine direction (MD) and cross direction (CD) (Supplementary Fig. 13). The results show that the tensile Young's modulus and strength of our paperboard sheets in the MD direction is about 7.9 GPa and 57 MPa, respectively, roughly two times higher than the values obtained in the CD direction. The difference we observe can be attributed to the pulp and paper process which tends to align the network of cellulose fibers in the MD direction. Supplementary Fig. 13 Tensile test specimen dimensions and results of paperboard. Specimen geometry of base paperboard material along with their uniaxial stress-strain responses in the machine direction (MD) and cross (CD) direction.

S10.2. Response of periodic cellular material
Prior to carrying a series of experiments on our material system, we studied the properties of finite size specimens to ensure they are representative of those of their periodic counterparts. We study the role of in-plane and out-of-plane tessellation of our material building block with the goal of determining the smallest material system (minimum number of unit-cells in all three directions) that represents the response of an unbounded periodic material.
To perform our convergence analysis of tessellation level, we examine two representative units, N n and N n , and adopt the notation ( , ) to refer to an in-plane level of tessellation where the primitive unit is tessellated times in the in-plane direction e 1 and times in e 2 . Given the role of tessellation along the three directions, we denote our material system as N n , , and study two representative systems, N n 7,7 and N n 7,7 . The out-ot-plane tessellation is assessed by the number of layers n stacked in the third direction. Our experiments investigate: (A) In-plane tessellation. For a single layer sample N n , determine the minimum number of in-plane unit cells that is representative of the unbounded in-plane tessellation.
(B) Out-of-plane tessellation. For a multilayer specimen N n with minimum number of in-plane unit cells, assess the minimum number of stacking layers, n, that can parallel the response of the unbounded periodic domain.
(C) Comparison between of N n and N n at given levels of tessellation.
(A) Mechanical response of single-layer specimen. Supplementary Fig. 14a shows the primitive unit-cell of the locked unit N n along with the effective area of four tessellation levels (1,1), (2,2), (3,3) and (4,4), each shown by a color ranging from brown to yellow. In Supplementary Figs. 14b, c, the statistical values of Young's moduli and yield strength (as measured in the subset of Supplementary Fig. 14b) of a single layer sample show that by increasing the number of unit-cells, the convergence for both parameters occur at the tessellation level of (7,7). Beyond this level, negligible changes in properties can be observed. (7,7) can thus be considered as the minimum level of tessellation representing the periodic unbounded response.
(B) Mechanical response of multilayered specimens. The role of layer stacking (z-direction) was assessed for a specimen with (7,7) tessellation level ( Supplementary Fig. 14d). The results in Supplementary Figs. 14e, f show that adding n layers results in a convergence of the Young's modulus and strength at n = 6. Beyond this value, only marginal changes can be observed, indicating that the periodic response of N n with tessellation level (7,7) is attained with a minimum of 6 layers stacked in the zdirection.
(C) Comparison of the mechanical response of N n and N n at given levels of tessellation.
Supplementary Fig. 14g shows the primitive unit cell of the locked unit kinematic chain N n along with the manufactured paper sample with a tessellation level (7,7). Supplementary Fig. 14h, i show respectively the Young's modulus and the yield strength of the single layer N n and compare them with those of N n at given levels of tessellation. The results show that a single unit chain, i.e., tessellation level (1,1), has the highest stiffness, but the lowest strength. Upon increasing the tessellation level up to (7,7), the stiffness decreases for both units; on the other hand, for N n the strength monotonically increases, while for N n we observe first an increase and then a decrease. Beyond the tessellation level (7,7) the mechanical properties tested in our experiment do not show noticeable changes ( Supplementary Fig. 14h, i). Hence (7,7) is the minimum level of tessellation that we identify as representative of the unbounded periodic response. strength  * f for given numbers of stacking layers n. g Top view of representative fabricated N n sample with unit chain area in locked mode and tessellation level (7,7). Yellow hexagons represent the top area of our primitive unit cells used for the calculation of the stress values. Compressive Young's modulus h and yield strength i of N n and N n at given tessellation levels. For each tessellation level, five tests were conducted.

S11. Manufacturability
Most of existing reconfigurable metamaterial concepts using origami concepts require a laborious process of fabrication. The reason lies mainly in the non-developable nature of their patterns [8][9][10], which often adds an additional layer of manufacturing sophistication. In several cases, the fabrication of a single threedimensional unit cell requires several steps (cutting, folding, and gluing) [8,9]. The spatial assembly of multiple unit cells is another elaborate undertaking, as each cell needs to be attached one by one in space to its neighbor via edge-to-edge bonding [8][9][10][11]. While 3D printing has been effectively used to fabricate origami-based materials [10,[12][13][14], the high stiffness of the 3D printed hinges has been shown to sacrifice their foldability [10,12,13]. Additionally, 3D printing has been used mainly with soft materials [13,14], thereby bringing to realization prototypes with limited load-bearing capacity.
In contrast, our crease and cut patterns are fully developable. This enables us to use a relatively straightforward process entailing two basic steps of fabrication which can be easily automated. From a single layer of the base material, we introduce in one step: the tessellation of shaped voids and folding lines through laser cutting. After partially folding one kirigami layer along the origami crease pattern, we stack identical layers and bond (here, using glue) their triangular faces. This process does not require the assembly of 3D unit cells in space, rather three-dimensionality is imparted in the flat configuration where the constituent layers are stacked and bonded together. As a result, the flat multilayered crease pattern can be easily folded in space since each layer is geometrically constrained to lay between parallel planes (Supplementary Movie 1). Such characteristics enable to bypass the more laborious procedures currently pursued in the literature.