Zero modes activation to reconcile floppiness, rigidity, and multistability into an all-in-one class of reprogrammable metamaterials

Existing mechanical metamaterials are typically designed to either withstand loads as a stiff structure, shape morph as a floppy mechanism, or trap energy as a multistable matter, distinct behaviours that correspond to three primary classes of macroscopic solids. Their stiffness and stability are sealed permanently into their architecture, mostly remaining immutable post-fabrication due to the invariance of zero modes. Here, we introduce an all-in-one reprogrammable class of Kagome metamaterials that enable the in-situ reprogramming of zero modes to access the apparently conflicting properties of all classes. Through the selective activation of metahinges via self-contact, their architecture can be switched to acquire on-demand rigidity, floppiness, or global multistability, bridging the seemingly uncrossable gap between structures, mechanisms, and multistable matters. We showcase the versatile generalizations of the metahinge and remarkable reprogrammability of zero modes for a range of properties including stiffness, mechanical signal guiding, buckling modes, phonon spectra, and auxeticity, opening a plethora of opportunities for all-in-one materials and devices.

x − ō − ȳ. w(ȳ) and ϕ(ȳ) can be expressed as a superposition of a series of admissible functions: where M and N are the number of modes taken into account, A m and B n are the amplitude of each mode, and AG is the distance between vertices A and G in the undeformed configuration; the overbar convention will now be used in the following to denote the distance between two vertices in the undeformed configuration.In the local coordinate system x − ō − ȳ, the displacement of the vertex A + along the ȳ-axis is w ȳ=AG ; the projection to the x-axis is ϕ ȳ=AG * AD/2.
The absolute displacement of vertex A + in the global coordinate system x − o − y comprises the relative motion in the local coordinate system x − ō − ȳ and the translation associated with triangle EFG: where the rotation matrix R is given by Therefore, the position vector of vertex A + in the deformed configuration is The position vector of vertex C in the deformed configuration is where v C is the displacement of vertex C along the negative y−axis.The position vector of vertex A − in the deformed configuration can then be expressed as where θ ABC is the rotation of the rigid triangle ABC.
After fully defining the deformed configuration of the metahinge architecture, we now formulate its potential energy.The strain energy of the Timoshenko beam is given by where t 0 is the out-of-plane thickness (see caption of Supplementary Figure 1), E is the Young's modulus of the constituent material of the Timoshenko beam, G is the shearing modulus, and ζ is the cross-section coefficient which is selected as 0.856 [1].We introduce an energy term U A + A − to penalize the relative motion between vertices A + and A − , which can be expressed as where η is a sufficiently large parameter ensuring that the vertices A + and A − are kinematically constrained together.Similarly, we use U contact to penalize the rotation of the triangle ABC as where H(•) is the Heaviside function, and σ is a sufficiently large parameter ensuring the hard contact behavior along the normal direction.The external work can be given by Consequently, the total potential of the architecture is The equilibrium equations can then be expressed as where The equilibrium path and the energy landscape are obtained via the Pseudo Acr-length method (Algorithm 1).The above process is implemented in an in-house Matlab code.
Algorithm 1 Pseudo Arc-length method for computing the equilibrium path of the multibody architecture ds ← the arc-length increment.ε ← the numerical tolerance.E ← ∂Π/∂d.Let v C be the control parameter for the numerical continuation.
Let s current store the values of the symbolic vector s; it is initially a null vector.be obtained by examining the null space of its kinematic matrix [2].Our focus is not on the determination of the actual elasticity of the real lattice, but rather on providing a representative assessment which we can obtain by simply assuming that all NN and NNN bonds have an axial stiffness of 1, eliminating the need to distinguish between NN and NNN bonds.The position vector of metahinge i is X i , and the stretching force per unit length in bond j is denoted as b j .We define In the absence of any concentrated forces applied at each metahinge, the equilibrium equations for the metahinge i can be given by By assembling the equilibrium equations of each metahinge, we can derive the equilibrium equations of the entire lattice system as where A is the equilibrium matrix, and b =[b 1 , ..., b N NN +N NNN ] T .Next, the kinematic matrix B can be given by The null space of B is [c 1 , ..., c N 0 +3 ], which comprises N 0 zero modes and 3 rigid-body motions, represented by d i (i = 1, 2, 3).To extract the zero modes (internal mechanisms) from the null space, we apply the following orthogonalization The remaining non-zero c i constitutes the zero modes (internal mechanisms) of the lattice analogy.

tion of biaxial zero modes
To computationally investigate the biaxial zero modes illustrated in Fig. 3c, 3d, and 3e, we resorted to Abaqus 2020 (Dassault Systèmes) and created a Finite Element (FE) model comprising truss elements (T2D2) and axial connector elements (CONN2D2), which represent the nearest neighbor (NN) bonds and the next nearest neighbor (NNN) bonds respectively.Both the truss elements and connectors have an axial stiffness of 1.The deformation response is computed using the General Static solver with geometric nonlinearity.The lattice FE model, boundary conditions, and the corresponding physical specimens are illustrated in Supplementary Figure 2, where distinct colors are used to highlight the areas enclosed by the rotational triangles.In the FE model, the rotation angle of each triangle is evaluated as the average rotation angle of its three constituent edges.By comparing each enclosed area of the FE model and the physical experiment specimen, we can observe a good agreement, a result that validates that the deactivation of a metahinge is equivalent to adding a pair of NNN bonds to its adjacent rotation triangles.
Leveraging the FE model, we obtained the force-displacement relation as well as the evolution of the incremental stiffness versus the displacement, as shown in Supplementary Figure 3.The initial incremental stiffness in all three cases is zero but becomes positive as the displacement increases due to geometric nonlinearity, confirming that the corresponding biaxial zero mode is infinitesimal.
The lattice models in Supplementary Figure 3a and 3c manifest a stiffening behavior as the applied load amplifies; the lattice in Supplementary Figure 3b, however, exhibits a softening behavior once u 1 /l surpasses approximately 0.45.Comparing their incremental stiffness at a modest deformation level, e.g., when u 1 /l ≈ 0.2, it can be observed that the lattice model showing the smallest modal amplitude along the e 1 axis (Supplementary Figure 2a) has the highest stiffness (Supplementary Figure 3a), whereas the lattice model exhibiting the largest modal amplitude along the e 1 axis (Supplementary Figure 2c) is the softest in response to the applied force F 1 (Supplementary Figure 3c).

Supplementary Note 4: Minimum energy path of the periodic Kagome lattice with selectively added NNN bonds
The periodic Kagome lattice with selectively added NNN bonds is shown on the left of Supplementary Figure 4, where its primitive unit cell is highlighted in yellow with an initial area of A.
The nonlinear deformation of the periodic lattice exhibits translational symmetry, and hence we can extract its unit cell for analyzing the minimum energy path (MEP) during state transition.As illustrated in the middle and right of Supplementary Figure 4, the unit cell comprises independent metahinges (red), dependent metahinges (green), independent NN bonds (outlined in blue), and independent NNN bonds (outlined in red); the motion of dependent metahinges is dependent on the corresponding independent metahinges thanks to the translational symmetry.
The initial coordinates of each metahinge are defined in the Cartesian coordinate system x − o − y, and denoted by X i , where the subscript i is the label of each metahinge.The displacement and deformed position of each meta-hinge are u i and x i respectively, and hence we have To eliminate macroscopic rigid-body motions, we assume that the metahinge 5 is fixed at the origin point, and the metahinge 15 is only allowed to slide along the y−axis.We use a 15 × 15 matrix C NN to store the connectivity information of each NN bond and C NNN for each NNN bond.For example, as the metahinge 1 is connected to the metahinge 4 via an NN bond, the elements C NN are 0. The total strain energy of the unit cell U total incorporates the contribution of NN bonds and that of NNN bonds, and can be expressed as (S20) Algorithm 2 Iteratively update the MEP Given the initial guess of the MEP [q 2 , ..., q i , ..., q I−1 ], we first evaluate the driving force f current at the intermediate state points (see more details in Algorithm 3 to know how to evaluate the driving forces).f current ← [Force(q 2 ), ..., Force(q i ), ..., Force(q I−1 )] n max ← the maximum number of iterations.α 0 ← scale factor ε ← tolerance of the driving force for counter←1 to n max do if counter > 1 then if ∥f current ∥ < ε then return [q 2 , ..., q i , ..., ..., q i , ..., q I−1 ] ← [q 2 , ..., q i , ..., q I−1 ]+s current f current ← [Force(q 2 ), ..., Force(q i ), ..., Force(q I−1 )] end for return [q 2 , ..., q i , ..., q I−1 ] Algorithm 3 Evaluate the driving force at state point q i : Force(q i ) To evaluate the driving force at intermediate state q i , we need to evaluate the gradient of the potential energy at this point.F strain ← −∇U total (q i ) Also, we need to evaluate the tangent at state point q i (see Algorithm 3 for details).t ← Tangent(q i ).The projection of F strain with respect to t is F strain|| .F strain|| ← F strain • t The orthogonal component of F strain is F strain⊥ .F strain⊥ ← F strain − F strain|| We assume that there exists a fictitious elastic band connecting a sequence of intermediate states.This fictitious elastic band can prevent each state point from converging to the local energy minimum points during iterations.This fictitious elastic band has a tangential stiffness of k band .The driving force generated by the elastic band is F band .k band ← the tangential stiffness of the fictitious elastic band.F band ← k band * ∥q i+1 − q i ∥ • t − k band * ∥q i − q i−1 ∥ • t The total driving force comprises F strain⊥ and F band .F total ← F strain⊥ + F band return F total Algorithm 4 Evaluate the tangent vector at state point q i : Tangent(q i ) To evaluate the tangent at state point q i , we also need to know the energy information of its two adjacent state points q i−1 and q i+1 .U 1 ← U total (q i−1 ) U 2 ← U total (q i ) U 3 ← U total (q i+1 ) a max ← max( and (U 2 ≥ U 3 or U 2 ≤ U 1 ) then t ← (q i+1 − q i ) * a max + (q i − q i−1 ) * a min else if U 3 ≤ U 1 and (U 2 ≤ U 3 or U 2 ≥ U 1 ) then t ← (q i+1 − q i ) * a min + (q i − q i−1 ) * a max end if end if end if end if t ← t/∥t∥ return t For a two-dimensional lattice comprising N hinge metahinges connected by N NN nearest neighbor (NN) bonds and N NNN next nearest neighbor (NNN) bonds, the infinitesimal zero modes can