# Oligomodal metamaterials with multifunctional mechanics

See allHide authors and affiliations

Edited by Katia Bertoldi, Harvard University, Cambridge, MA, and accepted by Editorial Board Member John A. Rogers April 2, 2021 (received for review September 9, 2020)

## Significance

Mechanical metamaterials are man-made materials with extraordinary properties that come from their geometrical structure rather than their chemical composition. For instance, they can be engineered to be extremely light and stiff; to shrink sideways when compressed, instead of expanding; or to exhibit programmable shape changes. Such properties often rely on zero-energy modes. In this work, we created a class of mechanical metamaterials with zero-energy modes that can exhibit multiple properties at the same time within a single structure. In particular, we created a metamaterial that can either shrink or expand on the side when compressed, depending on how fast it is compressed. These metamaterials could lead to novel adaptable devices for, for example, robotics and energy absorption applications.

## Abstract

Mechanical metamaterials are artificial composites that exhibit a wide range of advanced functionalities such as negative Poisson’s ratio, shape shifting, topological protection, multistability, extreme strength-to-density ratio, and enhanced energy dissipation. In particular, flexible metamaterials often harness zero-energy deformation modes. To date, such flexible metamaterials have a single property, for example, a single shape change, or are pluripotent, that is, they can have many different responses, but typically require complex actuation protocols. Here, we introduce a class of oligomodal metamaterials that encode a few distinct properties that can be selectively controlled under uniaxial compression. To demonstrate this concept, we introduce a combinatorial design space containing various families of metamaterials. These families include monomodal (i.e., with a single zero-energy deformation mode); oligomodal (i.e., with a constant number of zero-energy deformation modes); and plurimodal (i.e., with many zero-energy deformation modes), whose number increases with system size. We then confirm the multifunctional nature of oligomodal metamaterials using both boundary textures and viscoelasticity. In particular, we realize a metamaterial that has a negative (positive) Poisson’s ratio for low (high) compression rate over a finite range of strains. The ability of our oligomodal metamaterials to host multiple mechanical responses within a single structure paves the way toward multifunctional materials and devices.

Flexible metamaterials use carefully designed arrangements of deformable building blocks to achieve unusual and tunable mechanical functionalities (1). Such mechanical responses rely on on-demand deformation pathways that cost a relatively low amount of elastic energy. A useful and widely applicable paradigm for the design of such pathways leverages the limit in which their elastic energy is zero—these pathways then become mechanisms or zero-energy modes. Flexible metamaterials based on such principle are, so far, either monomodal (Fig. 1*A*) or plurimodal (Fig. 1*C*). On one hand, monomodal metamaterials feature a single zero-energy mode and a single functionality (2⇓⇓⇓⇓⇓–8), which is typically robust and easy to control with a simple actuation protocol, that is, a protocol that requires a single actuator, for example, uniaxial compression. On the other hand, plurimodal metamaterials feature a large number of zero-energy modes, which increases with system size (9, 10). The presence of these multiple zero modes offers multiple possible functionalities in principle, but they are hard to control in practice; that is, they require complex actuation protocols—protocols that require more than one actuator—for successful execution (9). The challenge we address here is whether it is possible to find a middle ground between monomodal and plurimodal metamaterials. In other words, can we design and create metamaterials that have more than one zero-energy mode, but not a number that grows with system size? For convenience and clarity, we term such metamaterials *oligomodal* (Fig. 1*B*). Could oligomodal metamaterials be actuated in a robust fashion with a simple actuation protocol (Fig. 1*B*)? Could oligomodal metamaterials host distinct mechanical properties within a single structure?

## Combinatorial Analysis

To first design oligomodal architectures, we turn toward combinatorial metamaterials, which is a particularly fruitful paradigm for the design of advanced mechanical functionalities (4, 11, 12). In combinatorial metamaterials, the structural complexity is reduced to a discrete design space, typically, by controlling the orientation of the constitutive unit cells. Such discreteness makes the design space much easier to explore and has recently been leveraged to create nonperiodic metamaterials with shape-changing (4, 11) and topological properties (12), yet only for single zero-energy mode metamaterials so far. Here, we generalize this combinatorial approach to metamaterials with more than one zero-energy mode. To this end, we introduce a unit cell (Fig. 2*A*) with two zero-energy modes (Fig. 2 *B* and *C*). A central ingredient in combinatorial design is that the unit cells have fewer symmetries than that of the lattice they live on. In the present case, the unit cell does not have any planar rotational symmetry. Thereby, it can be tiled in four different orientations on a square lattice (Fig. 2*D*), in turn, allowing us to construct a wide range of combinatorial metamaterials.

We first analyze the combinatorial landscape of the metamaterials and notice that the design space is extremely large; that is, for an *E*), which is a very small fraction of the design space; we discuss an example of quasiperiodic tilings in *SI Appendix*. Any choice of four orientations defines a *F*). A numerical study of the kinematics (*Materials and Methods*) then reveals a rich spectrum of zero-energy modes. This spectrum contains the three classes of metamaterials defined in Fig. 1. For system sizes *G* and *H*, red) have a single zero-energy mode; oligomodal tilings have a constant zero-energy mode number larger than one (Fig. 2 *G* and *I*, green); and plurimodal tilings have a number of zero-energy modes that increases linearly with n (Fig. 2 *G* and *J*, blue).

All three classes of tilings share a common zero-energy mode, which is the rotating-squares mechanism (13), which spans the size of the system, and where the unit cells have the same deformation in an alternating manner with a lattice spacing periodicity of two (Fig. 2*H*). Monomodal tilings only have this zero-energy mode, while oligomodal and plurimodal tilings have additional zero-energy modes. In particular, plurimodal tilings include most periodic tilings and feature modes that are typically localized along lines (Fig. 2*J* and *SI Appendix*), which we term *line modes*. The existence of line modes in a periodic tiling implies that such a tiling is plurimodal, and explains why the number of modes grows with the size of the sample n. Such lines modes have, in fact, previously been identified in highly symmetric lattices (14, 15), and highly localized deformations related to line modes have been observed in foams under compression (16⇓⇓⇓⇓⇓–22). We discuss a particularly intriguing example of nonperiodic plurimodal tilings, in *SI Appendix*, namely, a quasi-crystalline tiling that gains a single zero-energy line mode when its side length is doubled, thus leading to a number of modes that scales logarithmically with the size of the sample n.

## Oligomodal Tilings

In contrast to the highly localized line modes, the additional zero-energy modes present in oligomodal tilings are markedly different: They typically exhibit deformations covering the whole sample (Fig. 2*I*). By definition, their number remains constant with increasing system size. There is, in principle, no upper bound on the number of modes in an oligomodal tiling, but, by construction, this number of modes in the tiling cannot exceed the number of modes in the supercell. These spatially extended deformations, which, to the best of our knowledge, have not been reported before, offer the prospect of unconventional bulk properties.

In the following, we focus on the simplest oligomodal configuration that we uncovered: the bimodal tiling of Fig. 3*J*, which hosts two bulk zero-energy modes irrespective of system size. In order to understand the nature of these zero-energy modes (and, in fact, also to interpret and classify all three classes of tilings; see *SI Appendix*), we introduce a vertex representation that maps out angle deformations onto a directed graph, that is reminiscent of vertex models found in two-dimensional (2D) statistical physics (23⇓–25) (*Materials and Methods*). This representation allows for a simple graphical representation of single-cell deformations. Each vertex corresponds to a unit cell, each edge corresponds to a hinge, and the number of arrowheads on each edge quantifies the hinge deformation (Fig. 3 *A*, *D*, and *G*). In this representation, the geometric constraints can be formulated as follows: Only linear combinations of the two vertices shown in Fig. 3 *D* and *G* are allowed. These geometric constraints are propagated from one vertex to its neighbors via the edges, through a combinatorial search. Unlike the standard structural analysis computation that we used to predict the number of zero-energy modes shown in Fig. 2*G* (*Materials and Methods*), the vertex model provides direct insight into the spatial shape of the possible soft deformation modes. Performing such an analysis on the oligomodal tiling shown in Fig. 3*J* reveals two distinct system-spanning zero-energy modes. First, we find the rotating-squares mode (Fig. 3*K*) as expected, using the first unit cell mode of Fig. 3 *D*–*F* only. Second, we find a new bulk bidomain zero-energy mode, which features a symmetry axis going from southwest to northeast, splitting the deformation into two domains (Fig. 3*L*), combining both unit cell modes of Fig. 3 *D*–*I*. This mode spans the size of the system, with a constant gradient of deformation along the southeast to northwest axis. Therefore, for a tiling that would be twice as large, the gradient of deformation would the same, and the deformation would be twice as large. If the tiling were rectangular, the gradient of deformation would remain the same, oriented in the southeast to northwest axis. In the remainder of the paper, we use the oligomodal geometry with two distinct zero-energy modes to create flexible metamaterials with two distinct soft deformation modes. We demonstrate that we can selectively actuate these soft deformation modes—and thus obtain two distinct effective mechanical properties—by using either boundary texture or strain rate-dependent viscoelasticity.

## Selective Actuation by Textured Boundaries

First, we demonstrate that sets of suitably chosen textured boundary conditions under simple uniaxial compression allow preferential exciting of a zero-energy mode of choice while frustrating the others. Using multimaterial 3D printing, we create a metamaterial made of *M*), following the oligomodal tiling of Fig. 3*J*. The 3D-printed unit cells are made of rigid curved struts and compliant hinges (Fig. 3*B*; see *Materials and Methods* for details), such that their shapes coincide with the idealized unit cell (Fig. 3*C*). Crucially, the 3D-printed unit cells have two soft deformation modes (Fig. 3 *E* and *H*) that closely mimic the zero-energy modes of the idealized unit cell (Fig. 3 *F* and *I*). Following the deformation prescribed by the vertex representation (Fig. 3 *K* and *L*), we then compress our metamaterial using two sets of eight indenters on each side, that apply the load on every other unit cell and where the bottom set of indenters is either antialigned (Fig. 3*N*) or aligned (Fig. 3*O*) with the top set of indenters. To quantify the deformations of each unit cell, we use particle tracking (OPENCV [a computer vision library] and Python) and custom-made tracking algorithms to quantify the flattening f and orientation ϕ with respect to the horizontal of each pore and calculate the polarization (26) *N*) or the bidomain mode (Fig. 3*O*) by changing the location of load application (see also Movie S1). This example demonstrates that the zero-energy modes of oligomodal metamaterials can be selectively actuated by tuning the boundary texture. In the following, we build on this result to achieve mode selection under a simpler actuation protocol that does not resort to boundary texture.

## Selective Actuation by Dissipation

To achieve this goal, we leave the realm of static responses and harness viscoelastic dissipation to tune the dynamical response of our metamaterial (26⇓–28), previously applied to alter the direction of a mechanical mode (26, 28). Instead, we now harness viscoelasticity to switch between modes altogether. To this end, we use the oligomodal metamaterial design of Fig. 3*J*, which we rotate by *A*–*D*). With this rotation of the lattice, we expect the domain wall seen in Fig. 3*O* to rotate by *B* and *D*), and thus to become perpendicular to the loading axis as well as to provide mirror symmetry about the transverse axis.

The selective actuation of modes via viscoelasticity is possible because the two modes—the rotating-square mechanism (Fig. 4 *A* and *C*) and the bidomain mode (Fig. 4 *B* and *D*)—involve distinct sets of hinges. While, in the rotating-square mechanism, all hinges but the one at the center of the supercell are actuated (white hinges in Fig 4*A*), in the bidomain mode, all hinges including the one at the center of the supercell are actuated (white and black hinges in Fig 4*B*). Therefore, in order to favor the actuation of the rotating-squares mechanism (respectively, the bidomain mode), the central hinge should be stiffer (respectively, softer) than the other hinges. Conversely, in order to favor the actuation of the bidomain mode, the central hinge should be softer than the other hinges. To obtain a selective actuation of the modes, we tune the stiffness of hinges by making them viscoelastic. Specifically, we use 1) viscoelastic hinges—whose stiffness increases with the strain rate at which they are actuated—for the hinges that are actuated in the rotating-squares mechanism (white hinges in Fig 4 *A* and *B*) and 2) elastic hinges—whose stiffness does not depend on the strain rate at which they are actuated—for the hinge that is only actuated in the bidomain mode (black hinge in Fig 4*B*).

Next, we prove the feasibility of the concept described above with finite element simulations (Abaqus; see *Materials and Methods* and *SI Appendix*). In order to selectively actuate one mode or the other, we use a combination of stiff material for the rigid parts and soft material for the compliant hinges. In particular, we dress our metamaterial with two types of elastic and viscoelastic hinges (see Fig. 9*B*) whose properties correspond to (visco)elastic elastomers we will use in the experiments (Fig. 5 *A* and *B*; see *Materials and Methods* for details about the simulations). For slow actuation rates, both types of compliant hinges have a very similar torsional stiffness (*Materials and Methods*). When the actuation rate increases, the viscoelastic hinges stiffen dramatically: They become more than 5 times stiffer than the elastic hinges (*Materials and Methods*). As expected, we find that, on the one hand, compressing the metamaterial at a low strain rate with respect to the relaxation timescales of the viscoelastic material of *E* and *G*). On the other hand, compressing the metamaterial at a large strain rate with respect to the relaxation timescales of the viscoelastic material of *F* and *H*). A closer look at the central hinge of the supercell (Fig. 4 *E* and *F*, *Insets*) reveals that this hinge is not actuated for low strain rates (Fig. 4 *E*, *Inset*) and actuated for large strain rates (Fig. 4 *F*, *Inset*), which is consistent with the mechanism of mode selection described above.

Moreover, we confirm that the mode selection is applicable across a wide range of strains, including at large strains (10% strain). Both small and large strain rates initially buckle at a strain of *G–I*), as can be seen from the presence of two domains of opposite polarization Ω in the top and bottom of the sample (blue and red regions in Fig. 4 *G* and *H*), which then builds up an average polarization that remains close to zero (light red region in Fig. 4*I*). This first deformation mode is due to the fact that the sample has been designed with a small imperfection aimed at favoring the bidomain mode (see *SI Appendix*). At low strain rates, however, the sample ultimately deforms into the rotating-squares mode, as can be seen from the uniform distribution of polarization (red region in Fig 4*G*), which then builds up a large average polarization (red region in Fig 4*I*). A systematic scan of various strain rates (Fig 4*I*) reveals the following: 1) We find a sharp transition between the actuation of the rotating-squares mode (at low strain rates) and the actuation of the bidomain mode (at large strain rates). While the fact that strain rate allows selection of the mode was expected from our arguments above, the fact that the transition is sharp is more surprising. At this transition, the deformations abruptly snap from the bidomain mode toward the rotating-squares mode, as can be seen from the drastic change of polarization. 2) We find that the onset of buckling toward the rotating-squares mode is shifted to larger values of strains at larger strain rates.

This shift is the result of a relatively lower stiffness of the bidomain actuated elastic hinges with respect to the viscoelastic hinges. Hence, the imperfection favoring the bidomain mode is more difficult to overcome at higher strain rates, and buckling occurs at larger strains.

We use the numerical results described above for the design of our experiments, and we fabricate a metamaterial with plastic rigid parts, viscoelastic hinges, and elastic hinges (Fig. 5 *A* and *B*) using a combination of 3D-printing, molding, and casting techniques (see *Materials and Methods*). We compress the metamaterial slowly at a strain rate of *C*). The polarization Ω, overlaid on the experimental pictures, allows clear visualization of the overall negative polarization corresponding to the rotating-squares mode (see also Movie S2). We also compress the metamaterial fast at a strain rate of *D*). Indeed, the overlaid polarization clearly splits into two domains separated by a relatively undeformed region. Tracking the polarization in Fig. 5 *E* and *F* for small and large strain rates, respectively, shows that the simulations and experiments are in qualitative agreement. First, the applied bidomain mode imperfection (see *Materials and Methods* and *SI Appendix*) always induces a bidomain mode deformation, with a widening polarization range, centered around *E*, while the symmetry of the bidomain mode keeps the average polarization close to zero at all strain values in Fig. 5*F*. However, in both cases, it appears that buckling occurs earlier numerically than experimentally. This could be because the numerical model neglects the hinge bonds’ flexibility and is therefore more prone to buckling, or the experimental imperfections favor the bidomain deformation less.

In conclusion, the combination of oligomodal architectures and viscoelastic hinges provides a robust strategy to achieve mode selection by strain rate, which can be designed with the help of finite element methods. Employing an oligomodal tiling is crucial to achieve such controlled mode selection; indeed, the many additional localized modes present in plurimodal tilings could introduce unwanted degrees of freedom that perturb the mechanical response.

## Strain Rate-Dependent Poisson’s Ratio

In this last section, we ask whether the selection of distinct modes can also lead to selection of distinct mechanical properties. Specifically, we probe whether the rotating-squares mechanism and the bidomain modes lead to distinct values of the Poisson’s ratio, and, therefore, whether one can obtain distinct values of the Poisson’s ratio using strain rate. Intuitively, we expect the presence of the domain wall to frustrate the overall lateral deformation and therefore to induce a larger Poisson’s ratio for the bidomain mode than for the rotating-squares mechanism—see *Materials and Methods*.

We measure, both experimentally and numerically, the Poisson’s ratio as a function of the strain ε for the same different strain rates as before. While we always find a positive Poisson’s ratio in the prebuckling phase, we find that compressing at a small strain rate leads to a response with a negative nonlinear Poisson’s ratio at strains *A*, blue), while compressing at large strain rate induces larger values of the Poisson’s ratio (Fig. 6*B*, blue). To obtain more insight on these results, we measure the Poisson’s ratio as a function of the strain ε and the strain rate *A* and *B*, orange and Fig. 6*C*) and find that the Poisson’s ratio remains positive in the prebuckled region, while it becomes negative in the buckled region. The rotating-squares mechanism and the bidomain mode both exhibit negative Poisson’s ratios, yet the Poisson’s ratio is larger and therefore less negative in the case of the bidomain mode, which is consistent with the experiments. Notably, the Poisson’s ratio in the rotating-square mode shows larger spatial dispersion; this is a signature of the snapping transition that occurs when the deformation snaps from the bidomain mode into the rotating-square mode. Interestingly, as the buckling strain is pushed to larger values for larger strain rates, there is a range of strains

## Discussion

To conclude, we have introduced a class of materials, oligomodal metamaterials, that strike a delicate kinematic balance that allows them to exhibit a limited number of zero-energy modes. Based on these zero modes, we have created flexible metamaterials, whose deformations can then be selectively and robustly actuated using relatively simple loading protocols such as uniaxial compression, thus providing an example of multifunctional mechanics.

Our approach is fundamentally different from and is complementary to recent developments based on poroelasticity (29) and magnetoelasticity (30). In these works, multifunctionality is achieved through elaborate actuation multiphysics couplings and via loading in the bulk: The loading is exerted on the bulk by external fields, such as a chemical potential or a magnetic field. This is different from our case, where the deformation modes are selected in a purely mechanical fashion, and where the loading is exerted at the boundaries. These approaches could, in fact, benefit from employing oligomodal architectures, in order to reduce the number of involved degrees of freedom to the minimum required by the desired functions.

Our approach is broadly applicable for the design of flexible metamaterials with low-energy strain pathways. In particular, our approach could be further generalized to other types of unit cells, to other lattices, and to 3D tilings and applied to more types of functionalities such as chiral responses (31⇓–33), shape morphing (4, 11, 34), enhanced energy absorption (26, 35) and nonlinear wave propagation (36, 37). In particular, we expect the switchable nature of oligomodal metamaterials to enable the control of deformation pathways, which could be potentially useful for dissipation of energy. Moreover, there is no reason to assume that a zero-energy mode-based method as we use is the only way to design flexible oligomodal metamaterials. Alternative analysis and design methods that can induce clear modal separation, such as topology optimization, may be used instead. Also, our design space, mostly restricted to periodic tilings of *SI Appendix* for an example), or by harnessing defects (12, 38). Finally, data-driven approaches such as machine learning could provide a powerful alternative to better understand the structure–property relationship (33) and to rationally design oligomodal tilings. Additional open questions that remain are how to effectively design large deformation responses, which variety of multifunctionality can be achieved, and how to produce these materials on a large scale effectively.

## Materials and Methods

### Kinematics of Multimode Combinatorial Metamaterials.

Repeating our bimode cell with arbitrary orientations on a square lattice yields a large variety of configurations. To understand how the balance between constraints and degrees of freedom gives rise to zero-energy modes and states of self-stress, we use the Maxwell–Calladine index theorem (14). It provides a relation between the number of states of self-stress *A*). Constraining displacements to the plane, this yields**2** implies a lower bound for the number of states of self-stress of **2** provides a negative lower bound for *G*, we calculated the nullity of the compatibility matrix numerically using the LAPACK (Linear Algebra PACKage) driver xGESDD (a divide-and-conquer singular-value decomposition method).

### Vertex Model.

The approach described above has, however, a drawback: It yields arbitrary superpositions of the zero-energy modes, potentially providing very little insight on the nature of the involved deformations. To complement this standard method, we therefore introduced an approach in *Oligomodal Tilings*, based on graphs with directed edges. We now derive this approach in detail, starting from the fully nonlinear geometric constraints of the primitive cell. They can be encoded in three trigonometric equations, namely,*B*. To construct the directed graph model, we perform a few algebraic manipulations and linearization. First, we change variables and consider

Then, we linearize Eqs. **4** and **5** in α, β, γ, δ, and ε. This linearization around the rest position of our primitive cell yields*A*, *D*, and *G*. We draw the cell as a vertex with five edges, one for each hinge, and make use of the fact that only integer coefficients appear in Eq. **7** to draw arrows on the edges. Setting *D* (Fig. 3*G*). Any linear combination of those two configurations also produces a compatible vertex. The final ingredient in our model consists in remarking that every hinge is also subject to angle conservation, allowing us to concatenate single-cell graphs to obtain a directed graph describing the full tiling (Fig. 3 *K* and *L*).

Importantly, by an angle-fixing argument, the model could be trivially extended to six, seven, and eight edges per vertex by taking our existing vertex and adding a diagonal edge in either of the three available corners. If we then fix the value of one of these diagonal edges, we reduce the problem to the pentagonal case. Since this yields three, four, or five independent vertices, it exhausts the number of available zero-energy modes. In principle, this linearization procedure is also applicable to other primitive cells, and provides a convenient graphical tool in the cases for which the resulting coefficients are integers.

### Finite Mechanisms and Infinitesimal Zero-Energy Modes.

We discuss here the nature of the zero-energy modes determined above. Two classes of zero-energy modes can exist (14): 1) the finite mechanisms, in which the zero-energy deformations can persist in the nonlinear regime and lead to large rotations without stretching or compressing any bars; 2) infinitesimal zero-energy modes, in which the deformations are only zero-energy in the linear range. Further deformations in the nonlinear regime typically cost elastic energy through elongation of the rods. As demonstrated in our experiments, such infinitesimal modes translate to soft deformation modes in practice.

A visual inspection of the kinematics of the elementary modes in *SI Appendix*, Table S1 reveals that only the rotating-squares mode is finite. All of the other zero-energy modes require length change of the rods to be extended into the nonlinear regime and, as such, are infinitesimal.

### Sample Design and Fabrication.

The 16*M*–*O* are designed following the tilings shown in Fig. 3*J*. The samples’ dimensions are *J*, rotated by

For these two types of samples, special care has been devoted to the design of the details of the unit cells: Imperfect compliant hinges, experiencing shear and tension instead of pure bending/torsion, can affect the linear and nonlinear deformation of mechanisms significantly (39). For the samples of Figs. 3 and 5, we negate these effects in two ways. First of all, the central unit cells are made from a much more rigid material than the elastic hinges, such that the central unit cells do not deform during loading. Next, we use tapered rigid parts as seen in Fig. 8. A central narrow section between the rigid sections greatly increases resistance to shear and tension while offering only a negligible resistance to bending.

When the actuation is aligned with the primitive lattice, the bidomain mode leads to shear–compression coupling. However, such a mode in a *SI Appendix*). Therefore, we choose to rotate the sample by *C*. However, this preference is tuned by designing an imperfection (40) favoring the bidomain mode of Fig. 5*D*, which has been designed using nonlinear viscoelastic finite element methods (Abaqus; see *SI Appendix*). Furthermore, as the samples are tested horizontally, the bottom of the sample features small pockets with 5-mm steel ball bearings to reduce the effects of friction with the bottom surface.

We produce the samples in Fig. 3 using additive manufacturing using a Stratasys Objet500 Connex3 3D printer. The hinges are produced from a highly viscoelastic rubber-like PolyJet Photopolymer (Stratasys Agilus 30), which is cocured with the central parts, which are made from a much stiffer Photopolymer (Stratasys VeroWhitePlus). The square hinges and central parts of the sample of Fig. 5 are produced in the same way. However, the triangular hinges of the sample of Fig. 5 are cast from an elastic two-component silicon-based rubber (Zhermack Elite Double 32), and are then glued to the sample using Wacker Elastosil E43 glue. We anticipate that advances in additive manufacturing, using, for instance, silicone resins, will allow single-step manufacturing of the elastic–viscoelastic samples of Fig. 5 in the near future.

### Experimental Methods.

We compress the periodic bimode sample in Fig. 3 *M*–*O* using a uniaxial testing device (Instron 5943 with a 500-N load cell), with laser-cut Plexiglas indenters, denoted by the white arrows in the figure, at a strain rate of *SI Appendix*, Fig. S1 is produced and tested in a similar way. However, the sample is compressed, instead, in the direction from the upper left to the bottom right.

The sample in Fig. 5 is compressed using the same uniaxial testing machine with straight surfaces on the top and bottom, at various strain rates. The sample is also positioned horizontally to enhance overall mechanical stability, with gravity also preventing out-of-plane buckling. For Fig. 3 (Fig. 5), the tests are recorded using a high-resolution

### Hinge Material Properties.

We measure the mechanical properties of a single compliant hinge of the sample in Fig. 5 using the same uniaxial testing device. We perform three different experiments: bending, tension, and shear. For bending, we pull on the hinge up to 0.30-rad deflection at rates of 0.002 and 20 mm/s. For tension and shear, we pull on the hinge up to 1.0-mm deflection at rates of 0.002 and 20 mm/s. We report the average of the measured stiffnesses in Table 1.

We observe that the stiffnesses at long timescales are similar for both types of hinges. However, at short timescales, the stiffnesses of the elastic (Elite Double 32) hinges increases only by approximately 20 to 30% compared to long timescales, while the stiffnesses of the viscoelastic (Agilus 30) hinges can increase by more than 1,300%. This difference in loading-rate dependency allows us to obtain the switchable response observed in Fig. 5.

### Simulating Oligomodal Metamaterials.

We use dynamic explicit nonlinear finite element methods (Abaqus 2018, Dassault Systèmes) to predict and design the dynamic response of oligomodal metamaterials, as well as to study the effects of changing material and geometry parameters. We model the geometry of the sample of Figs. 4 and 5 as seen in Fig. 9 using 2D triangular quadratic elements (CPE6). Plane strain elements are used as the rigid material to prevent the hinges from expanding and contracting in reality (as in Fig. 5). A mesh density featuring at least four elements through the hinge thickness and at least eight elements through the hinge thickness at the center is used, after performing a mesh convergence study showing that mesh refinement leads to little to no differences in global deformation.

We model the stiff plastic components as a linear elastic material with a Young’s modulus of *B*) have *B*) is accounted for using a three-term Prony series, with relaxation strengths

Reference points are created in the middle of the top and bottom stiff square segments (yellow in Fig. 9*A*), which are rigidly coupled to all nodes on these square elements, with the boundary conditions visualized in white. The bottom reference nodes are restricted from moving in a vertical direction, and one is also constrained horizontally. The top reference nodes are homogeneously compressed vertically using displacement control. All nodes are free to rotate.

To prevent generating jerky stress waves, the load application follows a smooth step in the time domain, where the first and second derivatives of the displacement with respect to time are zero, with average strain rates between

After running the analyses, we compute the polarization of the holes, as defined in Fig. 5, by tracking the nodes corresponding to the points highlighted in red in Fig. 9*A*.

## Data Availability

Experimental and numerical data (codes, raw images, and postprocessed images) have been deposited in Zenodo (https://doi.org/10.5281/zenodo.4699990) (42) .

## Acknowledgments

We thank D. Giesen and S. Koot for skillful technical assistance and S. Janbaz, R. van Mastrigt, D. Bonn, and M. van Hecke for insightful discussions. C.C. acknowledges funding from European Research Council Grant ERC-StG-Coulais-852587-Extr3Me.

## Footnotes

↵

^{1}A.B. and D.M.J.D. contributed equally to this work.- ↵
^{2}To whom correspondence may be addressed. Email: coulais{at}uva.nl.

Author contributions: A.B., D.M.J.D., and C.C. designed research; A.B., D.M.J.D., J.v.d.L., and C.C. performed research; A.B., D.M.J.D., and C.C. analyzed data; and A.B., D.M.J.D., and C.C. wrote the paper.

The authors declare no competing interest.

This article is a PNAS Direct Submission. K.B. is a guest editor invited by the Editorial Board.

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2018610118/-/DCSupplemental.

- Copyright © 2021 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- ↵
- K. Bertoldi,
- V. Vitelli,
- J. Christensen,
- M. van Hecke

- ↵
- K. Bertoldi,
- P. Reis,
- S. Wilshaw,
- T. Mullin

- ↵
- J. Shim,
- C. Perdigou,
- E. Chen,
- K. Bertoldi,
- P. Reis

- ↵
- ↵
- ↵
- ↵
- ↵
- G. P. T. Choi,
- L. H. Dudte,
- L. Mahadevan

- ↵
- Y. Cho et al.

- ↵
- D. Liarte,
- O. Stenull,
- T. Lubensky

- ↵
- P. Dieleman,
- N. Vasmel,
- S. Waitukaitis,
- M. van Hecke

- ↵
- A. Meeussen,
- E. Oguz,
- Y. Shokef,
- M. van Hecke

- ↵
- J. Grima,
- K. Evans

- ↵
- ↵
- S. D. Guest,
- J. W. Hutchinson

- ↵
- S. D. Papka,
- S. Kyriakides

- ↵
- ↵
- M. F. Ashby et al.

- ↵
- ↵
- W. Y. Jang,
- S. Kyriakides

- ↵
- S. Gaitanaros,
- S. Kyriakides

- ↵
- S. Gaitanaros,
- S. Kyriakides,
- A. M. Kraynik

- ↵
- ↵
- E. Lieb

- ↵
- R. Baxter

- ↵
- D. Dykstra,
- J. Busink,
- B. Ennis,
- C. Coulais

- ↵
- M. Stern,
- V. Jayaram,
- A. Murugan

- ↵
- S. Janbaz,
- K. Narooei,
- T. van Manen,
- A. A. Zadpoor

- ↵
- H. Zhang,
- X. Guo,
- J. Wu,
- D. Fang,
- Y. Zhang

- ↵
- S. M. Montgomery et al.

- ↵
- T. Frenzel,
- M. Kadic,
- M. Wegener

- ↵
- ↵
- M. A. Bessa,
- P. Glowacki,
- M. Houlder

- ↵
- E. Siefert,
- E. Reyssat,
- J. Bico,
- B. Roman

- ↵
- T. Frenzel,
- C. Findeisen,
- M. Kadic,
- P. Gumbsch,
- M. Wegener

- ↵
- L. Jin et al.

- ↵
- A. Zareei,
- B. Deng,
- K. Bertoldi

- ↵
- B. Pisanty,
- E. C. Oguz,
- C. Nisoli,
- Y. Shokef

- ↵
- C. Coulais,
- C. Kettenis,
- M. van Hecke

- ↵
- C. Coulais,
- A. Sabbadini,
- F. Vink,
- M. van Hecke

- ↵
- Dassault Systèmes

- ↵
- A. Bossart,
- D. M. J. Dykstra,
- J. van der Laan,
- C. Coulais

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Physical Sciences