Temperature-Dependent Configurational Entropy Calculations for Refractory High-Entropy Alloys

The cluster expansion formalism for alloys is used to construct surrogate models for three refractory high-entropy alloys (NbTiVZr, HfNbTaTiZr, and AlHfNbTaTiZr). These cluster expansion models are then used along with Monte Carlo methods and thermodynamic integration to calculate the configurational entropy of these refractory high-entropy alloys as a function of temperature. Many solid solution alloy design guidelines are based on the ideal entropy of mixing, which increases monotonically with N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document}, the number of elements in the alloy. However, our results show that at low temperatures, the configurational entropy of these materials is largely independent of N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document}, and the assumption described above only holds in the high-temperature limit. This suggests that alloy design guidelines based on the ideal entropy of mixing require further examination.

High entropy alloys (HEAs) are a promising class of multicomponent alloys in which all the constituent elements are present in near equal proportions and the solid solution phase is expected to be favorable compared to ordered phases. [4,16,17] This is very different from traditional alloys, in which solid solution phases are typically formed at the terminal side of a phase diagram (i.e. typically one element has a much higher concentration than the rest) and intermetallics or ordered phases are formed when the compositions of the constituent elements are high. [17] Given the large number of possible high entropy alloys, methods are needed to narrow the scope of the search. Traditionally, such methods use the Hume-Rothery rules as a guideline for determining a potential alloy's stability; that is, they examine the atomic size mismatch, electronegativities and enthalpy of mixing between constituent elements [10] While atomic size mismatch in many HEAs are within the 15% upper limit of the Hume-Rothery rules, HEAs such as NbVTiZr, WNbMoTa, and WNbMoTaV [7] have a large mismatch between atomic sizes and yet exist as solid solution alloys and exhibit very high ductility and strength. Indeed, as Troparevsky et al. have noted, '' [...] consideration of the enthalpy of formation of the HEA alone is not a predictor of its stability as a single-phase solid solution and [...] the discovery of new single-phase compounds requires more than Hume-Rothery arguments''. [10] Thus, HEAs are quite different from traditional alloys and their particular microstructure and properties offer many potential applications, such as tools, molds, dies, and mechanical and furnace parts that require high strength, thermal stability, and wear and oxidation resistance. [16] It is conjectured that the high entropy of mixing, DS mix , in HEAs can significantly lower the free energy of mixing, meaning that the propensity to form clusters or to phase segregate is diminished. [18] DS mix is often approximated by the ideal entropy of mixing (Eq 1), where N is the number of elements and x i is the atomic fraction of species i in the alloy.
Here, we note that DS ideal is an upper bound and often DS mix DS ideal . Additionally, it is apparent from this equation that DS ideal is independent of the chemical species and only dependent on the relative fractions of each constituent species.
Since DS ideal increases with an increase in N, Yeh [17] conjectured that solid solution alloys can be stabilized by simply increasing the number of alloying elements. However, by analyzing the properties of high entropy alloys published in the literature, Senkov et al. [8] concluded that ''the likelihood of forming equimolar alloys that contain only disordered solid solution phases decreases with an increase in the number of constituents''. Given that a large number of HEAs are known to exhibit phase segregation, the authors also suggested that configurational entropy alone does not govern the formation of solid solution alloys. This is also evident from the experiments done by Li et al. [3] . The authors prepared multiple HEAs, each containing 5 alloying elements, that exhibit drastically different microstructures: a few of these alloys are solid solutions, but the remaining alloys contain a mixture of BCC and FCC phases, or intermetallics. Since the ideal mixing entropy is same for these alloys, the presence of a diverse set of microstructures clearly suggests that either the enthalpy, vibrational or electronic entropy contributions can play a dominant role. However, it is also possible that actual entropy of mixing, DS mix , for these alloys can be significantly different from DS ideal . Yeh and co-workers have conjectured that configurational entropy contributions can be greater than contributions from vibration, electronic or magnetic entropy. [16,17] This is in line with the analysis of 4 element alloys done by Gao et al. using a combination of CALPHAD and ab-initio density functional theory (DFT) based calculations. However, the small system sizes accessible in DFT calculations can drastically affect the many-body correlations and the entropy of mixing. [2] On the other hand, in the case of stainless steel (consisting of Cr, Fe, Mn, Mo, and Ni), Wang et al. [14] reported that vibrational entropy contributions are higher than TDS ideal at temperatures above 600 K. At temperatures less than 600 K, even though the configurational entropy calculated is higher than the other contributions, it is not clear if this is valid for the actual configurational entropy contribution, which is supposed to be lower than DS ideal . Thus, a systematic analysis of the temperature dependence of the short range order and the configurational entropy is also important because many solid solution HEA design guidelines that use DS ideal have been proposed in the literature. For example, based on their analysis of about 100 multi-component alloys for which the enthalpies of mixing lie in the range of À15 DH mix 5 kJ =mol, Yang and Zhang [15] suggested that the parameter, X ¼ T m DS ideal DH mix can be used to identify solid-solution HEAs. But, Tsai et al. [11] found that predictions made using these parametric models can be incorrect and suggested using the enthalpy of mixing and a parameter that captures local distortion. Similarly, experiments done by Singh and Subramaniam, Pradeep et al., and by Otto et al. also point to the limitations of alloy design guidelines simply based on the ideal entropy of mixing for the identification of solid solution alloy compositions. [5,6,9] To help clarify the situation, we use first principles density functional theory and cluster expansion calculations to perform a systematic analysis of the configurational entropy in three refractory high entropy alloys (NbTiVZr, HfNbTaTiZr, and AlHfNbTaTiZr) and two conventional alloys (stainless steel consisting of Cr, Fe, Mn, Mo, and Ni and Ni 3 Al with Co, Fe, and Ti impurities). The refractory HEAs selected here represent a sampling of the more promising HEAs, and others will be studied in the future. In this study, we are interested in understanding the following issues: (a) Since DS ideal increases with N, does this mean that DS mix also increases monotonically with N? (b) Is the entropy of mixing, DS mix sensitive to the temperature? (c) Is the short range ordering of elements important in these alloys? If so, then at what temperature do they become important?
These questions will be answered in the concluding sections of this paper. A closely related question is whether the short range order is indeed zero at the melting point. Answering these questions can have important ramifications on the various solid-solution design principles based on DS ideal .
To investigate these issues, we need to accurately capture the effect of chemistry. Therefore, ab-initio density functional theory (DFT) calculations are well suited for this purpose. However, there are two competing effects which combine to make brute-force DFT calculations impractical. On the one hand, proper equilibration of the system requires sampling many configurations. However, decreasing the supercell size enough to make this feasible leads to a supercell too small to sufficiently capture shortrange order. Thus, to overcome this length-scale limitation, we have used the cluster expansion formalism for our analysis.
The cluster expansion formalism for alloys (Eq 2) constructs a linear model which expands the energy of a crystal structure as a series of interactions between different clusters of atoms a. m a is the multiplicity of a particular cluster, h Q i . . .i is a cluster function that satisfies certain properties, and J a is the effective cluster interaction (ECI) coefficient. [ In order to construct a cluster expansion, the coefficients J a must be determined. Typically, this is done by calculating energies of a number of structures, constructing a linear system of equations based on Eq 2, and solving it for the coefficients. [13] The cluster expansions for NbTiVZr and AlHfNbTaTiZr are generated by fitting the coefficients J a in Eq 2 using least-angle regression (LAR) (the quality of fit is determined by calculating the R 2 score for a subset of data set aside as validation data), while the cluster expansion for HfNbTaTiZr is generated through manual optimization of a cluster expansion generated by mmaps, part of the Alloy Theoretic Automated Toolkit (ATAT), where the quality of fit is determined by the leave-one-out cross-validation (LOOCV) score. The final cluster expansions used in the following work are described both by the distribution of included clusters and the relative importance of each cluster, denoted by the magnitude of the ECI (Fig. 1). The final NbTiVZr expansion contains 2-body clusters up to 11.3 Å , 3-body clusters up to 6.3 Å , and 4-body clusters up to 3.8 Å . The HfNbTaTiZr expansion contains 2-body clusters up to 8.2 Å , 3-body clusters up to 5.4 Å , and 4-body clusters up to 5.4 Å . Finally, the AlHfNbTaTiZr expansion contains 2-body clusters up to 7.6 Å , 3-body clusters up to 5.4 Å , and 4-body clusters up to 3.8 Å .
The short-range order parameters used here are based on the Warren-Cowley Short-Range Order (SRO) [1] parameters that track the variations in probability of finding different elements pairs within a specified cut-off distance. Short-range order parameters like these facilitate studying the evolution of local correlations as a function of temperature. An SRO parameter of 0 signifies perfect disorder, while a positive SRO parameter (up to a maximum of 1) signifies repulsion and a negative parameter signifies attraction. In this work, to analyze the effect of temperature on the short range ordering of alloying elements, SRO parameters are calculated by using the first nearest neighbors of an atom.
Before determining SRO parameters, one must first determine the temperature and composition range where the system exhibits a single-phase thermodynamic equilibrium. To this effect, grand-canonical Monte Carlo simulations are run at various temperatures between 1000 and 2500 K. In order to run Monte Carlo simulations under the grand-canonical ensemble, the concentration of each element is no longer fixed, and the chemical potential of one or more elements is varied during the course of the trajectory. This leads to changes in the chemical composition of the supercell, and sudden jumps in composition indicate a phase transition (assuming the step in chemical potential is small enough). The range of chemical potentials of interest is determined by doing a wide sweep at 2500 K until the end members are encountered and narrowing the range until the equiatomic composition is bracketed in the chemical potential grid. For NbTiVZr, the chemical potentials range from -0. . It is clear that the equiatomic phase is quite unstable at 1000 K in both alloys. For example, there is a phase transition from a roughly equiatomic alloy phase to a Ta-rich phase in AlHfNbTa-TiZr at 1000 K as the chemical potential of Ta is changed, which can be seen in the jumps in the ''Ta'' curve in Fig. 2c-e. Those jumps in composition are not seen in Fig. 2h-j, showing that the equiatomic phase is stable at that temperature. Figures like the those in Fig. 2 for intermediate temperatures can be found in the supplemental materials. The NbTiVZr equiatomic phase stabilizes at 2500 K, while the AlHfNbTaTiZr equiatomic phase stabilizes at 2000 K. This means that SRO parameters at 2500 K and above are meaningful, and the SRO parameters of these alloys at 2500 K as a function of cut-off distance (as mentioned in the preceding paragraph) are given in the supplemental materials. In order to understand the effect of temperature on local chemical ordering, the cluster expansions generated above are utilized to run Monte Carlo simulations in the canonical ensemble. Supercells ranging from 432 atoms (6 Â 6 Â 6 supercell) to 1024 atoms (8 Â 8 Â 8 supercell), large enough to capture the phase segregation and chemical ordering mentioned previously, are used to obtain a set of 4000-5000 equilibrium structures at each temperature using Monte Carlo. The temperatures explored here range from 250 to 3000 K for NbTiVZr and AlHfNbTaTiZr and 100 to 5000 K for HfNbTaTiZr. To ensure that the phase space corresponding to the distribution of alloying elements is properly sampled, multiple Monte Carlo trajectories from different initial structures are launched and allowed to equilibrate.
The structures obtained from the canonical Monte Carlo trajectories are analyzed and mean energies are calculated at different temperatures. Specific details of the number of structures used and the corresponding uncertainty in the mean energies for the different systems are shown in the supplementary material. Next, free energies are calculated using the thermodynamic relationship o bF ð Þ ob ¼ hEi. This relationship is valid when the number of atoms in the system and the volume are held constant. Here, k B is Boltzmann's constant and b ¼ k B T ð Þ À1 . For our purpose, this thermodynamic relationship takes the following form: However, in order to utilize Eq 3, the free energy at one of the end points must be known. Given that the entropy of a purely disordered alloy can be described by Braggs-Williams entropy (Eq 1) (where x i is the atomic concentration of a particular species in the alloy), the free energy at the high temperature limit can be computed and the integral in Eq 3 can be rewritten as Eq 4, where T 2 \T 1 and T 1 is the showing the V and Zr concentrations along the same trajectory. Jumps in concentrations signify a phase transition. For example, in (b), the sudden jumps from an Zr-rich phase to a Zr 0:75 V 0:25 alloy phase to a V-rich phase happen over minute changes in chemical potential, suggesting a highly unstable Zr 0:75 V 0:25 phase. Likewise, the sudden jump from a Nb-Zr rich phase to a Nb-V rich phase happens as the chemical potential of Nb is changed. In AlHfNbTaTiZr, there is a sudden jump from a phase rich in all 6 elements (slightly favouring Al and Zr) to a Ta-rich phase as the chemical potential of Ta is changed.
There is similarly a jump from a Nb-Ta rich phase to a Al-Hf-Nb-Ti-Zr rich phase as the chemical potential of Al is changed and a jump in the opposite direction as the chemical potential of Nb is changed. At 2500 K, both equiatomic alloys are stable, as shown in (f-j)-all composition curves are smooth, indicating a lack of phase transitions at this temperature. Further calculations at intermediate temperatures are provided in the supplementary materials temperature at which the short-range order parameters are 0.
This methodology therefore allows for calculating the free energy at a given temperature by discretizing the integral in Eq 4-the scheme used here is the canonical trapezoidal integration scheme. In the case of NbTiVZr and AlHfNb-TaTiZr, T 1 is taken as 3000 K; however, in the case of HfNbTaTiZr, the alloy has not yet become fully disordered at that temperature and only reaches the desired state around 3500 K. Thus, the configurational entropy per atom in units of Boltzmann's constant k B can be obtained as a function of temperature using Eq 5.
The curves in Fig. 3 are calculated for the configurational entropy of each of the high-entropy alloys in question. The configurational entropy of HfNbTaTiZr approaches the ideal entropy of 1:61 kB = atom around 3500 K, while the configurational entropies of NbTiVZr and AlHfNbTaTiZr approach the ideal entropies of 1:39 kB = atom and 1:79 kB = atom , respectively, around 3000 K. Interestingly, HfNbTaTiZr converges to the ideal entropy only around 3500 K, a temperature far greater than the material's melting temperature. The SRO parameters of this alloy also seem to confirm that at the melting temperature, the SRO parameters are still significantly far from zero in many cases. This is not the case for NbTiVZr and AlHfNbTaTiZr.
It is also instructive to compare the shape of the configurational entropy curve of these equiatomic HEAs to those of an Fe-and Ni-rich steel with the composition of Cr 0:11 Fe 0:48 Mn 0:09 Mo 0:12 Ni 0:21 and a dopant-enriched Ni 3 Al alloy. The configurational entropy curve of the steel approaches the ideal entropy of 1:38 kB = atom around 2900 K, and its ideal entropy is similar to that of NbTiVZr and far less than those of HfNbTaTiZr and AlHfNbTaTiZr. Meanwhile, the configurational entropy of the dopant-enriched Ni 3 Al alloy approaches the ideal entropy (0:91 kB= atom ) around 3500 K. The ideal entropies of all of the high-entropy alloys are far larger than that of the dopant-enriched Ni 3 Al due to the equal concentrations of all of the elements in the HEAs. Further, in NbTiVZr and AlHfNbTaTiZr, the ideal entropies are reached at lower temperatures than in the case of HfNbTaTiZr, which may point to the stability of the single-phase solid solution at lower temperatures in these systems.
All the temperatures reported above are for a cluster expansion model that does not include vibrational entropy contributions. The well-documented fact that including such effects tends to reduce the temperature scale [12] should be taken into consideration when interpreting these results.
This work calculates the configurational entropy of NbTiVZr, HfNbTaTiZr, and AlHfNbTaTiZr by combining the cluster expansion formalism, Monte Carlo simulations, and thermodynamic integration. The degree of disorder in the material at each temperature is quantified by examining the short-range order parameters. From these results, it is clear that the configurational entropy in these high-entropy alloys is sensitive to temperature, but only weakly  The data that support the findings of this study are available from the corresponding author upon reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/. For each alloy, the following are given: lattice parameter used here, enthalpies at two different temperatures, entropies at two different temperatures, and ideal entropy