Machine learning enables the discovery of 2D Invar and anti-Invar monolayers

Thermal expansion in solids arises from changes in lattice vibrations with temperature and is quantified by the linear coefficient of thermal expansion (α). Most materials show positive thermal expansion (PTE), while some exhibit negative thermal expansion (NTE), and very few show near-zero thermal expansion (ZTE, Invar effect). Since Guillaume’s 1897 discovery of Invar alloys, various Invar materials (ferromagnetic alloys, ferroelectric relaxors, silica-based glasses, and nanoporous crystals) have been identified and used in precision devices. With electronic devices scaling to nanometer and Ångstrom dimensions, 2D monolayers are key building blocks, and controlling their thermal expansion is critical to device performance and reliability. Reported 2D crystals exhibit either PTE (e.g., MoS₂) or NTE (e.g., graphene), but no 2D monolayer with ZTE had been discovered. The emergence of large 2D materials databases (C2DB, 2DMatPedia, JARVIS) enables data-driven discovery. Given that both PTE and NTE exist in 2D crystals, the authors hypothesize that ZTE counterparts should also exist and can be identified via machine learning using physically interpretable descriptors.

The study situates itself within prior work on Invar and thermal expansion: Invar behavior has been extensively studied in alloys and complex oxides and glasses, with applications in precision instruments. Traditional explanations involve compensation between lattice PTE and other mechanisms (e.g., magnetic or electronic transitions). For 2D materials, prior theory and experiments have documented PTE (MoS₂) and NTE (graphene, h-BN) behaviors and linked 2D NTE to out-of-plane flexural (ZA) phonons with negative Grüneisen parameters. Active learning and ML have been used recently to discover bulk Invar alloys. Databases like C2DB, 2DMatPedia, and JARVIS enable high-throughput computational screening and ML-guided property prediction. Prior works also address mechanics of 2D membranes, thermal fluctuations, finite-temperature elasticity, and relationships between bonding characteristics and thermal expansion of bonds, providing the foundation for the descriptors used here.

Overview: The authors use a combined data-driven and first-principles workflow to identify interpretable mechanical descriptors for 2D thermal expansion and to discover ZTE (Invar) and extreme thermal expansion (ELTE: EPTE and ENTE) monolayers. Feature selection and classification: Starting with physically motivated features (thickness h, volume per atom, average bond length, average atomic radius, bond density, in-plane stiffness E2D, out-of-plane bending stiffness D), a random forest ranked E2D and D as the most important for classifying thermal expansion behavior. Using the SISSO method, the authors trained a classification model to separate PTE and NTE 2D crystals based on E2D and D, yielding two descriptors S1 = 1/D and S2 = 1/E2D. A linear SVM separation line cleanly classified all training data into PTE vs NTE groups. Descriptor interpretation: Physical arguments relate S1 to in-plane contraction induced by out-of-plane thermal fluctuations (membrane statistical mechanics indicates α scales approximately as 1/D) and S2 to in-plane bond-driven thermal expansion, inversely related to E2D via empirical relations between bond valence, force constants, and thermal expansion of bonds. Symbolic regression for property prediction: Because D is absent in C2DB, D was predicted via symbolic regression using inputs {E2D, density, volume per atom, thickness h, average bond length, bond density}, achieving RMSE ≈ 1.34 eV. The LTEC at 500 K (α500K) was then modeled by symbolic regression using quantities including E2D, D, and isothermal compressibility KT, giving an expression α500K = (14.31 × KT^-0.038) × E2D + 14.27 × δ (units and constants as in Supplementary). Predicted α500K correlates well with QHA results with RMSE 1.67 × 10⁻⁶ K⁻¹ (training) and 1.35 × 10⁻⁶ K⁻¹ (testing), R² of 0.91 and 0.93, respectively. High-throughput screening workflow: From C2DB, 1224 stable 2D crystals were selected based on elastic stability (C > 0), dynamical stability (phonon frequencies positive at high-symmetry points; later full Brillouin zone for fine screening), and thermodynamic stability (ΔHf⁰ < 0, ΔH < 0.2 eV). ZTE candidates were initially screened using the regression model with |α500K| < 4 × 10⁻⁶ K⁻¹, then validated via Grüneisen theory (GT), dynamical stability in full q-space, and quasi-harmonic approximation (QHA). EPTE candidates were selected by D > 8 eV and E2D < 80 N/m; ENTE candidates by D < 1 eV; both were evaluated by GT and confirmed via QHA with thresholds |α| > 15 × 10⁻⁶ K⁻¹ (at relevant temperatures). First-principles calculations: DFT calculations used VASP with PAW and PBE functionals, plane-wave cutoff 600 eV, k-point mesh 17 × 17 × 1, and ~20 Å vacuum. Phonons were computed by finite differences with PHONOPY using 6 × 6 × 1 supercells and 4 × 4 × 1 k-point sampling. For QHA, small strains within ±1.5% were applied; equilibrium lattice constants at temperature T were obtained by minimizing Helmholtz free energy, and LTEC defined as α(T) = da/a. Grüneisen-theory-based LTECs were also computed for efficiency, with validation against QHA. Effects of vdW corrections (DFT-D3, vdW-DF) on LTECs were tested and found negligible. Validation against literature data indicated good agreement. Mechanical property extraction: Bending stiffness D was obtained by fitting the flexural ZA dispersion near Γ; in-plane elastic constants Cij were computed via small ±1% strains and used to derive E2D = (C11 C22 − C12 C21)/C22. Areal compressibility Ks was calculated as S11 + S12 from the elastic compliance matrix S = C⁻¹.

- Two interpretable mechanical descriptors, in-plane stiffness E2D and out-of-plane bending stiffness D, effectively classify PTE vs NTE 2D crystals. SISSO yields S1 = 1/D (capturing out-of-plane fluctuation-induced contraction) and S2 = 1/E2D (capturing bond-driven in-plane expansion). A linear SVM using these descriptors perfectly separates PTE and NTE in the training data. - Trends within materials families are consistent with descriptor physics: for 2H-MoX₂ (X = S, Se, Te) at 300–800 K, LTEC follows MoS₂ < MoSe₂ < MoTe₂; correspondingly, E2D decreases and D increases along S→Te, explaining larger PTE for MoTe₂. - Phonon mode analysis: The ZA-mode contribution to LTEC (α_ZA) is negatively correlated with D; fit yields α_ZA = β1 D with β1 ≈ −9.2 × 10⁻⁶ K⁻¹/eV at 500 K. The non-ZA contribution (α500K − α_ZA) anti-correlates with E2D; fit yields α500K − α_ZA = β2 E2D with β2 ≈ −117.49 × 10⁻⁶ K⁻¹/(N/m), supporting S1 and S2 interpretations. - Symbolic regression model for α500K achieves RMSE of 1.67 × 10⁻⁶ K⁻¹ (training) and 1.35 × 10⁻⁶ K⁻¹ (testing), R² = 0.91/0.93, sufficient for rough ZTE screening. - Discovery of ZTE (2D Invar) monolayers: ZrO₂ and HfO₂ exhibit near-zero LTECs over 300–600 K. At 500 K, α = 0.47 × 10⁻⁶ K⁻¹ (ZrO₂) and 1.80 × 10⁻⁶ K⁻¹ (HfO₂). Minimum LTECs are −6.90 × 10⁻⁶ K⁻¹ (ZrO₂) and −3.16 × 10⁻⁶ K⁻¹ (HfO₂), with NTE at low T transitioning to PTE near ~400 K (ZrO₂) and ~170 K (HfO₂). These TMOs have relatively large E2D (152–157 N/m) and moderate D (3.67–4.48 eV), and very negative ZA Grüneisen parameters near Γ (≈ −54 for ZrO₂; ≈ −81 for HfO₂), much lower than MoS₂ (≈ −10). - Discovery of EPTE (extreme PTE, anti-Invar) monolayers: SbSe₂, HfBi₂, HfSb₂. At 500 K, α = 24.49, 16.99, and 19.81 × 10⁻⁶ K⁻¹, respectively, all >15 × 10⁻⁶ K⁻¹. They share large thicknesses and bending stiffnesses (D ≈ 9.03–13.74 eV) and relatively low E2D (31.88–70.99 N/m), consistent with EPTE placement far below the SVM line (high D, low E2D). - Discovery of ENTE (extreme NTE, anti-Invar) monolayers at low temperatures: As₂ (α_min = −36.39 × 10⁻⁶ K⁻¹ at 40 K), Ge₂ (−30.46 × 10⁻⁶ K⁻¹ at 90 K), Sn₂ (−61.27 × 10⁻⁶ K⁻¹ at 80 K). These have small D (≈ 0.74–0.85 eV) with moderate-to-large E2D, consistent with placement above the SVM line (low D, high E2D). ENTE manifests at low T due to early activation of low-frequency ZA modes with negative Grüneisen parameters. - Additional observations: Graphene and h-BN exhibit NTE below 1000 K but not extreme enough for ENTE by the adopted threshold. Some NTE materials (e.g., Si₂, BP) can show ZTE at high temperatures, but are categorized as NTE due to strongly negative α at low T.

Findings support a unifying picture that thermal expansion in 2D monolayers emerges from competition between in-plane bond expansion (favoring PTE and scaling with 1/E2D) and out-of-plane fluctuation-induced in-plane contraction (favoring NTE and scaling with 1/D). The interpretable descriptors E2D and D provide both predictive and explanatory power, enabling classification and targeted discovery of ZTE and ELTE materials. The results highlight design strategies: pairing PTE and NTE layers in heterostructures can yield ZTE; increasing layer number n tends to shift NTE monolayers toward ZTE because E2D ∝ n while D scales faster (∝ n² or with exponent >1), as exemplified by graphite exhibiting smaller |α| than graphene. The work underscores the technological relevance of ZTE and ELTE 2D materials for nanoelectronics and functional architectures (e.g., exploiting mismatch strains for self-curving 3D nanostructures). It also emphasizes that though free-standing, infinite-sheet models capture key trends, boundary conditions, substrates, finite sizes, and entropy-driven waviness can modify phonon populations and effective E2D and D, altering thermal expansion responses.

The study identifies two physically interpretable mechanical descriptors—E2D (in-plane tensile stiffness) and D (out-of-plane bending stiffness)—that accurately classify PTE vs NTE in 2D crystals and illuminate the mechanisms governing thermal expansion. Using high-throughput DFT, SISSO-based classification, and symbolic regression, the authors predict and validate 2D monolayers with near-zero thermal expansion (ZrO₂, HfO₂) and extreme thermal expansion (EPTE: SbSe₂, HfBi₂, HfSb₂; ENTE: As₂, Ge₂, Sn₂). The approach demonstrates robust predictive accuracy for α500K and offers generalizable design rules linking E2D and D to thermal expansion behavior. Future directions include extending searches to other databases (2DMatPedia, JARVIS), developing theoretical frameworks that incorporate finite-size and boundary effects, computing temperature-dependent E2D and D, and engineering heterostructures or multilayers to achieve targeted ZTE or ELTE behavior.

- The quasi-harmonic approximation (QHA) models free-standing, infinite-size membranes and does not capture finite-size effects, boundary conditions, substrates, or entropy-driven wavy conformations, which can alter phonon populations and thermal expansion. - Descriptor values E2D and D are computed at 0 K for flat configurations; their temperature and configuration dependence is not included, which may impact predictions for real devices. - The symbolic regression model for α500K is trained within a limited domain and is not used to extrapolate to ELTE regimes; EPTE/ENTE identification relied on descriptor-based screening plus GT/QHA. - The prediction of D via symbolic regression carries an RMSE (~1.34 eV); improved accuracy would benefit from larger training datasets. - Initial stability screening uses phonon checks at high-symmetry points before full-q validation; some marginal cases may require more exhaustive stability analyses. - ENTE behavior identified occurs at low temperatures; applicability at operating temperatures may be limited for some candidates.