Composition-based phase stability model for multicomponent metal alloys

,


I. INTRODUCTION
Multicomponent metal alloys, including those labeled "high-entropy" that contain several elements in nondilute concentrations, constitute a vast space of possible structural materials. 1,2For example, there are approximately 3.5 billion possible compositions that contain six out of 55 metallic elements in multiples of 10% (e.g., Cr 0.1 Mn 0.2 Fe 0.3 Co 0.1 Ni 0.2 Os 0.1 ). 3 The vastness of this space, as compared to the smaller space of more traditional alloys containing fewer elements, presents a significant challenge to materials design efforts.][11][12][13] In this work we focus on assessing whether a given hypothetical alloy composition is realizable as a singlephase solid solution.Our target is alloys that are thermodynamically stable (i.e., rather than kinetically stable) given our interest in high-temperature applications for which metastable materials may be less robust.0][21][22][23][24] First-principles electronic structure calculations of alloy free energy [25][26][27][28] in principle could provide a rigorous, accurate, and broadly applicable alternative model that avoids errors relating to limited experimental training data, which is particularly important for application to unexplored parts of composition space.
Along these lines, Troparevsky et al. developed a model to predict the formation of a single-phase solid solution purely based on density functional theory (DFT) energetics. 29They reasoned that a single-phase solid solution will not form if there exists, for at least one pair of constituent elements, either (1) one or more highly stable (i.e., very negative formation enthalpy) binary intermetallic phases or (2) all highly unstable (i.e., very positive formation enthalpy) binary intermetallic phases.A multiphase alloy is expected for the first case due to precipitation of such intermetallic(s) and for the second case related to a tendency for elemental phase separation.The criteria were assessed by looking up, for each pair of constituent elements, the lowest binary intermetallic formation enthalpy (for an ordered structure of any stoichiometry) present in an existing DFT database.The free energy contribution corresponding to the ideal configurational entropy of a five-element equiatomic composition and a typical experimental annealing temperature was used to set the lower formation enthalpy limit, whereas the upper limit was chosen empirically.
Here, we develop a thermodynamic model that is based on the same well-reasoned physical intuition of Troparevsky et al. and retains its computational affordability, but has several key advantages.Since we explicitly model the stoichiometry-dependent Gibbs free energy and determine ground states via the convex hull construction, our model ( 1) is a proper stoichiometry-dependent model that takes into account, for example, the location of competing intermetallic phases in composition space, (2) has no heuristic and/or empirical formation enthalpy bounds, (3) provides more complete information, such as the complete set of phases expected to form for a multiphase alloy, and (4) is applicable to both equiatomic com-positions, the focus of the model of Troparevsky et al., and non-equiatomic compositions.Our model achieves an accuracy of 70-75% when validated on high-entropy alloy experiments and provides a platform upon which additional physical effects (e.g., short-range order, vibrational entropy, defect configurational entropy) could be incorporated for improved accuracy.

II. THERMODYNAMIC MODEL
We consider a single solid solution phase and a finite number of ordered intermetallic competing phases that are treated as line compounds.The intermetallic phases are among the approximately one million materials in the Open Quantum Materials Database (OQMD), 30,31 which contains DFT calculations of both experimentally known phases from the Inorganic Crystal Structure Database (ICSD) and hypothetical ones derived from decoration of crystal prototypes.Our model is a composition-only model in that we consider a single solid solution phase of unspecified crystal structure at each composition.
Our model for the Gibbs (formation) free energy for the solid solution ∆G alloy at temperature T consists of an enthalpy term ∆H alloy and a contribution from the entropy ∆S alloy : ∆G alloy = ∆H alloy − T ∆S alloy . ( Since we only consider dense solid phases at ambient pressure in this work, the pressure-volume term in the enthalpy is negligible and we take all formation enthalpies to equal formation energies.The entropy is approximated as the ideal configurational entropy ∆S alloy = −k B i x i ln(x i ), where x i is the atomic percentage of each element i in the composition (normalized with i x i = 1) and k B is the Boltzmann constant.The alloy entropy thus has no temperature dependence in our model.Other entropy contributions such as non-ideal configurational entropy and vibrational entropy are not included in this work to avoid a more expensive and complex model.
The alloy enthalpy, similarly treated as independent of temperature, is derived from (zero-temperature) DFT energetics.To evaluate the solid solution enthalpy term, we first construct the enthalpy of a binary solid solution ∆H alloy α,β (for distinct elements α and β).In principle, one can estimate the enthalpy of the solid solution for a particular lattice by averaging the enthalpies of many different ordered structures on that lattice, to the extent the enthalpy is a well-behaved function of local atomic configurations and the ordered structures appropriately sample configurations contained in the solid solution.As an illustration of this, the solid solution formation energy obtained via quasirandom structure calculations is commonly seen to take on an average of the formation energies of ordered structures considered when generating a cluster expansion (see, for example, Ref. 32).Taking inspiration from this property, we determine ∆H alloy α,β based on the formation energies of intermetallic phases in the OQMD, which includes various binary crystal structures that are orderings on the fcc (e.g., D0 22 , L1 0 , L1 1 , L1 2 ), bcc (e.g., B2, D0 3 , B 19 ), hcp (e.g., D0 19 ), and other (e.g., B1, B3, B4, B h ) lattices for every α-β pair, in addition to experimental crystal structures from the ICSD.This approach is similar in spirit to the small set of ordered structures (SSOS) method, 33 although we do not consider any weighting of different ordered structures based on their local structural correlations.Since intermetallic phases corresponding to multiple underlying lattices are present, we obtain a general (i.e., not associated with a specific lattice) enthalpy tendency suitable for our compositiononly model.
To obtain an average, we fit the ordered structure formation energies to the polynomial form where A is a scalar fitting parameter encapsulating the enthalpy behavior for each elemental pair.The fitting form in Eq. 2 is the zeroth-order Redlich-Kister polynomial 34 and corresponds to the enthalpy term in a regular solution model.However, since we are not fitting the formation energies of phases that all correspond the same parent lattice, we can view the x(1 − x) form as essentially a simple fitting function to capture the overall binary enthalpy trend.To avoid disruption of the fit by any outlier phases (e.g., stemming from structures that are extremely high energy for the composition), we ignore any phases whose formation energy is beyond the mean value, for the full binary space, by more than 1.5 times the standard deviation.We use the fitted binary formation enthalpies to estimate the formation enthalpy of the multicomponent solid solution (containing N elements) via the linear interpolation equation where the sum is over distinct elemental pairs.As in the Kohler method, 35 here the stoichiometry of the relevant binary composition is determined simply by the relative amount of the two elements, i.e., x α /(x α + x β ).The weighting factor w α,β = (x α + x β )/(N − 1) is chosen to preserve mass balance.As a simple example, in this method the formation enthalpy of NbTaW is equal to an average of ∆H alloy Nb,Ta (x = 1/2), ∆H alloy Ta,W (x = 1/2), and ∆H alloy W,Nb (x = 1/2).We note that alternative interpolation schemes such as Muggianu's approach exist, 36-40 but they are not explored in this work given we expect the binary enthalpy contributions themselves to be a more important source of uncertainty.
We consider two enthalpy models for the multicomponent alloy.In the first, we only include a contribution from ∆H alloy α,β if all intermetallic phases in the α-β binary space are unstable, i.e., w α,β is set to zero if there are any stable (negative formation energy) intermetallics in the space.This enthalpy model (labeled "only positive enthalpy contributions") is in the spirit of the work of Troparevsky et al. since in their model, with respect to the upper formation energy limit, it is only the presence of an elemental pair with all sufficiently unstable formation energies -regardless of possible stabilizing enthalpic contributions from other elemental pairs -that leads to the prediction of an unstable solid solution.By contrast, in our second enthalpy model (labeled "positive and negative enthalpy contributions"), we simply include contributions from all elemental pairs.Finally, for the Gibbs free energy of the intermetallic competing phases ∆G intermetallic , we only include the formation enthalpy part (obtained directly from the OQMD), i.e., they have no entropy.Configurational entropy associated with defects/off-stoichiometry and vibrational entropy are not included to avoid the need for expensive calculations.We note that since the Gibbs free energy of the intermetallic competing phases therefore has no temperature dependence, only intermetallics that appear on the zero-temperature convex hull have the possibility to appear on the Gibbs free energy convex hull (at any temperature) in our model.Figure 1 shows a schematic example of the Gibbs free energy behavior within our model.The solid lines are the convex hull of Gibbs free energy versus composition (for a particular temperature), which is determined from the solid solution Gibbs free energy ∆G alloy (dashed line) and those of the intermetallic competing phases ∆G intermetallic (points).The black segments correspond to multiphase cases in which at least one intermetallic phase is contained in the thermodynamic ground state.
In contrast, the red parts of the convex hull (whose limits are determined by, e.g., the common tangent construction) are regions in which the ground state is a singlephase solid solution.Only a binary case is shown for simplicity but the behavior for ternary and higher cases is analogous.We note that in the actual implementation of our thermodynamic model for the general multicomponent case, we employ a discretized stoichiometry grid for the solid solution phase to enable use of the existing convex hull construction capabilities in the OQMD based on the quick hull algorithm in Qhull, 41 as opposed to common tangent construction.
The schematic in Fig. 1 is drawn to emphasize the possibility of, considering solely the intermetallic phases, a convex hull that is appreciably asymmetric, i.e., not invariant to permutation of different elements.Given our model is capable of representing such an asymmetry via full inclusion of the intermetallic competing phases, it can capture the resulting asymmetry of the predicted region(s) of single-phase solid solution.This contrasts with approaches that only consider the value of an intermetallic formation enthalpy and not its corresponding stoichiometry.

III. RESULTS AND DISCUSSION
Figure 2 shows the fitting of ∆H alloy α,β (x α ) for four example elemental pairs: Nb with Ir, Al, Y, and Ce.Nb-Ir and Nb-Al generally have negative (stable) formation energies leading to a negative (favoring mixing) enthalpy contribution, in contrast to the phase separating tendency of Nb-Y and Nb-Ce.We observe that the fits capture the overall enthalpy trends reasonably well.The symmetric (about x α = 1/2) nature of the zerothorder Redlich-Kister polynomial does not appear problematic given the overall formation energy behavior (distinct from the enthalpy convex hull discussed in Sec.II) is not especially asymmetric.Additional examples of fits of ∆H alloy α,β (x α ) are shown in the Supplemental Material.A set of the fit A α,β parameters considered in this work are displayed in Fig. 3. Like the table of the single most stable intermetallic formation energy contained in Fig. 1 of Ref. 29, this matrix is determined via a large number of DFT calculations (derived from the OQMD in this work).In the present case, the matrix quantifies the overall binary miscibility trends for the full stoichiometry range.
For the set of elements shown in Fig. 3, we observe a region of strong phase separation for elemental pairs including one of the lanthanides (La or Ce) with Re or an element from group V or VI.In contrast, strongly negative enthalpy contributions are found for elemental pairs with (1) Rh or Ir and (2) one of several rare-earth elements (La, Ce, Sc, Y), refractory elements (Ti, Zr, Hf, Ta), or Al.We note that the presence of very stable binary intermetallic competing phases, which may prevent the formation of a single-phase solid solution, con- tributes to the negative binary formation enthalpy contribution.Some similar trends can be observed in the corresponding table from Troparevsky et al., but there are significant quantitative differences since Fig. 3 represents the overall enthalpy trend as opposed to the single lowest formation energy.Similarly, the corresponding A α,β values obtained via the semiempirical Miedema model (commonly used to estimate multicomponent enthalpy; see, e.g., Refs.17,21,42) have similar broad trends but are quantitatively very different, as shown in the Supplemental Material.
To validate our thermodynamic model, we assess how accurately it predicts whether a given multicomponent alloy composition corresponds to a single-phase solid solution (classified here as "positive") as opposed to a combination of multiple phases (classified as "negative") in experiment.The positive (negative) case corresponds to the red (black) region in Fig. 1.In attempting to compare to experiment, several challenges include (1) uncertainty of whether an alloy is truly thermodynamically (rather than kinetically) stable, (2) dependence on synthesis method and processing conditions (not generally captured in our model), (3) experimental limitations in detecting small amounts of secondary phases, and (4) limited amounts of experimental data.We apply our model to 90 cases selected from a database of experimental measurements of high-entropy alloys. 43We only consider annealed samples and apply our thermodynamic model at the annealing temperature, with the assumption that the annealing treatment is sufficient to achieve thermodynamic equilibrium and therefore the observed single-phase or multiphase nature of the alloy is dictated by the thermodynamic phase stability at that temperature.We note that the presence of multiple distinct solid solution phases (e.g., with different lattices or rich in different elements) can occur in experiment but is not observed in our model due in part to its composition-only nature, which represents one source of error.
For comparison, we also test the model of Troparevsky et al.While their model used DFT data primarily from the AFLOW database, 29 here we use data in the OQMD to enable a fair comparison to our thermodynamic model.We do not find significant differences in predictions stemming from using the OQMD instead of AFLOW, which can be understood in part by the good agreement between the OQMD and AFLOW for intermetallic formation energies (in particular, median absolute difference of 2.7%). 44Additional details of the validation are included in the Supplemental Material.
Table I shows the validation results in the form of a truth table.For our thermodynamic model, we test both the version that only includes ∆H alloy α,β for elemental pairs if all intermetallic phases in the α-β binary space are unstable and the version with enthalpy contributions from all elemental pairs, as described in Sec.II.The best model in terms of accuracy is the model with only positive enthalpy contributions, which has an accuracy of 74%.If we include both positive and negative enthalpy contributions, the accuracy is slightly lower at 69%.However, given including both positive and neg- The improved performance of our thermodynamic model indicates the utility of a rigorous, non-heuristic approach to phase stability predictions for multicomponent alloys.Furthermore, the nonempirical nature of the method offers the potential for excellent transferability given the absence of any parameters fit to experiments on a specific chemistry.While the increased accuracy is desirable, we view the broadened applicability (to arbitrary temperature and stoichiometry) and more complete information (e.g., full set of competing phases) to be a key benefit of our model.For instance, our approach can be straightforwardly employed to help search for alloy compositions amenable to precipitation hardening, in addition to single-phase alloys.
With correct predictions observed approximately three out of four times, the thermodynamic model developed in this work appears to be moderately accurate but far from perfect, which is understandable given (1) the various neglected physical effects and significant approximations involved in the model construction (in part to ensure low computational cost) and (2) uncertainty in the experimental results.Therefore, we view the developed model as most applicable as part of broad, high-throughput searches (e.g., as an initial downselect), rather than for making highly quantitatively accurate predictions for a small number of compositions of interest.With the ability to utilize standard information existing in DFT databases, our approach can enable the acceleration of discovery and development within the vast space of multicomponent alloys.

IV. CONCLUSIONS
We develop a first-principles phase stability model for multicomponent metal alloys that considers an alloy composition of unspecified crystal structure.Taking inspiration from the more heuristic model of Troparevsky et al., we build a rigorous model by explicitly constructing the Gibbs free energy of the solid solution phase, with binary enthalpy contributions and ideal configurational entropy, and that of intermetallic competing phases, for which only enthalpy is considered.Our model captures the tendencies for intermetallic formation and elemental phase separation that can prevent formation of a singlephase solid solution.Since our model is computationally inexpensive (only requiring energies from DFT calculations that already exist in a database) and is reasonably accurate at classifying an alloy as single-phase or multiphase, it can aid high-throughput efforts for the discovery and development of multicomponent alloys.

FIG. 1 .
FIG.1.Schematic of the convex hull of Gibbs free energy in our model for a binary case.The blue points correspond to intermetallic competing phases (and the pure elements at the compositions x = 0 and x = 1) and the dashed purple line corresponds to the solid solution.The convex hull consists of (1) regions with a single-phase solid solution ground state (red curves) and (2) multiphase regions (black lines).

FIG. 2 .
FIG. 2. Binary formation enthalpy ∆H alloy α,β (xα) = A α,β xα(1 − xα) fit (blue line) to formation energies of ordered intermetallic structures (blue points) in the OQMD for four elemental pairs: Nb-Ir, Nb-Al, Nb-Y, and Nb-Ce.The abscissa is 1−x Nb .Red points are considered outliers and not included in the fit.The orange point indicates the region of enthalpy values that contributes to the multicomponent alloy enthalpy ∆H alloy for a Nb-rich alloy.

FIG. 3 .
FIG. 3. A α,β parameters representing the enthalpy behavior for a given pair of metal elements.Only the top-right of the table is shown as the bottom-left contains redundant information (A α,β = A β,α ).
Material for "Composition-based phase stability model for multicomponent metal alloys" Jay C. Spendlove, * Bryan H. Fong, John H. Martin, Mark R. O'Masta, Andrew Pan, Tobias A. Schaedler, and Eric B. Isaacs † HRL Laboratories, LLC, 3011 Malibu Canyon Road, Malibu, California 90265, USA (Dated: September 22, 2023) I. ADDITIONAL EXAMPLES OF BINARY ENTHALPY FITS In Fig. S1, we show additional examples of binary enthalpy fits for elemental pairs including Nb.As discussed in the main text, the fits capture the overall enthalpy trends reasonably well.
FIG. S1.Binary formation enthalpy ∆H alloy α,β (xα) = A α,β xα(1 − xα) fit (blue line) to formation energies of ordered intermetallic structures (blue points) in the OQMD for various elemental pairs including Nb.The abscissa is 1 − x Nb .Red points are considered outliers and not included in the fit.The orange point indicates the region of enthalpy values that contributes to the multicomponent alloy enthalpy ∆H alloy for a Nb-rich alloy.
FIG. S2.A α,β values from the semiempirical Miedema model.Only the top-right of the table is shown as the bottom-left contains redundant information (A α,β = A β,α ).