Nano-architected GaN spinodoid metamaterials with tailorable anisotropic piezoelectric properties

Abstract

Piezoelectric properties of nanomaterials are often constrained by their intrinsic crystallographic structures. Inspired by spinodal phase separation, this study develops gallium nitride (GaN) spinodoid metamaterials with enhanced and anisotropic piezoelectric properties. Molecular dynamics simulations reveal that these metamaterials exhibit significantly improved piezoelectric stress and strain constants (e.g., d₃₃ enhanced by up to 12 times) and increased piezoelectric anisotropy (e.g., d₃₁ ≠ d₃₂) compared to bulk GaN. These enhancements in piezoelectric performance are strongly affected by their underlying nano-architecture, which is governed by the evolutionary time during spinodal decomposition. Due to the asymmetric topology designs, GaN spinodoid metamaterials can possess more independent non-zero piezoelectric stress/strain constants as well as elastic constants compared to the bulk piezoelectric GaN. Relative density is found to further modulate the piezoelectric properties and anisotropy of the nano-architected piezoelectric materials through the contribution of surface effects and tuning the surface-to-volume ratio. This work underscores the potential of topology engineering to overcome crystallographic constraints in piezoelectric nanomaterials, opening avenues for their applications in nano-energy harvesters and three-dimensional pressure mapping/sensing nano-devices.

---

New concepts

We demonstrate a spinodoid topology engineering strategy to construct nano-architected gallium nitride (GaN) metamaterials with programmable anisotropic piezoelectric properties. Inspired by spinodal decomposition, our approach generates GaN metamaterials with spontaneously formed, interconnected topologies. In contrast to conventional metamaterial design that relies on computer-aided geometric modeling, our construction process is both automatic and rapid, yet highly controllable: the topology is governed by evolution time, and the relative density is dictated by the phase ratio. A key breakthrough of this work is the demonstration that spinodoid nanostructures can transcend crystallographic symmetry constraints, enabling anisotropic piezoelectric responses (e.g., d₃₁ ≠ d₃₂) and achieving up to a twelve-fold enhancement in piezoelectric constants compared to bulk GaN. This topology-driven electromechanical coupling establishes a new design strategy for reprogramming intrinsic material properties via nanoscale architecture rather than through conventional methods such as doping, crystallographic growth control, and grain/domain alignment. This study not only advances the fundamental understanding of piezoelectricity in nano-engineered materials but also provides an efficient and generalizable strategy for designing metamaterials with tunable piezoelectric properties, opening opportunities for nano-energy harvesting, tactile sensing, and adaptive nanorobotic systems.

Introduction

Piezoelectric materials, which convert mechanical energy to electric energy and vice versa, have been widely used in energy harvesting, sensing, actuation, nanorobotics, and implantable biomedical devices.^1–10 However, the piezoelectric constants of these materials are intrinsically bounded by the crystallographic structure of their constituent material, leading to directionally limited piezoelectric responses.¹¹ This limitation poses challenges for applications requiring high directional pressure sensitivity.¹¹ Motivated by these needs, alternative strategies have been proposed for achieving piezoelectric anisotropy and directional responses. The first strategy, which is crystal engineering, involves modifying crystallographic structures and symmetry through doping or controlled growth techniques.^12–14 However, this technique encounters with challenges in achieving a precise control over crystallographic orientation and dopant distribution as well as limitations in selecting suitable doping agents that affect the anisotropic piezoelectric design space.¹⁵ The second approach, which is texture engineering, aligns grains and domains to enhance anisotropic piezoelectric properties.^16–18 In return, this strategy also negatively impacts phase stability¹⁹ and mechanical properties.²⁰ The metamaterial design as a third and more recent strategy involves engineering underlying structures to achieve tailored anisotropic piezoelectric responses, enhanced piezoelectric constants, and novel functionalities.^5,11,21–26 Although promising, this approach relies heavily on the complex design of topological materials and precise manufacturing.^11,22–24

Building upon these strategies, the advent of additive manufacturing has revolutionized the design of piezoelectric metamaterials,^{11,22,23,27,28} offering possibilities for achieving exceptional piezoelectric properties and multifunctionalities through the manipulation of topologies. Different from crystal engineering and texturing approach, additive manufacturing enables precise control over metamaterial architectures, expanding the design space for piezoelectric responses. For example, three-dimensional (3D) printed lead zirconate titanate (PZT) piezoceramic metamaterials exhibit customizable directional responses through structural topology manipulation.¹¹ This capability allows for the enhancement or even reversal of voltage outputs along specific directions in response to mechanical stress. In addition, 3D printed piezoelectric metamaterials have demonstrated high specific piezoelectric coefficients and electromechanical sensitivity.^11,21–23 For example, piezoelectric constant (d₃₃) of 3D printed PZT piezoceramic can reach 583 pC N⁻¹, comparable to pristine ceramic, while achieving significantly lower material densities.²² Similarly, 3D printed lead-free barium titanate (BaTiO₃) octet truss lattice has shown a d₃₃ as high as 849 pC N⁻¹, exceeding the corresponding values for BaTiO₃ ceramics achieved through material processing and sintering optimization.⁵ Beyond enhanced piezoelectric properties and piezoelectric anisotropy, these architected metamaterials offer additional functional advantages, with potential applications in self-sensing robotic systems with feedback control,²³ wireless self-sensing sporting equipment,²¹ and miniaturized ultrasound transducers for biomedical devices.²²

While considerable progress has been made in developing piezoelectric materials at the micro- and macroscale, the nanoscale domain remains relatively unexplored. Recent advances in nanoscale additive manufacturing have made it possible to fabricate nano-architected metamaterials,^29–31 unlocking new opportunities for tailoring piezoelectric properties at previously inaccessible length scales. At the nanoscale, unique mechanisms such as size and surface effects become increasingly significant, leading to substantial differences in both mechanical and piezoelectric behavior compared to their larger-scale counterparts.^32–41 For example, with the increased high surface-to-volume ratios, wurtzite zinc oxide (ZnO) and gallium nitride (GaN) nanowires (NWs) exhibit enhanced piezoelectric stress constant (e₃₃) and reduced relative dielectric constant (k₃₃), attributed to surface effects that are negligible at the micro- and macroscale.^32,33,36 In addition, size-dependent behavior has been observed in GaN nano-architected metamaterials: as the characteristic unit size decreases, d₃₃ decreases while d₁₅ increases.¹⁰ These nanoscale-specific behaviors necessitate further exploration of topology engineering, where the interplay between surface effects and topology design may lead to different piezoelectric phenomena, such as enhanced piezoelectric effect and directional anisotropy, that are not achievable in larger-scale systems.

In this study, we focus on GaN as a representative example of nanoscale piezoelectric materials and investigate its anisotropic piezoelectric properties through topological design using molecular dynamics (MD) simulation. Inspired by the natural self-assembly process, we construct different GaN spinodoid metamaterials based on the Cahn–Hilliard model of spinodal phase separation. Previous studies on cubic, octahedron, and triply periodic minimal surface (TPMS) GaN metamaterials revealed enhanced piezoelectric properties compared to bulk GaN, while maintaining symmetry-driven identical responses in [011̄0] and [2̄110] directions.¹⁰ Given the reported anisotropic mechanical properties of spinodoid metamaterials,^42,43 it is hypothesized that these metamaterials may exhibit directionally tailored piezoelectric properties. Our results indicate that topology designs can increase the number of independent piezoelectric (i.e., e_ij and d_ij) and elastic (C_ij) constants. The proposed GaN spinodoid metamaterials demonstrate anisotropic piezoelectric responses (e.g., e₃₁ ≠ e₃₂, d₃₁ ≠ d₃₂), which are absent in bulk GaN. Moreover, the piezoelectric properties of GaN spinodoid metamaterials are significantly enhanced compared to bulk GaN, with relative density playing a crucial role in modulating the piezoelectric anisotropy. This study highlights the potential of GaN spinodoid metamaterials for next-generation nano-engineered systems, enabling multifunctional applications in energy harvesting, tactile sensing, and adaptive nanodevices with directionally tailored responses.

Results and discussion

We utilize wurtzite GaN to construct nanoscale architectures with different topologies inspired by spinodal decomposition (Fig. 1). The atomic structure of wurtzite GaN is characterized by lattice parameters a = 3.21 Å, c = 5.20 Å, b = 5.55 Å, and u = uc/c = 0.378¹⁰ (Fig. 1a). Spinodal decomposition refers to a spontaneous phase separation process in binary fluid mixtures, where two components differentiate into two distinctive phases (phase a and phase b) through diffusion-driven mechanisms.^42,44 This process occurs when the system becomes thermodynamically unstable, often upon cooling below a critical temperature, leading to interconnected bicontinuous structures. After the phase separation, selective removal of one phase results in a porous network that retains the characteristic spinodal topology (Fig. 1a). The spinodal decomposition process is described by the Cahn–Hilliard equation^44,45where u(x, y, z, t) (u ∈ [−1, 1]) represents the concentration difference between the two phases, t is the evolutionary time, and θ = 0.5 defines the width of transition region between the phases. The free energy function f(u) is modeled as a double-well potential function f(u) = (u² − 1)²/4 in this study.⁴⁴

Fig. 1 Atomistic models of nano-architected GaN piezoelectric spinodoid metamaterials. (a) Topology generation through spinodal decomposition simulation and the lattice structure of wurtzite GaN. (b) Distribution patterns of phase concentration u in the top section (x–y plane) of the simulation box from t = 5000τ to 30 000τ (τ = 0.015). (c) Atomistic models of nano-architected GaN spinodoid metamaterials constructed based on Cahn–Hilliard equation at various evolutionary times t. The schematics depict nanoarchitectures characterized by L ≈ 20 nm and [small rho, Greek, macron] = ρ/ρ₀ ≈ 0.5.

The 3D spinodal topologies (or spinodoid metamaterials) are generated by solving eqn (1) in a cuboid simulation box, discretized into an N_x × N_y × N_z lattice, where N_x, N_y, N_z = 36, 62, and 39, respectively. Discretizing eqn (1) using a central approximation of spatial derivatives yields

where u^m_ijk represents the discrete values of u(ia, ja, ka, τ) at lattice node point (ia, ja, ka), with a as the mesh size and τ = 0.015 as the integration time step.⁴⁴ A discrete Laplacian operator, is utilized to approximate the continuous Laplacian in the lattice.⁴⁴ During the spinodal decomposition, phase separation occurs uniformly, producing an interconnected pattern that evolves over time. The Cahn–Hilliard equation is used here purely as a topology generation tool to obtain bicontinuous architectures that resemble spinodal morphologies. This is a mathematical construction process, not a physical simulation of spinodal decomposition in GaN. Periodic boundary conditions are applied in all three directions of the simulation box. Fig. 1b shows representative patterns of u in the top region (x–y plane) of the simulation box at different evolutionary times (from t = 5000τ to 30 000τ). To differentiate the phases (phase a and phase b) to generate 3D spinodal topologies, a standard step function is introduced as:

G^m_ijk = H(u^m_ijk − u^m_c) (3)

where u^m_c is a threshold value that distinguishes the phases. When u^m_ijk > u^m_c, G^m_ijk = 1, indicating that the lattice node point (ia, ja, ka) belongs to phase a, otherwise G^m_ijk = 0, corresponding to phase b. To achieve a specific relative density ([small rho, Greek, macron]) of a spinodal topology (phase a), u^m_c is adjusted to satisfyFig. 1c displays the nano-architected GaN spinodoid metamaterials at different evolutionary times characterized by L ≈ 20 nm and a relative density ρ/ρ₀ ≈ 0.5 (ρ₀ = 6.1 g cm⁻³ is the density of bulk GaN), which are utilized to study the topology effect on piezoelectric properties. The two phases generated by the Cahn–Hilliard equation represent mathematical partitions of space used to define bicontinuous architectures. Only one phase is retained and assigned GaN material properties, while the other is removed. Therefore, phase a and phase b do not coexist or contribute simultaneously to the piezoelectric response; all electromechanical behavior arises from the geometry of the retained GaN topology. As shown in Fig. 1b and c, increasing the evolutionary time leads to more continuous regions of separated phases with reduced intermixing. As a result, the surface area between the phases, or the surface-to-volume ratio, is reduced over evolutionary times (see Fig. S1a in SI).

After constructing the nano-architected GaN spinodoid metamaterials, we proceed to investigate their piezoelectric characteristics. To begin with, the piezoelectric stress constants of bulk GaN are determined through MD simulation (see Computational Methods section). Wurtzite GaN exhibits three independent nonzero piezoelectric stress constants: e₃₃, e₃₁ (= e₃₂), and e₁₅ (= e₂₄).⁴⁶ In this study, we focus specifically on e₃₃, e₃₁, and e₃₂. The computed values for bulk GaN are e₃₃ = 0.895 C m⁻², e₃₁ = −0.45 C m⁻², and e₃₂ = −0.457 C m⁻², aligning well with existing experimental and computational reports (see Table S1 in SI).^{10,36,39,47–49} Next, we evaluate the piezoelectric stress constants for the GaN spinodoid metamaterials. Fig. 2a shows the polarization change (ΔP₃) along the [0001] crystallographic direction in response to uniaxial strain applied in the same direction, comparing the results for bulk GaN and GaN spinodoid metamaterials at various evolutionary times. A schematic of the external strain used to induce the polarization change is also presented in Fig. 2a. In addition, the polarization change (ΔP₃) under the strain applied in the [011̄0] or [2̄110] directions, used to compute e₃₁ and e₃₂, respectively, is shown in Fig. S2 (SI). The finding reveals that ΔP₃ increases with applied strain and is influenced by evolutionary times.

Fig. 2 Piezoelectric properties of nano-architected GaN spinodoid metamaterials. (a) The change of the polarization for GaN spinodoid metamaterials (phase a) and bulk GaN as a function of mechanical strain, where external normal strains are applied along the [0001] direction. (b) The piezoelectric stress constant e₃₃ of bulk GaN and GaN spinodoid metamaterials at different evolutionary times (t₁, t₂, t₃, t₄, t₅, and t₆ correspond to 2500τ, 5000τ, 10 000τ, 15 000τ, 20 000τ, and 25 000τ, respectively. τ = 0.015). (c) Piezoelectric stress constants (e₃₁ and e₃₂) of bulk GaN and GaN spinodoid metamaterials. (d) The ratio between e₃₁ and e₃₂ for bulk GaN and GaN spinodoid metamaterials at different evolutionary times. (e) Elastic constants (C₁₁, C₂₂, C₁₃, and C₂₃) of bulk GaN and GaN spinodoid metamaterials.

To further characterize the piezoelectric performance of GaN spinodoid metamaterials, we calculate the piezoelectric stress constants (e₃₃, e₃₁, and e₃₂). It is important to note that piezoelectric stress constants reported in Fig. 2b and c are proper piezoelectric stress constants⁴⁶ (see Computational methods section). As shown in Fig. 2b and c, e₃₃ is significantly enhanced in GaN spinodoid metamaterials compared to bulk GaN, with certain spinodoid topologies achieving an increase of at least 65%. In contrast, the magnitudes of e₃₁ and e₃₂ are decreased in GaN spinodoid metamaterials compared to bulk GaN (Fig. 2c). In addition, we find that the piezoelectric stress constants of GaN spinodoid metamaterials depend on evolutionary time of spinodal decomposition. As shown in Fig. 1c, the spinodoid topology evolves over time, with the surface-to-volume ratio decreasing as the evolutionary time increases (Fig. S1 in SI). At the nanoscale, the surface effect introduced by the free surface of the nanostructure significantly influence piezoelectric behavior.^10,50 Specifically, surface atoms carry lower partial charges compared to bulk atoms, leading to reduced surface polarization and negatively impacting piezoelectric properties.^10,51 Fig. S10a–c in the SI compare the piezoelectric stress constants of GaN metamaterials with and without accounting for reduced surface partial charges, showing that consideration of surface charge reduction consistently lowers the predicted piezoelectric response. On the other hand, the presence of free surfaces also induces notable surface contraction at the nanoscale, reducing the overall volume and thereby enhancing piezoelectricity. The interplay between these two competing mechanisms, i.e., reduced surface polarization and volume contraction, ultimately governs the piezoelectric response of GaN spinodoid metamaterials at the nanoscale. For example, in Fig. 2b, the maximum e₃₃ for phase a is obtained at t = 5000τ (surface-to-volume ratio = 0.491) rather than at t = 2500τ (surface-to-volume ratio = 0.631), highlighting the trade-off between these competing effects.

It is also interesting to note the anisotropic piezoelectric behaviors observed in GaN spinodoid metamaterials. While bulk GaN has three independent piezoelectric stress constants (e₃₃, e₃₁ = e₃₂, and e₁₅ = e₂₄), GaN spinodoid metamaterials, through topology design, possess more independent piezoelectric stress constants. Fig. 2c and d compare e₃₁ and e₃₂ of bulk GaN and GaN spinodoid metamaterials at different evolutionary times. Different from bulk GaN, where e₃₁ and e₃₂ are almost identical, these constants are different and independent in GaN spinodoid metamaterials. Previous study on GaN metamaterials with topologies from cubic, octahedron, and TPMS families showed that e₃₁ = e₃₂, same to bulk GaN.¹⁰ We also find that while surface effects influence the magnitude of piezoelectric constants, they do not significantly affect the anisotropy of the piezoelectric response (Fig. S10d).

With the increased independent piezoelectric stress constants through topology design, the number of independent elastic constants is also expected to be increased in GaN spinodoid metamaterials. For Wurtzite crystal structures, five independent nonzero elastic constants are typically present: C₁₁, C₁₂, C₁₃ = C₂₃, C₃₃, and C₄₄ = C₅₅. These constants must satisfy specific mechanical stability criteria: C₁₁ > 0, C₁₁ − C₁₂> 0, C₄₄ > 0, and (C₁₁ + C₁₁) C₃₃ > 2 C²₁₃.⁵² In this study, MD simulations yield the following elastic constants for bulk GaN: C₁₁ = 353.58 GPa, C₁₂ = 139.93 GPa, C₁₃ = 123.89 GPa, C₃₃ = 369.53 GPa, and C₄₄ = 97.18 GPa, which satisfy the stability criteria and are in good agreement with previously reported results from DFT,^53,54 MD,^36,55 and experiment⁵² (Table S3). However, GaN spinodoid metamaterials are found to have nine independent elastic constants, i.e., C₁₁, C₁₂, C₁₃, C₂₃, C₃₃, C₄₄, C₅₅, and C₆₆ (see Tables S4 and S5 in SI). Fig. 2e demonstrates that, in GaN spinodoid metamaterials, C₁₁ and C₁₂ (C₁₃ and C₂₃) are different, whereas they are identical in bulk GaN, highlighting the potential of topology design to engineer mechanical anisotropy.

Having understood the effect of nanoscale topology and evolutionary time on the piezoelectric stress constants, we now explore another essential parameter: piezoelectric strain constants (d₃₃, d₃₁, and d₃₂). These piezoelectric constants quantify the mechanical strain generated in response to an applied electric field.¹¹ As shown in Fig. 3a, when an electric field is applied along the [0001] direction, the GaN spinodoid metamaterial (t = 5000τ, [small rho, Greek, macron] = 0.5) undergoes elongation or contraction along the same direction. The resulting electric field-induced mechanical strains are captured through MD simulations, revealing the strain response as a function of the applied electric field strength. Fig. 3b and d demonstrate that GaN spinodoid metamaterials stretch along the field direction (i.e., [0001]) while contracting along the transverse directions (i.e., [011̄0] or [2̄110]). By analyzing the relationship between the induced strain (ε_zz, ε_xx, and ε_yy) and the applied electric field strength (E_z), the corresponding piezoelectric strain constants, namely d₃₃, d₃₁, and d₃₂, can be determined.

Fig. 3 Piezoelectric strain constants of GaN spinodoid metamaterials at different evolutionary times. (a) An electric field applied to GaN spinodoid metamaterial (phase a, t = 2500τ) along the [0001] direction, with corresponding atomic configurations of GaN metamaterial showing compression and stretching. (b) Electric field-induced strains (ε_zz) in bulk GaN and GaN spinodoid metamaterials at different evolutionary times as a function of the applied electric field along the [0001] direction. The inset shows a schematic of the applied external electric field and the recorded strain from MD simulations used to calculate the piezoelectric strain constants (d₃₃). (c) Piezoelectric strain constants (d₃₃) of bulk GaN and GaN spinodoid metamaterials obtained from theoretical analysis and MD simulations. (d) Electric field-induced strains (ε_xx and ε_yy) in bulk GaN and GaN spinodoid metamaterials as a function of the applied electric field along the [0001] direction. The inset shows a schematic of the applied electric field and the recorded strains from MD simulation used to calculate the piezoelectric strain constants (d₃₁ and d₃₂). (e) Piezoelectric strain constants (d₃₁ and d₃₂) of bulk GaN and GaN spinodoid metamaterials obtained from theoretical analysis and MD simulations. (f) Comparison of d₃₁ and d₃₂ among GaN spinodoid metamaterials, BCC GaN,¹⁰ octet truss GaN,¹⁰ gyroid GaN,¹⁰ and bulk GaN. The data highlights the anisotropic piezoelectric properties of GaN spinodoid metamaterials due to asymmetric topology designs, as evidenced by deviations from the diagonal line (d₃₁ = d₃₂).

We employ two approaches to obtain the piezoelectric strain constants: direct MD simulations and theoretical calculations based on the elastic constants C_ij and piezoelectric stress constants e_ij⁵⁶

Detailed formulations for the computation of piezoelectric strain constants are provided in SI Section S2. The theoretical results obtained from eqn (5) can be utilized to verify our MD simulation results. It is important to note that the elastic constants of GaN spinodoid metamaterials are determined with respect to their cubic volume V₀, while polarization is computed based on the architectural volume (V = [small rho, Greek, macron]V₀).^10,32 Therefore, the piezoelectric stress constants need be appropriately adjusted (i.e., e_ij[small rho, Greek, macron]) prior to calculating the piezoelectric strain constants using eqn (5). For bulk GaN, MD simulation yields d₃₃ = 3.605 pm V⁻¹, d₃₁ = −1.799 pm V⁻¹, and d₃₂ = −1.793 pm V⁻¹ (d₃₁ ≈ d₃₂), aligning closely with the theoretical predictions (d₃₃ = 3.674 pm V⁻¹, d₃₁ = −1.784 pm V⁻¹, and d₃₂ = −1.959 pm V⁻¹) and previously reported experimental data (Table S2 in SI).^56–58 In addition, the expected relationship between d₃₁ and d₃₃, d₃₁ = −1.793 pm V⁻¹ ≈ −1/2 d₃₃ = −1.803 pm V⁻¹, is well preserved.⁵⁶ The strong consistency confirms the reliability of our MD simulation methodology for evaluating the piezoelectric properties of both bulk GaN and GaN spinodoid metamaterials.

In Fig. 3c and e, we report the piezoelectric strain constants (i.e., d₃₃, d₃₁, and d₃₂) for GaN spinodoid metamaterials, demonstrating excellent agreement between theoretical and MD results. Compared with piezoelectric strain constants of GaN, nanoarchitecture engineering significantly enhances piezoelectric strain constants (i.e., d₃₃, d₃₁, and d₃₂). For example, the lowest and highest d₃₃ obtained in GaN spinodoid metamaterials are 10.908 pm V⁻¹ (phase a, t₆) and 17.454 pm V⁻¹ (phase b, t₁), approximately 3× and 5× greater than bulk GaN (3.605 pm V⁻¹). Notably, d₃₃ of GaN spinodoid metamaterial (phase b, t₁) is the highest piezoelectric strain constant ever reported in GaN metamaterials, exceeding previous study (d₃₃ = 15.469 pm V⁻¹ for octet truss GaN metamaterials at [small rho, Greek, macron] ≈ 0.5).¹⁰ For d₃₁ and d₃₂, Fig. 3e and Fig. S3 in SI highlight two key findings. First, both d₃₁ and d₃₂ are enhanced in GaN spinodoid metamaterials. The highest |d₃₁| and |d₃₂| among the GaN spinodoid metamaterials are 4.344 pm V⁻¹ (phase b, t₁) and 6.426 pm V⁻¹ (phase a, t₂), respectively, which is much higher than bulk GaN (d₃₁ = −1.799 pm V⁻¹, and d₃₂ = −1.793 pm V⁻¹). Second, the asymmetric topology of GaN spinodoid metamaterials results in that d₃₁ ≠ d₃₂, different from bulk GaN and other symmetric nanoarchitectures (e.g., topologies from cubic, octahedron, and TPMS families),¹⁰ where d₃₁ = d₃₂. This study is the first report of anisotropic piezoelectric properties in GaN metamaterials achieved through asymmetric topology design at the nanoscale.

To better show the piezoelectric anisotropy in GaN metamaterials, we compare d₃₁ and d₃₂ of GaN spinodoid metamaterials, BCC GaN,¹⁰ octet truss GaN,¹⁰ gyroid GaN,¹⁰ and bulk GaN in Fig. 3f and Fig. S4 in SI. The relative densities of all nano-architected GaN metamaterials are identical ([small rho, Greek, macron] ≈ 0.5). As shown in Fig. 3f, only GaN spinodoid metamaterials exhibit points significantly deviating from the diagonal line, indicating anisotropic piezoelectric properties in [011̄0] and [2̄110] directions. These enhancements and anisotropy (d₃₁ ≠ d₃₂) suggest that GaN spinodoid metamaterials are not only suitable for nano-energy harvesting but also enable new functionalities, such as 3D pressure sensing, mapping, and directionality detection.¹¹

While we have found that GaN spinodoid metamaterials characterized by [small rho, Greek, macron] = 0.5 exhibit enhanced and anisotropic piezoelectric properties, we now investigate the potential for further enhancement and anisotropy by tuning the relative densities. Using the method introduced in Section 2.1, we construct GaN spinodoid metamaterials with varying [small rho, Greek, macron] at t = 5 × 10³τ (Fig. 4a). Previous studies on GaN NWs with different cross-sectional shapes have demonstrated the crucial role of the free surface in influencing the piezoelectric stress constant, with a focus on e₃₃.^{32–34,36,59} These studies employ a core–shell model to explain the size effect, revealing that e₃₃ increases as the cross-sectional size decreases (i.e., surface-to-volume ratio increases).^{32–34,36,59} For GaN spinodoid metamaterials at the same evolutionary time, tuning the relative density is an efficient method to adjust the surface-to-volume ratio (Fig. 4b), thereby modifying their piezoelectric properties. Fig. 4b shows that as [small rho, Greek, macron] decreases, the surface-to-volume ratio increases, which is expected to enhance piezoelectric properties.^32,33,36 At the nanoscale, the influence of increasing surface-to-volume ratio is not purely monotonic, because free surfaces introduce two competing effects on piezoelectricity.¹⁰ Reduced atomic coordination at free surfaces leads to deviations in bond lengths, bond angles, and charge distribution (Fig. S13),^10,47 which suppress local polarization and tend to reduce piezoelectric constants. In contrast, surface relaxation induces contraction of the spinodoid architecture, resulting in a reduced effective volume (Fig. S12). Because polarization scales inversely with volume, this contraction enhances the overall piezoelectric response. As a result, the net piezoelectric behavior of GaN spinodoid metamaterials reflects a competition between reduced surface polarization and topology-induced volume contraction.

Fig. 4 Effect of relative density on the piezoelectric properties of GaN spinodoid metamaterials. (a) Atomic configurations of GaN metamaterials (t = 5 × 10³τ, phase a) characterized by varying relative densities. (b) Ratio of surface atoms of GaN metamaterials as a function of [small rho, Greek, macron]. (c) Polarization changes as a function of mechanical strains for bulk GaN and GaN spinodoid metamaterials (phase a) characterized by different [small rho, Greek, macron], where the external normal strain is applied along the [0001] direction (see the inset). (d) Piezoelectric stress constants (e₃₃, e₃₁, and e₃₂) of bulk GaN and GaN metamaterials characterized by different [small rho, Greek, macron]. (e) Ratio between e₃₁ and e₃₂ of Bulk GaN and GaN metamaterials characterized by different [small rho, Greek, macron]. (f) Elastic constants of bulk GaN and GaN spinodoid metamaterials (phase a) characterized by different [small rho, Greek, macron].

Next, we examine the piezoelectric properties (i.e., e_ij and d_ij) of GaN spinodoid metamaterials with varying [small rho, Greek, macron]. The polarization changes (ΔP₃) of GaN spinodoid metamaterials subjected to strain along the [0001] direction indicate that the slope of the strain-induced ΔP₃ increase as [small rho, Greek, macron] decreases, signifying an enhancement in e₃₃ (Fig. 4c). The ΔP₃ of GaN spinodoid metamaterials characterized by different [small rho, Greek, macron] with strain applied along the [011̄0] or [2̄110] directions are presented in Fig. S5 in SI. Based on the results shown in Fig. 4c and Fig. S5, we calculate the piezoelectric stress constants (e₃₃, e₃₁, and e₃₂) for GaN spinodoid metamaterials (both phase a and phase b) with different [small rho, Greek, macron]. The results, presented in Fig. 4d, reveal that e₃₃, e₃₁, and e₃₂ all increase with the decreased [small rho, Greek, macron]. Furthermore, deviations between e₃₁ and e₃₂ are observed, highlighting the anisotropic nature of the piezoelectric properties in GaN spinodoid metamaterials. To better understand the effect of [small rho, Greek, macron] on piezoelectric stress constant anisotropy, Fig. 4e shows the ratio between e₃₁ and e₃₂ as a function of [small rho, Greek, macron]. The results indicate that anisotropy is enhanced with the decreased [small rho, Greek, macron], as the ratio between e₃₁ and e₃₂ increases. While at very high relative densities, where the surface-to-volume ratio approaches that of bulk GaN, both the enhancement in piezoelectric magnitude and the degree of anisotropy are reduced, consistent with diminished topological and surface contributions. In addition to piezoelectric stress constants, we explore the effect of [small rho, Greek, macron] on the elastic constants. Fig. 4f presents the elastic constants C₁₁ and C₂₂ (C₁₃ and C₁₂) for GaN spinodoid metamaterials with varying [small rho, Greek, macron]. Different from bulk GaN, where C₁₁ = C₂₂ and C₁₃ = C₁₂, the different mechanical responses are observed in the x- and y-directions in GaN spinodoid metamaterials. The elastic constants of GaN spinodoid metamaterials with varying [small rho, Greek, macron] are detailed in Tables S6 and S7 in SI.

Finally, in Fig. 5, we report the piezoelectric strain constants of GaN spinodoid metamaterials characterized by different [small rho, Greek, macron]. Under the same electric field applied in [0001] direction, the electric field-induced strain in GaN spinodoid metamaterials is increased with the decreased [small rho, Greek, macron] (Fig. 5a). In addition to the d₃₃ obtained directly from MD simulation, we also calculate d₃₃ based on eqn (5); the theoretical predictions and MD results show excellent agreement for GaN spinodoid metamaterials with different [small rho, Greek, macron] (Fig. 5b). A remarkable enhancement in d₃₃ is found with the decreased [small rho, Greek, macron]. The highest d₃₃ found in phase a and phase b GaN spinodoid metamaterials are 26.76 pm V⁻¹ and 44.39 pm V⁻¹, respectively, at [small rho, Greek, macron] = 0.3, which are approximately over seven and twelve times of bulk GaN (d₃₃ = 3.605 pm V⁻¹). d₃₃ = 44.39 pm V⁻¹ is the highest piezoelectric strain constants ever obtained in nano-architected GaN metamaterials with [small rho, Greek, macron] = 0.3, while the d₃₃ obtain in previous study is 19.10 pm V⁻¹, 34.87 pm V⁻¹, and 22.59 pm V⁻¹ in BCC GaN, octet truss GaN, and Gyroid GaN metamaterials ([small rho, Greek, macron] ≈ 0.3), respectively.¹⁰ Piezoelectric strain constants (i.e., d₃₃, d₃₁, and d₃₂) of phase b GaN spinodoid metamaterials with different [small rho, Greek, macron] can be found in Fig. S6 in SI.

Fig. 5 Piezoelectric strain constants of GaN spinodoid metamaterials characterized by different relative densities. (a) Electric field-induced strains (ε_zz) in bulk GaN and GaN spinodoid metamaterials with different [small rho, Greek, macron] as a function of applied electric field along the [0001] direction. The inset shows a schematic of the applied electric field and the recorded strain from MD simulation used to calculate the piezoelectric strain constants (d₃₃) (b) Piezoelectric strain constants (d₃₃) of bulk GaN and GaN spinodoid metamaterials with different [small rho, Greek, macron] obtained from theoretical analysis and MD simulations. (c) Electric field-induced strains (ε_xx and ε_yy) in bulk GaN and GaN spinodoid metamaterials as a function of the applied electric field along the [0001] direction. The inset shows a schematic of the applied electric field and the recorded strain from MD simulation used to calculate the piezoelectric strain constants (d₃₁ and d₃₂). (d) Piezoelectric strain constants (d₃₁ and d₃₂) of bulk GaN and GaN spinodoid metamaterials characterized by different [small rho, Greek, macron] obtained from theoretical analysis and MD simulations. (e) Comparison of d₃₁ and d₃₂ among GaN spinodoid metamaterials, BCC GaN,¹⁰ octet truss GaN,¹⁰ gyroid GaN,¹⁰ and bulk GaN. The data highlights the anisotropic piezoelectric properties of GaN spinodoid metamaterials due to asymmetric topology designs, as evidenced by deviations from the diagonal line.

The similar enhancement resulting from the decreased [small rho, Greek, macron] is also found in d₃₁ and d₃₂ (Fig. 5c and d). Both |d₃₁| and |d₃₂| are found significantly increased with the decreased [small rho, Greek, macron]. For example, as shown in Fig. 5d, the minimum |d₃₁| and |d₃₂| are 12.79 pm V⁻¹ and 7.72 pm V⁻¹, respectively, when [small rho, Greek, macron] = 0.3. Interestingly, the difference between d₃₁ and d₃₂ in GaN spinodoid metamaterials is found increased with the decreased [small rho, Greek, macron], indicating the enhanced piezoelectric anisotropy (Fig. 5d and Fig. S7 in SI). In Fig. 5e, we compare such anisotropy found in GaN spinodoid metamaterials at different [small rho, Greek, macron] with bulk GaN and other nano-architected GaN metamaterials (i.e., BCC GaN,¹⁰ octet truss GaN,¹⁰ gyroid GaN¹⁰ at different [small rho, Greek, macron]). It is clear to found that only GaN spinodoid metamaterials exhibit points significantly deviating from the diagonal line, indicating anisotropic piezoelectric properties, which are enhanced with the decreased [small rho, Greek, macron].

Note that the anisotropy in the piezoelectric properties of GaN spinodoid metamaterials arises from their asymmetric topology designs. To further investigate this asymmetry, we then examine the mechanical properties of GaN spinodoid metamaterials, as these are intrinsically linked to their piezoelectric behavior. Fig. 6a compares the stress–strain curves for bulk GaN (black curve) and GaN spinodoid metamaterial (t = 5000τ, [small rho, Greek, macron] = 0.5) (blue curve) under tension applied along the [011̄0], [2̄110], and [0001] directions, respectively. The results show that the stiffness, ultimate tensile strength, and failure strain of the GaN spinodoid metamaterials are significantly decreased compared to bulk GaN. In bulk GaN, stress is uniformly distributed under tension, whereas in spinodoid metamaterials, stress concentration is observed, which leads to decreased mechanical properties (Fig. 6b). Interestingly, the stress–strain curves reveal different mechanical responses in the [011̄0] and [2̄110] directions. For bulk GaN, the mechanical responses are nearly identical in the [011̄0] and [2̄110]. directions, with E_[011̄0] = 271.69 GPa and E_[2̄110] = 262.09 GPa, σ_u[011̄0] = 50.25 GPa and σ_u[2̄110] = 49.19 GPa. However, for GaN spinodoid metamaterial, significant differences are observed: E_[011̄0] = 62.59 GPa and E_[2̄110] = 22.71 GPa, σ_u[011̄0] = 5.34 GPa and σ_u[2̄110] = 2.33 GPa.

Fig. 6 Mechanical properties of GaN spinodoid metamaterials. (a) Stress–strain curves of GaN spinodoid metamaterial (phase a, t = 5000τ, and [small rho, Greek, macron] = 0.5) and bulk GaN subjected to tension in different directions. The black and blue curves represent bulk GaN and GaN spinodoid metamaterial, respectively. (b) von Mises (VM) stress distributions in bulk GaN and GaN spinodoid metamaterial (phase a, t = 5000τ, and [small rho, Greek, macron] = 0.5) at a tensile strain ε_zz = 0.1. Stress–strain curves of GaN spinodoid metamaterials characterized by different [small rho, Greek, macron] with tensile strain applied in the (c) [0001] and (d) [011̄0] or [2̄110] directions. (e) 2D polar plots of stiffness for GaN spinodoid metamaterials ([small rho, Greek, macron] = 0.3 and [small rho, Greek, macron] = 0.7) and bulk GaN. (f) Universal anisotropy index (A_U) of bulk GaN and GaN spinodoid metamaterials characterized by different [small rho, Greek, macron]. If a material is isotropic, A_U = 0. (g) Normalized stiffness and (h) normalized ultimate tensile stress of GaN spinodoid metamaterials as a function of [small rho, Greek, macron]. E_B and σ_uB represent the stiffness and ultimate tensile stress of bulk GaN in specific directions. (i) Young's modulus Ashby map comparing different scale architected metamaterials as well as commercial bulk materials. PC represents pyrolytic carbon.

Since relative density significantly affects piezoelectric properties (Fig. 4 and 5), we next explore its effect on their mechanical properties (Fig. 6c and d). As expected, stiffness and ultimate tensile strength decrease with the decreased relative density. The mechanical differences between the [011̄0] and [2̄110] directions persist across all densities studied (Fig. 6d). Detailed data on stiffness and strength are provided in Tables S8 and S9 in SI. To clarify elastic anisotropy, Fig. 6e presents 2D polar plots of stiffness for GaN spinodoid metamaterials ([small rho, Greek, macron] = 0.3 and [small rho, Greek, macron] = 0.7) and bulk GaN. In bulk GaN, the plot forms a perfect circle in the x–y plane (yellow curve), indicating isotropy. For GaN spinodoid metamaterial with [small rho, Greek, macron] = 0.7, the plot becomes elliptical, reflecting anisotropy with maximum and minimum stiffness along the ellipse's major and minor axes. When [small rho, Greek, macron] decreases further to 0.3, the polar plot assumes a propeller-like shape in the x–y plane, with pronounced difference in stiffness between the [011̄0] and [2̄110] directions, indicating extreme anisotropy. This increasing mechanical anisotropy with decreasing density parallels the piezoelectric anisotropy observed in Fig. 4e and Fig. S7 in SI. To quantify elastic anisotropy, we calculate the universal anisotropy (A_U) for GaN spinodoid metamaterials with different [small rho, Greek, macron] (Fig. 6f and SI Section S3). As shown, A_U values deviate from zero in GaN spinodoid metamaterials, confirming anisotropic elasticity. We also find that lower relative densities correspond to greater deviations, indicating increased anisotropy, which is consistent with piezoelectric anisotropy observed in Fig. 4e.

Fig. 6g and h show the variation of Young's modulus and ultimate tensile strength, normalized by their bulk counterparts (E_B and σ_uB), for different [small rho, Greek, macron]. We find that the normalized Young's modulus and strength of GaN spinodoid metamaterials follow a power-law scaling with relative density (E ∝ E_B[small rho, Greek, macron]^m and σ_u ∝ σ_uB[small rho, Greek, macron]ⁿ) with components (m and n) ranging between 2 and 3. While the performance is inferior to ideal stretching-dominated structures, it aligns with traditional lightweight and bending-dominated materials. Finally, Fig. 6i places these material properties on a Young's modulus versus density material selection chart, together with other commercially available bulk materials and reported ultralight metamaterials.^60–71 While carbon-based metamaterials (e.g., carbon Neovius shells⁷¹ and pyrolytic carbon (PC) octet truss⁶⁸) achieve comparable stiffness at lower densities, GaN spinodoid metamaterials outperform metal-based architected materials (e.g., Cu, Ni, and Au)^62,65,67 in Young's modulus. The observed reduction in stiffness and ultimate strength compared to bulk GaN arises from the porous spinodoid architecture and associated free-surface effects, which concentrate stress within slender ligaments.⁷² While these mechanisms govern the overall mechanical degradation, the asymmetric spinodoid topology further induces direction-dependent stress redistribution, leading to elastic anisotropy that is critical for the anisotropic piezoelectric response.

Computational methods

MD simulation is conducted using the large-scale atomic molecular massively parallel simulator (LAMMPS)⁷³ to study the mechanical and piezoelectric properties of bulk GaN and GaN spinodoid metamaterials. The Open Visualization Tool (OVITO)⁷⁴ is employed to visualize the evolution of atomic structures. The interactions between gallium and nitrogen atoms are described by the Stillinger–Weber (SW) potential⁷⁵ in this work, which has been widely employed to evaluate the mechanical and piezoelectric properties of wurtzite GaN.^10,34,36 The potential E contains a two-body term φ₂ and a three-body term φ₃ as follows:where subscripts i, j, k represent different atoms; r_ij is the distance between atoms i and j; θ_ijk is the angle between the atomic bonds ji and jk; a is the cut-off distance; δ is the cohesive energy of the bond; d is a length-unit parameter; A, B, C, and γ are dimensionless parameters. The potential parameters used in this study can be found in ref. 76. Periodic boundary conditions are applied in all three directions. Energy minimization is performed by the conjugate gradient algorithm. The integration time step is set as 0.5 fs. The bulk GaN and the constructed GaN spinodoid metamaterials are fully relaxed at 300 K in the isothermal–isobaric ensemble (NPT) for 500 000 steps (250 ps), during which the potential energy converged smoothly to a plateau, confirming equilibration (Fig. S8). To further assess stability across a broader thermal range, a representative GaN spinodoid metamaterial (t = 10 000τ, L ≈ 20 nm, and [small rho, Greek, macron] = ρ/ρ₀ ≈ 0.5) is relaxed at 100, 200, 300, 400, and 500 K. In all cases, the potential energy stabilized after relaxation (Fig. S8), and comparison of the architectures before and after relaxation revealed no structural collapse or topological distortion, indicating that the GaN metamaterials preserve their topology across the examined temperatures (Fig. S9). After reaching the equilibrium state, different approaches are applied to study their mechanical and piezoelectric properties.

To calculate the elastic constants, a small deformation strain (±0.1%) is applied. Then, the corresponding stress tensor is obtained after another relaxation process where the deformation strains are fixed.⁷² The elastic constant can be calculated from the applied strain and stress tensor. To obtain the mechanical deformation, GaN metamaterials are first fully relaxed at 300 K in NPT ensemble. Then, GaN metamaterials are subjected to a homogeneous tensile deformation by rescaling the coordinates of all atoms in the simulation box in the specific direction.⁷⁷ The stress along the directions that normal to the tension is controlled by using NPT ensemble to ensure the uniaxial loading condition.^78–80 The strain rate is set at 1 × 10⁹ s⁻¹. The axial tensile stress is calculated by averaging the virial stress of all atoms in the simulation system.

To study the piezoelectric properties, we need to first calculate the polarization of the system. Taking the piezoelectric response in the [0001] direction as an example, the polarization P₃ is defined as

P₃ = P₃(u(ε₃),ε₃) = P^elec₃(ε₃) + P^dis₃(u(ε₃)) (7)

where P^elec₃(ε₃) and P^dis₃(u(ε₃)) represent the polarization induced by separation between the core and the electrons, and the changes in the fractional atomic coordinate u (i.e., the displacement between Ga and N atoms after applying mechanical or electric fields), and ε₃ represents the strain field along the [0001] direction. In this study, we only consider the polarization induced by the fractional atomic coordinate changes (i.e., P^dis₃(u(ε₃))), because previous studies have demonstrated that the polarization between nucleus and electron cloud (i.e., P^elec₃(ε₃)) is negligible compared to the polarization induced by the displacement between atoms (i.e., P^dis₃(u(ε₃))).^34–36,81 Hence, the polarization along the [0001] direction is then given bywhere q_i is the electric charge of atom i; z_i is the coordinate of atom i in the [0001] direction; n is the total number of the atoms; V = [small rho, Greek, macron]V₀ ([small rho, Greek, macron] is the relative density of nano-architected metamaterials and V₀ = L_xL_yL_z is the volume of the architected cube at a certain mechanical strain) represents the volume of the solid parts of nano-architected metamaterials. In this study, we reduce the partial charges of these surface atoms to 75% of their bulk values, considering that the surface Ga and N atoms have 75% of the neighbors of the corresponding bulk atoms.^10,50,51

Without considering P^elec₃(ε₃), the piezoelectric stress constants (i.e., e₃₃, e₃₁, and e₃₂) can be obtained as⁴⁶

To calculate the piezoelectric stress constants (e₃₃, e₃₁, and e₃₂) of GaN metamaterials and bulk GaN, a small deformation strain (±1%) is applied by rescaling the coordinates of all atoms in the simulation box. It is worthwhile to distinguish the proper and improper piezoelectric stress constants here.^10,82 The improper piezoelectric stress tensor is calculated e_ijk = ∂P_i/∂ε_ijk, and the improper piezoelectric stress coefficients vary depending on the terminal surface in crystal lattice in calculation.⁴⁶ To calculate the invariant piezoelectric stress constants, three proper piezoelectric stress coefficients (i.e., ẽ₃₃, ẽ₃₁, and ẽ₃₂) are utilized in this study:^10,82

ẽ₃₃ = e₃₃ (10a)

ẽ₃₁ = e₃₁ + P₃ (10b)

ẽ₃₂ = e₃₂ + P₃ (10c)

The piezoelectric strain constants (d_ij) measure the strain induced by an external electric field, which can be described by⁸³

ε_i = S_ijσ_j + d_ikE_k (11a)

D_m = d_mjσ_j + κ_mkE_k (11b)

where S_ij is the elastic compliance matrix. σ_i and D_m are the stress and electric displacement, respectively. ε_i and E_k represent the strain and electric fields, respectively. κ_mk represents the dielectric constant. In this study, the piezoelectric strain constants are obtained by recording the strain changes with different external electric fields while the stress is fixed, givingwhere ε^E_j represent the strain induced by electric field (E). After reaching the equilibrium state, an external electric field with a strength from −0.08 V Å⁻¹ to 0.08 V Å⁻¹ is applied at a rate of 0.0016 V Å⁻¹ per picosecond. The electric field-induced strain is recorded to calculate the piezoelectric strain constants (i.e., d₃₃, d₃₁, and d₃₂) for GaN metamaterials and bulk GaN.

Conclusions

This study demonstrates the potential of topologically-engineered GaN spinodoid metamaterials for overcoming the intrinsic crystallographic limitations by achieving significantly enhanced piezoelectric properties and an increased number of independent non-zero piezoelectric stress/strain constants. Inspired by the spinodal decomposition, MD simulation results demonstrate that the GaN spinodoid metamaterials exhibit topology-dependent enhancements in piezoelectric characteristics (i.e., piezoelectric stress and strain constants) and anisotropy (e.g., d₃₁ ≠ d₃₂), distinguishing them from bulk GaN and other symmetric nano-architectures (e.g., face centred cubic, octet truss, and gyroid).¹⁰ These advancements are attributed to the role of asymmetric topologies and surface effects at the nanoscale, which are significantly influenced by the evolutionary times during the spinodal decomposition and relative density of spinodoid topologies. An important aspect of this work is the competing influence of free surfaces on the piezoelectric properties of GaN metamaterials. Under-coordinated surface atoms reduce local polarization,^50,51 while surface-induced contraction decreases the effective volume and enhances polarization.^10,84 The balance between these opposing mechanisms determines the overall piezoelectric response and provides a powerful additional degree of freedom for property tuning. Although our simulations assume pristine surfaces, we note that surface oxidation or hydroxylation, known to reconstruct the surface atoms, would also affect the piezoelectric behavior of GaN metamaterials.⁸⁵ Incorporating realistic surface chemistry therefore represents an important direction for future exploration. In addition to piezoelectric characteristics, the GaN spinodoid metamaterials reveal increased independent elastic constants compared to bulk GaN, with the number of independent constants rising from five to nine. While this study focuses on GaN, the design principles and methodologies presented here can be potentially extended to other piezoelectric materials such as ZnO, BaTiO₃, and PZT. Regarding the fabrication of the nanoscale piezoelectric metamaterials developed in this study, although experimentally challenging, recent advances in 3D nanofabrication techniques suggest feasible pathways toward realizing spinodoid piezoelectric architectures.^86,87 This progress motivates future experimental exploration of these topologically engineered, nano-architected piezoelectric metamaterials (see Section S5 in the SI). This study provides a pathway for developing next-generation intelligent materials and nano-devices through nanoscale topological engineering.

Author contributions

J. C. performed the computational studies and wrote the paper draft. J. C., A. S., B. S., and Z. Y. analyzed the data. J. C. and A.H. A. conceptualized the research. P. X. performed the computational studies. The research was conducted under the direction and supervision of A.H. A. All authors contributed to the review and editing of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5mh02127h.

Acknowledgements

This research was supported by the Canada Research Chairs program (Programmable Multifunctional Metamaterials), Natural Sciences and Engineering Research Council of Canada through NSERC Discovery Grant (RGPIN-2022-04493), Canada Foundation for Innovation (CFI) through John R. Evans Leaders Fund, and Quebec Research Fund Nature and Technologies (Team Grant 326563, and Doctoral Awards B2X). The authors gratefully acknowledge the Calcul Québec and Compute Canada for supporting the computational resources.

References

1. J. A. Paradiso and T. Starner, IEEE Pervasive Computing, 2005, 4, 18–27.
2. Z. L. Wang, MRS Bull., 2012, 37, 814–827.
3. J. H. He, C. L. Hsin, J. Liu, L. J. Chen and Z. L. Wang, Adv. Mater., 2007, 19, 781–784.
4. Z. L. Wang, MRS Bull., 2023, 48, 1014–1025.
5. J. Shi, K. Ju, H. Chen, V. Orsat, A. P. Sasmito, A. Ahmadi and A. Akbarzadeh, Adv. Funct. Mater., 2024, 35, 2417618.
6. W. Xiao, Z. Chen, X. Liu, Z. Zhou, Z. Fu, Y. Tang and R. Liang, Mater. Horiz., 2024, 11, 5285–5294.
7. A. Li, J. Yang, Y. He, J. Wen and X. Jiang, Nanoscale Horiz., 2024, 9, 365–383.
8. M. Gill, M. T. Mathe, P. Salamon, J. T. Gleeson and A. Jakli, Mater. Horiz., 2025, 12, 8920–8942.
9. J. Zhang, Q. Xu, Y. Zhang, W. Guo, H. Zeng, Y. He, J. Wu, L. Guo, K. Zhou and D. Zhang, Mater. Horiz., 2025, 12, 3494–3504.
10. J. Cai, L. Yan, A. Seyedkanani, V. Orsat and A. Akbarzadeh, Nano Energy, 2024, 129, 109990.
11. H. Cui, R. Hensleigh, D. Yao, D. Maurya, P. Kumar, M. G. Kang, S. Priya and X. R. Zheng, Nat. Mater., 2019, 18, 234–241.
12. K. Nadaud, G. F. Nataf, N. Jaber, B. Negulescu, F. Giovannelli, P. Andreazza, P. Birnal and J. Wolfman, ACS Appl. Electron. Mater., 2024, 6, 7392–7401.
13. X. Liu, M. Jiang, Y. Zeng, Y. Ouyang, Y. Xu, S. Cao, J. Song, L. Li, S. Cheng and G. Rao, J. Mater. Sci.: Mater. Electron., 2024, 35, 780.
14. Y. Zeng, M. Jiang, X. Liu, T. Wang, Y. OuYang, Y. Xu, S. Cao, J. Song and G. Rao, J. Mater. Sci.: Mater. Electron., 2024, 35, 1510.
15. D. T. Harris, M. J. Burch, E. J. Mily, E. C. Dickey and J.-P. Maria, J. Mater. Res., 2016, 31, 1018–1026.
16. Y. Yan, L. D. Geng, H. Liu, H. Leng, X. Li, Y. U. Wang and S. Priya, Nat. Commun., 2022, 13, 3565.
17. C. J. Ramos-Cano, M. Miki-Yoshida, R. Herrera-Basurto, F. Mercader-Trejo, L. Fuentes-Cobas, O. Auciello and A. Hurtado-Macías, Appl. Phys. A, 2023, 129, 113.
18. A. B. Haugen, M. I. Morozov, M. Johnsson, T. Grande and M.-A. Einarsrud, J. Appl. Phys., 2014, 116, 134102.
19. J. Schultheiß, O. Clemens, S. Zhukov, H. von Seggern, W. Sakamoto and J. Koruza, J. Am. Ceram. Soc., 2017, 100, 2098–2107.
20. H. Liu, Y. Yan, H. Leng, A. Heitmann, J. B. Blottman and S. Priya, Appl. Phys. Lett., 2020, 116, 252901.
21. D. Yao, H. Cui, R. Hensleigh, P. Smith, S. Alford, D. Bernero, S. Bush, K. Mann, H. F. Wu, M. Chin-Nieh, G. Youmans and X. Zheng, Adv. Funct. Mater., 2019, 29, 1903866.
22. H. Lu, H. Cui, G. Lu, L. Jiang, R. Hensleigh, Y. Zeng, A. Rayes, M. K. Panduranga, M. Acharya, Z. Wang, A. Irimia, F. Wu, G. P. Carman, J. M. Morales, S. Putterman, L. W. Martin, Q. Zhou and X. R. Zheng, Nat. Commun., 2023, 14, 2418.
23. H. Cui, D. Yao, R. Hensleigh, H. Lu, A. Calderon, Z. Xu, S. Davaria, Z. Wang, P. Mercier, P. Tarazaga and X. R. Zheng, Science, 2022, 376, 1287–1293.
24. J. Shi, K. Ju, H. Chen, A. Mirabolghasemi, S. Akhtar, A. Sasmito and A. Akbarzadeh, Nano Energy, 2024, 123, 109385.
25. J. He, Y. Wang, Z. Shen, L. Xia and Y. Xiong, Mater. Horiz., 2024, 11, 6371–6380.
26. J. Cai, Y. Chen, A. Seyedkanani, G. Shen, M. Cerruti and A. Akbarzadeh, Adv. Sci., 2025, e14597,  DOI:10.1002/advs.202514597.
27. A. Seyedkanani and A. Akbarzadeh, Adv. Funct. Mater., 2022, 32, 2207581.
28. B. Shahryari, H. Mofatteh, A. Sargazi, A. Mirabolghasemi, D. Meger and A. Akbarzadeh, Adv. Funct. Mater., 2024, 35, 2407651.
29. D. W. Yee, M. L. Lifson, B. W. Edwards and J. R. Greer, Adv. Mater., 2019, 31, e1901345.
30. C. M. Portela, B. W. Edwards, D. Veysset, Y. Sun, K. A. Nelson, D. M. Kochmann and J. R. Greer, Nat. Mater., 2021, 20, 1491–1497.
31. C. M. Portela, A. Vidyasagar, S. Krodel, T. Weissenbach, D. W. Yee, J. R. Greer and D. M. Kochmann, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 5686–5693.
32. H. Jiang, Y. Su, J. Zhu, H. Lu and X. Meng, Nano Energy, 2018, 45, 359–367.
33. R. J. Wang, C. Y. Wang, Y. T. Feng and C. Tang, Nano Energy, 2018, 53, 906–915.
34. J. Zhang, Nano Energy, 2021, 79, 105489.
35. J. Zhang, Nano Energy, 2021, 86, 106125.
36. J. Zhang and S. A. Meguid, Nano Energy, 2015, 12, 322–330.
37. Y. Zhang, J. Cai, Q. Cai, L. Wang and X. Gou, Compos. Struct., 2025, 119036,  DOI:10.1016/j.compstruct.2025.119036.
38. Y. Zhang, J. Cai, C. Mi and A. Akbarzadeh, Adv. Theory Simul., 2022, 5, 2200339.
39. J. Cai, B. Shahryari, A. Seyedkanani, A. P. Sasmito and A. Akbarzadeh, Nano Lett., 2025, 25, 7603–7610.
40. Y. Zhang, J. Xie, X. Gou, Q. Cai and J. Cai, Compos. Sci. Technol., 2025, 271.
41. Y. Zhang, J. Cai, C. Mi, F. Wang and A. H. Akbarzadeh, Acta Mech., 2022, 233, 233–257.
42. S. Kumar, S. Tan, L. Zheng and D. M. Kochmann, npj Comput. Mater., 2020, 6, 73.
43. W. Deng, S. Kumar, A. Vallone, D. M. Kochmann and J. R. Greer, Adv. Mater., 2024, 36, e2308149.
44. X. Y. Sun, G. K. Xu, X. Y. Li, X. Q. Feng and H. J. Gao, J. Appl. Phys., 2013, 113, 023505.
45. J. W. Cahn and J. E. Hilliard, J. Chem. Phys., 1958, 28, 258–267.
46. S. Dai, M. L. Dunn and H. S. Park, Nanotechnology, 2010, 21, 445707.
47. J. Zhang, Nanotechnology, 2020, 31, 095407.
48. F. Bernardini, V. Fiorentini and D. Vanderbilt, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 56, R10024–R10027.
49. J. Cai, A. Seyedkanani, B. Shahryari, H.-C. Lin and A. Akbarzadeh, Int. J. Heat Mass Transfer, 2026, 256, 127911.
50. S. Dai and H. S. Park, J. Mech. Phys. Solids, 2013, 61, 385–397.
51. D. Tan, Y. Xiang, Y. Leng and Y. Leng, Nano Energy, 2018, 50, 291–297.
52. A. Polian, M. Grimsditch and I. Grzegory, J. Appl. Phys., 1996, 79, 3343–3344.
53. K. Kim, W. R. Lambrecht and B. Segall, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 1502–1505.
54. A. F. Wright, J. Appl. Phys., 1997, 82, 2833–2839.
55. J. Zhang, Appl. Phys. Lett., 2014, 104, 253110.
56. I. L. Guy, S. Muensit and E. M. Goldys, Appl. Phys. Lett., 1999, 75, 4133–4135.
57. S. Muensit, E. M. Goldys and I. L. Guy, Appl. Phys. Lett., 1999, 75, 3965–3967.
58. I. Vurgaftman and J. R. Meyer, J. Appl. Phys., 2003, 94, 3675–3696.
59. J. Zhang, C. Wang, R. Chowdhury and S. Adhikari, Scr. Mater., 2013, 68, 627–630.
60. J. Bauer, C. Crook, A. Guell Izard, Z. C. Eckel, N. Ruvalcaba, T. A. Schaedler and L. Valdevit, Matter, 2019, 1, 1547–1556.
61. J. Bauer, A. Schroer, R. Schwaiger and O. Kraft, Nat. Mater., 2016, 15, 438–443.
62. J. Biener, A. M. Hodge, A. V. Hamza, L. M. Hsiung and J. H. Satcher, J. Appl. Phys., 2005, 97, 024301.
63. L. Dong, V. Deshpande and H. Wadley, Int. J. Solids Struct., 2015, 60–61, 107–124.
64. A. J. Jacobsen, S. Mahoney, W. B. Carter and S. Nutt, Carbon, 2011, 49, 1025–1032.
65. S. N. Khaderi, M. R. J. Scherer, C. E. Hall, U. Steiner, U. Ramamurty, N. A. Fleck and V. S. Deshpande, Extreme Mech. Lett., 2017, 10, 15–23.
66. L. R. Meza, S. Das and J. R. Greer, Science, 2014, 345, 1322–1326.
67. X. Wendy Gu and J. R. Greer, Extreme Mech. Lett., 2015, 2, 7–14.
68. X. Zhang, A. Vyatskikh, H. Gao, J. R. Greer and X. Li, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 6665–6672.
69. X. Zheng, H. Lee, T. H. Weisgraber, M. Shusteff, J. DeOtte, E. B. Duoss, J. D. Kuntz, M. M. Biener, Q. Ge, J. A. Jackson, S. O. Kucheyev, N. X. Fang and C. M. Spadaccini, Science, 2014, 344, 1373–1377.
70. C. Crook, J. Bauer, A. Guell Izard, C. Santos de Oliveira, E. S. J. Martins de Souza, J. B. Berger and L. Valdevit, Nat. Commun., 2020, 11, 1579.
71. Y. Wang, X. Zhang, Z. Li, H. Gao and X. Li, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2119536119.
72. C. T. Chen, D. C. Chrzan and G. X. Gu, Nat. Commun., 2020, 11, 3745.
73. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
74. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.
75. F. H. Stillinger and T. A. Weber, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 5262–5271.
76. A. Béré and A. Serra, Philos. Mag., 2006, 86, 2159–2192.
77. J. Cai, C. Mi, Q. Deng and C. Zheng, Mech. Mater., 2019, 139, 103205.
78. J. Cai, B. Yang and A. Akbarzadeh, ACS Nano, 2024, 18, 894–908.
79. J. Cai, E. Estakhrianhaghighi and A. Akbarzadeh, Carbon, 2022, 191, 610–624.
80. J. Cai and A. Akbarzadeh, Mater. Des., 2021, 206, 109811.
81. T. Dan, M. Willatzen and Z. L. Wang, Nano Energy, 2019, 56, 512–515.
82. D. Vanderbilt, J. Phys. Chem. Solids, 2000, 61, 147–151.
83. I. Chopra, AIAA J., 2002, 40, 2145–2187.
84. R. Agrawal and H. D. Espinosa, Nano Lett., 2011, 11, 786–790.
85. J. Zhang and J. Zhou, Nano Energy, 2018, 50, 298–307.
86. I. Bita, J. K. Yang, Y. S. Jung, C. A. Ross, E. L. Thomas and K. K. Berggren, Science, 2008, 321, 939–943.
87. H. S. Suh, D. H. Kim, P. Moni, S. Xiong, L. E. Ocola, N. J. Zaluzec, K. K. Gleason and P. F. Nealey, Nat. Nanotechnol., 2017, 12, 575–581.

---