## Abstract

Modern fabrication techniques, such as additive manufacturing, can be used to create materials with complex custom internal structures. These engineered materials exhibit a much broader range of bulk properties than their base materials and are typically referred to as metamaterials or microstructures. Although metamaterials with extraordinary properties have many applications, designing them is very difficult and is generally done by hand. We propose a computational approach to discover families of microstructures with extremal macroscale properties automatically. Using efficient simulation and sampling techniques, we compute the space of mechanical properties covered by physically realizable microstructures. Our system then clusters microstructures with common topologies into families. Parameterized templates are eventually extracted from families to generate new microstructure designs. We demonstrate these capabilities on the computational design of mechanical metamaterials and present five auxetic microstructure families with extremal elastic material properties. Our study opens the way for the completely automated discovery of extremal microstructures across multiple domains of physics, including applications reliant on thermal, electrical, and magnetic properties.

## INTRODUCTION

Microstructures can exhibit remarkable physical properties that extend beyond the properties of their constituent materials. Many microstructure types have been developed to demonstrate applications in mechanics (*1*–*6*), acoustics (*7*, *8*), and electromagnetics (*9*–*11*). These microstructures are typically designed by domain experts using time- and labor-intensive manual processes. These designs are often programmable in the sense that they have a small number of control parameters for generating a family of geometries. A given microstructure family can be tested by performing simulations or experimental measurements on samples with different sets of parameter values. This testing process allows the mapping of the control parameters to physical properties, which in turn helps to reveal the underlying design principles that explain the behavior of the structures. In practical applications, this mapping also enables the selection of a family member that has a desired trade-off between physical properties, allowing it to achieve some target performance objective (*12*). Unfortunately, it is rare for manually designed microstructure families to reach extremal material properties, for the reason that the space of possible microstructure designs is combinatorial and therefore impossible to explore exhaustively. One standard approach to bypass this design challenge is to use computational methods, such as topology optimization (*13*–*15*), with a computer simulation in their inner loop to find a microstructure with a specified trade-off of physical properties. However, although topology optimization is a powerful tool for finding individual structures with target material parameters, constructing parametric models from these optimized structures requires further expertise and manual design effort (*16*). In contrast to previous work, we present the first computational method that automatically explores the space of microstructure designs and generates parametric microstructure representations, and we apply this technique to the discovery of parametric families optimized for competing properties.

Our methodology is designed to be modular, allowing users to choose suitable components for different physical properties. Here, we apply our method to the design of mechanical metamaterials (*17*). Specifically, we ask our algorithm to construct templates for a particularly interesting type of mechanical microstructure—auxetic materials—which have a negative Poisson’s ratio. These materials have the unusual property of becoming laterally thinner under axial compression. Two-dimensional (2D) auxetic structures are well understood because of their relatively simple shape and geometric features, with common classes including reentrant structures (*13*, *18*), chiral structures (*19*, *20*), and rotating mechanisms (*21*–*23*). Unfortunately, generalizing existing 2D structures to 3D is challenging because a simple arrangement of 2D mechanisms often results in orthotropic or other anisotropic structures with low shear resistance. These structures will easily shear when the load is not well aligned with the auxetic direction. In addition, because Poisson’s ratio for orthotropic structures is unbounded, orthotropic auxetic structures are much easier to find than isotropic ones (*24*). In a pioneering work, Lakes (*18*) fabricated and tested the first isotropic 3D auxetic structure. Lakes’ work proved the feasibility of 3D auxetic structures but did not aim to provide any means to control elastic material parameters such as Young’s modulus, Poisson’s ratio, and shear modulus. Designing manufacturable 3D auxetic structures with desired elastic parameters remains a complex task, and only a handful of 3D designs have been fabricated and measured (*24*, *25*). Our results show that it is possible to automatically discover common design patterns of auxetic structures that can be parameterized by a small set of control parameters. We demonstrate this capability by presenting five discovered auxetic families that cover a wider range of material parameters, which have been shown in previous works.

At a high level, our system takes a set of base materials as input and outputs parametric microstructure templates. By turning voxel-based microstructure representations into parametric shapes, structures of arbitrary resolution can be obtained. This allows the generation of microstructures with a resolution adapted to the chosen fabrication method. These newly generated structures achieve optimal trade-offs of physical properties. Instead of finding individual optimal structures, we aim to understand the common mechanisms that make the structures work, which allows us to come up with new structures using these very mechanisms. We propose a four-step software pipeline to find these patterns in a data-driven approach (Fig. 1). First, we generate a database of microstructures covering the range of achievable material parameters. We call the range of realizable material parameters the material property gamut. With the gamut at hand, we then select a subset of microstructures near the gamut boundary that share similar material properties. Our algorithm analyzes the geometric similarity among these structures. Because many distinct structures can achieve similar properties, we separate them into families based on their geometric features. Once the families are separated, our algorithm extracts their topological structures to form templates. Thus, these templates reveal the underlying common structure of the microstructures within each family. They can also be used to generate new structures by varying a small number of their parameters. To create new structures in the vicinity of the material gamut boundary (that is, with nearly extremal properties), we compute directions along which to vary the template parameters using a regression method. The final output of our system is a set of templates controlled by reduced parameters for generating new resolution-independent structures along the gamut boundary.

The details of the four steps of our discovery pipeline are presented and discussed in what follows. The first step of the pipeline estimates the material property gamut given a set of base materials. We used two base materials in our experiments: a stiff material and a soft or void material. This step computes a set of microstructure samples that covers as much of the material property gamut as possible within our computational budget. We define a microstructure as a set of material assignments in a regular grid. This definition is more general than alternative definitions such as networks of trusses or geometries generated from a set of analytic functions. However, our definition also brings new challenges, because a broader design space means many more designs to explore. For example, if we use two base materials and a 64 × 64 × 64 grid, then the number of possible material combinations is 2 to the power of 64^{3}. Ideally, this step would try all possible combinations, remove the ones unsuitable for manufacturing, and compute the gamut achieved by the remaining structures. In practice, with a limited computational budget, we need an efficient sampling algorithm that allocates more computing resource to the promising structures near the gamut boundary. To this end, we rely on a sampling algorithm that alternates between continuous optimization and discrete stochastic search to progressively expand the gamut (*26*). Both steps randomly pick structures near the gamut boundary as promising initial guesses. The continuous optimization stage pushes structures past the explored gamut boundary along gradient directions, whereas the discrete stochastic stage escapes local minima by making jumps without following the gradient.

Our study on mechanical microstructures uses a topology optimization formulation inspired by the work of Andreassen *et al.* (*24*) for the continuous optimization stage. The disadvantage of the continuous optimization step is that the objective function for the microstructure parameters is highly nonconvex and has many local optima (*27*). The discrete stochastic stage overcomes the limitation by introducing discrete structural changes, which allows us to explore different local optima. Note that the discrete sampling is not as efficient as the continuous optimization algorithm in approaching local optima. A discrete change often makes a structure suboptimal. These trade-offs between discrete and continuous sampling algorithms lead us to believe that the interleaving of the two algorithms is necessary for general nonconvex material optimization problems. For different objectives and physical phenomena, the details of each step should be modified to incorporate domain knowledge.

In the second step, we identify common geometric traits among microstructures near the gamut boundary. As mentioned above, the microstructure design problem is nonconvex, with many different microstructures achieving similar material properties. To find meaningful similarities among structures, we first need to separate structures with distinct geometric features into different families. One naive approach would be to treat the material assignment of each microstructure as a feature vector and use an off-the-shelf clustering algorithm. However, as our preliminary study showed us, the resulting clusters are unsatisfactory because they include structures with different topological features while excluding structures that should be included. This is expected because structures with the same topology may seem very distant when compared to each other using a per-voxel material difference. Therefore, we use a nonlinear dimensionality reduction (NLDR) to plot the microstructures in a low-dimensional embedding space where topologically similar structures are arranged more closely. We chose Isomap (*28*) as the reduction method because it can discover long sequences of related structures while keeping distant points well separated. The effectiveness of NLDR depends on the distance metric that measures the geometric difference. A smoothed Euclidean norm is chosen for robustness (fig. S1). NLDR outputs an embedding of the microstructures in a low-dimensional space where similar structures are tightly packed. Microstructures in the embedding space are clustered using a Gaussian mixture model (*29*) where each cluster corresponds to a family. Families with a significant number of members (>200) are extracted for further analysis.

The third step constructs templates for each microstructure family. A template represents a class of structures that can be generated by varying continuous template parameters instead of editing per-voxel material assignments. Because templates are analytical expressions, they can be voxelized at different grid resolutions. A template for a family should be able to reproduce all the structures in the family it belongs to by varying its template parameters. In addition, a compact representation can reveal common design patterns within a family. Theoretically, one could use volumes bounded by polynomial surfaces as a very general type of building block for the templates. However, for this problem, such a general formulation adds unnecessary complexity to finding simple template expressions. Because most of the extremal structures appear to be composed of beams, plates, and blocks, we chose cuboids with different side lengths as building blocks for the microstructure templates. We note that some structures near the boundary do not fit the cuboid representation (fig. S2), but this is not a fundamental limitation because our modular software pipeline can be easily extended to use other types of geometric building blocks such as generalized cylinders and curved cuboids. To find a template from a family representative, we compute its topology using a morphological skeleton (*30*). The morphological skeleton is a set of connected voxels that largely preserves topological and branching characteristics of the structures (fig. S3). The skeleton is converted into a graph to represent a template. A cuboid is placed on each edge of the graph with optimized sizing and orientation to best match the representative structure. More details of the process are available in the Supplementary Text.

Finally, our algorithm computes reduced parameters to allow an intuitive navigation in the material property space. The output from the previous step is a template parameterized by tens of parameters. These parameters control geometric features, such as beam thicknesses and orientations, but they do not directly control the material properties of the structures. Therefore, it is still difficult to understand the fundamental design principles that affect the material properties. This step computes reduced parameters that allow direct edits of the material properties by varying all the relevant geometric features simultaneously. This step proceeds by computing template parameters for each structure in each family. To avoid outliers, we exclude microstructures with too large fitting errors (>5% voxel difference). Principal component regression (PCR) is then performed on the set of fitted template parameters to find principal directions in the template parameter space. Varying the template parameters in a principal direction corresponds to moving on the gamut boundary in a specific direction. A reduced parameter is assigned to each direction to control the amount of change in that direction. Using this representation, engineers or automatic algorithms can select extremal structures by only varying the reduced parameters. The final output of our system is a set of templates tunable by reduced parameters. The templates produced in our experiments are included in the Supplementary Materials in the form of functions that take in reduced parameters and output microstructures at desired resolutions.

## RESULTS AND DISCUSSION

The results of this study focus on elastic material properties: Young’s modulus, Poisson’s ratio, and shear modulus. The elastic material property gamut is estimated from 15,000 3D cubic-symmetric microstructures at a voxel resolution of 64^{3}. The voxel resolution is a power of 2 because that is necessary to achieve optimal performance of our multigrid finite element method (FEM) simulation. The specific resolution 64^{3} is chosen because it is sufficient for discovering auxetic structures with a wide range of relative shear modulus, whereas 32^{3} structures cannot achieve comparable complexity or property ranges. The macroscopic elastic parameters of each microstructure are computed using homogenization theory (*31*, *32*) assuming periodic boundary conditions (that is, the structure is repeated infinitely). Each microstructure consists of a per-voxel binary material assignment. Because of manufacturing limits on the minimum feature sizes, sensitivity filtering (*27*) is applied in the gamut sampling step to prevent the structures from having overly thin features. Sensitivity filtering works by applying a low-pass smoothing filter to the gradient of the 3D binary image of the microstructure to force its material distribution to have smooth variations. This results in a structure with coarser geometric features such as thicker beams instead of many thin beams or checkerboard patterns.

The computational bottleneck of the gamut sampling step is the simulation of microstructures. In the current implementation, it takes an average of 10 s to simulate the deformation of a multimaterial microstructure on a single graphics processing unit (GPU). Homogenization takes twice the amount of time because it requires two simulations. For other types of physical properties, such as thermal or electromagnetic properties, simulation is as costly as elastic properties. Because of the computational cost of simulation, we will still need efficient sampling algorithms to explore material property gamut of other physical properties.

Here, we report a deeper analysis on the auxetic structures present in the gamut, which are significantly more complex than families with positive Poisson’s ratios (fig. S4). Five families with significant numbers of members (Fig. 2B) were discovered using three Isomap embedding dimensions. We confirmed that Isomap associates seemingly distant structures through intermediate structures. For example, structures 5-1 and 5-3 from family 5 have very different beam thicknesses resulting in large geometric distances. However, the embedding reveals that there is a sequence of structures like 5-2 that make the connection between them.

Our system constructed parametric templates of all the five families. The initial topology of the templates is extracted from morphological skeletons (Fig. 3B). Although the topologies are visually complex, they are generated by mirroring a small number of beams (highlighted in red) reflected according to cubic symmetry (table S1). The most complex template 5 contains only six control beams. The five families cover similar ranges of Young’s moduli and Poisson’s ratios while spanning different ranges of shear moduli. Inspired by classical linear elasticity theory, we analyzed the shear modulus ratios of the microstructures defined aswhere *G* is the shear modulus, *E* is the Young’s modulus, and ν is the Poisson’s ratio of the metamaterial. For traditional isotropic materials, the theoretical ratio is one. A low ratio indicates low resistance to shear deformation. For auxetic materials, lower ratios are much easier to obtain than higher ones. Even with foam structures assumed to be isotropic, experimental data from previous work indicate that the ratio is less than one (*33*).

Template 1 shown in fig. S5A resembles the conceptual sketch by Lakes (*18*) and belongs to the reentrant class of geometry. The difference with our template is that the latter has only two beams mirrored by cubic symmetry, whereas Lakes’ sketch contains three beams (fig. S5B). In addition, using Lakes’ conceptual template as an initial configuration, our topology optimization algorithm could not optimize it toward a structure with a shear modulus ratio higher than 0.12. We note that Lakes’ template is intended as a conceptual design for explaining auxetic structures. His early work does not aim to provide full control over other elasticity parameters. Template 1 is the most straightforward auxetic template that we identified, because our microstructure database does not include any single-beam auxetic structure. The shear modulus ratio of this family falls between 0.07 and 0.24, which is the lowest range among all five families. Templates 2 and 3 are similar to each other and differ by a diagonal beam located in the face center (highlighted in green in template 3). Because their geometric difference is small, they are adjacent in the Isomap embedding space. The central beam is responsible for increasing the shear modulus of the structures. For structures with ν around −0.5, the additional beam increases the maximum shear modulus ratio from 0.34 to 0.90. Templates 4 and 5 also differ by a single beam. Even the most complex template 5 is optimized from a simple cube frame through our continuous optimization step, as shown in movie S1. The additional beam in template 5 makes the family stiffer overall. Both families can achieve shear modulus ratios greater than 1 for ν < −0.5.

For each family, principal directions of template parameters are extracted using PCR. The template definitions are included in the Supplementary Materials. Each template is represented as a function that maps from two reduced parameters to a microstructure. We kept two principal directions to tune the Young’s modulus and the shear modulus. We made some observations by varying the parameters along these directions. For families 2 and 3, the thickness of the slanted column (Fig. 4A, highlighted in red) is crucial for controlling the Poisson’s ratio. More specifically, the Poisson’s ratio increases quickly when increasing the beam thickness. For families 4 and 5, the thickness of the rotating triangle affects the trade-off between the Young’s modulus and the shear modulus (fig. S6). The effect of varying the reduced parameters for family 4 is shown in movie S2.

Although our cuboid-based templates have very few parameters, they are sufficient for replicating the auxetic behavior of the corresponding families. We validated the auxetic properties of the fitted microstructures using numerical simulation. The templates provide an efficient way for sampling new structures by varying the reduced parameters. We sampled 300 new structures from each family along the two main PCR coordinate directions. The coverage of the templates in the microstructure gamut (Fig. 2B) shows that they can effectively generate microstructures on the gamut boundary.

So far, all our simulations were carried out assuming linear elasticity, which is only accurate for infinitesimal deformations. We also made the common assumption that there is no self-collision, which imposes a limit on the maximum compressive strain applied to the structures before self-collision occurs. Representative structures from families 4 and 5 have the lowest limit at 7% compressive strain. In practice, nonlinear deformations, such as bending and rotation, are prevalent in our auxetic structures. These deformations can cause linear elasticity to incorrectly predict significant volume expansion of rotating parts (up to 20% expansion in our test cases). Thus, we tested our structures using a nonlinear deformation model to understand their behavior under large deformations. We simulated the nonlinear deformation behavior using a neo-Hookean material model (*34*). Up to a maximum allowed strain of 7%, the linear elasticity and the neo-Hookean model are still in acceptable agreement with an average error of 16% in the computed Poisson’s ratios. In addition to performing simulations, we also manufactured three example structures from each family with different Young’s moduli and material ratios. All the structures were printed using a single elastic material. Our structures demonstrated consistent auxetic behaviors (movie S3), although they were optimized using the linear elasticity assumption. It is noteworthy that our structures do not rely on structural instabilities (*35*) to gain their auxetic behavior and that they shrink uniformly as the applied load increases. This means that their deformations consistently follow the same pattern for different trials.

Our process automatically discovered two types of auxetic mechanisms: slanted columns and rotating triangles (Fig. 4). The slanted column transforms vertical compression to horizontal motions. The rotating triangles transform vertical compression into a winding deformation that pulls the right end of the mechanism toward the center of the microstructure. Their motions are shown in movie S3. The rotating triangle mechanism bears resemblance to existing 2D structures (*36*) known as chiral structures (fig. S6D). 3D chiral structures have been shown to exhibit extraordinary effects such as negative effective compressibility and twisting under compression (*37*, *38*). Extension of chiral structures to 3D cubic structure with large shear modulus has never been reported before. Unlike traditional design approaches, our algorithm discovers these mechanisms entirely automatically without imposing any artificial geometric restrictions—all microstructures are built from hexahedral voxels. To inspire future applications of these mechanisms, we report their loading behavior. These auxetic mechanisms are the most active parts in the microstructures. They act like joints that connect the more rigid scaffolding in microstructures. Because of this, they undergo the most deformation and concentrate a large amount of stress. For the rotating triangles, the stress is concentrated on the connections around the triangle. We computed the maximum principal strain in the structure with respect to the vertical compressive loading to provide insights into the strength of the block. At the maximal compressive loading (7%), the maximum principal strain in the structure is 7%. Calculation using a reported Young’s modulus of 80 MPa yields a von Mises stress of 6.72 MPa (Fig. 4E), whereas our print material has a reported strength of 8.5 MPa. The printed structures are approaching the strength limit under the load. Because the available material is relatively weak even compared to conventional materials such as acrylonitrile butadiene styrene (ABS) plastics and rubber, we believe that structural strength can be improved significantly with future manufacturing materials.

## CONCLUSION

We have shown a computational method that combines discrete sampling, continuous optimization, and dimensionality reduction methods for automatic discovery of new microstructure families and mechanisms that would have been challenging to design manually. The discovered structures are suitable for manufacturing because they avoid thin features and distribute deformation over beams instead. They also span a wide range of shear moduli, allowing engineers to balance between different macroscopic properties. Although our case study focuses on elastic material properties, the technique may be applied to other physical properties whenever predictive simulation exists. Our computational pipeline paves the way to the discovery of structures that balance mechanical, thermal, optical, acoustic, and electromagnetic properties. Moreover, it advances the understanding of underlying mechanisms that are key to extremal properties.

## MATERIALS AND METHODS

### Simulation of elasticity

We used a FEM with linear hexahedral elements to simulate the elastic behavior of the microstructures. The macroscopic behavior of a structure was computed using homogenization theory (*22*) assuming an infinite tiling of identical structures (periodical boundary condition). The material parameters of the microstructures were computed in the context of linear elasticity. Despite only being valid for small deformations, the linear elasticity assumption allows more efficient computation of the parameters. It also prevents issues related to structural instability such as buckling. For validation, we used a neo-Hookean material model to simulate and compare our microstructures to real experiments. For both linear and nonlinear simulations, most of the computational power was spent solving a linear system of the form **Ku** = **f**, where **K** is the stiffness matrix modified to account for the periodic boundary conditions. The linear system was solved using a GPU implementation of a geometric multigrid preconditioned conjugate gradient algorithm inspired by previous work (*39*). At a very high level, geometric multigrid algorithms work by merging neighboring elements into coarse voxels to quickly estimate solutions. They work best when the grid resolution is a power of 2. We chose 64^{3} as our grid resolution because it is the minimum resolution required to arrive at the auxetic structures reported here. Our simulation took an average of 12 s on a Titan Z GPU. Two simulations, including stretching and shearing, were needed to compute the macroscopic elastic parameters of a cubic symmetric microstructure (Young’s modulus *E*, Poisson’s ratio ν, and shear modulus *G*). We use a soft material to simulate empty space. The Young’s modulus of the soft material is chosen to be 5 × 10^{−5} relative to the Young’s modulus of the solid material.

### Continuous optimization of microstructure

Continuous optimization is one of the two methods used for sampling the microstructure gamut. It assigns a continuous “density” variable to each cell that represents the ratio of stiffness material. For extension to multiple materials, we can assign multiple variables representing the mixture ratio of different materials. We used the SIMP (solid isotropic material with penalization) scheme to interpolate the linear elasticity material properties (*24*). For a cell with density *x*, its stiffness matrix **K**(*x*) is computed by interpolating the stiffness matrices of the soft and the stiff material using a constant exponent *c*.

We set *c* = 3 in our experiments. The linear tensor of elasticity **C**(**x**) of a microstructure with a continuous material assignment **x** can be computed using homogenization theory (*33*). Here, we computed six vertex displacements **u**_{i} under periodic boundary conditions (stretch and shear along three axes). For cubic-symmetric structures, only two displacement fields, a stretch and a shear, need to be simulated. The rest can be obtained by permuting spatial coordinates. The elasticity tensor **C**(**x**) is given by

The material properties **p**(**x**) include Young’s modulus, Poisson’s ratio, shear modulus, and total density. The first three properties are functions of the elasticity of **C**, and the last property is the sum of **x** across all cells in a microstructure grid. Because we are interested in softer structures, Young’s modulus and shear modulus are represented in the log space.

Given a discrete microstructure and its material properties **p**, the level set gamut tells us the direction toward the exterior of the known gamut (*26*). We move along the direction for a fixed distance to arrive at a target material property **p̂**. Our objective is to move the current microstructure toward the target by changing its continuous material distribution. We use a quadratic objective function

Note that the density variable on each cell has a bound constraint between 0 and 1. This bound-constrained nonlinear least squares problem can be solved using gradient-based optimization algorithms. We derived the analytical gradient (sensitivity) of the objective function with respect to **x**. Each term in the gradient takes the form

The elasticity parameters depend on displacement vectors **u**_{i} computed from the homogenization step. For a material property *p*(**x**) related to the linear elasticity tensor **C**(**x**), it is further expanded as

We use the method of moving asymptotes (MMA) (*40*) to optimize the objective to reach local minima. MMA usually converges within 30 steps. We run MMA until either convergences or reaching a maximum of 50 steps. Because our experiments used cubic-symmetric structures, the continuous optimization effectively optimized material distribution inside 5984 cells instead of 64^{3} cells. Note that the objective function applies to general domains without any symmetry assumptions. Removing the symmetry constraint would allow us to sample more classes of structures at the cost of more computation time. At the end of the continuous optimization, the structure was converted back to a discrete structure by thresholding. Because our goal was to explore more structures, we thresholded at several density values instead of a single value to obtain multiple structures from an optimization. In the experiments, the thresholds were 0.1 to 0.5 with an interval of 0.1.

### Fabrication and measurement

Three structures with different material properties from each family were printed to verify our simulation (table S2). All structures were printed using an EOS SLS printer with a PEBA2301 elastic material. The printer required a minimum feature size of 0.9 mm for wire diameters and 0.8 mm for wall thicknesses. To satisfy the printing constraints, we scaled each cell to a side length of 2.54 mm. Simpler structures from families 1 to 3 were printed using a 2 × 2 × 2 grid arrangement, whereas more complex families are printed using a 3 × 3 × 1 arrangement to allow support material to escape. The base printing material was measured using an Instron 5944 with tensile tests instead of compression tests, because a solid block of base material was too stiff for our equipment. The Poisson’s ratio was measured using a video camera fixed on the test machine (fig. S7). Although all samples were printed using the same printer and the same materials, the Young’s modulus of the prints was highly variable. More specifically, the stiffest sample had a Young’s modulus that was twice as high as the one of the softest samples. On the other hand, the Poisson’s ratio of the base print material was measured to be 0.344 and had a much lower variance of 0.02 in the sample set. The 3D structures were measured using a compression test at a speed of 2 mm/min with a maximum strain of 6 mm. The compression plates were lubricated with oil to reduce friction. Significant variance in Poisson’s ratio was observed due to several factors in manufacturing and measurement. An additional challenge was that the printer did not reliably reproduce the geometry specified by input files. In practice, the printed models were thickened by 0.1 to 0.4 mm, which is significant compared to the thinnest feature size in our microstructures. The thickening stiffened the joints and reduced Poisson’s ratios. The effect was exacerbated by incomplete support removal. The support material was the same as the print material in powder form, which stuck to the print easily especially around hard-to-reach internal corners. We believe that the discrepancy can be reduced in the future by using more precise printing technologies with soluble support material.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/1/eaao7005/DC1

Supplementary Text

fig. S1. Smoothing for geometric difference metric.

fig. S2. Examples of auxetic structures in addition to the five selected families in the main paper.

fig. S3. Computing a microstructure template from a representative structure.

fig. S4. Structures with large Poisson’s ratios (ν > 0.3).

fig. S5. Microstructures that resemble designs from previous works.

fig. S6. Reduced parameters for family 4.

fig. S7. Test apparatus for measuring Young’s modulus and Poisson’s ratio.

table S1. Auxetic families and templates.

table S2. Simulated and measured Poisson’s ratios of example structures.

movie S1. Continuous search of a microstructure using topology optimization.

movie S2. Shape variation of structures from family 4 along principal directions.

movie S3. Compression testing and simulation of example from families 1, 3, and 5.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank C. Dick for providing the code of his GPU multigrid FEM solver that inspired us. We also thank D. Rus for providing access to the mechanical test machines that we used.

**Funding:**The research is funded by Defense Advanced Research Projects Agency (DARPA) Simplifying Complexity in Scientific Discovery (SIMPLEX) N66001-15-C-4030.

**Author contributions:**D.C. contributed to the software development and the physical experiments. M.S. contributed to the software development and the original ideas. B.Z. contributed to the software development and the original ideas. W.M. contributed to the original ideas and the method development. All authors contributed to editing the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).