Scale-rich network-based metamaterials

We begin by determining the mechanical characteristics of SR microstructures, comparing them against three single-scale microstructures frequently used in engineering applications: the square lattice, the hexagonal lattice and the irregular Voronoi tessellation. By generating finite domains of each microstructure with a sufficiently large number of characteristic features, we confirm that the P(λ), P(ℓ) and P(k) for square and hexagonal lattices result in delta functions, whereas the Voronoi tessellation is characterized by narrow distributions (Fig. 1i–k).

To extend the comparison beyond theoretical differences between architectures, we experimentally and numerically determined the mechanical properties of these metamaterial microstructures, using vat photopolymerization 3D-printing (Methods) to fabricate centimetre-scale specimens out of elastic–plastic polymer resin (Fig. 1d, left), with a feature resolution of approximately 100 μm. We focused on SR systems with T = 500—a point where the homogenized mechanical properties show convergence and boundary effects are minimized (Supplementary Information, section 2.3). Given this feature resolution and the 500-ligament constraint, the relative densities achievable with vat photopolymerization 3D-printing ranged between \(20 \% \le \bar{\rho }\le 75 \% \) (Supplementary Information, section 2.5); we can also demonstrate lower fill fractions within the tunable-density phase (\(\bar{\rho } < 20 \% \); Fig. 1c) through the use of two-photon lithography (Fig. 1d, right) with submicrometre minimum feature sizes. Next, we characterized the linear and nonlinear mechanical responses in the quasi-static and dynamic regimes of these SR microstructures.

Highly tunable linear-elastic responses

To determine the linear-elastic response of SR metamaterials, we performed computational homogenization using the finite-element method. Imposing affine boundary conditions, we extracted a directional modulus E_d in each direction e_d and calculated a non-dimensional modulus \(\bar{E}\equiv {E}_{{\rm{d}}}/{E}_{{\rm{s}}}\), where E_s is Young’s modulus of the constituent material (Methods).

At a relative density of \(\bar{\rho }=40 \% \), we performed this analysis on the three single-scale systems (square, hexagonal and Voronoi) and SR geometries, demonstrated by the representative realizations in Fig. 2a, via polar representations of \(\overline{E}\) as a function of the probing direction. While the hexagonal and Voronoi structures show a near-constant \(\overline{E}\) in all orientations (that is, an isotropic stiffness response, as expected⁴⁶; Fig. 2a(i),(ii)), the square lattice is highly anisotropic, displaying high stiffness parallel to its struts and minimal stiffness in any other direction (Fig. 2a(iii)). In contrast to these near extremes, we find the behaviour of the SR structures to be modulated by α: for α ≈ 0, the stiffness is nearly isotropic (Fig. 2a(iv)); and as α → 1, a principal orientation of maximum stiffness emerges (Fig. 2a(v),(vi)). Altogether, the full set of SR architectures exhibited highly tunable (an)isotropy exceeding that of the three conventional geometries combined.

a, Directional stiffness plots of the reference hexagonal (i), Voronoi (ii) and square (iii) lattices from finite-element simulation, compared with SR samples with α = 0 (iv), α = 0.3 (v) and α = 0.9 (vi) at the same relative density. The right panels show the experimentally measured directional stiffness for SR materials with α = 0 (vii),  α = 0.1 (viii) and α = 0.3 (ix), validating the theoretically predicted stiffness anisotropy. Scale bars, 10 mm. b, Anisotropy–density plot, illustrating the broad range of anisotropy achievable with the SR framework compared with the narrow range of values attained by the reference samples. The black filled circles (vii)–(ix) show the experimentally measured anisotropy. c, As the decay parameter α increases, anisotropy arises because the maximum specific directional stiffness increases, whereas the minimum specific directional stiffness is constant. The solid lines represent the mean across all samples, and the shaded regions correspond to a ±1 standard deviation from the mean value. d, Phase diagram in (α, λ₀) space of the anisotropy, which functions as a design map for a tunable anisotropic response across relative densities \(0.2\le \bar{\rho }\le 0.75\). e, Temporal evolution of ligament strain energy during the SR generation process. f, Snapshots at fixed time t are shown for two example ligaments in an SR network, corresponding to the markers in e. g, Universal scaling of the ligament energy U_i as a function of the degree, k_i, thickness, λ_i and angle, θ_i. For each α, data are averaged over n = 5 independent realizations, each with T = 5,000 steps.

To quantify the mechanical implications of this observed range in anisotropy across the full space of SR architectures, we plot the elastic anisotropy ratio \({\eta }_{{\rm{E}}}\equiv {E}_{\min }/{E}_{\max }\in (0,1]\), where \({E}_{\text{min}}\equiv \text{min}(\bar{E})\) and \({E}_{\text{max}}\equiv \text{max}(\bar{E})\), as a function of relative density for the reference geometries and 2,600 distinct SR microstructures with different α and λ₀ (Fig. 2b). We confirm that the hexagonal lattice (0.87 ≤ η_E ≤ 0.94) and the Voronoi tessellation (0.88 ≤ η_E ≤ 0.97) remain nearly isotropic across a wide range of relative densities. For the square lattice, the anisotropy decreases from η_E = 0.03 at \(\bar{\rho }=20 \% \) to η_E = 0.55 at \(\bar{\rho }=70 \% \). By contrast, we find that the SR materials cover 77% of the available area on the density–anisotropy plot (Fig. 2b), reaching values not accessible to the reference systems. Indeed, by changing α and λ₀, we can attain highly anisotropic systems (η_E ≤ 0.2) across a wide range of densities (\(17 \% \le \bar{\rho }\le 64 \% \)) as well as nearly isotropic systems (η_E ≥ 0.85) in the range \(33 \% \le \bar{\rho }\le 72 \% \). Simultaneously, we can generate systems of widely varying anisotropies at any fixed relative density: at \(\bar{\rho }=20 \% \), we span 0.04 ≤ η_E ≤ 0.78; at \(\bar{\rho }=40 \% \), we span 0.12 ≤ η_E ≤ 0.88; and at \(\bar{\rho }=70 \% \), we span 0.26 ≤ η_E ≤ 0.93. Directed growth approaches (Supplementary Information, section 2.15) further extend the achievable coverage to more than 96% of the density–anisotropy space.

To elucidate the origin of this tunable anisotropy, we plot the mass-normalized (that is, specific) minimum and maximum stiffness values (\({E}_{\min }/\bar{\rho },{E}_{\max }/\bar{\rho }\)) as a function of the decay parameter α across the entire sample set (Fig. 2c). Although the minimum specific stiffness is largely independent of α, increasing α while keeping density constant enables an increase in the maximum achievable specific stiffness, leading to robust tunability in the architecture space. This tunable anisotropy is best captured by the (α, λ₀) configuration space (Fig. 2d), which shows that anisotropy increases with α and decreases with λ₀. To validate our finite-element model predictions and confirm the isotropic and anisotropic responses, we performed uniaxial compression experiments across various SR geometries and orientations, confirming that at constant density increasing α directly leads to greater anisotropy (Fig. 2a(vii)–(ix), Fig. 2b, filled black circles, and Supplementary Information, section 2.7).

A scalar metric that accounts for the alignment of material along a certain loading direction indicates that the tunable directional stiffness of SR materials is the result of the inherent heterogeneity in ligament orientations and sizes. This is best captured by an orientation tensor M

where A_t is the area of ligament t and d_t is a unit vector along its centre line. M allows us to compute a geometric anisotropy η_G ≡ μ₁/μ₂ ∈ (0, 1], where μ₁ ≤ μ₂ are the eigenvalues of M. We find that η_G agrees with the elastic anisotropy η_E (Supplementary Information, section 2.6), confirming that the anisotropy of the SR system is a direct consequence of its length-scale heterogeneity.

Tunable energy-density scaling

To understand the role of multiscale features on the measured properties, we follow two ligaments, one introduced at step t = 20 and one at t = 100, through the full generation process. As Fig. 2e shows (see Fig. 2f for the network snapshots), these ligament-specific trajectories span more than eight orders of magnitude and are irregular and rather unpredictable. Yet, as we show in Supplementary Information, section 3, the strain energy distribution across all ligaments scales with the introduction time as i^{−(3/4+α/2)}, with the degree as \({k}_{i}^{3/2+\alpha }\), with the length as \({{\ell }}_{i}^{3/2+\alpha }\), and with the thickness as \({\lambda }_{i}^{3/(4\alpha )+1/2}\) (Supplementary Fig. 46). Combining these scaling laws, we find that the elastic strain energy follows a single per-ligament law (Fig. 2g and Supplementary Fig. 47)

where f(θ_i) is an orientation factor, with θ_i being the ligament orientation. This law is, remarkably, independent of α, confirming that all scales contribute to the total energy. Plotting the per-ligament energy against the combination in equation (2) collapses all ligaments onto a single line of slope approximately 1, independent of α (Fig. 2g), confirming that all ligaments follow the scaling law equation (2).

Integrating the per-ligament law, equation (2), over the SR distributions P(k) ≈ k⁻³ and P(λ) ≈ λ^−(1+1/α) allows us to analytically derive the scaling relation between the total elastic energy and relative density, \({U}_{\mathrm{tot}}\approx {\bar{\rho }}^{m(\alpha )}\) (Supplementary Fig. 49c), finding that the scaling exponent m depends on the SR exponent α as

where \(p(\alpha )\) is a fitted exponent describing the residual time dependence of the per-ligament energy. This expression matches the directly fitted slopes within ±0.2 (more details in Supplementary Information, section 3). The SR construction therefore interpolates continuously across the scaling map, with α as the only control parameter. These results show that the multiscale mechanics of SR networks can be reduced to a universal ligament-level energy law, which turns the growth exponent α into a macroscopic control parameter that continuously tunes the energy-density scaling response, with lower m exponents indicating more efficient mass utilization.

Delocalized nonlinear deformations

Next, we characterized the nonlinear response of SR metamaterials by conducting large-deformation uniaxial compression experiments (Methods) on the 3 reference geometries and SR systems of α = 0, chosen to minimize the impact of elastic anisotropy and emphasize the effects of the heterogeneity of ligament lengths (Fig. 3a; 3 specimens each, all at nominal \(\bar{\rho }=42 \% \); experimental results for SR systems with α > 0 are discussed in Supplementary Information, section 2.11.). Experiments on the square, hexagonal and Voronoi tessellations showed qualitative differences between each architecture (for example, in the densification strain); repeated samples of each architecture showed minimal variation (Supplementary Information, section 2.8). By contrast, all realizations of SR architectures (Fig. 3a, orange dotted lines) showed a monotonically increasing stress–strain response, coexisting with a wide envelope of responses, elucidating the tunability enabled by the SR framework and breaking the structure–response constraints observed in the single-scale reference geometries.

a, Experimental stress–strain curves of SR and reference samples, for \({\bar{\rho }}_{\mathrm{nom}}=42 \% \), under uniaxial compression. For the reference samples, the mean behaviour of three replicate specimens is plotted. For the SR geometries, the solid line represents the mean of all samples with α = 0, whereas the dashed curves are individual samples. The response of SR samples with α > 0 is given in Supplementary Fig. 33. b, Experimental toughness–density Ashby plot, showing the energy absorption capacity across all relative densities. c, Evolution of localized deformation with increasing macroscopic strain: still images together with computed areal Jacobians, J = A/A₀, for each cell. The geometries dominated by a characteristic length scale (square lattice, (i); hexagonal lattice, (ii); Voronoi tessellation, (iii)) all develop regions of compaction where localization is observed; only the SR sample (iv) avoids the formation of compaction bands. We tracked images up to a macroscopic strain of 25%, after which compaction precludes the accurate computation of Jacobian values. Scale bars, 10 mm. d, Evolution of the localization parameter \({\mathcal{L}}\) in each geometry together with snapshots of the deformation evolution as recovered from still images. Here, the mean and standard deviation of three unique SR specimens are plotted.

From each stress–strain response, we computed the toughness, measured as the integral under the stress–strain curve to an engineering strain of 40%—before the onset of complete densification. These toughness metrics, when plotted against relative density (Fig. 3b), showcase the ability of SR architectures to match or surpass the mechanical resilience of the three reference architectures.

To understand the deformation mechanisms that lead to these benefits, we extracted still frames from the experiments and computed the areal Jacobian J_i ≡ A_i/A_i,0 of each ith cell (the relative change in void area; Fig. 3c) where A_i,0 and A_i denote the undeformed and deformed cell areas, respectively⁴⁷. The square lattice was prone to buckling of vertically aligned ligaments, resulting in the observed local maxima in its stress–strain response (Fig. 3a, blue curve). Buckling events began at a single ligament, followed by a horizontal cascade of buckling events that spanned the whole length of the specimen (Fig. 3c(i)), leading to a sudden drop in its load capacity. This process repeated in a layer-by-layer fashion until the entire specimen buckled. Similarly, the hexagonal lattice (teal; Fig. 3c(ii)) also developed wide compaction bands that spanned the width of the lattice. Although the Voronoi tessellation did not succumb to well-defined shear or compaction bands owing to its lacking periodicity, consistent with previous experimental findings⁴⁸, it exhibited localized compaction near the bottom, ultimately forming a domain of compacted cells that spanned the width of the specimen (purple; Fig. 3c(iii)), whereas unit cells near the top remained largely undeformed.

In contrast, SR architectures lack compaction or shear bands, which can be attributed to the local coexistence of ligaments of different lengths and thicknesses; hence, the failure of the weakest ligament in one domain will not trigger the propagation of contiguous failure pathways (Fig. 3c(iv)). To understand the origin of this diffuse damage, we computationally evaluated the susceptibility of each unsupported ligament section to buckling under a macroscopic deformation (Supplementary Information, section 2.13). We find that the network consists of ligaments characterized by a broad distribution of critical buckling loads, and hence a local buckling event does not trigger a global cascade of instability, explaining the stable macroscopic response. Moreover, selectively removing a range of length scales degrades the load-bearing capacity of the network (Supplementary Information, section 2.14), suggesting that length-scale heterogeneity underlies the robustness of the system.

To further understand this delocalization, we introduce a localization parameter

where \({\mu }_{i}\equiv \frac{1}{M}{\sum }_{j\in {{\mathcal{B}}}_{r}(i)}{J}_{j}\) is the local average Jacobian of cell i, defined over a neighbourhood \({{\mathcal{B}}}_{r}(i)\) of i containing M cells (Supplementary Information, section 2.10), and N is the total number of cells. As a result, \({\mathcal{L}}\in [0,1]\), where \({\mathcal{L}}=0\) implies that the deformation in every neighbourhood is identical—the ideal case—and \({\mathcal{L}}=1\) indicates that the deformation variation between neighbourhoods is comparable to the deformation variation among all individual cells, corresponding to extreme localization.

For all samples, \({\mathcal{L}}\) initially increases with increasing strain (Fig. 3d); the square and hexagonal lattices exhibit rapidly increasing and high asymptotic \({\mathcal{L}}\) owing to the presence of compaction bands, whereas the disordered Voronoi lattice has lower \({\mathcal{L}}\). Consistently, the SR samples show the lowest value of \({\mathcal{L}}\)—with the localization within \({\mathcal{L}}\le 0.2\) through 25% macroscopic strain, which is 37% lower than the Voronoi geometry—confirming the material’s ability to uniformly distribute damage and offer a macroscopically stable load-bearing capacity. We can further reduce \({\mathcal{L}}\) by selectively nucleating the ligaments in an SR network (Supplementary Information, section 2.12) using a ‘checkerboarding’ approach, allowing us to constrain \({\mathcal{L}}\le 0.08\) through 25% strain, representing a 74% reduction over the Voronoi value.

Efficient dynamic impact absorption

Looking beyond the quasi-static regime, an ideal impact-absorbing material maximizes the energy dissipated during the impact event while minimizing the peak force⁴⁹. We performed impact experiments on SR and reference architectures at relative density \(\bar{\rho }=40 \% \), using a rigid impactor released with controlled initial energy corresponding to an instantaneous strain rate of \(\mathop{\varepsilon }\limits^{\cdot }=100\,{{\rm{s}}}^{-1}\) at impact.

Using a high-speed camera, we analysed frames during the impact event (Fig. 4a), finding that similar to the quasi-static experiments, the square lattice fails in a layer-by-layer fashion with fractured ligaments localized within sequentially forming horizontal crushing bands. The square lattice showed a sharp initial peak force that dropped to nearly zero followed by subsequent peaks corresponding to local compactions (Fig. 4b, blue curve). The hexagonal lattice also exhibited extreme localization owing to the fracture of diagonally oriented planes in as little as 2 ms post-impact. Although the Voronoi tessellation diffused damage more effectively, it also showed clear localized crush bands, leading to a high maximum force (Fig. 4b, purple curve). The SR geometry, by contrast, lacks not only well-defined crushing bands but also diffused damage through the structure, leading to the most stable force–time behaviour (Fig. 4b, orange curve), with low initial peak force and a lower compaction force than the Voronoi geometry.

a, The SR system efficiently distributes damage under dynamic impact loading, captured by high-speed photography. Scale bars, 10 mm. b, Force–time data reveal a flat, plateau-like force–time response in the SR materials, compared with the reference architectures (square, hexagonal, triangle, sea-sponge⁵³, and Kagome lattices; Gilbert⁵⁵, Voronoi, and E-Voronoi⁴⁸ tessellations; and isotropic (stretch h = 1) k-nearest foams⁵⁴). c, This corresponds to simultaneous minimization of peak force and maximization of specific energy dissipation during impact. d, Anisotropy and localization under quasi-static loading. All systems have relative density \(\bar{\rho }=40 \% \). Localization corresponds to the experimentally observed value at 25% macroscopic strain and anisotropy is computed from finite-element simulations.

We extracted the overall peak force, F_max, and the change in kinetic energy of the impactor before and after the event, representing the energy dissipated by the architecture, allowing us to measure the specific energy dissipation as

where m_impactor is the mass of the impactor, m_specimen is the mass of the specimen, and v_i and v_f represent the impactor velocities before and after the impact event, respectively. In Fig. 4c, we plot W_sp against F_max, finding that the SR system occupies the optimal space with maximum specific energy dissipation W_sp and minimum peak force F_max among all reference geometries considered in this work, indicating that the advantageous quasi-static properties of SR materials offer clear benefits under dynamic loading as well.

Steering elastic waves using SR materials

SR metamaterials also support programmable elastic-wave propagation, owing to their widely tunable elastic anisotropy, fully parameterized with only two parameters (α, λ₀). To quantify direction-dependent wave propagation in SR metamaterials, we performed Bloch wave analysis on SR microstructures across the (α, λ₀) design space—here reinterpreted as representative volume elements in an infinite tiling (Fig. 5a and Supplementary Information, section 2.16). Leveraging the tunability of SR metamaterials, we designed acoustic metamaterials out of quasi-isotropic SR architectures. Using the effective refractive index concept from optics n_eff = c₀/v_s, which compares the wave velocity in the constituent material c₀ to the velocity v_s of a wave travelling through an SR metamaterial made of the same constituent material, we find that for SR metamaterials, n_eff spans 1.5 ≤ n_eff ≤ 8.3 in the (α, λ₀) configuration space, demonstrating a substantial range in achievable refractive indices compared with the hexagonal and square lattices, which span 1.0 ≤ n_eff ≤ 1.9 within the same relative-density limits (\(20 \% \le \bar{\rho }\le 75 \% \); Supplementary Information, section 2.17).

a, SR samples with their corresponding directional velocity polar plots, demonstrating a wide range of achievable space of wave velocities and anisotropy. Curves, matching the colour of their respective geometries, denote the quasi-longitudinal wave velocities, and the green curve represents the wave velocity of the monolithic material. Samples in the grey–black colour scale show the breadth of wave anisotropy achieved and the samples in the rainbow scale demonstrate isotropic SR unit cells used to construct the elastic lens in b. We use each sample as a unit cell, adding a frame with small gaps to provide boundary conditions for building the superstructure. The effective refractive index, n_eff, is defined as the ratio of the green to blue curves along the x direction. b, One-quarter of a Luneburg-inspired lens constructed from SR samples with varying refractive indices, highlighted by colour (blue, high refractive index; orange, low refractive index). c, Schematic illustrations of the target ray-focusing behaviour of an ideal Luneburg lens (top) and the designed refractive index profile of the full pixelated lens (bottom). d, Steady-state wavefield of the SR elastic Luneburg-like lens showing the elastic perturbation plane source (white) and focal point (red). The elastic wave is focused on the side opposite the source, as evidenced by the radial wave pattern emanating from the focal point. e, The normalized displacement intensity along the y-axis cross-section through the lens focal point (white dashed line in d). The Luneburg-like lens (solid blue curve) exhibits a sharp focal peak, whereas the randomized lens (inset, green dashed curve) shows no noticeable focusing. The residual intensity peak in the randomized lens—three orders of magnitude smaller—arises from weak energy trapping at the elastic circular boundary. A suitable excitation frequency was found in the range f ≈ 646.7 Hz, well within the long-wavelength regime of each pixel in the lens, where the refractive indices are still valid.

The refractive index tunability of SR microstructures enables multiple elastodynamic functionalities such as wave guiding or lensing, with applications in subwavelength imaging and energy harvesting^50,51. We designed a Luneburg-like cylindrical lens of diameter d in a pixelated gradient-index (GRIN) approach⁵² (Fig. 5b and in Methods), and validated it via steady-state dynamic finite-element method, focusing elastic perturbations from an incident plane wave (Fig. 5c). Displacement fields reveal that waves traversing the lens focus on the opposite side, producing semi-circular wave fronts at the focal point and indicating that SR lenses qualitatively follow Luneburg patterns⁵⁰. The displacement intensity (squared displacement magnitude) contour along the vertical line through the focal point shows a peak that is an order of magnitude larger than the side peaks, indicating energy concentration into a region approximately 0.1d wide (Fig. 5d). For comparison, a randomized reference lens constructed by randomly selecting SR cells at each radius produced a focal intensity roughly three orders of magnitude weaker, with a broader peak and side lobes of comparable magnitude, confirming that focusing arises from the GRIN arrangement of the SR pixels. Spanning such a breadth of n_eff enables the generation of graded materials with smoothly changing acoustic properties over a much larger range than any individual lattice family, achieved by tuning (α, λ₀).