Kirigami-based metastructures with programmable multistability

Significance Different from most existing multistable structures whose multiple stable states are achieved through the combinational effect of bistable units, we invent a generic tristable kirigami cuboid. The three stable states have fundamentally distinct geometric configurations and chirality, and the transformation among them can be realized by tension/compression or clockwise/counterclockwise twist. Tessellating the units in series, a family of multistable metamaterials can be constructed, the mechanical behaviors of which are programmable by the unit geometry, the material of the elastic joints, the number of units, and the loading conditions. As a demonstration of the potential applications, a frequency reconfigurable antenna for 5G triple-band communication is developed based on a tristable unit, and the frequency tunability is verified by experiments.


Kinematic Analysis of the Kirigami-inspired Foldable Cuboid
The kirigami is shown in Fig.S1A with four rotationally symmetrical limbs and each limb consists of five hinges. As the mechanism is a parallel mechanism, the truss method (1) is used to determine the number of mobility. For a truss without external forces and consisting of b bars and j joints, the equilibrium equation is At = 0, where A is the 3j × b equilibrium matrix, t is a b × 1 vector of bar axial forces. If r is the rank of the matrix A, the number of mobility is 3 6 m j r    , (S1) The truss form of the mechanism in Fig.S1A can be derived by using bars and nodes to replace edges of facets (the black solid lines in Fig.S1A). For triangular facets, three bars will make a facet rigid and for quadrilateral facets, arbitrary points V1 and V2 out of the facets ABCD and EFGH are introduced. By counting, the truss form in Fig.S1A contains j = 18 nodes and b = 46 bars. When taking the configuration in system 0 ( Fig.S1B), the coordinate of nodes can be obtained (Dataset S3, truss method data) and the equilibrium matrix of the mechanism can be established according to (1). Then, we obtain the rank of the matrix as r = 45. The number of mobility is m = 3. As the mechanism needs three inputs to determine its configuration, it is not easy to analyze and control the motion of the mechanism directly. Here, we consider the motion with rotational symmetry in the mechanism. A deployable configuration shown in Fig.S1A is constructed where all the limbs are symmetric. Hence, their rotated angles are correspondingly equal. Here, one of the limbs with the top and bottom facets is picked up and its corresponding mechanism diagram (Fig.S1B) is also raised. Coordinate systems are   in a limb. The mechanism with four limbs has three independent loops according to the Euler's equation 1 l g n    . (S2) As any two limbs can form a loop, we choose the three loops, L1L2, L2L3 and L3L4 to carry on the further study.
For a single closed loop linkage consisting of n links, the product of the transformation matrices is equal to the 4×4 identity matrix I4, as 21 32 ( 1) 1 4 n n n  Hence, we can obtain the following relationships from the three loops,

T T T T T T T T T T T T T T (S4)
where limb i = 1, 2, 3 and the transformation matrix ( T  T  R  R  T  T   T  T  R  R  T  T  R   T  T  R  R  T  T Here, rotation matrices rotate vectors by an angle  about the x-, y-, or z-axis of a and the translation matrix along vector v is   For the Eq. S4, we can rewrite it as 1 The homogeneous coordinate of the vector 0 0 Besides, we can obtain the homogeneous coordinates of points D in system 5, (S16) Then, we can obtain 2  2  2  2  2  1  1  3  3  2   2  2  2  2  2  1  2  3  1  1  2  1  3   1  1  1  1  cos  cos  sin  cos cos  cos sin  2  2  2  2  1  1  cos  sin  cos  cos sin  cos cos sin sin  2 and its corresponding equation with dihedral angles is 2  2  2  2  2  1  1  3  3  2   2  2  2  2  2  1  2  3  1  1  2  1  3   1  1  1  1  cos  cos  sin  cos  cos  cos sin  2  2  2  2  1  1  cos  sin  cos  cos sin  cos cos sin sin  2 According to Eqs. S12 and S18, relationships of θω vs. φ 3 in Fig.1D and G are derived. The z-axial coordinate of point J in the coordinate system 0 is the height of h, which is   for Fig.1E and H.

The Stored Energy of the System
A. The stored energy of foldable cuboid with torsional springs For a spring hinge with a torsional spring, the potential energy can be derived from where 0  and  are the rest angle and real-time angle of the spring hinge, respectively, K is the stiffness of the torsional spring which is determined by modulus of elasticity of the material E, wire diameter dw, mean spring diameter Dm, number of active coils N and can be calculated from Here, the torsional spring is manufactured by 65Mn Spring Steel with E = 206000 MPa, dw = 1 mm, Dm = 3.5 mm, and N = 5. Its spring constant is supposed to be K = 3.21 N.mm/rad. For the torsional springs of stiffness K 2,4 = K for creases 2 and 4 of four limbs in the kirigami cuboid with the rest angle of the springs φ 20 = φ 40 = 90°, the stored energy of each spring hinge can be derived from   2 2,4 2 20 The total energy of the system with eight hinges can be derived from U 2,4 = 8U whose curve is plotted in Fig.2A, where constant K s0 = 3.21 N.mm/rad. Similarly, for the structure only with torsional springs (stiffness K 1,5 = K) on creases 1 and 5 of four limbs with rest angle, φ 10 = φ 50 = 180°, the total energy can be derived from U 1,5 = 8×½ K 1,5 (φ 1 -φ 10 ) 2 with curve shown in Fig.2B. When these two sets of creases 1, 5 and 2, 4 are replaced by the abovementioned torsional springs, respectively and simultaneously, we will only obtain monostable configuration ④ because the total energy of the system U 1,5 + U 2,4 , as shown in Fig. 2C.

B. The stored energy of thick-panel foldable cuboid with elastic joints
In Fig. 2E, each elastic joint 1 (and 5) can be stretched by An elastic joint 2 (and 4) in Fig. 2E As there are eight elastic joints on 1(and 5) and eight elastic joints 2 (and 4), the whole stored energy of the system can be derived from In the tristable structure with elastic joints in Fig. 2F

Tension Experiment of the Elastic Joint for Stiffness
The elastic joints with elastic sheets of the foldable cuboid are constructed in Fig. 2E. Since the sheet is very thin, its bending stiffness can be negligible. Here, we carried out an experiment of an elastic joint with an elastic sheet under tension. The specimen consists of panels, wedges and elastic sheets with dimensions shown in Fig. S2A. The manufacturing process includes three steps, as shown in Fig. S2B. First, two panels and two wedges are fabricated by a Stratasys Dimension Elite 3D printer with ABS (acrylonitrile-butadiene-styrene), and the elastic sheet (natural rubber latex film, thickness 0.3mm) was cut by a Trotec Speedy 300 laser cutter (power: 70%, speed: 2.5%, Hz: 2000). Second, two panels P t1 and P t2 are symmetrically placed and glued together with the elastic sheet (LOCTITE 401 glue). Finally, the wedges are bolted to the panels.
The setup of the tension experiment is shown in Fig. S2C. An Instron Universal Testing Machine 5982 with a load cell of 100N is adopted for the experiment. The specimen is connected to the two fixtures on the machine through PE (polyethylene) wires. The loading speed is 1mm/min and the tension displacement is 8mm. The displacement and forces are recorded by the data collection system. The deformation process is recorded by a digital camera (Canon 70D).
Under the experiment, we can obtain the stiffness of the joint with the elastic sheet from whose relationships are shown in Fig. S2D (black solid polyline). This black polyline can be fitted by the red line with the equation Fe = 3.2187 y  with the coefficient of determination R 2 = 0.9942. Hence, the stiffness of elastic joints (length 31.04 mm) with latex film (thickness te0.3 = 0.3 mm) is 3.2187 N/mm, and the per length (mm) stiffness is Kep0.3 = 0.1037N/mm 2 . As the stiffness can also be derived from e ee e e where Ee is the modulus of elasticity of the material, te is the thickness of the film, be is the valid length of film on the crease, le is the rest width of the film in the joint, the per length stiffness of latex film with thickness te0.09=0.09 can be derived from K ep0.09 = K ep0.3 te0.09/te0.3 = 0.3K ep0. 3 . The per length stiffness of latex film with thickness te0.4=0.4 mm can be derived from K ep0.4 = K ep0.3 te0.4/te0.3 = 4K ep0.3 /3. Figure. 2E has shown the changes in hinges' positions during the motion. During the motion between ① and ④, hinges 1 and 3 are kept while hinge 2 becomes the hinge 2.

Kinematic Analysis of Foldable Cuboid with Thick Panels
When the motion is in the range between ④ and ⑥, hinge 1 is activated and hinges 2 and 3 are kept. Similar to Fig. S1, the thick-panel structures of the two moving parts are constructed in Fig. S3 with their corresponding mechanisms of a limb for the kinematic analysis. Deployable configuration of the structure and its corresponding mechanism diagram with a limb, which is the thick-panel theoretical model for the range between ① and ④; (C-D) Deployable configuration of the structure and its corresponding mechanism diagram with a limb, which is the theoretical model for the range between ④ and ⑥.
For the motion between ① and ④, its thick-panel structure and corresponding mechanism are shown in Fig. S3A and B. Similar to Eqs.S5a and S5b, the transformation matrix ( 1) i j j  T about the j+1th coordinate system and the jth coordinate system are Considering the symmetric condition 1 from the following elements of the matrix Q1     According to the relationships between dihedral angles and rotated angles For the twist angle of the structure, we can derive it from Eqs. S15 and S16 based on the homogeneous coordinates of points D in system 5, and its corresponding equation with dihedral angles is For the motion between ④ and ⑥, its thick-panel structure and corresponding mechanism are shown in Fig. S3C and D. The transformation matrices are For height h based on point J, we get    

Analysis of Energy Barriers of Tristable Structures
According to Eqs. S24 to S26, the energy barriers of the system between states ① and ④, and   Here,  (Fig. 2E) in state ⑤. From Eq. S32, we can obtain the angle 2,3  is only determined by t 2 and μ in state ③, when the section length a is defined, i.e., U S3 can be programmed by length ratio μ and thickness of panels t 2 without considering the stiffness of elastic joints. As the relationships of angles are implicit functions (Eqs. S32 and S40), 2,3  is the maximum 2  which can be derived from the kinematic analysis of the foldable cuboid under MATLAB function 'max', the energy barrier U S3 can be derived from the code 'F2J_EnergyUs3VSmu_thickness_80.m' in Dataset S3 under varying length ratio μ and thickness of panels t 2 , as shown in Fig. 2J. Similarly, from Eq. S40, we can obtain the angle 1,5  is only determined by t 1 and μ in stable state ⑤, when the section length a is defined, i.e., U S5 can be programmed by length ratio μ and thickness of panels t 1 without considering the stiffness of elastic joints. Similarly, the energy barrier U S5 can be derived from the code 'F2K_EnergyUs5VSmu_thickness_80.m' in Dataset S3 under varying length ratio μ and thickness of panels t 1 , as shown in Fig. 2K. Here, all the elastic joints are supposed to be with latex film (thickness 0.3 mm), the length a = 80 mm, and the real length of the elastic sheets in a joint is supposed to be a 0 = 70 mm, and the stiffness is supposed to be e2, e4 e1, e5 ep0.3 0 Figures 1C and F show φ3 gradually increases from 0° to 180°during the deployed sequence from state ⑥ to state ①, which indicates an arbitrary state of the folded cuboid corresponds to a unique value of angle φ3. Therefore, the angle φ3 can be used to represent the relative positions of the critical states ① to ⑥. As the state ②, the configuration with the maximum h, and state ③, the one with maximum dihedral angles φ2 (= φ4) which corresponds to the maximum energy state, are in the range between state ① and ④, and do not always appear in order, angles φ3's in the state ② and state ③ are chosen to analyze their relative positions under the variation of μ and thickness of panels t2. The relationships among φ32, φ33, μ and t2 are derived from the codes 'F2L_muVSpositionOFstate2and3_3D.m' in Dataset S3 with 'IFDoTThi2.m' in Dataset S3, as shown in Fig. 2L and Fig. S4A-B. From the relationships among φ3 in state ② (the blue surface), φ3 in state ③ (the black surface), and t2 in Fig. S4 A-B, we find the intersection of the two surfaces approaches a straight line which shows the influence of thickness on the relative position of two stable states ② and ③ can be negligible. Then, a foldable cuboid with the panel thickness 5 mm is chosen to show the relative position of state ② and state ③, the state ② is between ① and ③ with μ > 0.388, between ③ and ④ with μ < 0.388, and coincident with ③ with μ = 0.388, as shown in Fig.  S4C ('FS4C_muVSpositionOFstate2and3.m' in Dataset S3).

Torques and Forces of Thick-Panel Tristable Structures with Elastic Joints
According to the derivative of φ3 to φ2 in Eq. S32 with Eq. S47, we obtained The derivative of h to φ2 in Eq. S36 is

A. Fabrication of the foldable cuboids and tristable structures
Based on the geometric parameters in Fig. S5A, the panels of the specimen are fabricated by a Stratasys Dimension Elite 3D printer with ABS, whereas the elastic sheets are cut by a Trotec Speedy 300 laser cutter. Then the specimen is assembled following the process in Fig. S5B. First, one Ptop panel is connected to four P1 panels through four elastic sheets and PE wires for hinges 1 or 5 with LOCTITE 401 glue. Second, four P2 panels are respectively connected to the four P1 panels through four elastic sheets for hinges 2 and 4. Third, four P2 and four P3 panels are respectively joined by Tyvek paper. Fourth, the subassembly of P1, P2, P3, Ptop panels, and that of Pbot and P4 panels are assembled with elastic sheets for hinges 2 and 4. Finally, sixteen wedges are installed at hinges 2 and 4.
The tristable wood structure in Fig. 2F is assembled similar to Fig. S5B with natural rubber latex film adhesive to plywood sheets. The assemblies of two units in Fig. 4 are assembled the same as that in Fig. S5B. For the assemblies of Ut1&Ut2 (Fig. 4C) and Ut1&Ut3 (Fig. 4F and I), the elastic joints are made of latex films (thickness 0.3mm), where the latex films for Ut2 and Ut3 are cut with holes (Fig. S5A) to obtain , / , = , / , =2 and , / , = , / , =2, respectively.
As the parameters μ = 0, 0.2, 0.388, 0.417, 0.5, 0.706, 1, 1.5, 2 are specially selected, there is no fundamental difference in the curves of Fig. 1F-H. For the ideal case with μ = 0 (b = 0, a ≠ 0), axes of hinges 1, 2 and 3 on each limb will intersect at a point, and axes of hinges 3, 4 and 5 on each limb will intersect at another point and there is physical interference of the corner joints at states ① and ⑥. Hence, the SolidWorks® model is constructed, as shown in Fig. S6A. When μ = 2 a prototype with a = 80 mm was fabricated according to the process in Fig. S5B, its stable states are shown in Fig. S6B. Here, a prototype of a smaller scale tristable unit with a = 8mm and μ = 1 was fabricated according to the process in Fig. S5B, in which the panels of the tristable unit are made from carbon fiber board (thickness 0.2mm) and the elastic joints are made from Tyvek paper (thickness 0.142mm). Three stable states of this unit are shown in Fig. S6C.

B. Uniaxial tension experiment for a tristable structure
For the tristable structure with a = 80mm and μ = 0.706, a specimen is fabricated (Fig.  S5B) and a uniaxial tension experiment is carried out. A horizontal testing machine shown in Fig. S7A is utilized to avoid the influence of gravity. The load cell is 50N with an accuracy of 0.5%. Two fixtures are used to connect the specimen to the machine, where the left one is fixed to the load cell and the right one has a rotational degree of freedom about the z-axis (Fig. 1A) of the specimen. A two-step loading process is adopted, one is tensioning the specimen by 136.72 mm from stable state ⑥ to the max height state ② and the other is compressing the specimen by 24.4 mm from state ② to states ①. The loading rate is selected to be 0.2mm/s. The deformation process is recorded by a digital camera (Canon 70D).  Fig. 4 Units of the assembly with a = 80mm and μ = 0.706 are fabricated according to Fig. S5B. The setup of the experiment on a horizontal testing machine is shown in Fig. S7B. The fixed supporter is used to fix the left end of Ut1 and the right supporter is used to support the specimen through a metal shaft to avoid the influence of gravity. For the tension experiment, the load cell was connected to the right end of Ut3 through a rotational fixture by a PE wire (diameter 0.181 mm). The load cell is 50N with an accuracy of 0.5%. The loading rate is 0.2mm/s. For the tension process, displacement of 220 mm was set up and the assembly realise the transformation from ⑥-◇ 6 to ④-◇ 2 via ⑥-◇ 4 and ④-◇ 4 , and overcome the state ◇ 3 of Ut3. Then, for the release process, displacement of -40 mm was setup, the assembly went to state ④-◇ 1 . By tensioning 80 mm, Ut1 overcame the state ③, and the assembly got to state ①-◇ 1 by releasing displacement of -40 mm. The displacement and forces were recorded by the data collection system. The deformation process was recorded by a digital camera (Canon 70D). For the twisting process, the setup is similar to that in Fig. S7B, where the rotation is applied to the right end of the assembly through a cylinder with aramid wires and metal shafts connected to the assembly (Movie S5).

Simulation for the Tristable Structures and Their Assemblies
SolidWorks® (Dassault Systèmes S.A, France) (3) is used to simulate the single unit and the assemblies. For a single unit, the panels are modeled as rigid bodies. Joints 1, 2, 4 and 5 are modeled as wire springs with zero original length and identical stiffness to that of the elastic sheet, whereas joint 3 is an ideal revolute joint. For the two-layer assembly, a cuboid with a thickness of 15 mm is used to connect the two units to avoid physical interference in the motion process. Solid contact is applied with setting Dry and penetration 0.001mm. During deformation, the Pbot is fixed and a linear or rotary motor is applied to the Ptop to realize tension or twist. In addition, a loading rate of 0.2mm/s or 0.0115 rad/s is used to ensure a quasi-static process.
The single unit with μ = 0.706 and a = 80mm is first simulated, and the numerical results are presented together with the theoretical and experimental ones in Fig. 3E. A very good match is obtained, thus validating the theory and numerical model. Subsequently, the assemblies in Fig. 4 with different geometries and loading conditions are analyzed numerically using the same approach. Comparison of the theoretical, numerical, and experimental results are respectively presented in Fig. S8A-F. For each assembly, the numerically obtained transformation of stable states is identical to those obtained theoretically and experimentally. And the numerically predicted peak forces or torques at the transformation points match with the theoretical values very well. Therefore it can be concluded that the theoretical model can accurately predict the transformation process and required peak loads.

A. The tessellation of two units with different length ratio
To discuss the variation of length ratio μ, we choose an assembly of two tristable units Ut7 (right-handed, μ = 1) and Ut1 (left-handed, μ = 0.706), where the thickness of latex film on the hinges are identical to 0.3mm. The schematic diagram of the assembly and curves of forces vs. h/a are shown in Fig. S9 which indicates 0 . Here, the transformation from stable states ◇ 6 -⑥ is analyzed with tension ( Fig. S9B). Because the length ratio is different in the Ut7 and Ut1, the energy barriers of the two units are different which leads to the transformation path in the order of states ◇ 6 -⑥, ◇ 4 -⑥, ◇ 4 -④, ◇ 4 -① and ◇ 1 -① (Fig. S9C, Movie S6). Here, by changing the length ratio μ, the stable state ◇ 4 -⑥ is obtained, which is different from the transformation path in Fig. 4I.

B. A programmed tessellation of three units with multiple stable states.
A multistable tessellation of three tristable units (left-handed Ut4 and right-handed Ut5 with μ = 0.706, right-handed Ut6 with μ = 1) is carefully programmed through parameters in length ratio μ, elastic joint stiffness, cuboid chirality, loading modes, as shown in Fig. 5. For hinges 1 and 5, latex films with thickness 0.4 mm, 0.3 mm (cutting holes, see Fig. S5A), and 0.3 mm are used to Ut4, Ut5, Ut6, respectively. For hinges 2 and 4, films with thickness 0.3 mm, 0.4 mm, and 0.3 mm (cutting holes, see Fig. S5A) are used to Ut4, Ut5, Ut6, respectively. Then the stiffness in Fig. 5A is obtained, which determine the magnitude of critical forces/torques with length ratio μ according to SI appendix, section 7, and relationships between the stable states and the regions of normalized forces/torques during the transformation from an arbitrary stable state in the assembly are obtained, as shown in Fig. 5B. The transformation path of 12 stable states (Fig. 5C) is obtained by controlling input rotations or tensions (Movie S7). The setup of the rotation and tension process is the same as that shown in Fig. S7B. The fixed supporter is used to fix the left end of Ut4 and the right supporter is used to support the specimen through the metal shaft to avoid the influence of gravity.

Frequency Reconfigurable Antenna A. Design and fabrication
The antenna is designed using ANSYS High Frequency Structure Simulator (HFSS) according to the designed dimensions shown in Fig. S10A. The reconfigurable antenna includes two parts, antenna substrates with feeding network and radiator and a tristable metastructure, which is fed by an SMA connector for a coaxial cable. The FR4 with relative permittivity (εr) of 4.4, the thickness of 2 mm and loss tangent of 0.02 is selected as the antenna substrate. The antenna was fabricated by LPKF printed circuit broad (PCB) prototyping machine. The tristable metastructure was fabricated by Stratasys Dimension Elite 3D printer with ABS based on the dimensions in Fig. S10B, and the elastic sheet (length 24.2 mm and 50 mm, width 9 mm, thickness 0.3 mm) was cut by a Trotec Speedy 300 laser cutter. The assembly of the tristable metastructure is similar to that in Fig. S5B. The antenna and tristable metastructure are glued by 3M VHB double sticky tape.

B. Simulation
The antenna is simulated in HFSS. During the simulation, an airbox is employed and all six faces of the airbox are set as radiation boundaries. The distance between the antenna and radiation boundaries is around half wavelength at each resonant frequency, as shown in Fig. S10C.

C. Experiment
The setup of the electromagnetic experiment is shown in Fig. S10D. The antenna was measured in a near-field antenna pattern measurement chamber and connected to Rohde & Schwarz ZNA24 by the coaxial cable (Fig. S10E). By changing the state of the specimen and connecting the adjacent radiators with copper tape (3M 1181), frequencies can be measured and results are shown in Fig. 6D. Experiment and simulation radiation patterns at xz-and yz-planes of the three stable states resonant frequencies are exhibited in Fig.  S10F. Measurement agrees with the simulations well.  Movie S1. Folding sequence of the right-handed kirigami cuboid.
Movie S4. The tension experiment of a tristable structure with μ = 0.706.
Movie S5. Transformation path of stable states about assemblies of two units under twist in Fig.  4C and Fig. 4F.
Movie S6. Transformation path of stable states about assemblies of two units under tension in Fig. 4I and Fig. S9C.
Movie S7. Transformation path of stable states about the assembly of three units.
Dataset S1. The data about the experiment of an elastic joint in Fig. S2D Dataset S2. Multistable_structure_Datas_0.706 The data about the experiment and simulation of multistable structures with μ = 0.706 in Fig. 3E, Fig. S8. The data about the experiment and simulation of the antenna in Fig. 6.