Universal framework for efficient estimation of stability in multi-principal element alloys

Domain: chemistry
Published: 2026-02-23
Editorial outcome: major_revision

Abstract
Predicting the synthetic accessibility of multi-principal element alloys (MPEAs) across the global chemical space remains a challenge. In this study, we show that the synthesizability of MPEAs across broad compositional and structural spaces can be predicted using a physical model that expresses the total energy of any MPEA as a linear combination of energies from lower-dimensional subsystems. The model is validated with a large computational dataset and supported by the experimental synthesis of multiple MPEAs, achieving mean absolute errors near or below 7 meV/atom on a density functional theory dataset of 135,791 MPEAs spanning 28 metals and up to ten components. Its accuracy is comparable to state-of-the-art deep learning models while maintaining interpretability through cluster-expansion theory. Moreover, we show that the stability of high-entropy alloys can be predicted using a linear combination of energies from lower-dimensional systems with low errors, indicating a flatter energy landscape at high compositional complexity.

Introduction
Multi-principal element alloys (MPEAs)1–12, including medium- and high-entropy alloys, have emerged as promising materials in various fields due to their extensive compositional space and potential for diverse applications9–19. However, predicting the synthesizability of MPEAs remains a challenge, which has slowed down the synthesis explorations. Earlier work examined high entropy alloy formation using binary formation energies as an initial screening approach20. Subsequent studies have explored a variety of approaches using machine learning (ML) models and heuristic theories3,21–26. However, these models often struggle with limited accuracy and/or interpretability27,28. Despite notable progress in predicting synthesizability of MPEA3,21–26 and high entropy materials29,30, interpretable physical models for predicting phase stability across global chemical space of MPEAs remain lacking. Three key questions persist: (1) Is there a universal rule governing MPEA stability? (2) What is the crucial data to reveal MPEA stability across different compositions and dimensionality? and (3) What factors limit the accuracy of phase stability predictions for MPEAs?

In this work, we address the above questions by combining theoretical and experimental approaches. We present a physical model that can serve as universal model for predicting the synthesizability of MPEAs with any composition and dimensionality. This model is validated on a dataset of 135,791 density functional theory (DFT)-computed MPEA structures covering 28 typical metals, showing a mean absolute error (MAE) as low as 6.956 meV/atom for predicting synthesizability of MPEAs. Such performance matches the state-of-the-art graph neural network models31 trained on the same dataset32. The reliability of our computational dataset and stability metric can be supported by successful synthesis of nine previous unreported MPEAs and accurate prediction of 30 out of 32 experimentally reported MPEAs synthesized below 1273 K. It also supports the stability rule for MPEAs: the energy and hence stability of a high-dimensional MPEA can be reasonably approximated by a weighted combination of the energies of its lower-dimensional constituent subsystems. More importantly, we uncover two unprecedented phenomena: (1) accurate determination of the synthesizability of any MPEA compositional space can be achieved with small datasets; Specifically, when focusing on a subset of MPEA including N metals, as few as C3 N DFT calculations are sufficient to predict the energy of all corresponding MPEAs; and (2) there is a decaying of decomposition energies of higher dimensional MPEA to lower dimensional MPEA with the trend being more pronounced when dimensionality goes up.

Results
Formation energy and energy above hull
We create a comprehensive database to investigate the stability of MPEAs with increasing compositional complexity. This database includes MPEAs with two, three, four, and five equimolar compositions (denoted as 2M, 3M, 4M, and 5M) from 28 metals (8 noble metals and 20 non-noble metals, Fig. 1a). The 2M, 3M, and 4M MPEAs were generated by enumerating across this 28-metal space, yielding 24,129 (=C2 28 +C3 28 +C4 28 ) MPEAs with distinct compositions. We also enumerate all combinatorial space among quaternary equimolar MPEAs that can be entropy-stabilized below 1273 K with the leftover 24 metals, which yielded another 18,014 equimolar MPEAs compositions with five metal elements. For each composition, three different crystal structures, FCC, HCP, and BCC, are initialized for DFT calculations, resulting in a total of 126,429 DFT-optimized MPEA structures.

The stability of MPEAs is commonly characterized by two quantities: formation energy (E form) and energy above hull (E hull). The value E form is defined as the energy released when forming an MPEA from pure elemental phases and E hull is a generalized measure of stability that assesses the largest driving force for an MPEA to decompose into any potential ground state (metal or intermetallic) in its high-dimensional phase diagram33–35. We establish the E hull values by constructing phase diagrams that include all experimentally reported metals and intermetallic phases, as well as all hypothetical low-energy structures hosted in both Materials Project36 and the Open Quantum Materials Database (OQMD)37,38. This approach results in a competing phase pool containing 6723 metals and intermetallic phases. Particularly, 53, 754, and 11,996 quaternary MPEAs are capable of being stabilized below 500 K, 1000 K, and 2000 K, respectively (Fig. 1b). Both E form and E hull present in Fig. 1b are calculated at 0 K, with the black line (E form = E hull) indicating cases where decomposition leads to elemental phases. Structural-dependent plots are shown in Supplementary Fig. 1; the energy distributions across FCC, BCC, and HCP structures are similar, and a summary of the number of 4-metal MPEAs stabilized under different temperatures and phases is provided in Supplementary Fig. 2. To examine the correlation between the two stability descriptors and commonly used physical parameters, we adopt the same plotting style as in Fig. 1b, coloring the data points by valence electron concentration (VEC), average electronegativity (χ ave), and average off-lattice displacement (ε ave) (see Supplementary Fig. 3). These analyses reveal no clear or systematic trends, although compositions with large off-lattice displacements generally correspond to instability. These observations further highlight the complexity of stability prediction in high-entropy systems, emphasizing the need for models that are both predictive and physically interpretable.

Experimental validation of predicted MPEAs
Among all MPEA compositions, noble-metal-containing MPEAs emerge as a distinct and underexplored group in terms of compositional coverage. As shown in Fig. 2a, this compositional space remains relatively underexplored, with only about one hundred reported compositions published39–46 from our survey (Supplementary Fig. 4 and Supplementary Table 1), in contrast to over several thousand papers and compositions reported for other MPEAs4,47. Furthermore, the ease of reducing noble-metal precursors in a hydrogen atmosphere facilitates their synthesis using techniques such as scanning probe block copolymer lithography (SPBCL). Notably, noble-metal-containing MPEAs have also garnered significant attention due to their great potential in catalysis and other applications1,16,17,39,40,48. Therefore, our experimental validation efforts will primarily focus on noble-metal-containing MPEAs. In this work, using high-throughput phase diagram calculations, we find 647 ternary, 2178 quaternary, and 5478 quinary noble-metal-containing MPEAs that can be stabilized below 1273 K with ideal mixing entropy (Fig. 2b).

To verify the reliability of using E hull values as a stability descriptor, we select six quaternary and three quinary MPEAs with stable temperatures below 1273 K for synthesis (e.g., E hull < S ideal × 1273 K). The compounds are: FeNiCuPd, MnFeNiPt, NiCuPtAu, FeNiPdAu, FeNiCuPt, FeNiPdPt, CoNiPdPtRh, RuIrPtPdRh, and FeCoNiRhPd (Supplementary Table 2, and E hull, E form value shown in Supplementary Fig. 5). An example of how synthesizability is evaluated can be found in Supplementary Note 1. These nine MPEAs are carefully selected based on synthesis feasibility and theoretical predictions. We focus on noble-metal-rich compositions due to their enhanced reducibility under H2, excluding early transition metals (e.g., Ti, V, Cr), other active metals such as Mg and Al, and toxic or volatile elements (e.g., Os, Zn, Sn) for safety and experimental limitations. Selection prioritized alloys with low predicted stabilization temperatures, high noble metal content, and compositions that are never reported before (with reported compounds listed in Supplementary Table 1). It should be noted that even though those materials are synthesized as nanoparticles (NPs), the size is controlled to be rather large (>50 nm) on purpose, which goes beyond the region of typical size-dependent phase stability, which usually limit to less than and close to 10 nm49,50. We thus consider the experiment to be comparable to our DFT calculations.

The MEPANPs were prepared by SPBCL, a technique enabled by Dip Pen Nanolithography51 and Polymer Pen Lithography52. (see Methods “Multimetallic nanoparticle (NP) synthesis”)53–57. In this process, various metal precursors with equimolar ratio were dissolved in a block copolymer solution, which was then delivered by an atomic force microscope tip to form polymer dome arrays on the substrate. After thermal annealing treatments, the metal precursors were reduced to form one NP in each dome (Supplementary Fig. 6). In SPBCL, the amount of metal precursors in the polymer inks can be adjusted, allowing precise control of the chemical compositions of single NPs. More importantly, the SPBCL-prepared NPs can be heated at high temperatures (up to 1000 °C) without particle aggregation, which is necessary to evaluate the thermal stabilities of MPEAs. The scanning transmission electron microscope (STEM) images and energy dispersive X-ray (EDX) elemental mappings of 4M (Fig. 2c) and 5M (Fig. 2d) MPEAs confirmed the successful synthesis of all predicted stable MPEAs with uniform element distribution and equal element concentrations (Supplementary Fig. 7 and Supplementary Table 3). The experimental results validated that E hull can be a useful metric to assess MPEA stability.

Elemental stability rules
In our established dataset of 126,429 DFT-computed MPEAs up to 5 metals, 8,303 MPEAs (6.56% of all compounds) with three to five elements are predicted to have stabilization temperatures below 1273 K. We consider these MPEAs as synthetically accessible MPEAs. A chord diagram (Fig. 3a), which charts the overall elemental compatibility between metals, was used to categorize the 28 metals considered in this study into noble metals (colored) and non-noble metals (gray). The thickness of each segment in the diagram is proportional to the occurrence (n M stable) of each metal (M) appearing in a stable MPEA that can be entropy stabilized below 1273 K. The connections between segments indicate the occurrence (n M1 ∩ M2 stable) when both metals (M1 and M2) coexist in an MPEA that can be entropy stabilized below 1273 K. The connections between non-noble metals are shown in Fig. 3c. Given the substantial enhancement in catalytic performance brought by the inclusion of noble metals, their impact is specifically emphasized in Fig. 3b, d–k.

The elemental stability mapping in Fig. 3 provides insights into the compatibilities and preferences of metal species within MPEAs. Specifically, we find that Rh and Pt stand out as the most compatible noble metals (Fig. 3a, b, d, e), while Ni and Fe are the most compatible non-noble metals (Fig. 3a, c). In contrast, Ag and Au are the most selective noble metals when forming MPEAs (Fig. 3a, b, j, k). Sc, Hf, Zr, and Ta are found to form stable MPEA with a limited range of other metals (Sc: 203, Hf: 370, Zr: 313, Ta: 553). In terms of elemental pairing preferences, Rh, Pt, and Ru generally show broad compatibility with most elements, as highlighted by the extensive pairing data (Rh: 3175, Pt: 2777, Ru: 2391). In comparison, stable MPEAs containing Ag are selectively formed, primarily with Au, Pt, Pd, and Zn.

The observed trends align with general periodic table groupings and proximities; metals from the same group (e.g., Pd-Pt, Ru-Os, Ag-Au) or adjacent positions in a row (e.g., Ir-Pt, Ru-Rh, Rh-Pd) often exhibit higher pairing probabilities. These trends are supported by their close atomic properties, such as electronic structure and atomic radius, in accordance with the Hume-Rothery rules58. The latter group FCC metals, including noble metals such as Rh, Pt, Pd, and Ir, as well as non-noble metals such as Ni and Cu, present enhanced compatibility in MPEAs across different structures. This compatibility extends to Fe, Co, Mn, and Cr, aligning with the discovery of the Cantor alloy (CrMnFeCoNi)59 and its various single-phase modifications1,60–62.

The elemental stability rules we have uncovered provide insights into previously reported experimental efforts on MPEAs (detailed in Supplementary Fig. 4 and Supplementary Table 1). The MPEAs containing Ru, Pd, Ag, and Pt have been frequently investigated, aligning with our findings that noble-metal-containing MPEAs tend to be stable (6928 of the 8303 stable MPEAs contain noble metals). We can also link the high frequency of specific elemental pairs, such as Pd-Pt, Rh-Pt, and Rh-Pd, in the literature to our findings (Fig. 3b). Prior research on the high probability of Ag-Zn in binary alloys supports our prediction63,64. In this work, the successful synthesis of nine MPEAs (Fe, Ni, Cu, Pd, Mn, Pt, Au, Rh, Ru, Ir, Co) further verifies the reliability of our computational dataset and elemental stability rules.

Universal framework across extended dimensionalities and compositions
Beyond elemental compatibility, the stability of MPEAs is also influenced by their dimensionality—the number of distinct elements in a single-phase alloy. With growing interest in the synthesis and application of high-dimensional MPEAs, the challenge lies in the ability to predict their stability accurately and efficiently at extended dimensionality. Here, we demonstrate that a simple linear model can be used to effectively estimate the stability of MPEAs across various dimensionalities.

The core hypothesis of our model is that the total energy of an MPEA can be expressed as a linear summation of contributions from various atomic clusters. These clusters refer to local atomic arrangements present in the MPEA structure, defined in the same manner as in the typical cluster expansion theory65. The construction of total energy from energies of different local clusters has been widely used in cluster expansion theory34,65–71, from simple alloys to high-entropy solid solutions. This concept can be expressed as:
E = Σ i J i n i + Σ i,j J ij n ij + Σ i,j,k J ijk n ijk + ... (1)
where J i, J ij, and J ijk represent the energy contributions from singlets (single atom), pairs, and triplets, respectively, while n i, n ij, and n ijk refer to the concentrations or counts of these clusters. This formalism has inspired the development of our model, which proposes that the total energy of higher-dimensional MPEAs can be approximated by summing the energies of lower-dimensional MPEAs, thus closely matching the averaged cluster probability of the target MPEA. Specifically, the model can be expressed as:
E N ≈ Σ i∈(1,CM N) εi E M (2)
where E N and E M represent the total energy of an MPEA with N and M elements, respectively, and εi is the weighting factor in the summation. To represent an N-element MPEA as a combination of its M-component subspaces, a convenient and physically meaningful choice is εi = 1 / CM N (see Supplementary Note 2), where CM N is the binomial coefficient denoting the number of possible M-component combinations from an N-component system. The uniform choice of εi ensures that the energy contributed by the singlet term cancels exactly when averaging over all subspaces, thereby reducing the overall cluster probability difference. It also helps conserve the elemental concentration across the ensemble of subspaces, as our model focuses on equimolar compositions, which lie at the center of the composition space and maximize configurational entropy. Numerically, we can use special quasi-random (SQS) structures to serve as the basis of MPEAs, as the cluster probabilities will closely approach the random limit72,73 defined by the algorithms used to generate these structures. Utilizing SQS structures allows us to quantify the cluster concentrations for each MPEA structure.

As a result of setting εi = 1 / CM N while the energy of each MPEA taken from SQS structures, the probabilities of pairs, triplets, and higher-order clusters are closely aligned. The differences in cluster concentration can be calculated as:
Δn q p = [ N 1 ] p [ C C M N M ] q(M 1 ) p (3)
N
where Δn q p represents the cluster concentration difference of N-dimensional MPEAs as a linear basis (full derivation in Supplementary Note 3). Utilizing Eq. (2), we can predict the energy of an MPEA of any dimension using lower-dimensional MPEAs with the analytic difference of cluster concentration determined by Eq. (3). To show the quantitative trends, the analytic energy difference contributed by self and cross terms of pair interaction can be represented by the top two panels of Fig. 4a, while the energy difference contributed by cross term of triplet interactions can be represented by the bottom panel in Fig. 4a. For all three clusters shown in Fig. 4a, the energy difference decays and approaches zero as M (x value in Fig. 4a) approaches N (indicated as color in Fig. 4a). In other words, when fitting high dimensional MPEA energy (N→∞) with basis with similar but lower dimensionality (M→N-1), the fitting error for all clusters will approximate infinitely to zero. This also indicates that the prediction of high-dimensional MPEA stability will be much more accurate as dimensionality goes up. To validate our model with higher-dimensional DFT-computed MPEAs, we generate an additional dataset comprising 9,362 DFT-calculated equimolar MPEAs ranging from 6M to 10M MPEAs. This includes: (i) a BCC refractory dataset containing 6M, 8M, and 9M MPEAs composed of Ti, V, Cr, Zr, Nb, Mo, Hf, Ta, and W; (ii) two FCC datasets covering 6M, 8M, 9M, and 10M MPEAs, one constructed from the compositional space of Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, and Al, and the other from V, Cr, Mn, Fe, Co, Ni, Cu combined with at least one metal from Ru, Rh, Pd, Ag, Ir, Pt, and Au. The fitting errors using lower-dimensional basis to predict on the higher-dimensional data are shown in Fig. 4b. Our results confirm that low-dimensional MPEAs can reliably estimate the energies of higher-dimensional systems. The same trend of MAE evolution (Fig. 4a) is observed when fitting lower-dimensional DFT-computed SQS structures to higher dimensions (Fig. 4b), where mean absolute error (MAE) is calculated as the average of the absolute errors between the predicted and DFT total energies. Notably, prediction accuracy improves as the dimensionality of the basis approaches that of the target system, and tends to converge at higher dimensions, with MAEs as low as 6.956 meV/atom when using 8M MPEA basis to predict 9M compositions. The predicted total energies of higher-dimensional alloys can then be used to construct the convex hull and determine the corresponding E hull values, which can rapidly estimate synthesizability. When comparing with typical Hume-Rothery descriptors, such as the average and standard deviation of atomic radius and electronegativity, the relationship with formation energy is weak, with the maximum absolute Pearson correlation coefficient (|PCC|) being only 0.380, as shown in Supplementary Fig. 8. This further underscores the strength of our physical model in providing both quantitative accuracy and physical interpretability.

It should be noted that the model is constructed from zero-temperature total energies; it can serve as a baseline for estimating entropy-stabilization temperatures but more accurate estimation of all types of entropies will greatly enhance the accuracy of the stability prediction74. Our DFT fitting analysis not only supports the viability of previously reported theoretical attempts, which use binary formation energy or pairwise energy to predict promising single-phase high-dimensional alloys20,60,75, but also explains why the reported fitting errors are much larger than what we obtained and how to improve.

We also compare the accuracy of our physical model with the state-of-the-art convolutional neural network (e.g., crystal graph convolutional neural network (CGCNN))31. In recent work32, we showed that CGCNN stands out among several graph neural network models for predicting the stability of MPEAs, so we use it for this comparison. By fitting both models with the same dataset, our physical model shows similar, and in most cases, better performance compared to CGCNN (Fig. 4c), and MAEs are shown in Supplementary Table 4. For target dimensionalities greater than five, the training data for CGCNN consists of all lower-dimensional compositions sampled within the same expanded compositional space, rather than the entire dataset of lower-dimensional systems. More importantly, our approach requires no training on GPUs or high-performance CPUs as all we need to do is to determine coefficient εi once the energy basis is computed, while the easiest yet effective εi is 1 / CM N (the comparison with us coefficient is provided in Supplementary Fig. 9). In particular, when using materials with mixed dimensionalities to predict higher-dimensional MPEAs, the prediction error can be further optimized as the cluster concentration from low-dimensional basis can match better with higher-dimensional MPEAs (detailed in Supplementary Note 4, Supplementary Table 5 and Supplementary Fig. 10).

Methods
Multi-principal element alloys (MPEA) structures and phase stability
All disordered MPEAs investigated in this work were generated using the special quasi-random structure (SQS) approach72. The pymatgen package (v2024.5.1)88 and ATAT codes89 were used to generate the SQS structures. Three typical structures, BCC, FCC, and HCP, were generated for every composition3. A supercell consisting of 96 atoms was used for all binary, ternary and quaternary structures, a supercell consisting of 120 atoms was used for quinary structures. BCC and FCC systems with six or eight metals were modeled using 120- or 144-atom supercells; 9-metal systems used 144-atom supercells, and 10-metal FCC systems used 180-atom supercells. We have generated a complete dataset for MPEAs consisting of two, three, or four metal elements selected from 28 elements3, and a customized dataset for high-entropy alloys with 5-10 metals. The structures of competing phases were extracted from the Materials project36 and OQMD datasets37. We used pymatgen to construct the phase diagram33,90 and calculate the energy above hull (denoted as E hull, in units of meV/atom) for all MPEAs.

Electronic structure calculations
All electronic structure calculations were performed using the Vienna Ab-initio Simulation Package (VASP, version 6.2.1)91–93 with PBE functionals94 and PAW pseudopotentials95,96. A reciprocal space discretization of 25 k-points per Å-1 was used to sample the Brillouin zone97. The kinetic energy cutoff was 520 eV. The Methfessel-Paxton scheme was used to smooth the partial occupancies for each orbital98, with a smearing width of 0.2 eV. The convergence criterion for the self-consistent field (SCF) step was 10−5 eV/atom. Geometric optimizations were executed until the force on each atom was less than 0.05 eV.

Experimental materials
Poly(ethylene oxide)-b-poly(2-vinylpyridine) (PEO-b-P2VP, M n = 2.8-b-1.5 kg/mol, M w/M n = 1.11) and polystyrene (M n = 10 kg/mol, M w/M n = 1.09) were purchased from Polymer Source, Inc. All other chemicals were purchased from Sigma-Aldrich and used as received, including manganese(II) nitrate tetrahydrate (Mn(NO 3) 2 · 4H 2O, 99.9%), iron(III) nitrate nonahydrate (Fe(NO 3) 3 · 9H 2O, 99.999%), cobalt(II) nitrate hexahydrate (Co(NO 3) 2 · 6H 2O, 99.999%), nickel(II) nitrate hexahydrate (Ni(NO 3) 2 · 6H 2O, 99.999%), copper(II) nitrate hemi(pentahydrate), (Cu(NO 3) 2 · 2.5H 2O, ≥99.99%), palladium(II) nitrate dihydrate (Pd(NO 3) 2 · 2H 2O, ~40% Pd basis), chloroplatinic acid hydrate (H 2PtCl 6 · xH 2O, ≥99.9%), gold(III) chloride trihydrate (HAuCl 4 · 3H 2O, ≥99.9%), iridium(III) chloride hydrate (IrCl 3 · xH 2O, 99.9%), ruthenium(III) chloride (RuCl 3, 99.99%), rhodium(III) chloride hydrate (RhCl 3 · xH 2O, ≥99.9%), hexamethyldisilazane (HMDS), toluene, and hexane. Pen arrays (type-M, without gold coating) used for dip-pen nanolithography (DPN) were purchased from Advanced Creative Solutions Technology, Inc. Silicon nitride Transmission Electron Microscopy (TEM) window grids (100 nm with 0.1 × 0.1 mm windows) and carbon films (15 nm to 25 nm) coated molybdenum TEM grids were purchased from Ted Pella, Inc.

Preparation of polymer solutions
Polymer ink solutions were prepared by dissolving PEO-b-P2VP and various metal salts in deionized water, with a polymer concentration of 5 mg/mL and total metal salt concentrations of 4.0 mM. The molar ratio of different metal precursors was controlled to be equal. The pHs of the ink solutions were adjusted between 3 and 4 using HNO 3, followed by overnight stirring to ensure complete dissolution before use39.

Multimetallic nanoparticle (NP) synthesis
Scanning probe block copolymer lithography (SPBCL) technology was used to prepare multimetallic NPs. In a typical experiment, a pen array was dip-coated with the as-prepared block copolymer ink and air dried at room temperature. The pen array was then loaded into a modified atomic force microscope instrument (XE-150; Park Systems) in a chamber with humidity at 80%. TEM grids were first washed with isopropanol and dried under a stream of N₂, then rendered hydrophobic by either vapor-coating with HMDS for 12 h in a desiccator or spin-coating with a 5% (w/v) polystyrene solution in toluene at 2500 rpm, followed by baking at 110 °C for 2 min on a hot plate. The modified TEM grids were used as the patterning substrates and transferred into a tube furnace after patterning. The thermal treatment conditions were as follows: ramp to 150 °C under H 2 in 10 min, hold at 150 °C for 10 h, ramp to 500 °C in 1 h, hold at 500 °C for 10 h, ramp to 1000 °C in 1 h, hold at 1000 °C for 2 h, and quench to room temperature.

Characterization
Scanning transmission electron microscopy (STEM) images were taken with a Hitachi HD-2300 STEM at an acceleration voltage of 200 kV. The composition of various NPs was determined with the dual energy-dispersive X-ray spectroscopy (EDS) detectors on the same instrument. Elemental mapping and composition quantification were conducted using the Lα peaks for Au, Pt, Pd, Ir, Rh, and Ru, and the Kα peaks for Mn, Fe, Co, Ni, and Cu from the EDS spectra. Atomic compositions of the NPs were calculated based on EDS spectral analysis with the Cliff-Lorimer correction method, with an inherent error of less than 5% due to X-ray absorption and fluorescence. Background subtraction in the EDS maps was performed using Thermo Scientific NSS software.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors aim to provide a physical model for the synthesizability of multi-principal element alloys (MPEAs), which\napproximates the energy of a high-dimensional alloy (where the dimension refers to the number of elements) using an\naveraged, cluster-expanded, sum of low-dimensional subsets. Their approximation has the attractive feature of a seemingly\nrapid convergence with respect to the subset dimensionality. To validate their model, the authors generated a large density\nfunctional theory (DFT) dataset comprising 126,429 MPEAs based on 28 metals and synthesized nine new MPEAs. The\nmodel’s performance was also compared to a state-of-the-art deep learning model and was found to be comparable, even\nslightly outperforming it in some cases, with the added advantage of interpretability. Key takeaways from the paper include\nan analysis of the compatibility and selectivity of the 28 metals based on their occurrence in the predicted stable MPEAs, as\nwell as the discovery that as the dimension of an MPEA increases, predicting its stability based on lower subset alloys\nbecomes more reliable due to the diminishing influence of mixing enthalpy.\n\nA major revision is required before the paper can be considered further. The concerns raised significantly impact the clarity\nand reliability of the study. Additionally, there may be other issues that become evident once these revisions are made. A\nthorough reassessment of the manuscript after revision is necessary to ensure that all concerns have been adequately\naddressed.\n\nA) Comments on Physics\n1. The authors base the stability model of an MPEA on the relationship E_hull < S_ideal * 1273 K (line 103). This is a well\nknown but rough approximation. While it can serve as a useful heuristic predictor, it does not provide a reliable foundation\nfor a predictive model. Firstly, it neglects important factors such as vibrational free energy, and chemical short-range order.\nFurther, it does not account for the possibility of phase separation into a mixture of medium entropy alloys.\n2.The key novelty of this work lies in systematizing the cluster expansion approach in a manner that easily generalizes to a\nbroad family of alloy systems. Therefore, the title of the paper, \"Universal Stability Rules for Multi-Principal Element Alloys,\"\nas well as the overall narrative, should be revised to better reflect the actual contribution of this study. Specifically, the paper\nis NOT about a \"stability rule\", rather it is about their approximation method.\n\nB) Comments on Results and Discussion\n1. The \"Results\" section primarily describes data preparation rather than presenting actual findings; the key results are\ninstead discussed in the \"Discussion\" section.\n2.The authors synthesized nine new MPEAs. Were these selected at random, or was there a specific rationale behind their\nchoice? Additionally, were any unreported attempted cases unsuccessful?\n3. The model should be validated against MPEAs with more than five elements. For example, in Fig. 4b and 4c, do the errors\ncontinue to decrease beyond 5M MPEAs, or do they plateau? We think the needed data is available because some aspects\nare in 4a.\n\nC) Comments on the Supplementary Notes 1 and 2\n1. Leading up to Eq. 2, the authors state: “As shown in Eq. 3, when using M-dimensional SQS to predict the energy of N-\ndimensional SQS structures (M < N), the formula is:” However, it is unclear how Eq. 2 is derived from Eq. 3, and additionally\nEq. 2 should logically be placed following Eq. 3 if it is derived from Eq.3.\n2. In Eq. 2, the authors do not explicitly write the summation over M, whereas it is included in Eq. 1. This inconsistency\nshould be clarified. Additionally, it should be explicitly stated that an M-subset alloy is composed of clusters.\n3. The replacement of \\epsilon with 1/C_N^M​ in Eq. 2 should be clarified. Is this merely a convenient choice (as mentioned\nin line 194 of the main text) or can it be derived on some basis? In the main text it would be good to explain how this choice\n“help(s) conserve the concentration”.\n4. In Eq. 4, how does the singlet cluster cancel out? In fact, why is it present in Eq. 3? It seems that M=1 terms vanish for\nenthalpies of formation (hence M>=2 in Eq. 1) so why is it present in (3) but not (4)?\n5. In Eq. 7, it does not appear that \\Delta n is nonzero for a singlet cluster with p = q = 1, which contradicts what Eq. 4\nsuggests. Could the authors explain this discrepancy?\n6. The statement in line 64—“It is possible ... concentration difference.”—is unclear and needs rewording for better\ncomprehension.\n7. Additionally, it seems that Eqs. 9 and 10 are missing 1/C_N^M​ factors.\n\nD) Comments on Writing and Clarity\nThe abstract should be rewritten for clarity and conciseness. Specifically:\n1. The physical model is not explicitly defined, and the fact that it is based on cluster expansion and lower dimensional\nsubsystems is not mentioned.\n2. The last sentence is unclear. The phrase “with much fewer DFT data than traditionally expected” is vague and lacks a\nclear definition or quantification. Additionally, “the unprecedented diminishing of mixing enthalpy” should be rewritten for\nclarity.\nSpecific sentences:\n1. \"The universal rule ... by simply averaging the stability ...\" (line 45) \"Averaging the stability\" is not a well-defined term and\nshould be rephrased more precisely.\n2. \"We uncovered ... can be achieved with small datasets.\" (line 48) The term \"small datasets\" should be clearly defined or\nquantified to provide meaningful context.\n3. \"Noble-metal-containing MPEAs ....\" (line 85) This sentence does not fit well as the opening sentence of the Experimental\nValidation of Predicted MPEAs subsection. The sudden shift in focus to noble-metal-containing MPEAs lacks prior context\nand should be introduced more smoothly.\n4. \"The value of C_N^M​ represents the total types of M-dimensional compositional spaces among the N-dimensional\nsystem.\" (line 191) This is misleading; it is not the total types but rather the number of possible choices (i.e., the binomial\ncoefficient). Please mention that C_N^M is a binomial coefficient.\n5. Line 34 in the supplementary notes: \"can be expressed in a spacious way.\" This should be reworded as \"can be\nexpanded.\"\n\nE) Comments on Figures\n1. Fig. 1A and b presents widespread knowledge and may not be necessary. Consider removing them. Instead be sure to\nrefer to 1c in the text and discuss its significance.\n2. Fig. 2a, which summarizes published papers on noble-metal-containing MPEAs, is more suited for a review paper rather\nthan this study. It does not fit well within the scope of this work and should be reconsidered.\n3. Fig. 3 is difficult to interpret, particularly due to the large number of elements. Additionally, in line 163, the authors state:\n\"We can also link the high frequency of specific elemental pairs, such as Pd-Pt, Rh-Pt, and Rh-Pd reported in the literature,\nto their high probability (thick segment, Fig. 3b).\" However, using visual metrics like \"size\" and \"thickness\" makes direct\ncomparisons challenging. We are never told how the thicknesses were chosen. We suggest replacing Fig. 3 with a clearer\nfigure that directly reflects the key insights highlighted in the text, such as the most compatible and most selective elements.\nThe central figure 3a is very hard to interpret.\nThe remaining data could be moved to tables in the supplementary materials for better clarity.\n4. Fig. 4a presents data up to an 8M MPEA, whereas Fig. 4b only includes data for 4M and 5M MPEAs. Would it be possible\nto add more data to Fig. 4b for consistency?\n5. Fig. S1 is confusing and unclear. How can the diagonal entries be understood? i.e. what does the pair Ru-Ru mean in a\nbinary?", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_3": "Authors used two stability criteria to predict MPEA synthesizabilty, Ehull and Eform. The first was obtained by constructing a\nphase diagram from already available structure databases (pymatgen and OQMD) and an in-house high-throughput DFT\nenergy calculation. The energy above hull is used as a criterion for entropic stabilization by assuming an ideal entropy\ncontribution to the free energy. This stability criterion is then tested with experimental synthesis of 9 equimolar MPEAs.\nAs a shortcut for total energy calculation without the need of DFT calculations a linear model based on Cluster Expansion is\nproposed and trained on the previously generated DFT dataset. The model is based on extrapolation from lower to higher\ndimensionality MPEAs by approximating the total energy to weighted contributions of lower dimensionality MPEAs. A\nmathematical analysis is carried out in the supplementary material to theoretically assess the error accompanied by the\nmethod. The authors state that theoretically the error should be a negligible fraction of the pair, three-body, and the rest of\ncluster interactions.\n\nThe contribution from the authors to MPEAs design is substantial: a quick estimation of stability through a linear\napproximation of the free energy of MPEAs. The work is comprised by two main parts: stability rule and DFT total energy\nmodel.\n\nI would advise the following major revisions for the manuscript to be ready for publication.\n◦ While the Hume-Rothery rules are mentioned, a proper comparison with previously stablished rules and composition\ntrends is not fleshed out. Specially how this new stability criteria squares against previously shown soft rules in HEA design.\n◦ The claim is bold but only one scatter plot is given with little information on what is being shown here. The x-axis and\ncolor-scale communicate both the same information. Good properties to compare could be averaged mixing enthalpy, ideal\nmixing entropy, VEC, electronegativity, lattice misfit, etc.\n◦ Figure 1c. Specify the energies are calculated at zero temperature (zero Kelvin). It is recommended to provide subplots by\ncrystal structures: FCC, BCC, and HCP. what is the significance of the black solid line (E_form = E_hull)? Phase stability is\nall about comparing the energies among different crystal structures. What matters the most here is the E_hull that is\ncalculated based on E_form. I expect more meaningful plots of E_hull vs crystal structure vs number of component vs critical\ntemperature above which entropy would favor forming solid solution (T * S >= E_hull).\n◦ Mind that you could use the simple expression: T * S >= E_hull to estimate the entropy stabilizing effect, however, T * S >=\nE_hull does not always guarantee forming solid solution at temperatures above the critical temperature. This is because\ntemperature brings in the contribution of vibrational free energy, and competing compound phases could exhibit\ncompositional homogeneity range so that configurational entropy will play as well for compounds.\n◦ The authors wrote: “In our newly established dataset of 126,429 DFT-computed MPEAs, 8,303 MPEAs (6.56% 127 of all\ncompounds) with three to five elements are predicted to be entropy-stabilizable below 1273 K.” Do you really mean above\n◦ The authors wrote: “The accuracy of our physical model is comparable to, or better than, that of state-of-the-art deep\nlearning models trained on the same dataset, while retaining full interpretability.” This appears to be overstatement without\nproviding evidence to support their claim. The authors wrote: “Such performance matches the state-of-the-art graph neural\nnetwork [43] models trained on the same dataset [28].” However, reference [28] does not address high entropy alloys nor on\nthe same database.\n◦ Synthesizability of high entropy materials is reported and modeled in the literature, and citing them is very appropriate, to\nlist a few: (1) Sarker et al. High-entropy high-hardness metal carbides discovered by entropy descriptors. Nat Commun 9,\n4980 (2018), (2) Eckert et al. \"A formula to predict the synthesizability of high-entropy materials.\" NATURE (2024).\n◦ For the following model to be based on the cluster expansion method, the configurational entropy is glossed over in the\nstability rule definition. Previous works have shown entropy is not trivial at temperatures discussed here. Especially, when\nreferring to entropic stability as low as in 500K. I think a distinction must be made for the model to have a comfortable\napplication temperature threshold and a discussion on the entropy exceptions could be added.\n◦ The high-throughput effort is impressive here, and experimental validation is added and meaningful to the material. But the\nscope of the project is too big to be assessed by nine MPEAs. I recommend a comparison with previously reported literature\n(a MPEA experimental database).\n◦ An analysis on non-equimolar compositions could enhance the model extrapolation capabilities. To the best of my\nunderstanding the model is a truncated cluster expansion (only point clusters) trained on only M-sized equimolar\nsubsystems. The point cluster interaction should still be viable to use for non-equimolar systems but at a higher error cost.\nGiven the compositional design having an important role in MPEA design I believe this may be worth discussing.\n◦ I highly encourage the authors to share all the compositions vs energies data to benefit the entire community.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_5": "The authors address a very important and urgent issue: predicting the synthesizability of MPEAs. They propose a physical\nmodel to predict MPEA stability, achieving an accuracy of 13.38 meV/atom. The model has been validated by fabricating\nnine novel MPEAs based on its predictions. The results obtained are potentially significant and interesting for the materials\nscience community. However, there are two serious concerns that prevent me from recommending the article for publication\nin its current form.\n1. The authors compare their results with those obtained using a state-of-the-art data-driven CGCNN model. While this\ncomparison suggests the reliability of their proposed model, it raises the question of whether the computational resources\nrequired for their approach are justified, given that its accuracy is only slightly better than that of an existing model. The\nauthors should discuss this point in more detail, demonstrating whether their approach offers superior accuracy,\nperformance, or transferability compared to other models like CGCNN.\n2. The proposed approach relies on a large dataset of 126,429 DFT-computed MPEAs, implying significant computational\nresources were expended to generate this data. there is a state-of-the art approach to perform any atomistic simulations\n(particularly the same performed be the authors) which is based on Universal Interatomic Potentials (or Foundation Models).\nThese models, based on Graph Neural Networks and trained on extensive ab initio datasets, demonstrate high accuracy\nand transferability. The best-performing models achieve accuracy comparable to that reported by the authors (see\nhttps://matbench-discovery.materialsproject.org/ for reference). Moreover, their computational efficiency is orders of\nmagnitude higher than that of ab initio calculations. Thus, using such models might yield similar accuracy to the one\npresented in the manuscript, which could significantly diminish the work's novelty and impact. The authors should address\nthis concern and discuss their findings in this context."}, "author_response": "Response Letter for\nUniversal Framework for Efficient Estimation of Stability in Multi-\nPrincipal Element Alloys\nLin Wang1†, Bo Shen2,3†, Zheng-Da He1, Zihao Ye2,3, Yan Zeng1, Chad Mirkin2,3, and Bin\nOuyang1,*\n1Department of Chemistry and Biochemistry, Florida State University, Tallahassee, FL\n32304, USA.\n3International Institute for Nanotechnology, Northwestern University, Evanston, Illinois\n60208, USA.\n†These authors contributed equally to this work\n*To whom correspondence may be addressed. E-mail: bouyang@fsu.edu\nmodified/inserted texts in the manuscript and supporting information are colored in red. We believe\nthis work is greatly improved to fit better for publishing in Nature Communications.\nREVIEWER COMMENTS:\nReviewer #1 (Remarks to the Author):\nThe authors aim to provide a physical model for the synthesizability of multi-principal element\nalloys (MPEAs), which approximates the energy of a high-dimensional alloy (where the dimension\nrefers to the number of elements) using an averaged, cluster-expanded, sum of low-dimensional\nsubsets. Their approximation has the attractive feature of a seemingly rapid convergence with\nrespect to the subset dimensionality. To validate their model, the authors generated a large density\nfunctional theory (DFT) dataset comprising 126,429 MPEAs based on 28 metals and synthesized\nnine new MPEAs. The model’s performance was also compared to a state-of-the-art deep learning\nmodel and was found to be comparable, even slightly outperforming it in some cases, with the\nadded advantage of interpretability. Key takeaways from the paper include an analysis of the\ncompatibility and selectivity of the 28 metals based on their occurrence in the predicted stable\nMPEAs, as well as the discovery that as the dimension of an MPEA increases, predicting its\nstability based on lower subset alloys becomes more reliable due to the diminishing influence of\nmixing enthalpy.\nA major revision is required before the paper can be considered further. The concerns raised\nsignificantly impact the clarity and reliability of the study. Additionally, there may be other issues\nthat become evident once these revisions are made. A thorough reassessment of the manuscript\nafter revision is necessary to ensure that all concerns have been adequately addressed.\nResponse: We would like to thank the reviewer for their positive comments on our work. We have\naddressed all questions and suggestions and provided a point-by-point response below. The revised\nmanuscript has been improved in terms of clarity and reliability.\nA) Comments on Physics:\n1. The authors base the stability model of an MPEA on the relationship E_hull < S_ideal * 1273\nK (line 103). This is a well-known but rough approximation. While it can serve as a useful heuristic\npredictor, it does not provide a reliable foundation for a predictive model. Firstly, it neglects\nimportant factors such as vibrational free energy, and chemical short-range order. Further, it does\nnot account for the possibility of phase separation into a mixture of medium entropy alloys.\nResponse: We thank the reviewer for the insightful comments. We agree that the stability analysis\ncriterion based on E < S * 1273 is indeed a well-known rough approximation. However, it\nhull ideal\nremains a valuable heuristic for predicting material stability and has facilitated numerous\nsignificant discoveries in compositionally complex and high-entropy materials, as demonstrated\nby the following recent studies with DOIs being: 10.1126/sciadv.1600225, 10.1038/s41563-020-\n00816-0, 10.1038/s41467-024-45309-9, 10.1002/aenm.202404982 and 10.1038/s41467-021-\n26006-3. Moreover, in this revised manuscript, we further demonstrate that, in addition to\npredicting the 9 newly synthesized MPEAs, this stability approximation can successfully predict\n30 out of 32 experimentally reported MPEAs synthesized at temperatures below 1273 K (Tab. S2).\ncompositions for experimental validation.\nThe reviewer has also provided suggestions to consider multiple factors, such as vibrational\nentropy, chemical short-range order, and phase separation to low entropy disordered phases in\nstability evaluation. While these are all useful physical considerations, accounting for them would\nrequire millions to billions of additional calculations, which is currently computationally\nunfeasible. If such factors are to be estimated using new methodologies, the effort and complexity\ninvolved would merit a dedicated study of their own. We have elaborated on these points in the\nfollowing paragraphs, supported by relevant literatures and back-of-the-envelope estimates of the\ncomputational workload. Finally, we have also provided a clear articulation of both the limitations\nand the novel contributions of our work to ensure transparency and clarity.\nVibrational Entropy: We acknowledge that including vibrational free energy would enhance the\naccuracy of stability predictions. However, the computational cost of estimating vibrational\nentropy is prohibitively high and currently impractical, particularly for the >100,000 materials\nexamined in this study. To illustrate, estimating the vibrational entropy of a single MPEA structure\ncontaining 160 atoms, such as those used in our quinary MPEA calculations, requires generating\napproximately 960 displaced structures for static DFT calculations. Scaling this process to over\n100,000 materials would demand more than 90 million DFT calculations. While we are actively\ndeveloping more efficient methods for estimating vibrational entropy, the complexity and scope\nof this task represent a significant departure from the main focus of the present work. As such, we\nbelieve this effort is better suited for a dedicated future publication.\nChemical short-range order (CSRO): Similar to vibrational entropy, the inclusion of chemical\nshort-range order (CSRO) would indeed be valuable for enhancing the accuracy of stability\npredictions. However, such inclusion is again computationally impractical, and more importantly,\nless relevant to the main objective of this work. To elaborate, we draw on our eight years of\nexperience in modeling CSRO (DOIs: 10.1002/aenm.201903240, 10.1016/j.chempr.2019.10.001,\n10.1038/s41563-020-00816-0, 10.1103/PhysRevB.106.144202) and calibrate modelled CSRO\nwith experimental studies (DOIs: 10.1038/s41467-023-43356-2, 10.1021/acs.chemmater.2c03827)\nto estimate the computational cost of incorporating CSRO in the present study. These prior works\nindicate that approximately 1,000 DFT calculations are required per compositional space to obtain\na reasonable estimation of CSRO. A compositional space refers to, for instance, the quinary system\nCo–Cr–Fe–Mn–Ni. Given that this study aims to cover 28 elements across all possible quinary\ncombinations, the total number of required calculations would be:\n1000 structures per composition space × 𝐶5 compositional spaces × 3 crystal structures (BCC,\nFurthermore, our theoretical framework for estimating the energy of multi-principal element alloys\n(MPEAs) is not designed to just predict synthesizability or CSRO decorated systems, but all\nMPEAs across broad chemical space. We employ a large pool of special quasi-random structures\n(SQS) because they are widely regarded as a reasonable and efficient approach for approximating\nthe energetics of disordered alloy systems. While CSRO may exist in MPEAs, it is generally\nbelieved to be less pronounced than in low-entropy systems (DOIs: 10.1038/s41563-020-00816-\n0, 10.1038/s41563-023-01570-9, 10.1038/s41586-020-2275-z, 10.1038/s41467-020-16083-1),\nmaking the SQS approach particularly well-suited for high-entropy alloy modeling.\nestimation is simplified by not accounting for the potential decomposition into lower-dimensional\ndisordered phases. This is a well-known and unresolved challenge in MPEA simulations.\nAccurately capturing the competition with such disordered phases would require knowledge of the\ncomplete free energy landscape across all phases and compositions within a given alloy system.\nPerforming this level of analysis for all high-entropy alloys involving the 28 elements considered\nin this study is currently unfeasible due to the immense computational cost. As detailed in the\nprevious two paragraphs, such an undertaking would require both accurate vibrational entropy and\nCSRO estimations and inclusion of all possible phases. The total computational cost would exceed\nthe product of the individual costs associated with estimating vibrational entropy and CSRO,\nmaking it far beyond current practical capabilities.\nStrategies to strengthen this work and clarify the novel point:\nIn light of the many outstanding challenges of MPEA theory that remain unsolved, we believe the\nmost appropriate way to present our work is by avoiding unnecessary hype, openly acknowledging\nthe limitations, and emphasizing the novelty and scientific advances it offers. In this work, despite\nusing a relatively simple yet widely accepted approximation for phase stability (i.e., E ), our\nhull\napproach has demonstrated practical utility in guiding experimental MPEA design, validated by\nthe successful synthesis of nine predicted compounds and previous experimentally reported\nequimolar materials. We evaluate the E values at 1273 K (E -1273 × S ) for experimentally\nhull hull Ideal\nreported equimolar compounds synthesized below 1273 K, and found that they exhibit E values\nhull\naround or below 0 meV/atom at 1273 K. This observation supports the use of the E < S × T\nhull ideal\nas a reasonable stability approximation. More importantly, our theoretical framework enables rapid\nestimation of both total energy and phase stability for MPEAs with minimal computational effort,\nwhile achieving accuracy comparable to that of state-of-the-art graph neural networks. We believe\nthese advances provide a valuable and efficient toolset for the community to address broader\nchallenges, including those highlighted by the reviewer, in future studies.\nThe revisions have been detailed as below and colored in red in the revised manuscript:\nTo weaken our tone and acknowledge the limitation of the current stability estimation of MPEA,\nwe have added the following statement in the revised manuscript on pages 6-7, which is also\nattached below:\nPage 6: “Although the criterion E < S × T is a rough approximation, it has been widely\nhull ideal\nadopted as an effective indicator for predicting stability in compositionally complex and high-\nentropy materials. Incorporating vibrational free energy and chemical short-range order (CSRO),\nwould certainly enhance accuracy, but doing so remains computationally unfeasible for large-scale\ndatasets, especially when considering hundreds of thousands of compositions across a 28-element\nspace. It would also be beneficial to account for potential competition with lower-dimensional\ndisordered phases. However, such efforts are currently impractical, as they require mapping the\ncomplete free energy landscape across all relevant phases and compositions, including accurate\nassessments of both vibrational entropy and CSRO. Given these constraints, the approximation\nadopted in this study remains the most practical and computationally efficient approach for guiding\nthe discovery among global space of MPEAs, which is also validated by the successful synthesis\nof many predicted compounds.”\nhull\nof 32) or below (30 out of 32) zero, which are listed in Tab. S2. This observation supports the use\nof the E < S × T criterion as a reasonable approximation for thermodynamic stability.”\nhull ideal\nTo strengthen the major focus, novelty, and advancement of the field of this work, we have added\nthe following statement in the revised manuscript on pages 15 and 17, which is also attached below:\nPage 16: “More importantly, our approach requires no training on GPUs or high-performance\nCPUs as all we need to do is to determine coefficient 𝜀𝑖 once the energy basis is computed, while\nthe easiest yet effective 𝜀𝑖 is 1 .”\n𝐶𝑀\n𝑁\nPages 18: “Building on this thermodynamic interpretation, our physical linear model predicts the\nenergy and hence stability of high-dimensional MPEAs using weighted linear average of the lower-\ndimensional MPEA as basis. Fitting such model eliminates the need for training on GPUs or high-\nperformance CPUs. As determination of coefficient is through solving a set of linear equations or\nbe as simple as setting weighting factor 𝜀𝑖 = 1 . For specific subspace of MPEA of interest, there\n𝐶𝑀\n𝑁\nis no need to develop a comprehensive training data, but simply lower entropy equimolar basis at\na certain dimension, e.g., ternary or quaternary MPEAs. Moreover, the application of our model\ndoes not need further relaxation using Hamiltonian such as machine learning force field, which are\nrequired by the state-of-the-art foundational models. We thus believe this discovery will facilitate\nthe rapid exploration of novel MPEA space.”\n2.The key novelty of this work lies in systematizing the cluster expansion approach in a manner\nthat easily generalizes to a broad family of alloy systems. Therefore, the title of the paper,\n\"Universal Stability Rules for Multi-Principal Element Alloys,\" as well as the overall narrative,\nshould be revised to better reflect the actual contribution of this study. Specifically, the paper is\nNOT about a \"stability rule\", rather it is about their approximation method.\nResponse: Thanks a lot for the reviewer to acknowledge our novelty. We greatly appreciate the\nsuggestion for revising the title, the title has been changed to “Universal Framework for Efficient\nEstimation of Stability in Multi-Principal Element Alloys” in our revised manuscript and\nsupplementary information.\nB) Comments on Results and Discussion\n1. The \"Results\" section primarily describes data preparation rather than presenting actual findings;\nthe key results are instead discussed in the \"Discussion\" section.\nResponse: Thanks for the constructive feedback. We have revised the manuscript by merging all\nkey findings into the new 'Results' section and reserving the 'Discussion' for interpretation and\npractical implications of our models. The change can be reflected by the updated section titles on\npages 9 and 11 and the contents on pages 9-17.\n2.The authors synthesized nine new MPEAs. Were these selected at random, or was there a specific\nrationale behind their choice? Additionally, were any unreported attempted cases unsuccessful?\nhas been provided below:\n1. They were systematically selected from theoretically predicted candidates based on the lab\nsynthesis limitations. Specifically, we focused on noble-metal-based compositions, since\nseveral non-noble metal precursors (especially early transition metals such as Ti, V, and Cr)\nare significantly more challenging to reduce under an H atmosphere due to their unfavorable\n2\nreaction Gibbs free energies. Although entropy can facilitate metal reduction in multi-element\nalloys, the presence of noble metals will significantly enhance the overall reducibility and\nphase formation.\n2. Additionally, we excluded certain combinations incorporating highly toxic metals (e.g., Os)\nfor safety considerations, even though many Os-containing MPEAs show reasonable stability\nas indicated in Fig. 3. Furthermore, our annealing furnace has a maximum operating\ntemperature of 1000 °C, further limiting the accessible compositional space and eliminating\nelements that are easily melt or boiled (e.g., Zn and Sn).\n3. Finally, it is important to note that this technique has been attempted synthesis of nanoparticles\nincorporating early transition metals and other active metals (Ti, V, Cr, Mg, Al), but only metal\noxide can be obtained instead of metals (DOI: 10.1002/adma.202500088), hence we did not\nattempt such cases here.\n4. To conclude, we prioritize MPEAs that 1) exhibit low predicted entropy stabilization\ntemperatures (T ); 2) contain a high fraction of noble metals to ensure synthetic feasibility,\ns\nand; 3) represent compositions that have either not been reported previously or remain\nunderrepresented in the literature (Tab. S1); 4) Contained no Ti, V, Cr, Mg, Al, Sn and Zn.\nTo reflect the above-mentioned changes, we have modified page 8 of the revised manuscript. The\nchanged statement is also attached below for quick reference:\n“These nine MPEAs were carefully selected based on synthesis feasibility and theoretical\npredictions. We focused on noble-metal-rich compositions due to their enhanced reducibility under\nH , excluding early transition metals (e.g., Ti, V, Cr), other active metals such as Mg and Al, and\n2\ntoxic or volatile elements (e.g., Os, Zn, Sn) for safety and experimental limitations. Selection\nprioritized alloys with low predicted stabilization temperatures, high noble metal content, and\ncompositions that are novel and never reported before (with reported compounds list in Tab. S1).\n3. The model should be validated against MPEAs with more than five elements. For example, in\nFig. 4b and 4c, do the errors continue to decrease beyond 5M MPEAs, or do they plateau? We\nthink the needed data is available because some aspects are in 4a.\nResponse: We want to thank the reviewer for this insightful question. We fully agree that\nincorporating MPEA data with more than five elements would further strengthen our work. In\nresponse, we have extended our dataset to include compositions with up to 10 constituent metals,\nenabling a broader validation of our model. Due to the exponential growth of the compositional\nspace, the number of equimolar 6-element combinations from 28 elements reaches 376,740 (𝐶6 ),\n28\nwhile the number for 10-element combinations increases to 13,123,110 (𝐶10). Recognizing that it\n28\nis computationally unfeasible to exhaustively evaluate all possible high-order combinations as we\ndid for 3M and 4M MPEAs, we instead focused on a representative subset of high-dimensional\nMPEA compositions. Specifically, we targeted FCC and BCC MPEAs composed of more\ncommonly studied elements, and constructed the following three extended datasets:\n6\n1) For BCC-MPEAs, we generated all possible equimolar combinations among nine typical\nelements commonly found in refractory high-entropy alloys, e.g., Ti, V, Cr, Zr, Nb, Mo, Hf, Ta,\nand W. We extended the compositional dimensionality to 6M, 8M, and 9M, resulting in a total of\n94 (𝐶6+𝐶7+𝐶9) new BCC-MPEAs.\n9 9 9\n2) We considered 10 elements typically observed in FCC MPEAs: Ti, V, Cr, Mn, Fe, Co, Ni, Cu,\nZn, and Al. By generating all possible combinations in 6M, 8M, 9M, and 10M configurations, this\ndataset includes 266 (𝐶6 +𝐶8 +𝐶9 +𝐶10) additional FCC-MPEAs.\n10 10 10 10\n3) To strengthen the connection with both our experimental work (Fig. 2) and prior literature (Tab.\nS1), we created a noble-metal-enriched dataset. This includes all 6M, 8M, 9M, and 10M\ncombinations composed of at least one noble metal (Ru, Rh, Pd, Ag, Ir, Pt or Au) along with non-\nnoble metals (V, Cr, Mn, Fe, Co, Ni, Cu). This approach yielded 9,002 additional FCC-MPEAs.\nObservations:\nAs a result, we have added 9,362 new DFT-computed MPEAs, significantly expanding our dataset\nto include compositions from 6M up to 10M. This extended dataset enables a more robust\nvalidation of our predictive model and reinforces several key conclusions from our original\nmanuscript:\n1) Transferability across dimensions: Low-dimensional MPEAs can serve as an effective basis\nfor estimating the energies of higher-dimensional systems. For example, using 8M-MPEAs to\npredict 9M-MPEAs yields a minimum mean absolute error (MAE) as low as 6.956 meV/atom.\n2) Dimensional proximity improves accuracy: The prediction error generally decreases as the\ndimensional gap between the basis and target systems narrows. For instance, using 6M-MPEAs to\npredict 9M-MPEAs results in a lower MAE than using 4M-MPEAs as the basis.\nUpdates on comparison with CGCNN model: We have also trained additional CGCNN models\nusing the same training and testing datasets as our physical model in Fig. 4b. This guarantees that\nthe comparison in Fig. 4c reflects models evaluated under identical training and testing conditions.\nWe have updated Figs 4b-c to reflect the trend from additional dataset. Corresponding revisions\nhave also been made to the manuscript text on pages 14–16 to clearly present these new findings.\nPages 14-15: “To validate our model with higher-dimensional DFT-computed MPEAs, we\ngenerated an additional dataset comprising 9,362 DFT-calculated equimolar MPEAs ranging from\n6M to 10M MPEAs. This includes: (i) a BCC refractory dataset containing 6M, 8M, and 9M\nMPEAs composed of Ti, V, Cr, Zr, Nb, Mo, Hf, Ta, and W; (ii) two FCC datasets covering 6M,\n8M, 9M, and 10M MPEAs, one constructed from the compositional space of Ti, V, Cr, Mn, Fe, Co,\nNi, Cu, Zn, and Al, and the other from V, Cr, Mn, Fe, Co, Ni, Cu combined with at least one metal\nfrom Ru, Rh, Pd, Ag, Ir, Pt, and Au. The fitting errors using lower-dimensional basis to predict on\nthe higher-dimensional data are shown in Fig. 4b. Our results confirm that low-dimensional\nMPEAs can reliably estimate the energies of higher-dimensional systems.”\nPage 15: “Notably, prediction accuracy improves as the dimensionality of the basis approaches\nthat of the target system, and tends to converge at higher dimensions, with MAEs as low as\n6.956 meV/atom when using 8M MPEA basis to predict 9M compositions.”\nPages 15-16: “By fitting both models with the same dataset, our physical model shows similar,\nand in most cases, better performance compared to CGCNN (Fig. 4c), and MAEs are shown in\nTab. S5. For target dimensionalities greater than five, the training data for CGCNN consists of all\n7\nlower-dimensional compositions sampled within the same expanded compositional space, rather\nthan the entire dataset of lower-dimensional systems.”\nFig. 4. Performance of our predictive model. a, The theoretically derived MAE of predicting energy\ncontribution from selective cluster 4M-10M MPEA using 2M-9M MPEA as basis. J is used to represent the unit of\n0\nenergy, with the physical meaning being the energy contribution of certain cluster, e.g., 𝑱 , 𝑱 , 𝑱 respectively. b,\n𝒊𝒊 𝒊𝒋 𝒊𝒋𝒌\nThe MAE for using DFT-computed total energy of 2M-9M MPEA to predict 4M-10M MPEAs as a contrast to the\ntheoretical derivation. c, The comparison between our physical model with a crystal graph convolutional neural\nnetwork (CGCNN).\nC) Comments on the Supplementary Notes 1 and 2\nResponse: We sincerely thank the reviewer for careful reading through our equations and for\nidentifying both typographical errors and potential issues of ambiguity. We have addressed all\nseven equation-related questions individually, as detailed in the following responses. Additionally,\nwe have thoroughly reviewed all equations in the manuscript and supplementary information to\nensure that no similar typographical issues remain.\n1. Leading up to Eq. 2, the authors state: “As shown in Eq. 3, when using M-dimensional SQS to\npredict the energy of N-dimensional SQS structures (M < N), the formula is:” However, it is\nunclear how Eq. 2 is derived from Eq. 3, and additionally Eq. 2 should logically be placed\nfollowing Eq. 3 if it is derived from Eq.3.\nResponse: We would like to clarify that equation 2 is not derived from equation 3. We realized it\nwas a typo in our previous statement. It has been changed to “As shown in Eq.1” We added 𝜀𝑖 =\n1\nto help understand the Eq. 2. The changes are updated in the Supplementary Information, page\n𝐶𝑀\n𝑁\n2.\n8\n“One straightforward approach to implement such a linear combination is to assign equal weights\nto each basis composition while ensuring that the elemental fractions collectively match those of\nthe target higher-dimensional MPEA. This results in: 𝜀𝑖 = 1 ”\n𝐶𝑀\n𝑁\n2. In Eq. 2, the authors do not explicitly write the summation over M, whereas it is included in Eq.\n1. This inconsistency should be clarified. Additionally, it should be explicitly stated that an M-\nsubset alloy is composed of clusters.\nResponse: We appreciate the reviewer’s careful effort to validate and interpret our equations. To\nclarify, Eq. 1 should not include a summation over M, as the enumeration of all possible M-\ncomponent alloys is already captured by the index 𝑖 ∈ (1,𝐶𝑀), which refers to all equimolar M-\n𝑁\ncomponent MPEAs within a pool of N elements. We have updated the equation accordingly on\nSupplementary Page 2, as shown below:\n𝐸 ≈ ∑ 𝜀𝑖𝐸𝑖 (1)\n𝑁 𝑀\n𝑖∈(1,𝐶𝑀) 𝑀\n𝑁\n3. The replacement of \\epsilon with 1/C_N^M in Eq. 2 should be clarified. Is this merely a\nconvenient choice (as mentioned in line 194 of the main text) or can it be derived on some basis?\nIn the main text it would be good to explain how this choice “help(s) conserve the concentration”.\nResponse: We appreciate the reviewer’s thoughtful question. The setup of 𝜀𝑖 = 1 is a convenient\n𝐶𝑀\n𝑁\nchoice; in such a case we can conserve the concentration of each element in both left and right of\nEq. 2. We have updated our supplementary information in page 2 to avoid further confusion. It is\nalso shown as below:\n“One straightforward approach to implement such a linear combination is to assign equal weights\nto each basis composition while ensuring that the elemental fractions collectively match those of\nthe target higher-dimensional MPEA. This results in: 𝜀𝑖 = 1 ”\n𝐶𝑀\n𝑁\n4. In Eq. 4, how does the singlet cluster cancel out? In fact, why is it present in Eq. 3? It seems\nthat M=1 terms vanish for enthalpies of formation (hence M>=2 in Eq. 1) so why is it present in\n(3) but not (4)?\nResponse: With our clarification in response to questions C1-C3 above, the reviewer may already\nunderstand the rationale behind: the setup 𝜀𝑖 = 1 is primarily intended to ensure cancellation of\n𝐶𝑀\n𝑁\nthe singlet term (i.e., M=1 terms). For a given singlet i, its concentration in an M-component\nequimolar MPEA is 1/M, whereas in an N-component equimolar MPEA, it is 1/N. There are 𝐶𝑀−1\n𝑁−1\nM-component subsets that contain element i. When summing over all M-component equimolar\nMPEAs under N-dimensional space with weighing factor 𝜀𝑖 = 1 , the total contribution of singlet\n𝐶𝑀\n𝑁\ni becomes:\n9\n1 1 1\n×𝐶𝑀−1 × =\n𝐶𝑀 𝑁−1 𝑀 𝑁\n𝑁\nwhich exactly matches the concentration of i in N-component system. As a result, the singlet\ncontribution effectively cancels out when computing the energy differences between 𝐸 (𝑆𝑄𝑆) and\n𝑁\n∑ 𝐶 𝑁 𝑀 1 𝐸𝑖 (𝑆𝑄𝑆) . We have added the following sentence in the revised supplementary\n𝑖=1𝐶𝑀 𝑀\n𝑁\ninformation in page 3 to clarify this point.\n“The singlet contributions cancel out because the weighted sum over all M-component subsets\nexactly reproduces the concentration in the N-component system.”\n5. In Eq. 7, it does not appear that \\Delta n is nonzero for a singlet cluster with p = q = 1, which\ncontradicts what Eq. 4 suggests. Could the authors explain this discrepancy?\nResponse: Such outcome is actually expected and further validate the consistency of our\ntheoretical framework. When p=q=1, ∆𝑛 represents the difference in the concentration of a singlet\nbetween the N-component system and the average over its M-component subsets. As shown in our\nderivation and clarifications to questions C1-C4, this difference evaluates to zero:\n1 p C M−q 1 p 1 CM−1 1 1 M 1 1 1\nN−q N−1\nΔn q = ( ) − ( ) = ( )− ( ) = ( )− ( ) = ( )−( ) = 0\nc p N CM M N CM M N N M N N\nN N\n6. The statement in line 64, “It is possible ... concentration difference.”, is unclear and needs\nrewording for better comprehension.\nResponse: Thank the reviewer for raising this question. Based on our responses to questions C1-\nC5, a natural follow-up question is: “why should we stick with 𝜀𝑖 = 1 ? Is there a better way of\n𝐶𝑀\n𝑁\nchoosing 𝜀𝑖?”. The referred statement by the reviewer here is actually to address such potential\nquestion. We want to demonstrate that there are numerous ways of determining linear combination\ncoefficient 𝜀𝑖, the general principle is to cancel as many clusters as we can by choosing a good\nenough 𝜀𝑖 as shown in Eq. 1, so that the clusters defined in Eq. 3 can be cancelled as much as\npossible in a way in analog to Eq. 4.\nWe have revised the sentence for clarity. The updated sentence is on pages 18 in the revised\nsupplementary information, which is also attached below:\n“Differences in higher-order clusters, such as pairs, triplets, and quadruplets, cannot be fully\neliminated when using only single-dimensional MPEA SQS structures as the basis. However, some\nof these cluster contributions can be effectively cancelled by combining MPEAs of different\ndimensionalities. In general, the inclusion of a broader range of dimensionalities in the basis allows\nmore cluster terms to cancel out, potentially leading to improved prediction accuracy. In this work,\nwe demonstrate two examples for determining the linear combination coefficients when using\nMPEAs with two lower dimensions to predict the properties of higher-dimensional MPEA.”\n7. Additionally, it seems that Eqs. 9 and 10 are missing 1/C_N^M factors.\nand 10 and also shown as below:\n𝐶𝑀1 𝐶𝑀2\n𝑁 𝑁\n𝑎0−𝑎2 1 𝑎1−𝑎0 1\n𝐸 (𝑆𝑄𝑆) ≈ × × ∑𝐸𝑖 (𝑆𝑄𝑆)+ × × ∑𝐸𝑖 (𝑆𝑄𝑆) (9)\n𝑁 𝑎1−𝑎2 𝐶𝑀1 𝑀1 𝑎1−𝑎2 𝐶𝑀2 𝑀2\n𝑁 𝑁\n𝑖=1 𝑖=1\n𝐶𝑀1 𝐶𝑀2\n𝑁 𝑁\n𝑏0−𝑏2 1 𝑏1−𝑏0 1\n𝐸 (𝑆𝑄𝑆) ≈ × × ∑𝐸𝑖 (𝑆𝑄𝑆)+ × × ∑𝐸𝑖 (𝑆𝑄𝑆) (10)\n𝑁 𝑏1−𝑏2 𝐶𝑀1 𝑀1 𝑏1−𝑏2 𝐶𝑀2 𝑀2\n𝑁 𝑁\n𝑖=1 𝑖=1\nD) Comments on Writing and Clarity\nThe abstract should be rewritten for clarity and conciseness. Specifically:\n1. The physical model is not explicitly defined, and the fact that it is based on cluster expansion\nand lower dimensional subsystems is not mentioned.\nResponse: We would like to thank the reviewer for this suggestion. To clarify, our physical model\nis: the energy of MPEAs can be estimated by a linear combination of lower dimension MPEAs\n(basis) energy in a way that cancel-off the cluster contribution to total energy. The analytical\nvalidation of this model is done by cluster expansion formalism while the numerical validation has\nbeen done by fitting > 100,000 DFT data up to 10M.\nIn the revised manuscript, we have added the following sentences in the abstract on page 2, which\nis also shown below:\n“In this study, we propose that the synthesizability of MPEAs across diverse compositional and\nstructural spaces can be reliably predicted using a physical model that constructs the total energy\nof any MPEA as a linear combination of the energies of lower-dimensional subsystems.”\n2. The last sentence is unclear. The phrase “with much fewer DFT data than traditionally expected”\nis vague and lacks a clear definition or quantification. Additionally, “the unprecedented\ndiminishing of mixing enthalpy” should be rewritten for clarity.\nResponse: Thank you for the suggestion of further clarification. “with much fewer DFT data than\ntraditionally expected” refers to two aspects of data efficiency:\n1) Efficiency for predicting complete MPEA space: To give an example, with all 𝐶3 = 3,276\n28\nternary MPEAs as basis, we can predict any equal-molar MPEAs up to 10M with MAE being\naround 20 meV/atom, which covers a space of 24,817,650 (𝐶4 +𝐶5 +𝐶6 +𝐶7 +𝐶8 +𝐶9 +𝐶10)\n28 28 28 28 28 28 28\nequal molar MPEA compositions. Moreover, the MAE can be reduced to 10 meV/atom as long as\none uses a quaternary basis with 𝐶4 = 20,475 basis. More specifically, given that our work has\n28\nalready enumerated all ternary and quaternary dataset, we are already readily available to predict\nthe stability of any MPEA with BCC/FCC/HCP structures with reasonable efficiency and accuracy.\n2) Efficiency for predicting specific MPEA subspace: In case one needs to predict specific\nsubspace MPEA, our theoretical model suggests only lower-dimensional counter alloys are needed.\nFor example, to predict any refractory MPEA made from Ti, V, Cr, Zr, Nb, Mo, Hf, Ta, and W,\n11\nonly 𝐶3 = 84 calculations will be needed to establish the ternary MPEA basis to estimate\n9\nenergetics up to 9M. It is counterintuitive that a machine learning model or other complex model\ntrained with 84 data points can predict any refractory MPEAs. Therefore, we regard that our model\ncan greatly aid sample efficiency for predicting specific MPEA composition subspace\n“the unprecedented diminishing of mixing enthalpy” refers to the interpretation of the trend\nobserved in Fig. 4. The MAE in our model effectively represents a generalized mixing enthalpy of\nthe target MPEA formed from the basis MPEAs. The observed monotonic decrease and eventual\nconvergence of the MAE suggest that the mixing enthalpy of high-dimensional MPEAs gradually\ndeclines and approaches zero. This phenomenon fundamentally implies that the thermodynamic\nenergy landscape of MPEAs becomes increasingly flat as the compositional dimensionality\nincreases.\nTo further clarify the point, we have revised our statement in the abstract on page 2 of the revised\nmanuscript, which is also attached below:\n“Moreover, we reveal a universal trend that the stability of high-entropy alloys can be predicted\nusing a linear combination of the energies of lower-dimensional systems, such as ternaries or\nquaternaries, with low prediction errors. This observation further implies that the energy landscape\nbecomes increasingly flat in ultra-high entropy compositional spaces.”\nSpecific sentences:\n1. \"The universal rule ... by simply averaging the stability ...\" (line 45) \"Averaging the stability\"\nis not a well-defined term and should be rephrased more precisely.\nResponse: We appreciate the reviewer’s insightful comment. We agree that the phrase “averaging\nthe stability” is vague. To improve clarity, we have revised the sentence on page 3 to more precisely\ndescribe the approximation method used in the predictions of MPEAs:\n“It also supports the stability rule for MPEAs: the energy and hence stability of a high-dimensional\nMPEA can be reasonably approximated by a weighted combination of the energies of its lower-\ndimensional constituent subsystems.”\n2. \"We uncovered ... can be achieved with small datasets.\" (line 48) The term \"small datasets\"\nshould be clearly defined or quantified to provide meaningful context.\nResponse: We appreciate the reviewer’s suggestion to clarify the term “small datasets”. We have\nadded the following sentence in the revised manuscript on page 4 to clarify it.\n“Specifically, when focusing on a subset of MPEA including N metals, as few as 𝐶3 DFT\n𝑁\ncalculations are sufficient to predict the energy of all corresponding MPEAs.”\n3. \"Noble-metal-containing MPEAs ....\" (line 85) This sentence does not fit well as the opening\nsentence of the Experimental Validation of Predicted MPEAs subsection. The sudden shift in focus\nto noble-metal-containing MPEAs lacks prior context and should be introduced more smoothly.\nmajor reason for focusing on noble-metal contained MPEAs is due to 1) The experimental\nlimitation, as noble-metal-containing MPEAs are more readily accessible during SPBCL method\ndemonstrated, which is also detailed in response to B2 of the reviewer’s question for rationalizing\nthe choice of experimental compounds; 2) The noble-metal-containing MPEAs are less explored,\nas shown by Fig. 2(a), only 168 compositions have been reported, and 148 papers have been found\nin the past 25 years.\nTo improve the flow and provide a reasonable transition, we changed it into the following sentence\nto introduce the noble-metal containing MPEAs more smoothly on page 6, and also show below:\n“Among all MPEA compositions, noble-metal-containing MPEAs emerge as a distinct and\nunderexplored group. As shown in Fig. 2a, this compositional space is underexplored, with only\n168 compositions reported and 148 papers published (Fig. S4 and Tab. S1). Furthermore, the ease\nof reducing noble-metal precursors in a hydrogen atmosphere facilitates their synthesis using\ntechniques such as scanning probe block copolymer lithography. Notably, Noble-metal-containing\nMPEAs have also garnered significant attention due to their usefulness in catalysis and other\napplications. Therefore, our experimental validation efforts will primarily focus on noble-metal-\ncontaining MPEAs.”\n4. \"The value of C_N^M represents the total types of M-dimensional compositional spaces among\nthe N-dimensional system.\" (line 191) This is misleading; it is not the total types but rather the\nnumber of possible choices (i.e., the binomial coefficient). Please mention that C_N^M is a\nbinomial coefficient.\nResponse: We thank the reviewer for pointing out the ambiguity in our description. We have\nrevised the sentence to clarify that 𝐶𝑁 refers to the binomial coefficient, which denotes the number\n𝑀\nof possible combinations rather than distinct “types” of compositional spaces. The sentence has\nbeen changed in the revised manuscript on page 12, also as below:\n\"A convenient and physically meaningful choice is 𝜀𝑖 = 1 , where 𝐶𝑀 is the binomial coefficient\n𝐶𝑀 𝑁\n𝑁\ndenoting the number of possible M-component combinations from an N-component system.\"\n5. Line 34 in the supplementary notes: \"can be expressed in a spacious way.\" This should be\nreworded as \"can be expanded.\"\nResponse: We thank the reviewer for the suggestion. The statement has been rewritten to use “be\nexpanded” rather than “spacious” in page 3 of the supplementary notes, which is also shown below\nas:\n“We can then further write ∑ Δ𝑛 𝐽 such that the cluster concentration for both the self-\n𝑖𝑗 𝑖𝑗 𝑖𝑗\ninteraction term 𝐽 and cross-interaction term 𝐽 is explicitly expanded as below:”\n𝑖𝑖 𝑖𝑗\nE) Comments on Figures\n1. Fig. 1A and b presents widespread knowledge and may not be necessary. Consider removing\nthem. Instead be sure to refer to 1c in the text and discuss its significance.\na more meaningful periodic table style. In the revised Fig. 1a, all studied metals are highlighted\nwith colored triangles. Triangle in the bottom right represents the room temperature phase of\ncorresponding elemental metal, while noble metals are highlighted with red top-right triangle.\nFollowing the reviewer’s suggestion, we have removed the original Fig. 1b.\nFurthermore, appropriate revisions have been made to ensure that the original Fig. 1c (now Fig.\n1b) is clearly referenced and more discussions (structural dependent study and comparison\nbetween stability with descriptors, such as VEC, electronegativity and off-lattice displacement)\nare updated in the revised manuscript page 5 and also show below.\n“Both E and E present in Fig. 1b are calculated at 0 K, with the black line (E = E )\nform hull form hull\nindicating cases where decomposition leads to elemental phases. Structural-dependent plots are\nshown in Fig. S1, the energy distributions across FCC, BCC, and HCP structures are similar, and\na summary of the number of 4-metal MPEAs stabilized under different temperatures and phases is\nprovided in Fig. S2. To expand more on the stability, we have systematically plotted E and E\nform hull\nagainst other descriptors such as valence electron concentration (VEC), average electronegativity,\nand off-lattice displacement (see Fig. S3). These analyses show no clear systematic trends,\nalthough large off-lattice displacements generally correspond to instability.”\nFig. 1. Chemical space and data distribution of the MPEA dataset. a, The bottom left triangle represents\nthe room temperature phases as elemental metal, Mn and Sn marked as “other” because the room temperature phases\nare none of the three typical phases (BCC/FCC/HCP). The top-right triangle distinguishes each metal as either a noble\nor a non-noble metal; b, Scatter plot of formation energy (𝑬 ) and energy above hull (𝑬 ) for all 4M MPEAs\n𝐟𝐨𝐫𝐦 𝐡𝐮𝐥𝐥\nthat can be entropically stabilized below 2000 K (the stabilization temperature is calculated from 𝑬 /𝑺 , where\n𝐡𝐮𝐥𝐥 𝐢𝐝𝐞𝐚𝐥\n𝑺 the ideal mixing entropy). Blue and red represent low and high temperatures, respectively. The solid black line\n𝐢𝐝𝐞𝐚𝐥\nrepresents the line 𝑬 =𝑬 . Three dashed lines represent the maximum 𝑬 values of 4M MPEAs that can be\n𝐟𝐨𝐫𝐦 𝐡𝐮𝐥𝐥 𝐡𝐮𝐥𝐥\nstabilized at 500 K, 1000 K, and 2000 K, respectively, assuming ideal mixing entropy. Colormap represents a. VEC; b. Average electronegativity and c. average off-lattice\ndisplacement.\nFig. S1. Structural dependent scatter plot of formation energy (𝐸 ) and energy above hull (𝐸 ) for all a.\nform hull\nBCC 4M MPEAs; b. FCC 4M MPEAs; and c. HCP 4M MPEAs.\nand 2000 K.\nFig. S3. Scatter plot of formation energy (𝐸 ) and energy above hull (𝐸 ) for all 4M MPEAs that can\nform hull\nbe entropically stabilized below 2000 K (the stabilization temperature is calculated from 𝐸 /𝑆 , where 𝑆 the\nhull ideal ideal\nideal mixing entropy). The solid black line represents the line 𝐸 =𝐸 . Three dashed lines represent the\nform hull\nmaximum 𝐸 values of 4M MPEAs that can be stabilized at 500 K, 1000 K, and 2000 K, respectively, assuming\nhull\nideal mixing entropy. Colormap represents a. VEC; b. Average electronegativity and c. average off-lattice\ndisplacement.\n2. Fig. 2a, which summarizes published papers on noble-metal-containing MPEAs, is more suited\nfor a review paper rather than this study. It does not fit well within the scope of this work and\nshould be reconsidered.\nResponse: We would like to thank the reviewer for the suggestion. The intention of Fig. 2a is not\nto provide a comprehensive literature review, but rather to form a contrast of the number of new\nnoble-metal-based MPEAs discovered in this work versus what has been discovered in the past 25\nyears. It can help highlight the significant advancement of the data coverage on MPEAs. We have\nrevised the manuscript on pages 6-7 in order to avoid potential confusion. The revised contents are\nalso attached below:\n“Among all MPEA compositions, noble-metal-containing MPEAs emerge as a distinct and\nunderexplored group. As shown in Fig. 2a, this compositional space is underexplored, with only\n168 compositions reported and 148 papers published (Fig. S4 and Tab. S1).”\n15\n“Therefore, our experimental validation efforts will primarily focus on noble-metal-containing\nMPEAs. Among all the compounds screened in Tab. S1, equimolar compounds synthesized at\nhull\nout of 32) or below (30 out of 32) zero, which are listed in Tab. S2. This observation supports the\nuse of the E < S × T criterion as a reasonable approximation for thermodynamic stability. As\nhull ideal\na contrast to the less explored noble-metal-containing MPEAs, in this work, using high-throughput\nphase diagram calculations, we found 647 ternary, 2,178 quaternary, and 5,478 quinary noble-\nmetal-containing MPEAs that can be stabilized below 1273 K with ideal mixing entropy (Fig. 2b).”\n3. Fig. 3 is difficult to interpret, particularly due to the large number of elements. Additionally, in\nline 163, the authors state: \"We can also link the high frequency of specific elemental pairs, such\nas Pd-Pt, Rh-Pt, and Rh-Pd reported in the literature, to their high probability (thick segment, Fig.\n3b).\" However, using visual metrics like \"size\" and \"thickness\" makes direct comparisons\nchallenging. We are never told how the thicknesses were chosen. We suggest replacing Fig. 3 with\na clearer figure that directly reflects the key insights highlighted in the text, such as the most\ncompatible and most selective elements. The central figure 3a is very hard to interpret.\nThe remaining data could be moved to tables in the supplementary materials for better clarity.\nResponse: We thank the reviewer for this thoughtful feedback. In response, we have revised Fig.\n3 to improve clarity and interpretability. Specifically, we have updated the figure caption to clarify\nthat the thickness of the connection lines represents the relative probability of forming MPEAs\nbelow 1273 K. Therefore, the “thickness” is quantified by proportion to how many noble MPEA\nwith E < S * 1273 are there for given metal pairs. Due to the high-dimensional nature of\nhull ideal\nelemental compatibility, this chord diagram is regarded as the best way to illustrate the\ncompatibility when combining different metals as the audience can use thickness to directly\ncompare the relative stability between say M1-M2 pair in comparison to M3-M4.\nTo highlight such insights and avoid confusion, we have implemented the following changes:\nCaption in Fig. 3: “Thickness is set in proportional to the probability of finding MPEA with given\ntwo metals with stabilization temperatures below 1273 K.”\nIn the revised manuscript, page 9: “The thickness of each segment in the diagram is proportional\nto the occurrence (𝑛𝑀 ) of each metal (M) appearing in a stable MPEA that can be entropy\n𝑠𝑡𝑎𝑏𝑙𝑒\nstabilized below 1273K. The connections between segments indicate the occurrence (𝑛𝑀1−𝑀2)\n𝑠𝑡𝑎𝑏𝑙𝑒\nwhen both metals (M1 and M2) coexist in an MPEA that can be entropy stabilized below 1273 K.\nThe connections between non-noble metals are shown in Fig. 3c. Given the substantial\nenhancement in catalytic performance brought by the inclusion of noble metals, their impact is\nspecifically emphasized in Fig. 3(b, d–k).”\n4. Fig. 4a presents data up to an 8M MPEA, whereas Fig. 4b only includes data for 4M and 5M\nMPEAs. Would it be possible to add more data to Fig. 4b for consistency?\nResponse: We agree with the proposal that adding MPEAs data with more than five elements will\nbe helpful. This request closely aligns with Question B3 from the same reviewer, so we have\nprovided a similar response.\nof our model. The theoretical model is also updated to be up to 10 metals to be consistent (Fig. 4a\nin revised manuscript). Due to the exponential growth of the compositional space, the number of\nequimolar 6-element combinations from 28 elements reaches 376,740 (𝐶6 ), while the number for\n28\n10-element combinations increases to 13,123,110 (𝐶10). Recognizing that it is computationally\n28\nunfeasible to exhaustively evaluate all possible high-order combinations as we did for 3M and 4M\nMPEAs, we instead focused on a representative subset of high-dimensional MPEA compositions.\nSpecifically, we targeted FCC and BCC MPEAs composed of more commonly studied elements,\nand constructed the following three extended datasets:\n1) For BCC-MPEAs, we generated all possible combinations among nine typical elements\ncommonly found in refractory high-entropy alloys, e.g., Ti, V, Cr, Zr, Nb, Mo, Hf, Ta, and W. We\nextended the compositional dimensionality to 6M, 8M, and 9M, resulting in a total of 94\n(𝐶6+𝐶7+𝐶9) new BCC-MPEAs.\n9 9 9\n2) We considered 10 elements typically observed in FCC MPEAs: Ti, V, Cr, Mn, Fe, Co, Ni, Cu,\nZn, and Al. By generating all possible combinations in 6M, 8M, 9M, and 10M configurations, this\ndataset includes 266 (𝐶6 +𝐶8 +𝐶9 +𝐶10) additional FCC-MPEAs.\n10 10 10 10\n3) To strengthen the connection with both our experimental work (Fig. 2) and prior literature (Tab.\nS1), we created a noble-metal-enriched dataset. This includes all 6M, 8M, 9M, and 10M\ncombinations composed of at least one noble metal (Ru, Rh, Pd, Ag, Ir, Pt, Au) along with non-\nnoble metals (V, Cr, Mn, Fe, Co, Ni, Cu). This approach yielded 9,002 additional FCC-MPEAs.\nObservations:\nAs a result, we have added 9,362 new DFT-computed MPEAs, significantly expanding our dataset\nto include compositions from 6M up to 10M. This extended dataset enables a more robust\nvalidation of our predictive model and reinforces several key conclusions from our original\nmanuscript:\n1) Transferability across dimensions: Low-dimensional MPEAs can serve as an effective basis\nfor estimating the energies of higher-dimensional systems. For example, using 8M-MPEAs to\npredict 9M-MPEAs yields a minimum mean absolute error (MAE) as low as 6.956 meV/atom.\n2) Dimensional proximity improves accuracy: The prediction error generally decreases as the\ndimensional gap between the basis and target systems narrows. For instance, using 6M-MPEAs to\npredict 9M-MPEAs results in a lower MAE than using 4M-MPEAs as the basis.\nUpdates on comparison with CGCNN model: We have also trained additional CGCNN models\nusing the same training and testing datasets as our physical model in Fig. 4b. This guarantees that\nthe comparison in Fig. 4c reflects models evaluated under identical training and testing conditions.\nWe have updated Figs 4b-c to reflect the incorporation of this additional dataset. Corresponding\nrevisions have also been made to the manuscript text on pages 14–16 to clearly present these new\nfindings.\nPages 14-15: “To validate our model with higher-dimensional DFT-computed MPEAs, we\ngenerated an additional dataset comprising 9,362 DFT-calculated equimolar MPEAs ranging from\n6M to 10M MPEAs. This includes: (i) a BCC refractory dataset containing 6M, 8M, and 9M\nMPEAs composed of Ti, V, Cr, Zr, Nb, Mo, Hf, Ta, and W; (ii) two FCC datasets covering 6M,\n8M, 9M, and 10M MPEAs, one constructed from the compositional space of Ti, V, Cr, Mn, Fe, Co,\nIr, Pt, and Au. The fitting errors using lower-dimensional basis to predict on the higher-dimensional\ndata are shown in Fig. 4b. Our results confirm that low-dimensional MPEAs can reliably estimate\nthe energies of higher-dimensional systems.”\nPage 15: “Notably, prediction accuracy improves as the dimensionality of the basis approaches\nthat of the target system, and tends to converge at higher dimensions, with MAEs as low as\n6.956 meV/atom when using 8M MPEA basis to predict 9M compositions.”\nPages 15-16: “By fitting both models with the same dataset, our physical model shows similar,\nand in most cases, better performance compared to CGCNN (Fig. 4c), and MAEs are shown in\nTab. S5. For target dimensionalities greater than five, the training data for CGCNN consists of all\nlower-dimensional compositions sampled within the same expanded compositional space, rather\nthan the entire dataset of lower-dimensional systems.”\nFig. 4. Performance of our predictive model. a, The theoretically derived MAE of predicting energy\ncontribution from selective cluster 4M-10M MPEA using 2M-9M MPEA as basis. J is used to represent the unit of\n0\nenergy, with the physical meaning being the energy contribution of certain cluster, e.g., 𝑱 , 𝑱 , 𝑱 respectively. b,\n𝒊𝒊 𝒊𝒋 𝒊𝒋𝒌\nThe MAE for using DFT-computed total energy of 2M-9M MPEA to predict 4M-10M MPEAs as a contrast to the\ntheoretical derivation. c, The comparison between our physical model with a crystal graph convolutional neural\nnetwork (CGCNN).\n5. Fig. S1 is confusing and unclear. How can the diagonal entries be understood? i.e. what does\nthe pair Ru-Ru mean in a binary?\ndiagonal originally indicated MPEAs that contain only a single noble metal element. The bottom\nleft triangle indicates MPEAs that contain two noble metal elements, while the top right triangle\nindicates MPEAs that contain three noble metal elements. We understand that it may be confusing\nin the previous figure that the diagonal data was displayed using the colormap corresponding to\nMPEAs containing two noble metals. To address this, we have removed the entries containing only\none noble metal from this figure. Additionally, the figure numbering has been reorganized to match\nthe order of appearance in the manuscript, and the updated figure is now labeled as Fig. S4, which\nis also shown as below.\nFig. S4. Pair-wise analysis of the frequency of reported MPEAs containing noble metals, based on data from\nTab. S1. Blue bars represent alloys containing exactly two noble metals, while purple bars correspond to three or more\nnoble metals.\nReviewer #2 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We sincerely thank both reviewers for their 24 insightful questions, organized into five\nsections (A-E sections). We believe our manuscript has been significantly improved after\naddressing all their comments and making the corresponding revisions.\nReviewer #3 (Remarks to the Author):\nAuthors used two stability criteria to predict MPEA synthesizabilty, Ehull and Eform. The first\nwas obtained by constructing a phase diagram from already available structure databases\n(pymatgen and OQMD) and an in-house high-throughput DFT energy calculation. The energy\nabove hull is used as a criterion for entropic stabilization by assuming an ideal entropy contribution\nto the free energy. This stability criterion is then tested with experimental synthesis of 9 equimolar\nMPEAs.\nAs a shortcut for total energy calculation without the need of DFT calculations a linear model\nbased on Cluster Expansion is proposed and trained on the previously generated DFT dataset. The\nmodel is based on extrapolation from lower to higher dimensionality MPEAs by approximating\nthe total energy to weighted contributions of lower dimensionality MPEAs. A mathematical\nanalysis is carried out in the supplementary material to theoretically assess the error accompanied\nby the method. The authors state that theoretically the error should be a negligible fraction of the\npair, three-body, and the rest of cluster interactions.\nThe contribution from the authors to MPEAs design is substantial: a quick estimation of stability\nthrough a linear approximation of the free energy of MPEAs. The work is comprised by two main\nparts: stability rule and DFT total energy model.\nI would advise the following major revisions for the manuscript to be ready for publication.\nResponse: We would like to thank the reviewer for their recognition of novelty of our work. We\nhave addressed all questions and suggestions and provided a point-by-point response below.\n◦ While the Hume-Rothery rules are mentioned, a proper comparison with previously stablished\nrules and composition trends is not fleshed out. Specially how this new stability criteria squares\nagainst previously shown soft rules in HEA design.\nResponse: We thank the reviewer for the valuable suggestion. In the revised manuscript, we have\nclarified the relationship between our proposed physical linear model and traditional soft rules\nsuch as the Hume-Rothery guidelines. Specifically, we now include a comparison between\nformation energy and common descriptors used in Hume-Rothery-type rules, i.e., average radius\nof metals (𝑟 ), standard deviation of metal radius (𝑟 ), average electronegativity (𝜒 ) and\n𝑎𝑣𝑒 𝑠𝑡𝑑 𝑎𝑣𝑒\nstandard deviation of electronegativity (𝜒 ). As shown in the new analysis (Fig. S8), the highest\n𝑠𝑡𝑑\nPearson correlation coefficient (|PCC|) observed between formation energy and these parameters\nis only 0.380, which is 𝐸 versus 𝜒 . This weak correlation underscores the limitations of\n𝑓𝑜𝑟𝑚 𝑠𝑡𝑑\nHume-Rothery-type descriptors in quantitatively capturing the stability of multi-principal element\nalloys (MPEAs), which further highlights the effectiveness of our model.\nWe add the figure in supplementary information as Fig. S8 and the following statement in the\nrevised manuscript on page 15, and show below:\n“When comparing with typical Hume-Rothery descriptors, such as the average and standard\ndeviation of atomic radius and electronegativity, the relationship with formation energy is weak,\nwith the maximum absolute Pearson correlation coefficient (|PCC|) being only 0.380, as shown in\nFig. S8. This further underscores the strength of our physical model in providing both quantitative\naccuracy and physical interpretability.”\nFig. S8. The relations and Pearson correlation coefficients between E and a. average radius, b. standard\nform\ndeviation of radius, c. average electronegativity, and d. standard deviation of electronegativity.\n20\n◦ The claim is bold but only one scatter plot is given with little information on what is being shown\nhere. The x-axis and color-scale communicate both the same information. Good properties to\ncompare could be averaged mixing enthalpy, ideal mixing entropy, VEC, electronegativity, lattice\nmisfit, etc.\nResponse: We thank the reviewer for suggesting improvements to our figure. Based on the context,\nwe believe the reviewer is referring to Fig. 1c (now updated as Fig. 1b). Following the suggestion,\nwe explored using VEC, average electronegativity, and average off-lattice displacement as\nalternative color codes. The corresponding visualizations are presented in Fig. S3. As shown in\nFig. S3 as well as below, no clear correlation can be observed, which further underscores the\ncomplexity of stability in high-entropy systems. Therefore, we decided to retain the original plot.\nWe note that while the color coding in the original figure somewhat overlaps with the physical\nmeaning of the E values, it provides a more intuitive representation of the corresponding lowest\nhull\nentropy stabilization temperature of certain E values. To accommodate potential interest in VEC,\nhull\nelectronegativity, and off-lattice displacement distributions, we have also included additional\nstatements in the text to guide readers to Fig. S3.\nThe revision has been implemented in the revised manuscript on page 5-6, which is also attached\nbelow:\n“To examine the correlation between the two stability descriptors and commonly used physical\nparameters, we adopted the same plotting style as in Fig. S3, coloring the data points by valence\nelectron concentration (VEC), average electronegativity ( 𝜒 ), and average off-lattice\n𝑎𝑣𝑒\ndisplacement (𝜀 ) (see Fig. S3). These analyses reveal no clear or systematic trends, although\n𝑎𝑣𝑒\ncompositions with large off-lattice displacements generally correspond to instability. These\nobservations further highlight the complexity of stability prediction in high-entropy systems,\nemphasizing the need for models that are both predictive and physically interpretable.”\nFig. S3. Scatter plot of formation energy (𝐸 ) and energy above hull (𝐸 ) for all 4M MPEAs that can\nform hull\nbe entropically stabilized below 2000 K (the stabilization temperature is calculated from 𝐸 /𝑆 , where 𝑆 the\nhull ideal ideal\nideal mixing entropy). The solid black line represents the line 𝐸 =𝐸 . Three dashed lines represent the\nform hull\nmaximum 𝐸 values of 4M MPEAs that can be stabilized at 500 K, 1000 K, and 2000 K, respectively, assuming\nhull\nideal mixing entropy. Colormap represents a. VEC; b. Average electronegativity and c. average off-lattice\ndisplacement.\n◦ Figure 1c. Specify the energies are calculated at zero temperature (zero Kelvin). It is\nrecommended to provide subplots by crystal structures: FCC, BCC, and HCP. what is the\nsignificance of the black solid line (E_form = E_hull)? Phase stability is all about comparing the\nenergies among different crystal structures. What matters the most here is the E_hull that is\ncalculated based on E_form. I expect more meaningful plots of E_hull vs crystal structure vs\n21\nnumber of component vs critical temperature above which entropy would favor forming solid\nsolution (T * S >= E_hull).\nResponse: We have clarified in the caption that the formation energies in Fig. 1c (now Fig. 1b in\nthe revised manuscript) are calculated at 0 K. The black solid line (E = E ) represents cases\nform hull\nwhere the decomposition products are elemental phases. To provide additional insight, we have\nadded subplots in Fig. S1 showing the distribution of E and E for different crystal structures\nform hull\n(FCC, BCC, HCP). The energy distributions across these three phases do not show significant\ndifferences. Additionally, we include a supplementary table summarizing the number of 4-metal\nMPEAs stabilized under different temperatures and crystal structures to aid interpretation of\nentropy-stabilized phase formation, which is shown in Fig. S2. The following descriptions have\nbeen added to the revised manuscript on page 5.\nPage 5: “Both E and E present in Fig. 1b are calculated at 0 K, with the black line (E =\nform hull form\nE ) indicating cases where decomposition leads to elemental phases. Structural-dependent plots\nhull\nare shown in Fig. S1, the energy distributions across FCC, BCC, and HCP structures are similar,\nand a summary of the number of 4-metal MPEAs stabilized under different temperatures and\nphases is provided in Fig. S2.”\nFig. S1. Structural dependent scatter plot of formation energy (𝐸 ) and energy above hull (𝐸 ) for all a.\nform hull\nBCC 4M MPEAs; b. FCC 4M MPEAs; and c. HCP 4M MPEAs.\nFig. S2. The number of stable MPEA compositions with BCC, FCC, and HCP structures at 500 K, 1000 K,\nand 2000 K.\n◦ Mind that you could use the simple expression: T * S >= E_hull to estimate the entropy stabilizing\neffect, however, T * S >= E_hull does not always guarantee forming solid solution at temperatures\nabove the critical temperature. This is because temperature brings in the contribution of vibrational\n22\nfree energy, and competing compound phases could exhibit compositional homogeneity range so\nthat configurational entropy will play as well for compounds.\nResponse: We appreciate the reviewer’s comment. We agree that vibrational entropy plays an\nimportant role in determining stabilization temperature and, as such, the temperature we report,\nbased on the E < S × T, is intended only as a temperature reference for experiments. Even\nhull ideal\nthough it should be acknowledged that the inclusion of vibrational free energy would be useful.\nHowever, computational estimation of vibrational entropy will be very expensive and currently\nimpractical, especially for estimating vibrational entropy over 100,000 materials. We thus propose\nto exclude the consideration for this work, while a systematic and independent work would be\nneeded to tackle this challenge.\nIn our revised manuscript, we have clarified the limitation of our stability descriptor and added the\nfollowing sentences on page 6, which are also attached below:\n“Although the criterion E < S × T is a rough approximation, it has been widely adopted as a\nhull ideal\npractical indicator for predicting stability in compositionally complex and high-entropy materials.\nIncorporating vibrational free energy and chemical short-range order (CSRO) would certainly\nenhance accuracy, but doing so remains computationally unfeasible for large-scale datasets,\nespecially when considering hundreds of thousands of compositions across a 28-element space.”\n◦ The authors wrote: “In our newly established dataset of 126,429 DFT-computed MPEAs, 8,303\nMPEAs (6.56% 127 of all compounds) with three to five elements are predicted to be entropy-\nstabilizable below 1273 K.” Do you really mean above 1273 K?\nResponse: Thank you for the comment. We did mean below 1273 K, as these MPEAs are predicted\nto have stabilization temperatures lower than this threshold. To clarify, we have revised the\nsentence on page 9:\n“In our newly established dataset of 126,429 DFT-computed MPEAs up to 5 metals, 8,303 MPEAs\n(6.56% of all compounds) with three to five elements are predicted to have stabilization\ntemperatures below 1273 K.”\n◦ The authors wrote: “The accuracy of our physical model is comparable to, or better than, that of\nstate-of-the-art deep learning models trained on the same dataset, while retaining full\ninterpretability.” This appears to be overstatement without providing evidence to support their\nclaim. The authors wrote: “Such performance matches the state-of-the-art graph neural network\n[43] models trained on the same dataset [28].” However, reference [28] does not address high\nentropy alloys nor on the same database.\nResponse: We would like to thank the reviewer for the careful reading and thoughtful question.\nWe carefully check and realize that [43] in this comment may refer to the line number in the\ncorresponding page of the manuscript, not a literature citation. Reference [28] Phys. Rev. Lett. 120,\n145301, reports CGCNN, the model we use to compare with our physical model. We have done\nthe following revisions and clarifications to avoid potential confusion:\ntrained on the same dataset”, we will need to show a) The graph neural network we used, i.e.,\nCGCNN is a representative state-of-the-art graph neural network; b) We have made fair\ncomparisons with state-of-the-art graph neural network.\nTo clarify that CGCNN is a strong representative of state-of-the-art graph neural networks for\nstudying MPEAs, we have cited the original CGCNN paper as reference 30 in the revised\nmanuscript. Additionally, we have included a citation to our recent benchmarking study\n(ChemRxiv, 2024; doi:10.26434/chemrxiv-2024-fmq8p), which demonstrates that CGCNN\nachieves the best performance among several widely used graph neural network (GNN) models\nfor MPEA prediction. In that study, we systematically evaluated the performance of three GNN\nmodels with increasing architectural complexity: CGCNN, M3GNet (which incorporates three-\nbody interactions), and CHGNet (which includes both three-body interactions and magnetic\nmoments). All models were trained and tested on the same MPEA dataset, which contains\ncompositions with up to four metal elements, across multiple prediction tasks. The results showed\ncomparable performance among the models, with CGCNN outperforming M3GNet and CHGNet\nin some cases. These findings support the use of CGCNN as a representative state-of-the-art GNN\nfor MPEA modeling. This benchmarking study has been cited as reference 31 in the revised\nmanuscript.\nTo clarify that we have made fair comparison between CGCNN and our physical model, Fig. 4c\nshows the comparison between our physical model and CGCNN using the same training data to\npredict the same dataset of specific MPEAs. We observe that our physical model performs\ncomparably to, and in most cases outperforms, CGCNN. This serves as strong evidence for the\ngood performance of our physical linear model.\nThe citation for CGCNN is correctly included after “graph neural network model”, and the\nsentence has been updated on page 3 and also shown as below:\n“Such performance matches the state-of-the-art graph neural network model [30] trained on the\nsame dataset [31].”\n◦ Synthesizability of high entropy materials is reported and modeled in the literature, and citing\nthem is very appropriate, to list a few: (1) Sarker et al. High-entropy high-hardness metal carbides\ndiscovered by entropy descriptors. Nat Commun 9, 4980 (2018), (2) Eckert et al. \"A formula to\npredict the synthesizability of high-entropy materials.\" NATURE (2024).\nResponse: We thank the reviewer for providing two prior works on the synthesizability of high-\nentropy materials. Both are great papers on high entropy materials and are highly relevant. We\nhave now cited the suggested references in the revised manuscript as references 28 and 29:\n28. Sarker, P. et al. High-entropy high-hardness metal carbides discovered by entropy descriptors.\nNature communications 9, 4980 (2018).\n29. Eckert, H. & Curtarolo, S. A formula to predict the synthesizability of high-entropy materials.\nNature (2024).\n◦ For the following model to be based on the cluster expansion method, the configurational entropy\nis glossed over in the stability rule definition. Previous works have shown entropy is not trivial at\n24\ntemperatures discussed here. Especially, when referring to entropic stability as low as in 500K.\nI think a distinction must be made for the model to have a comfortable application temperature\nthreshold and a discussion on the entropy exceptions could be added.\nResponse: We appreciate this insightful comment. We agree that configurational entropy can\ncomplicate the prediction of stability, which is also suggested by DOI: 10.1038/s41467-021-\n25979-5 (now cited as Ref. 77). In terms of the concern of significant deviation when comes to\ntemperature below 500K, we want to clarify that the experimental session as well as well as major\nfocus is on synthesis temperature around 1000℃ or above (1273K), which is far away from the\ntemperature that the review worries about. Moreover, we have added the following clarifications\nin our revised manuscript on page 15 to highlight the potential challenge of vibrational entropy:\n“It should be noted that the model is constructed from zero-temperature total energies; it can serve\nas a baseline for estimating entropy-stabilization temperatures, but more accurate estimation of all\ntypes of entropies will greatly enhance the accuracy of the stability prediction.”\n◦ The high-throughput effort is impressive here, and experimental validation is added and\nmeaningful to the material. But the scope of the project is too big to be assessed by nine MPEAs.\nI recommend a comparison with previously reported literature (a MPEA experimental database).\nResponse: We agree that the space of MPEAs is vast, and that synthesizing only nine compositions\ncannot fully cover the entire design space. Following the reviewer’s suggestion, we have now\nfurther investigated compounds reported in our experimental literature databases, which is\nsummarized Tab. S1. We evaluated the E values at 1273 K (E -1273 × S ) for\nhull hull Ideal\nexperimentally reported equimolar compounds synthesized below 1273 K, and found that most\nexhibit E values around (2 out of 32) or below 0 meV/atom at 1273 K (30 out of 32). This\nhull\nobservation further supports the predictive capability of our models and illustrates how our work\ncomplements and extends the current experimental MPEA dataset.\nThe following statements are added in the revised manuscript on pages 6-7 to provide additional\nvalidation on the previously reported synthesized equal-molar compounds:\n“Among all the compounds screened in Tab. S1, equimolar compounds synthesized at 1273 K were\nexamined, and their corresponding E values at this temperature are all close to (2 out of 32) or\nhull\nbelow (30 out of 32) zero, which are listed in Tab. S2. This observation supports the use of the\nE < S × T criterion as a reasonable approximation for thermodynamic stability.”\nhull ideal\n◦ An analysis on non-equimolar compositions could enhance the model extrapolation capabilities.\nTo the best of my understanding the model is a truncated cluster expansion (only point clusters)\ntrained on only M-sized equimolar subsystems. The point cluster interaction should still be viable\nto use for non-equimolar systems but at a higher error cost. Given the compositional design having\nan important role in MPEA design I believe this may be worth discussing.\nResponse: We appreciate the reviewer’s insightful suggestion regarding the potential to enhance\nthe model’s extrapolation capability by incorporating non-equimolar compositions. In the current\n25\nstudy, we focus exclusively on equimolar compositions, which lie at the center of the\ncompositional space and exhibit the strongest high-entropy effect. Non-equimolar compositions\nwere not included in our dataset primarily due to the vastness of their compositional space, which\nlacks the constraint of equal atomic fractions and is therefore significantly larger and more\ncomplex to sample systematically. Nevertheless, our physical model is, in principle, applicable to\nnon-equimolar compositions. A more thorough exploration in that direction would require\nsubstantial effort in both data generation and model evaluation, which we believe is valuable and\ndeserving of a dedicated future study.\nThe following clarification has been added to the revised manuscript on page 17.\n“In the current validation dataset, we include exclusively equimolar MPEA compositions, which\nlie at the center of the compositional space and exhibit the strongest high-entropy effect. Non-\nequimolar compositions were not included in our dataset primarily due to the vastness of their\ncompositional space, which lacks the constraint of equal atomic fractions and is therefore\nsignificantly larger and more complex to sample systematically. Nevertheless, our physical model\nis, in principle, applicable to non-equimolar compositions, but further exploration of the linear\ncombination coefficients will be needed.”\n◦ I highly encourage the authors to share all the compositions vs energies data to benefit the entire\ncommunity.\nResponse: We would like to thank the reviewer for the suggestion. All data that enables the\nproduction of all figures in this work is provided. Meanwhile, we are in the process of open access\nto all raw DFT data, eventually in platforms like Materials Project. It may require some effort,\ngiven the large data size and time required for cleaning and transferring data. To clarify the\nsituation, we have also added in the data availability that early access to the raw data can be\nprovided to the corresponding author before the official data is made public in an open database\nsuch as Materials Project. The change has been made on page 28 of data and material availability,\nwhich is also attached below:\n“All data for reproducing this work have been included in manuscript and supplementary\ninformation. Raw DFT data will be made available upon reasonable request to the corresponding\nauthor. Moreover, the raw DFT data will be made available in public repository after processing\nand proper license.”\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We would like to thank the constructive questions from both reviewer 3 and 4, it is\nbelieved that our paper has become much stronger and clearer after implementing all these\nsuggestions.\nReviewer #5 (Remarks to the Author):\nThey propose a physical model to predict MPEA stability, achieving an accuracy of 13.38\nmeV/atom. The model has been validated by fabricating nine novel MPEAs based on its\npredictions. The results obtained are potentially significant and interesting for the materials science\ncommunity. However, there are two serious concerns that prevent me from recommending the\narticle for publication in its current form.\nResponse: We would like to thank the reviewer for their high recognition of our work. We have\naddressed all questions and suggestions and provided a point-by-point response below.\n1. The authors compare their results with those obtained using a state-of-the-art data-driven\nCGCNN model. While this comparison suggests the reliability of their proposed model, it raises\nthe question of whether the computational resources required for their approach are justified, given\nthat its accuracy is only slightly better than that of an existing model. The authors should discuss\nthis point in more detail, demonstrating whether their approach offers superior accuracy,\nperformance, or transferability compared to other models like CGCNN.\nResponse: We thank the reviewer for this thoughtful comment. The efficiency advantages of our\nmodel compared to deep learning approaches such as CGCNN primarily stem from the following\naspects:\n1. The \"training\" process for our model is straightforward and requires minimal computational\nresources, neither GPUs nor high-performance CPU are necessary. All that is needed is a dataset\nto serve as the basis, which can be as simple as equimolar ternary MPEAs. The “training process”\nessentially involves solving a set of linear equations to determine the coefficients, as described in\nSupplementary Notes 1 and 2. The easiest “training” can be setting 𝜀𝑖 = 1 as shown in Eq. 2 of\n𝐶𝑀\n𝑁\nsupplementary information\n2. Once trained, our model enables much faster energy and stability predictions compared to GNN-\nbased models. Specifically, it only requires a weighted linear average of known lower-entropy\nMPEAs, bypassing the need for structural relaxation using machine-learned force fields (MLFFs),\nwhich is typically required for GNN-based predictions. The completeness of our constructed\ndatabase shown in Fig. 1 and the reasonable prediction accuracy shown in Fig. 4 demonstrate that\nlinear averaging is already capable of producing reliable stability estimates for any MPEAs.\n3. For anyone who wants to fit our model for exploring a certain subspace of MPEA with their\nown DFT dataset, it can also be done very efficiently. To give an example, for predicting any\nrefractory MPEA made from Ti, V, Cr, Zr, Nb, Mo, Hf, Ta, and W, only 𝐶3 = 84 calculations will\n9\nbe needed to establish the ternary MPEA basis to estimate energetics up to 9M. It is\ncounterintuitive that a GNN trained with 84 data points can be reliable on predicting arbitrary\nrefractory MPEAs. Therefore, we regard that our model can greatly aid sample efficiency for\npredicting specific MPEA composition subspace\nTo further clarify the efficiency over deep learning models, we have added the following\nstatements in the revised manuscript on pages 16 and 18.\nPage 16 “More importantly, our approach requires no training on GPUs or high-performance CPUs\nas all we need to do is to determine coefficient 𝜀𝑖 once the energy basis is computed, while the\neasiest yet effective 𝜀𝑖 is 1 .”\n𝐶𝑀\n𝑁\nenergy and hence stability of high-dimensional MPEAs using weighted linear average of the lower-\ndimensional MPEA as basis. Fitting such model eliminates the need for training at GPU or high-\nperformance CPUs. As determination of coefficient is through solving a set of linear equations or\nbe as simple as setting weighting factor 𝜀𝑖 = 1 . For specific subspace of MPEA of interest, there\n𝐶𝑀\n𝑁\nis no need to develop a comprehensive training data but simply lower entropy equimolar basis at\ncertain dimensional, e.g., ternary or quaternary MPEAs. Moreover, the application of our model\ndoes not need further relaxation using Hamiltonian such as machine learning force field, which are\nrequired by the state-of-the-art foundational models. We thus believe this discovery will facilitate\nthe rapid exploration of novel MPEA space.”\n2. The proposed approach relies on a large dataset of 126,429 DFT-computed MPEAs, implying\nsignificant computational resources were expended to generate this data. there is a state-of-the art\napproach to perform any atomistic simulations (particularly the same performed be the authors)\nwhich is based on Universal Interatomic Potentials (or Foundation Models). These models, based\non Graph Neural Networks and trained on extensive ab initio datasets, demonstrate high accuracy\nand transferability. The best-performing models achieve accuracy comparable to that reported by\nthe authors (see for reference). Moreover, their computational efficiency is orders of magnitude\nhigher than that of ab initio calculations. Thus, using such models might yield similar accuracy to\nthe one presented in the manuscript, which could significantly diminish the work's novelty and\nimpact. The authors should address this concern and discuss their findings in this context.\nResponse: We appreciate the reviewer’s insightful comment and the opportunity to clarify the\nnovelty of our work, considering recent advances in universal interatomic potentials and\nfoundation models. We fully agree that these graph neural network (GNN)-based models, trained\non extensive ab initio datasets, represent a promising direction for accelerating atomistic\nsimulations and have demonstrated excellent accuracy and transferability in many cases.\nHowever, our approach differs fundamentally in both motivation and methodology. First, the core\ninnovation of our work lies in the discovery that the stability of high-dimensional MPEAs can be\naccurately predicted by leveraging a low-dimensional energy basis, specifically, from a relatively\nsmall set of DFT-computed ternary or (N-1) dimensional MPEAs. Once this basis is computed, no\nfurther training or iterative model optimization is needed. Predictions are made via simple,\ninterpretable operations (e.g., weighted averaging over basis energies), which offer significant\ntransparency and physical insight into how stability arises in complex alloys.\nMoreover, while foundation models offer dramatic speedups over DFT once trained, they still rely\nfinetune on additional data for specific subspace. In contrast, our model based on the\ncomprehensive dataset demonstrated in Fig. 1 is readily available for predicting any MPEAs. The\ninitial cost of generating the 126,429 MPEAs in our dataset was a one-time investment and is\nsubstantially amortized by the ability to predict the stability of tens of millions of higher-order\ncompositions without further DFT or ML model retraining. We have added a brief discussion in\nthe revised manuscript to clarify these distinctions and contextualize our contribution relative to\nmodern interatomic potential-based approaches.\nmodel once we have DFT data for MPEA systems of interest, which are detailed in response to the\nfirst question. Therefore, we believe such framework is not only accurate, interpretable, but also\nmuch easier to establish, rebuild, or apply.\nThe following statements have been added to the revised manuscript on page 18:\n“Building on this thermodynamic interpretation, our physical linear model predicts the energy and\nhence stability of high-dimensional MPEAs using weighted linear average of the lower-\ndimensional MPEA as basis. Fitting such model eliminates the need for training at GPU or high-\nperformance CPUs. As determination of coefficient is through solving a set of linear equations or\nbe as simple as setting weighting factor 𝜀𝑖 = 1 . For specific subspace of MPEA of interest, there\n𝐶𝑀\n𝑁\nis no need to develop a comprehensive training data but simply lower entropy equimolar basis at\ncertain dimensional, e.g., ternary or quaternary MPEAs. Moreover, the application of our model\ndoes not need further relaxation using Hamiltonian such as machine learning force field, which are\nrequired by the state-of-the-art foundational models. We thus believe this discovery will facilitate\nthe rapid exploration of novel MPEA space.”\nUniversal Framework for Efficient Estimation of Stability in Multi-\nPrincipal Element Alloys\nLin Wang1†, Bo Shen2,3†, Zheng-Da He1, Zihao Ye2,3, Yan Zeng1, Chad Mirkin2,3, and Bin\nOuyang1,*\n1Department of Chemistry and Biochemistry, Florida State University, Tallahassee, FL\n32304, USA.\n3International Institute for Nanotechnology, Northwestern University, Evanston, Illinois\n60208, USA.\n†These authors contributed equally to this work\n*To whom correspondence may be addressed. E-mail: bouyang@fsu.edu\nReviewer #1 (Remarks to the Author):\nWe thank the authors for their detailed response and acknowledge their attempt to substantially\nimprove the manuscript. Several major concerns still remain and should be carefully addressed\nand clarified in the text where appropriate.\nResponse: We would like to thank the reviewer for their positive comments and constructive\nadvice. In this revision, we have addressed all questions and suggestions and provided a point-by-\npoint response below. After addressing all 19 additional questions summarized with 16 bullet\npoints, we believe that the revised manuscript has been further improved in terms of clarity and\nreliability.\n1. This work aims to provide a framework to estimate the synthesizability of multi-principal\nelement alloys (MPEAs). Therefore, the central audience of the paper is a researcher asking: “Is\nthe MPEA X synthesizable or not?” However, the current layout and organization of the paper\nhinder the reader from clearly answering this question. For instance, the key results—Eqs. (1) and\n(2), which describe the model—appear late in the Results section, following a lengthy discussion\nof computational and experimental validation. The model should instead be presented earlier,\nbefore diving into its validation.\nResponse: We would like to thank the reviewer for the suggestion. We would like to clarify the\nlogic flow of this paper as follows:\n1. This work began with the development of a systematic dataset for stability assessment\nbased on DFT-calculated E values (Fig. 1), followed by experimental validation\nhull\ndemonstrating that E is a reasonable stability metric (Fig. 2; Fig. S5).\nhull\n2. With the reliability of E confirmed, we further established an understanding of elemental\nhull\ncompatibility using the developed dataset, as illustrated in Fig. 3.\n3. Building on the reliability demonstrated in Figs. 1–2 and the insights from Fig. 3, we\ndeveloped a universal model that can extend the understanding of MPEA stability to\nbroader compositional spaces (Fig. 4).\n4. We agree with the reviewer that the model itself represents a key novelty. However, the\nexperimental results are primarily used to demonstrate that DFT-calculated E values can\nhull\nreliably predict synthesizability, thereby showing that a model trained and validated on\nDFT-calculated energies can effectively guide synthesis. The experimental results do not\ndirectly verify the model. Consequently, if we were to begin with the universal model\nequations, it would be difficult to incorporate the experimental results into the workflow,\nand the overall narrative would be less effective in emphasizing other significances of this\nwork.\nWith this clarification, we hope to help the reviewer gain a better understanding of our work. At\nthe same time, we recognize that there may have been some confusion regarding the role of the\nglobal model presented in Fig. 4 and the corresponding equations. In response to questions in\nbullets 4–6, we have provided detailed explanations to clarify our model.\n2\n2. We strongly suggest guiding the reader through a concrete example—such as the successfully\nsynthesized FeNiCuPd—to show how the model is applied in practice. This should include\nnumerical values, e.g., the energies of relevant subspaces. A visual summary, such as a flowchart,\nwould also improve accessibility and help readers apply the framework to new systems.\nResponse: We thank the reviewer for this question. As discussed in our response to Question 1,\nthe experimental results are intended solely to demonstrate that the DFT-computed E serves as\nhull\na reliable metric for estimating experimental synthesizability, rather than being directly connected\nto the global model. The universal model is designed to predict a broader range of compositions\nbeyond those shown in Fig. 4 and is not specifically tied to FeNiCuPd.\nTo address the request regarding the evaluation of FeNiCuPd synthesizability: this can be\nevaluated by constructing a convex hull using DFT-calculated energies of all reported metals and\nintermetallics to identify the ground states. For the Fe–Ni–Cu–Pd system, the convex hull includes\n16 ground states (Cu, Cu Pd, CuPd , CuPd, Fe, Fe NiPd , FeCu Pd , FeNi , FeNi , FeNi, FePd ,\n3 3 2 6 4 5 3 8 2\nFePd , FePd , FePd , Ni, Pd). Our convex-hull analysis shows that FeNiCuPd lies ~99.59\n3 5 8\nmeV/atom above the ground state, competing primarily with FeCu Pd , FeNi, Cu Pd, and FeNi .\n4 5 3 3\nConsidering the ideal mixing entropy of FeNiCuPd (0.119 meV/atom/K), the entropy stabilization\ntemperature is estimated to be 836.89 K, which is below 1273 K. Therefore, we identify FeNiCuPd\nas synthetically accessible\nOur global model further enables estimation of the energies of additional MPEAs without the need\nfor new DFT calculations, providing an efficient framework for assessing synthesizability through\nE and phase diagram analysis.\nhull\nThe following discussion is updated in the revised Supplementary Note 1 to guide the audience on\nhow synthesizability is estimated:\n“To further clarify how the thermodynamic stability criterion is applied, we present an example\nfor a quaternary alloy FeNiCuPd. The convex hull was constructed using DFT energies of all\nrelevant reported metals and intermetallic ground states. For FeNiCuPd, the ground states on the\nconvex hull include Cu, Cu Pd, CuPd , CuPd, Fe, Fe NiPd , FeCu Pd , FeNi , FeNi , FeNi, FePd ,\n3 3 2 6 4 5 3 8 2\nFePd , FePd , FePd , Ni, and Pd. Among these, the competing phases identified for FeNiCuPd are\n3 5 8\nFeCu Pd , FeNi, Cu Pd, and FeNi . The calculated energy above hull (E ) is 99.59 meV/atom.\n4 5 3 3 hull\nConsidering the ideal configurational entropy of FeNiCuPd (S = 0.119 meV/(atom·K)), the\nideal\ncorresponding entropy stabilization temperature is 836.9 K. Since this value is lower than the\ntypical synthesis temperature of 1273 K, FeNiCuPd can be considered synthetically accessible at\n1273K. This example illustrates how the E < S × T criterion is practically employed to assess\nhull ideal\nthe thermodynamic feasibility of MPEAs.”\nWe also refer to Supplementary Note 1 on page 7 of revised manuscript as:\nExample of how synthesizability is evaluated can be found in Supplementary Note 1.\n3. To assess the robustness of the model, it is also important to consider negative cases. We\nrecommend evaluating noble-metal-containing alloys that the model predicts to lie near the\nsynthesizability boundary and attempting their synthesis, or alternatively, applying the model to\n3\nknown unsynthesizable MPEAs reported in the literature. These cases would provide a more\nbalanced validation of the model.\nResponse: We appreciate the reviewer’s suggestion to include negative cases as part of our\nvalidation. In the revised manuscript, we address both types of negative results proposed by the\nreviewer: (1) whether compounds predicted to be unstable at 1273 K by our DFT phase diagram\ncalculations are indeed not synthesizable in literature; and (2) whether there could be experimental\nfailures for compounds predicted to be synthesizable by our DFT phase diagram calculations,\nespecially those lie near the synthesizability boundary.\nTo address the first type of negative results, we have compiled 30 MPEAs that were attempted\nat equimolar composition but resulted in multiphase outcomes (Table S7). All of these cases show\npositive values of E – 1273K × S , indicating that our DFT phase diagram evaluation is also\nhull ideal\neffective in predicting the instability of MPEAs. Taken together with the successful predictions of\nsynthesizable materials (Fig. 2 and Table S6) and unsynthesizable ones (Table S7), we now\nprovide a more systematic justification of the reliability of our synthesizability evaluation scheme.\nFor the second type of negative results, we acknowledge that the DFT-computed E value may\nhull\nsometimes mispredict synthesizability, especially for those that lie near the synthesizability\nboundary. On one hand, failure to synthesize a compound after one or two experimental attempts\ndoes not necessarily imply that the compound is unsynthesizable. On the other hand, our DFT-\ncalculated E is not a perfect predictor of synthesizability, which has been discussed thoroughly\nhull\nin our response to Question 1 from the same reviewer in the first round of revision, and on page\n17 of the revised manuscript (highlighted in red and begins with “Incorporating vibrational free\nenergy and chemical short-range order (CSRO) would undoubtedly improve accuracy…”).\nWe believe that deliberately pursuing negative experimental outcomes would be very time-\nconsuming and would detract from the main focus of this work. Instead, we have strengthened the\nDiscussion section by emphasizing that DFT-computed energies and the E metric do not always\nhull\nguarantee successful synthesis. We have elaborated on the limitation of the current framework\nwhile suggesting more rigorous synthesis design when targeting compounds near the\nsynthesizability boundary.\nAccordingly, we have made the following revisions in the revised manuscript:\nPage 16: “In contrast, for the literature-reported equimolar compositions that failed to form a single\nphase (Tab. S7), all 30 collected compositions exhibit positive E - S × 1273K values. Taken\nhull ideal\ntogether, these results support the use of the criterion E < S × T as a reasonable approximation\nhull ideal\nfor synthesizability.”\nPage 17: “The E values and phase diagram calculations adopted in this work are certainly not\nhull\nperfect metric for estimating synthesizability. They may fail to predict the synthesizability and\nrequire a thorough synthesis design of some MPEAs, particularly for compounds near the\nsynthesizability boundary (e.g., E - S × 1273K being close to zero).”\nhull ideal\n4. The model relies on a linear combination of lower-dimensional subsystems. On page 12, the\nauthors write: \"…we assign weights \\epsilon^i to each subspace such that the overall cluster\nprobability difference is minimized. A convenient and physically meaningful choice is \\epsilon^i\n= 1/C_N^M\". However, the following points remain unclear and should be addressed:\nWhat exactly is meant by the “overall cluster probability difference”? How is this quantity\nminimized by the proposed choice of \\epsilon^i? We do not find a discussion in the main text or\nHow would one determine an optimal choice of weights?\nResponse: It appears that all the questions grouped under bullet point 4 stem from the reviewer’s\ndifficulty in understanding cluster expansion theory and our model. We therefore address them\ntogether, providing progressive explanations to clarify the concepts for the reviewer:\nWhat exactly is cluster probability difference:\nBy “overall cluster probability difference,” we refer to the deviation between the cluster\nprobabilities of the N-component system and those reconstructed from its M-component subspaces.\nAs a reminder, in cluster expansion theory, the total energy is defined as:\n𝐸 = ∑ 𝐽 𝑛 +∑ 𝐽 𝑛 +∑ 𝐽 𝑛 +⋯ (Eq 1 in manuscript)\n𝑖 𝑖 𝑖 𝑖,𝑗 𝑖𝑗 𝑖𝑗 𝑖,𝑗,𝑘 𝑖𝑗𝑘 𝑖𝑗𝑘\nwhere J represents the energy of each cluster, and n represents the probability of each cluster.\nAccordingly, the energy differences among different MPEA structures can be expressed as the\nproduct of cluster energies and cluster probability differences:\n∆𝐸 = ∑ 𝐽 ∆𝑛 +∑ 𝐽 ∆𝑛 +∑ 𝐽 ∆𝑛 +⋯\n𝑖 𝑖 𝑖 𝑖,𝑗 𝑖𝑗 𝑖𝑗 𝑖,𝑗,𝑘 𝑖𝑗𝑘 𝑖𝑗𝑘\nHere, the coefficients Δn here is the cluster probability difference for each cluster.\nHow is cluster probability difference minimized by the proposed choice of \\epsilon^i:\nWe propose that the energy fitting of an N-dimensional MPEA can be approximated by a linear\nsummation of the energies of its M-dimensional MPEA subspaces, expressed as:\n𝐸 ≈ ∑ 𝜀𝑖𝐸𝑖\n𝑁 𝑖∈(1,𝐶𝑀) 𝑀\n𝑁\nThe corresponding cluster probability difference of two sides of the equation can be written as:\np p\nΔn = ( 1 ) −𝜀𝑖C M−q ( 1 ) (from Eq. 3 of the manuscript)\nq\nc p N N−q M\nwhere p denotes the cluster size and q indicates the number of species in the cluster. Therefore,\nwhen we state that “such that the overall cluster probability difference is minimized”, we want to\nminimize |Δn | at specific p and q. For example, for a singlet cluster (p=q=1), the cluster\nq\nc\np\nprobability difference of singlet reduces to\n1 1 1 1\nΔn = ( ) −𝜀𝑖CM−1( )\nN N−1 M\nBy choosing 𝜀𝑖 = 1 , the |Δn| value is minimized to zero.\n𝐶𝑀\n𝑁\nAre there other reasonable choices for \\epsilon^i? How would they perform in comparison?\nHow would one determine an optimal choice of weights?\nAll those questions relate to well-established features of cluster expansion. We provide an in-depth\nwalkthrough below. More detailed explanations of cluster expansion theory can be found in DOIs:\n10.1103/PhysRevB.17.2926, 10.1016/0378-4371(84)90096-7, 10.1016/B978-0-12-525102-\n0.50013-X, and 10.1103/PhysRevB.106.144202.\n1) Choices of 𝜺𝒊: In principle, any value can be chosen for 𝜀𝑖.\n5\n2) Physical intuition for why 𝜺𝒊 = 𝟏/𝑪𝑴 is the best choice: To minimize the ∆𝐸, 𝜀𝑖 should be\n𝑵\nchosen such that the cluster probability difference is minimized. Setting 𝜀𝑖 = 1/𝐶𝑀 reduces the\n𝑁\nsinglet probability difference to zero. Since singlets are well established in cluster expansion\ntheory as the largest contributor to total energy to total energy (DOIs: 10.1103/PhysRevB.17.2926,\n10.1016/0378-4371(84)90096-7, 10.1016/B978-0-12-525102-0.50013-X, and\n10.1103/PhysRevB.106.144202), eliminating the singlet cluster difference ensures the minimum\npossible energy difference, e.g., the lowest mean absolute error (MAE).\n3) Numerical support of the fact that 𝜺𝒊 = 𝟏/𝑪𝑴 as optimal choice, and comparison with\n𝑵\nother weights: To provide direct numerical evidence for reviewer’s reference, we applied our\nmodel with different 𝜀𝑖 values when predicting quaternary MPEAs from ternary MPEAs.\nFor instance, when 𝜀𝑖 = 3/16, the pair clusters probability difference is minimized to zero. When 𝜀𝑖 =\n9/64, the triplet cluster probability difference is minimized to zero. Beyond these special cases,\nwe have evenly sampled 40 𝜀𝑖 values, as shown in Fig. S9 and below. The results clearly\ndemonstrate that choosing 𝜀𝑖 = 1/𝐶𝑀 = 1/𝐶3 = 1/4 yields the lowest prediction error.\n𝑁 4\nFig. S9. The relationship between different weights 𝜀 and predicting accuracy when using ternary alloys to predict\n𝑖\nthe energy of quaternary alloys, blue circles are arbitrary coefficients sampled between 1/3 and 1/7 and stars represent\nthe situation when the error generated by singlet, pair, and triplet cancel.\n4. A step-by-step guide on why 𝜺𝒊 = 𝟏/𝑪𝑴 minimizes singlet cluster probability difference to\n𝑵\nzero using a quaternary example:\nWe provide a simple example to clarify the reasoning. Consider predicting the quaternary system\nABCD using its ternary subspaces ABC, ABD, BCD and ACD. To find a linear combination of\nthe weights 𝜀 ABC, 𝜀 ABD, 𝜀 BCD, and 𝜀 ACD that reproduces the same elemental\n1 2 3 4\nconcentrations (singlet probability) as in ABCD, one needs to solve the equation:\n1 1 1 1\n𝜀 + 𝜀 + 𝜀 = (Conserve concentration of A)\n1 2 4\n3 3 3 4\n1 1 1 1\n𝜀 + 𝜀 + 𝜀 = (Conserve concentration of B)\n1 2 3\n3 3 3 4\n1 1 1 1\n𝜀 + 𝜀 + 𝜀 = (Conserve concentration of C)\n1 3 4\n3 3 3 4\n1 1 1 1\n𝜀 + 𝜀 + 𝜀 = (Conserve concentration of D)\n2 3 4\n3 3 3 4\nThe only solution is that 𝜀 = 𝜀 = 𝜀 = 𝜀 = 1/𝐶3 = 1 .\n1 2 3 4 4\n“The uniform choice of 𝜀𝑖 ensures that the energy contributed by the singlet term cancels exactly\nwhen averaging over all subspaces, thereby reducing the overall cluster probability difference”\nThe discussion on optimal choice of uniform weights 𝜀𝑖 = 1/𝐶𝑀 is updated in revised\n𝑁\nsupplementary Note 2 on page 10:\n“As shown in Eq. 8 in supplementary information, in principle, the weights can take any value,\nwhile it can also be designed to cancel the error contribution from clusters of any size. For example,\nwhen using ternary alloys to predict quaternary alloys, the uniform weights required to cancel the\nsinglet, pair, and triplet terms are 1/4, 3/16, and 9/64, respectively. As shown in Fig. S9 in yellow\nstars, the choice of 1/4 (corresponding to cancellation of the singlet term) shows the lowest\nprediction error. This behavior shows that the singlet term presents the dominant contribution to\nthe total energy of an alloy system. Failing to cancel this term naturally results in larger errors.”\nThe discussion on arbitrary choice of uniform weights 𝜀𝑖 = 1/𝐶𝑀 is updated in revised\n𝑁\nSupplementary Note 2 on page 10:\n“To better illustrate this point, we compared the performance of different weights by fitting all\nquaternary MPEAs using ternary MPEAs with arbitrary coefficients sampled between 1/3 and 1/7.\nAs shown by the blue circles in Fig. S9, these arbitrary choices lead to significantly higher mean\nabsolute errors (MAE) compared to the uniform choice of 1/𝐶3 = 1/4. Such observation confirms\n4\nthat cancelling the singlet contribution provides the most effective strategy for minimizing\nprediction error, consistent with the physical dominance of site energy terms.”\n5. While the authors refer to this as “our model,” they also acknowledge that it is built largely on\ncluster expansion ideas. The novelty introduced in this work should be clarified. Otherwise, It\nwould still be interesting to know that a simple, well-known approach can be used as a good\nestimator for the synthesizability of MPEAs, but the manuscript should be transparent about it.\nResponse: We hope that the clarifications provided in responses to questions 1–4 can help the\nreviewer better follow the logical flow of our paper and gain a clearer understanding of our model.\nWith such foundation, the novelty of our approach now becomes evident. We have also clarified\nas below:\nFirst, our model is the approximation of energy of any MPEA with 𝐸 ≈ ∑ 𝜀𝑖𝐸𝑖 . This\n𝑁 𝑖∈(1,𝐶𝑀) 𝑀\n𝑁\nmodel by itself is not cluster expansion.\nMoreover, the connection of our model to cluster expansion lies only in using cluster expansion\nformalism to justify that our model is both reasonable and physically interpretable. We did not\nperform a standard cluster expansion fitting.\nMore importantly, a cluster expansion fitting procedure is not feasible for global MPEA modeling.\nAs established in the cluster expansion literatures (DOIs: 10.1103/PhysRevB.17.2926,\n10.1103/PhysRevB.106.144202, and 10.1002/aenm.201903240), fitting even a single composition\nsubspace (e.g., A–B–C–D) typically requires hundreds to thousands of DFT calculations.\nExtending this to a classic cluster expansion involving 28 metals would demand millions or even\nbillions of calculations.\nTherefore, the central novelty of this work is that, guided by the physical intuition of cluster\nexpansion, we demonstrate for the first time that the stability of MPEAs can be estimated directly\nfrom a linear combination of low-dimensional basis, without performing conventional cluster\n7\nexpansion fitting. This approach drastically reduces the computational cost of modelling MPEAs,\nwhile its scientific validity remains firmly supported by cluster expansion theory.\nWe have also provided the following clarification statements in the revised manuscript on page 19:\n“While inspired by the cluster expansion formalism, our framework does not involve conventional\ncluster expansion fitting, which would require millions to billions of DFT calculations across a 28-\nelement space. Instead, we estimate the energy of an N-component equimolar alloy as a linear\ncombination of lower-dimensional subspaces. To our knowledge, this is the first demonstration\nthat Cluster Expansion concepts can be adapted to justify a more data-efficient and physically\ninterpretable model for predicting MPEA stability.”\n6. We also request clarification on how the mean absolute error (MAE), mentioned in Fig. 4 and\nelsewhere, is computed. Is it related to \\Delta n in Eq. (3)? How exactly does it relate to\nsynthesizability?\nResponse: We thank the reviewer for the question. The mean absolute error (MAE) is calculated\nas the average of the absolute differences between the energies predicted by our model and the\nDFT total energies. Since our model is based on 𝐸 ≈ ∑ 𝜀𝑖𝐸𝑖 , the MAE can be expressed\n𝑁 𝑖∈(1,𝐶𝑀) 𝑀\n𝑁\nas:\n𝑀𝐴𝐸 = 1/𝑁 ∙∑ |𝐸 −∑ 𝜀𝑖𝐸𝑖 |\n𝑀𝑃𝐸𝐴 𝑎𝑙𝑙 𝑀𝑃𝐸𝐴𝑠 𝑁 𝑖∈(1,𝐶𝑀) 𝑀\n𝑁\nAs clarified in response to question 1, the total energy of a system is the key for synthesizability,\nsince we use E derived from phase diagram calculations as our metric. As E is defined as:\nhull hull\n𝐸 = 𝐸 −𝐸\nℎ𝑢𝑙𝑙 𝑀𝑃𝐸𝐴 𝐶𝑜𝑚𝑝𝑒𝑡𝑖𝑛𝑔 𝑝ℎ𝑎𝑠𝑒𝑠\nwhile 𝐸 is fixed by the convex hull construction, the MAE of 𝐸 is thus\n𝐶𝑜𝑚𝑝𝑒𝑡𝑖𝑛𝑔 𝑝ℎ𝑎𝑠𝑒𝑠 𝑀𝑃𝐸𝐴\ndirectly corresponding to the MAE of E value, hence also the synthesizability. The capability\nhull\nof evaluating energies of MPEA through our model will thus enhance the efficiency of estimating\nsynthesizability.\nWe have also added the following clarifications to the revised manuscript:\nPage 14-15: “where mean absolute error (MAE) is calculated as the average of the absolute errors\nbetween the predicted and DFT total energies”\nPage 15: “The predicted total energies of higher-dimensional alloys can then be used to construct\nthe convex hull and determine the corresponding E values, which can rapidly estimate\nhull\nsynthesizability.”\n7. Beyond scientific and technical issues, the structure and flow of the manuscript need\nconsiderable improvement. The Results section should clearly present the findings, while\nexplanations, caveats, and interpretation should be reserved for the Discussion. Currently, the\nDiscussion section is short and underdeveloped, with only one subsection. In contrast, the\nResults section contains several paragraphs better suited to the Discussion. For example, the\nparagraph on page 6 starting with \"Although the criterion Ehull < Sideal × T is a rough approximation,\n…\" should be relocated to the Discussion.\nResponse: We appreciate the constructive comments from the reviewer. In the revised manuscript,\nwe have moved all discussions in the original Result section into Discussions as section 3.1.\nIn this section, we adopted the original context on page 6, starting with \"Although the criterion\nhull ideal\nthe capability and limitations of our synthesizability estimation framework.\nAdditional Comments:\n8. Section Titles: The section “2.4. Stability rules across extended dimensionality” introduces the\ncore model, yet the title does not reflect its central role. The title should be revised for clarity.\nResponse: We thank the reviewer for this suggestion. To better emphasize the importance of this\nsection, we have revised the title to “Universal framework across extended dimensionalities and\ncompositions”.\n9. Contradictory Claims: On page 6, the authors state that noble-metal-containing MPEAs are\n“underexplored,” then later note they have “garnered significant attention.” These statements\nconflict and should be reconciled.\nResponse: We thank the reviewer for this observation. To clarify our statement that noble-metal-\ncontaining MPEAs are underexplored, we want to emphasize that this compositional space remains\nrelatively less explored, with only 168 compositions reported across 148 papers (Fig. 2a). In\ncontrast, other MPEAs have been studied in several thousand papers, covering a few thousand\ncompositions, as indicated by the two references (DOIs: 10.1016/j.actamat.2016.08.081, and\n10.1038/s41598-019-50015-4).\nAt the same time, the noble-metal contained MPEAs have gained increasing significant attentions\nmore recently as promising catalysts, which can be represented by several recent reviews and\nresearch with DOIs: 10.1126/science.abn3103; 10.1126/sciadv.abg1600; 10.1126/sciadv.aaz0510,\n10.1038/s41578-019-0121-4, 10.1126/science.aav4302, 10.1126/science.aaf8402 and\n10.1126/sciadv.adn2877).\nThese two observations indicate that an increasing number of researchers are recognizing the\npotential of noble-metal-containing MPEAs for various applications, while the compositional\nspace remains relatively underexplored.\nTo further improve the clarify of our manuscript, we have revised the manuscript to highlight that\nnoble-metal-containing MPEAs are underexplored in terms of compositional coverage but have\nattracted increasing attentions for their catalytic applications.\nThe clarification is updated on pages 6-7 in the revised manuscript:\n“Among all MPEA compositions, noble-metal-containing MPEAs emerge as a distinct and\nunderexplored group in terms of compositional coverage. As shown in Fig. 2a, this compositional\nspace remains relatively underexplored, with only about one hundred reported compositions\npublished from our survey (Fig. S4 and Tab. S1), in contrast to over several thousand papers and\ncompositions reported for other MPEAs”\n“Notably, Noble-metal-containing MPEAs have also garnered significant attention due to their\ngreat potential in catalysis and other applications.”\n10. Disjointed Transitions: In the last sentence of page 6, the discussion abruptly shifts from noble-\nmetal-containing MPEAs to equimolar, non-noble-metal MPEAs. The logic of this transition is\nunclear. Furthermore, the validation based on previously reported MPEAs would be better placed\nin the Discussion section.\nand ensured that it now focuses exclusively on noble-metal-containing MPEAs. In addition, the\ncomparison between our DFT predictions and experimental reports for all MPEAs has been moved\nto the Discussion section, so that only noble-metal-containing MPEAs are discussed at the end of\npage 6 (now the first paragraph of Section 2.2 in the revised manuscript). The broader discussion\ncomparing DFT predictions with experimental reports for all MPEAs, including those without\nnoble metals, is now placed on page 16 as part of Section 3.1, which is shown below:\n“To further evaluate the capabilities and limitations of E beyond the nine experiments shown in\nhull\nFig. 2, we compiled 32 additional reported equimolar MPEA compositions (Tab. S6), all of which\nhave been experimentally synthesized below 1273 K. Among these 32 compositions, 30 are\npredicted to be synthesizable by our DFT-based phase diagram calculations. In contrast, for the\nliterature-reported equimolar compositions that failed to form a single phase (Tab. S7), all 30\ncollected compositions exhibit positive E - S × 1273K values. Taken together, these results\nhull ideal\nsupport the use of the criterion E < S × T as a reasonable approximation for synthesizability.\nhull ideal\nThis further implies that our universal framework (Fig. 4) can be used to rapidly predict the\nsynthesizability of a broad range of MPEAs.”\n11. Confusing Wording on Page 7: The paragraph beginning with “As a contrast to the less\nexplored noble-metal-containing MPEAs…” is difficult to follow. It is unclear whether it refers to\nnoble or non-noble alloys, and the connection to Fig. 2b should be made explicit.\nResponse: We appreciate the reviewer’s comment. We have revised the paragraph to avoid\npotential confusion. The revised sentence is highlighted in red on page 7 of revised manuscript as:\n“In this work, using high-throughput phase diagram calculations, we found 647 ternary, 2,178\nquaternary, and 5,478 quinary noble-metal-containing MPEAs that can be stabilized below 1273\nK with ideal mixing entropy (Fig. 2b).”\n12. Late Definitions: The quantity S_ideal is defined on page 8, despite being used earlier in Fig.\n1 (page 4). It should be introduced at first mention.\nResponse: We thank the reviewer for pointing this out. We have revised the manuscript to define\n𝑆 at its first mention (page 5 at the figure caption of Fig. 1). Specifically, we now update the\nideal\nfollowing description when it first appearance of 𝑆 :\nideal\n“Here, 𝑆 is the ideal mixing entropy, defined as 𝑆 = 𝑘 ln(𝑛 ) with 𝑘 as the Boltzmann\nideal ideal 𝐵 M 𝐵\nconstant and 𝑛 as the number of metal elements in MPEA”\nM\n13. Misplaced New Content: The red-colored text on page 14 beginning with “To validate our\nmodel with higher-dimensional DFT-computed MPEAs…” belongs with the methodology\npresented earlier in Section 2.1 on formation energy and energy above the hull.\nResponse: We thank the reviewer for the suggestion. We believe the paragraph beginning with\n“To validate our model with higher-dimensional DFT-computed MPEAs…” may not be best\nplaced in the methodology section. This paragraph describes how the higher-dimensional dataset\nof 9,362 equimolar MPEAs was constructed, with the primary goal of providing direct evidence\nfor the transferability of our model. The DFT calculation details are identical to those already\nreported in the methodology section. However, the purpose of including the composition\ninformation here is to directly support the evaluation of whether low-dimensional MPEA basis can\nbe used to predict the energies of higher-dimensional MPEAs beyond 5D. For this reason, we\n10\nbelieve the current placement is appropriate, as it highlights the significance of our validation\nanalysis.\n14. As emphasized, these are illustrative examples—further improvements to clarity, consistency,\nand structure are needed throughout the manuscript. We leave it to the authors' discretion to\naddress these issues thoroughly.\nResponse: We thank the reviewer for this valuable feedback. During the revision, we also carefully\nreviewed the entire manuscript to ensure clarity, consistency, and smooth flow.\n15. Finally, a few specific points on presentation:\nInterpretability: The abstract states: “...while retaining full interpretability.”\nThe claim of “full interpretability” is quite strong. The authors should emphasize in the text, e.g.,\nin the Discussion or Conclusion, what interpretability means in this context, and how their model\nsatisfies it.\nResponse: We thank the reviewer for raising this point. By “fully interpretable,” we mean that\nevery prediction made by our model can be decomposed into chemically meaningful clusters.\nEach coefficient J of the cluster as shown in Equations 1 reflects the contribution of a specific\ncluster, and the overall prediction is obtained as a transparent sum of these terms. This contrasts\nwith black-box machine learning models, where the internal decision process is often inaccessible.\nIn our framework, each contribution can be clearly traced to specific atomic interactions, so the\norigin is easy to understand. This provides not only accurate predictions but also clear physical\ninsight.\nTo avoid the absolute statement, the “fully interpretability” is changed into “while retaining full\ninterpretability through cluster expansion theory” on pages 2 and 17.\n16. Fig. 1b: The value of the scatter plot in Fig. 1b to the overall narrative of the manuscript is\nquestionable. It appears to consider only 4-element (4M) MPEAs and does not clearly\ncommunicate essential information, such as the number of stabilized phases or which data points\ncorrespond to specific alloys or elements. If the figure is to be retained, we suggest replacing it\nwith Fig. S5, which presents more informative data—including the formation energy and energy-\nabove-hull values for the nine synthesized phases. However, Fig. S5 also requires correction: its\ncaption refers to 4-element MPEAs, but it includes three 5-element MPEAs (marked with stars)\nthat were fabricated. This inconsistency should be addressed.\nResponse: We appreciate the reviewer’s suggestions regarding Fig. 1b. The purpose of this figure\nwas not to emphasize specific alloy compositions, but rather to illustrate the distribution of our\ndataset in terms of two commonly used synthesizability metrics, namely formation energy and\nenergy above hull. While we value the suggestion to highlight the nine experimentally attempted\ncompounds, doing so may disrupt the logical flow. At the point where we discuss Fig. 1b, we have\nnot yet introduced any experimental validations, nor have we transitioned to the rationale for\nfocusing on noble-metal-based MPEAs.\nRegarding Fig. S5, we thank the reviewer for pointing out caption inconsistency. The figure\ncaption has been changed into “Scatter plot of the formation energy (𝐸 ) and energy above hull\nform\n(𝐸 ) for the nine successfully synthesized MPEAs (yellow stars), while all computed MPEAs\nhull\nup to quinary systems are plotted with the same style as Fig. 1b.”\nReviewer #2 (Remarks to the Author):\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We would like to thank the constructive questions from both reviewers 1 and 2. It is\nbelieved that our paper has been further improved after implementing all these suggestions.\nReviewer #3 (Remarks to the Author):\nThe paper can be accepted as it is.\nResponse: We would like to thank the reviewer for taking the time to review our paper.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We would like to thank the reviewer for taking the time to review our paper.\nReviewer #5 (Remarks to the Author):\nThe authors have comprehensively addressed all the reviewer's comments. The paper is now\nclearer and better structured. I recommend publication of the paper in its current form.\nResponse: We would like to thank the reviewer for taking the time to review our paper.\nUniversal Framework for Efficient Estimation of Stability in Multi-\nPrincipal Element Alloys\nLin Wang1†, Bo Shen2,3†, Zheng-Da He1, Zihao Ye2,3, Yan Zeng1, Chad Mirkin2,3, and Bin\nOuyang1,*\n1Department of Chemistry and Biochemistry, Florida State University, Tallahassee, FL\n32304, USA.\n3International Institute for Nanotechnology, Northwestern University, Evanston, Illinois\n60208, USA.\n†These authors contributed equally to this work\n*To whom correspondence may be addressed. E-mail: bouyang@fsu.edu\nReviewer #1 (Remarks to the Author):\nPlease remove the statement on line 377 in the main text: \"To our knowledge, this is the first\ndemonstration that Cluster Expansion concepts can be adapted to justify a more data-efficient and\nphysically interpretable model for predicting MPEA stability.\" Many other's have used this notion,\nstarting with Troparevsky (2015).\nResponse:\nWe thank the reviewer for this helpful clarification. We agree that binary-interaction-based\napproaches, including Troparevsky et al. (2015), have previously shown that the energy of binaries\ncan guide screening of HEAs. We have now removed the statement implying the first\ndemonstration. To clarify prior context and position our contribution appropriately:\nIn the Result section, we revised the paragraph to acknowledge Troparevsky et al. and to describe\nour contribution as building on this physical foundation while providing a data-efficient and\ninterpretable formulation grounded in Cluster Expansion concepts. The following updates can be\nfound in our main text:\nOn pages 20-21: “Earlier work has demonstrated that binary interaction energetics provide\nguidance for identifying single-phase candidates in compositionally complex alloys20. Building\nupon this physical foundation, we demonstrate that Cluster Expansion concepts can be adapted to\njustify a more data-efficient and physically interpretable model for predicting MPEA stability.”\nIn the Introduction (page 3), we added text acknowledging earlier work that used binary formation\nenergies for initial screening.\n“Earlier work examined HEA formation using binary formation energies as an initial screening\napproach. Subsequent studies have explored a variety of approaches using machine learning (ML)\nmodels and heuristic theories”\nThese revisions ensure that the manuscript accurately credits prior work while clearly articulating\nthe distinct contribution of our approach.\nReviewer #2 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We would like to thank the comments from both reviewers 1 and 2. It is\nbelieved that our paper has been further improved after implementing the suggestion.\n2"}]

Decision letter