Engineering zero modes in transformable mechanical metamaterials

In the field of flexible metamaterial design, harnessing zero modes plays a key part in enabling reconfigurable elastic properties of the metamaterial with unconventional characteristics. However, only quantitative enhancement of certain properties succeeds in most cases rather than qualitative transformation of the metamaterials’ states or/and functionalities, due to the lack of systematic designs on the corresponding zero modes. Here, we propose a 3D metamaterial with engineered zero modes, and experimentally demonstrate its transformable static and dynamic properties. All seven types of extremal metamaterials ranging from null-mode (solid state) to hexa-mode (near-gaseous state) are reported to be reversibly transformed from one state to another, which is verified by the 3D-printed Thermoplastic Polyurethanes prototypes. Tunable wave manipulations are further investigated in 1D-, 2D- and 3D-systems. Our work sheds lights on the design of flexible mechanical metamaterials, which can be potentially extended from the mechanical to the electro-magnetite, the thermal or other types.


Supplementary Note 1. Homogenization and classification of the designed metamaterials A. Numerical-based homogenization method via the Cauchy-Born hypothesis
The homogenization theory [1][2][3]  Therefore, the microscopic displacement deformations at boundary can be expressed as a function of a macroscopic strain field [2,4] ( + ) = ( ) + • (1) where ( ) is the displacement at a node point located on the boundary of the meshed unit cell, and represents the translation vector between two nearby unit cells.
It is worth noting that these constraints automatically guarantee traction continuity on the unit cell's boundary. Next, the strain energy of the metamaterial unit cell 0 and that of an effective medium with the same area/volume can be calculated directly with the help of a commercial FE software. Finally, following the assumption of homogenization stating that the strain energy obtained from a metamaterial's unit cell is equal to that obtained from an effective medium, namely 0 = = 2 : : the components of the effective elasticity tensor can be determined [5].
By recognizing the symmetry of the elasticity tensor, 6 (or 21) independent components for the 2D (or 3D) case are unknown, and we need to specify at least 6 (or 21) different types of macroscopic strain fields to determine the effective elasticity tensor of a 2D (or 3D) metamaterial. Specifically, we specify 6 types of macroscopic strain fields as two uniaxial, one biaxial, one shear, and two combinations of uniaxial and shear strains for the 2D case, while we choose 21 types of macroscopic strain fields as three uniaxial, three shears, three biaxial, nine combinations of uniaxial and shear, and three combinations of shear strains [5].

B. Classification of the extremal metamaterials
It is well known that general Hooke's law states the linear relationship between stress and strain, as = : , where represents the fourth-order elasticity tensor which can be represented as a stiffness matrix c, a positive definite symmetric × matrix ( = 3 for 2D elasticity and = 6 for 3D elasticity). When elastic symmetry is considered [6], general Hooke's law can then be rewritten as where the stress ⃗ ⃗ and strain ⃗ are n-element column matrixes. Mathematically, the stiffness matrix c has real eigenvalues ( ) with their corresponding mutually orthogonal eigenvectors, ( ) . One can have where ( ( ) ) ( ) = , , = 1,2, … , . is the Kronecker delta. Comparing Eq. (3) with Eq. (4), the eigenvalues and eigenvectors of the stiffness matrix are related to the principal stiffnesses and strains, respectively.
Interestingly, if an eigenvalue ( ) goes to zero, the corresponding principal stiffness of the material also turns into zero, and therefore, no stress is generated when deformation is applied along the principal direction. As a result, a zero mode can be identified for the material. By counting the number of zero eigenvalues, the number of existing zero modes can be found accordingly and we can categorize all elastic materials into four types of extremal metamaterials from null-mode to tri-mode for 2D cases or seven types of extremal metamaterials from null-mode to hexa-mode for 3D cases [7].

Supplementary Note 2. Homogenization analysis of 2D metamaterials.
Here, a representative unit cell for 2D tessellations is shown in Supplementary Fig.   1. For each unit cell, the side length and the minimal angle of the rhombus are 0 = 5 mm and 0 = 12 ⁄ , respectively, and the radius at the hinge is 0 = 0.5 mm. For the first configuration, three eigenvalues of the effective stiffness matrix are larger than 1.0e-3 and therefore, no zero mode can be found, which indicates that it is a null-mode metamaterial. For the second configuration, only one eigenvalue (7.0e-5) is smaller than 1.0e-3 and therefore, a uni-mode metamaterial is identified with the only zero mode being shear deformation. For the third configuration, only one eigenvalue (1.4e-1) is larger than 1.0e-3 and therefore, a bi-mode metamaterial is identified with the shear deformation and the uniaxial deformation along the y-direction being related to the two zero modes, respectively. Moreover, this type of bi-mode metamaterial is different from the previously reported bi-mode metamaterial which only supports hydrostatic pressure. For the last configuration, all three eigenvalues are smaller larger than 1.0e-3 and therefore, a tri-mode metamaterial is identified.
8.0e-5 9.1e-4 9.1e-4 As shown in Supplementary Fig. 2a, the 3D unit cell consists of eight cubes and twelve rhomboids. From a mechanism design perspective, one can take the cubes and rhomboids as links, and take the rotation connections between the cube and rhomboid as revolute joints. Hence, the mechanical model of such unit cell in Supplementary  The planar 4-bar linkage is one of the simplest linkages with four links and four revolute joints to form a closed loop. Once one of four links is fixed to the ground, the linkage becomes one-degree-of-freedom (1-DoF) system, i.e., the rotation of any one revolute joint will decide the rotations of the rest three joints kinematically. As a result, the configuration of the linkage is known. For our unit cell, the two linkages in the same directional set always have the simultaneous motion as their four respective joints are collinear, i.e., A1a and A2a are collinear, etc., so these two linkages form a 1-DoF system. There are three such systems along three coordinate axes, respectively. Their motions are independent in each direction and therefore, a unit with decoupled motions of six 4-bar linkages in three orthogonal directions can be obtained subsequently. Hence, the unit has three degrees of freedom.
Since the unit cell is a 3-DoF system, the workspace of the unit cell is distributed in the entire three-dimensional space, and the boundary of the workspace is defined by the maximum rotation angle of the rhomboids in a single direction, which can be tuned by changing the angle of the rhomboids accordingly. Taking the linkage A2a -A2b-C2a-C2b as an example, we can find that when the small angle of the rhombus is /2, x can In order to reduce the anisotropy of the printed sample, the effect of the printing direction is considered and the 4×4×4 sample (in its initial all-CF configuration) is printed along its body-diagonal direction, as shown in Supplementary Fig. 4. Supplementary Fig. 4 The 3D printing setup for the 4×4×4 sample in all-CF configuration to be printed along its body-diagonal direction.
With the help of the numerical-based homogenization method, three configurations of the sample are investigated and are predicted to be tri-mode, tri-mode' (an alternative tri-mode configuration) and penta-mode, as shown in Supplementary Fig. 5. A three-step reconfiguration path is chosen from tri-mode to penta-mode, then to tri-  Additionally, both x-directional compression and tension tests are conducted for the sample in the tri-mode' configuration, as shown in Supplementary Fig. 8a. In this case, the sample is in the CF configuration along the x direction, where the cubes and rhomboids are bonded with adhesives which can be dissolved with degummant when necessary. As shown in the figure, under the two loading methods, the obtained forcedisplacement curves are relatively consistent and linear, and only small difference can be found.
To further evaluate the reversibility of the TPU material under thermal treatment, tensile tests have been conducted on a dog-bone specimen made with TPU, which was repeatedly thermally treated for four rounds. After each round of treatment, which is the same thermal treatment as that given to the 3D tessellation sample, a tensile test was performed and the results of the initially printed specimen and the specimen after four rounds of thermal treatments are shown in Supplementary Fig. 8b. It can be found that the obtained material modulus was softened in the first two rounds, but stabilized afterwards, which indicates that the TPU material is suitable for reversible thermal treatment. Supplementary Fig. 8

Supplementary Note 5. Homogenization analysis of 3D metamaterials.
We study the ten types of 3D metamaterials shown in Fig. 3c in the main text via the homogenization method. Noted that weaken connections are also considered instead of ideal revolute joints, so the corresponding eigenvalues of the effective stiffness matrix are not accurately zero. With the eigenvalues less than 1.0e-3 being recognized as the ones corresponding to the zero modes, the normalized effective stiffness matrix, the calculated eigenvalues, the corresponding zero mode and the classification of each 3D metamaterial are shown in Supplementary Fig. 9.
For the first configuration, i = 0 (i = x, y, z), all eigenvalues of the effective stiffness matrix are larger than 1.0e-3 and therefore, no zero mode can be found, which indicates that it is a null-mode metamaterial. For the second configuration, z = /4 while x = y = 0, only one eigenvalue is smaller than 1.0e-3 and therefore, a uni-mode metamaterial is identified with the only zero mode being shear deformation in y-z plane. For the third configuration, z = y = /4 while x = 0, two eigenvalues are smaller than 1.0e-3 and therefore, a bi-mode metamaterial is identified with the two zero modes being the two shear deformations in y-z and x-y planes. For the fourth configuration, x = y = 0 while z = /8, two eigenvalues are smaller than 1.0e-3 and therefore, an alternative bi-mode (bi-mode') metamaterial is identified with the two zero modes being the axial deformation along the z direction and the shear deformation in y-z plane. For the fifth configuration, i = /4 (i = x, y, z), three eigenvalues are smaller than 1.0e-3 and therefore, a tri-mode metamaterial is identified with the three zero modes being the shear deformations in x-y, y-z and x-z planes. For the sixth configuration, x = 0, y = /8 and z = /4, three eigenvalues are smaller than 1.0e-3 and therefore, an alternative tri-mode (tri-mode') metamaterial is identified with three zero modes being the axial deformation along the z direction and the two shear deformations in y-z and x-y planes.
For the seventh configuration, x = /8 while y = z = /4, four eigenvalues are smaller than 1.0e-3 and therefore, a quadra-mode metamaterial is identified with four zero modes being the three shear deformations and the axial deformation along the x direction. For the eighth configuration, x = 0 while y = z = /8, four eigenvalues are smaller than 1.0e-3 and therefore, an alternative quadra-mode (quadra-mode') metamaterial is identified with four zero modes being the two shear deformations in xy and y-z planes and the two axial deformations along the y and z directions. For the ninth configuration, x = y = /8 while z = /4, five eigenvalues are smaller than 1.0e-3 and therefore, a penta-mode metamaterial is identified with five zero modes being the three shear deformations and the two axial deformations along the x and y directions.
For the last configuration, i = /8 (i = x, y, z), all eigenvalues are smaller than 1.0e-3 and therefore, a hexa-mode metamaterial is identified with all axial and shear deformations being zero modes. It is the 3D version of the 2D tri-mode metamaterial.

Supplementary Note 6. Effective compressive and shear moduli of the ten configurations
We first transform the printed tessellation sample from the penta-mode configuration to the hexa-mode configuration with the fixture C. Sequentially, quadra-mode, tri-mode, bi-mode, uni-mode, null-mode, bi-mode', and quadra-mode' configurations can be obtained with eight fixtures D, E, F, G, H, I and J, respectively. These eight fixtures and their corresponding transformed configurations are shown in Supplementary Fig. 10.
Specifically, the penta-mode and tri-mode' configurations are obtained with fixtures A and B, as shown in Supplementary Fig. 6.
We also conduct compression and shear tests for each configuration and all ten where is the effective stiffness matrix obtained with the homogenization method.
Finally, the experimentally obtained moduli and numerically calculated moduli are normalized with the measured and calculated shear moduli of the tri-mode metamaterial, respectively.

Supplementary Note 7. Effects of different geometrical parameters
In this part, we first investigate the effects of folding angles on the normalized effective moduli along three orthogonal directions to validate the decoupled orthogonal mechanism design. In Supplementary Fig. 11, one can found that the elastic modulus changes only when the folding angle along the same direction changes, which validates the decoupled design. when / ⊥ = 1, the metamaterial is in CE configurations along the three orthogonal directions and therefore, it is a tri-mode with three shear zero modes. When / ⊥ decreases ( / ⊥ ≠ 0), the metamaterial is in PF configuration along z direction and therefore, the metamaterial's compressive deformation along the z direction becomes a zero mode. When / ⊥ = 0, the metamaterial is in CF configuration along z direction and therefore, two zero modes (the compressive deformation along the z direction and the shear deformation in the z-y plane) disappear, as shown in Supplementary Fig. 12a.
When we change the length of diamond bars 0 and the thickness of the hinge ℎ 0 ( / ⊥ = / ⊥ = / ⊥ = 1 being fixed), the tri-mode state (with only three shear zero modes) stays unchanged, as shown in Supplementary Fig. 12b and 12c, respectively. It can be concluded that the number of zero mode depends majorly on the folding angle, while both the thickness of the hinge ℎ 0 and the length of the rhombus 0 have little effect on the zero mode. Since the metamaterial's ability of polarized elastic wave manipulation is based on its reconfigurable effective elasticity, we demonstrate its frequency insensitivity by choosing different harmonic wave frequencies, as shown in Supplementary Fig. 16.

Supplementary Note 9. 2D and 3D wave analysis of transformable metamaterials
Here, we provide detailed analysis on the wave behaviors of the 2D and 3D transformable metamaterials based on the FE models and the equivalent effective metamaterial models.
For wave analysis in 2D metamaterials, the microstructures and the effective elastic stiffnesses of the uni-mode and bi-mode metamaterials shown in Supplementary Table 3 The normalized effective mass density. 0 represents the mass density of the constituent material.
To better understand the unique wave direction manipulations in the 2D and 3D metamaterials, the wave propagation characteristics in the effective homogenous media are investigated. Plane harmonic elastic wave propagations are determined by the Christoffel equation [9]: where is the wave vector, ⃗ is the wave polarization velocity vector, and is the angular frequency. Moreover, and represent the effective mass density and elastic tensor of a metamaterial, respectively. Therefore, an eigenvalue problem forms and its solution determines the dispersive nature of the metamaterial. By specifying a frequency, equi-frequency curves (EFCs) of uni-mode and bi-mode metamaterials can be drawn by calculating components of wave vectors along all directions.
EFCs of the uni-mode metamaterial are shown in Supplementary Fig. 17a. Red arrows normal to the EFCs indicate the wave energy propagation directions. First, when a wave propagates along the x or y direction, only one type of wave can be identified with its energy propagation direction being along the x or y direction. Second, when the propagation direction is away from the x or y direction, two types of waves can be identified with their energy propagation directions being always parallel with the x and y directions, respectively. Such strong anisotropic wave behavior in the uni-mode metamaterial explains the unique 90º wave-splitting results for the 45º L wave beam incidence in Fig. 4b in the main text.
EFCs of the bi-mode metamaterial with PF-y configuration are shown in Supplementary Fig. 17b. First, no wave propagates along the y direction. Second, when the propagation direction is away from the y direction, only one wave can be identified with its energy propagation direction being always parallel with the x direction, which also explains the only x-directional wave propagation result in Fig. 4c in the main text. Supplementary Fig. 17 Equi-frequency curves (EFCs) of the 2D transformable metamaterial.
a. The EFCs of the uni-mode metamaterial. b. The EFC of the bi-mode metamaterial. kx and ky represent two components of the wave vector, and k0 defined as k0 = /√ 0 / 0 , represents the wavenumber of the constituent material.
For wave analysis in 3D metamaterial, equivalent effective models of tri-mode (3) and penta-mode (5) are investigated. Each model is obtained by cutting a corner of a cube along the orange triangle area normal to the body diagonal of the cube (as illustrated in Fig. 4d-e in the main text). L wave incidence with a Gaussian distribution on the blue circle area is applied, and the operating frequency is 10 kHz. Absorbing layers are also applied to all three square ends to eliminate boundary reflections. Since the 3D tri-mode metamaterial is the 3D counterpart of the 2D uni-mode metamaterial, a three-wave-splitting phenomenon can be observed in Fig. 4d in the main text. On the other hand, the 3D penta-mode metamaterial (with CE-z, PF-x, and PF-y configuration) is the 3D counterpart of the 2D bi-mode metamaterial and therefore, the incident L wave only propagates along the z direction, as shown in Fig. 4e in the main text.