First-principles prediction of high-entropy-alloy stability

High entropy alloys (HEAs) are multicomponent compounds whose high configurational entropy allows them to solidify into a single phase, with a simple crystal lattice structure. Some HEA's exhibit desirable properties, such as high specific strength, ductility, and corrosion resistance, while challenging the scientist to make confident predictions in the face of multiple competing phases. We demonstrate phase stability in the multicomponent alloy system of Cr-Mo-Nb-V, for which some of its binary subsystems are subject to phase separation and complex intermetallic-phase formation. Our first-principles calculation of free energy predicts that the configurational entropy stabilizes a single body-centered cubic (BCC) phase from T = 1,700K up to melting, while precipitation of a complex intermetallic is favored at lower temperatures. We form the compound experimentally and confirm that it forms as a single BCC phase from the melt, but that it transforms reversibly at lower temperatures.


Introduction
Exploring new materials with outstanding properties is an eternal pursuit of scientists and engineers. High entropy alloys (HEAs) constitute a newly-emerging class of materials that offer enhanced mechanical and corrosion properties for practical applications, while presenting challenges to the scientist making predictions in the face of compositional complexity [1][2][3][4][5][6][7][8][9][10][11][12][13]. Alloys of two or three species often phase-separate, and when they form a compound, it is often a complex ordered structure. Thus, HEAs are remarkable for forming single phases that are, moreover, simple lattice structures. Trial-and-error experimental approaches to test the phase stability of HEAs are costly and time-consuming, while previously-reported phase formation rules are empirical and susceptible to the alloy systems [14,15], motivating a first-principles theoretical approach to accelerate the design of HEAs and control their properties through the prediction of phase stability.
Here we illustrate these ideas by demonstrating phase stability in the quaternary Cr-Mo-Nb-V alloy system. Some of its binary and ternary subsystems are subject to phase separation and complex intermetallic phase formation. Utilizing fully first-principles calculations of the free energies of the full quaternary and its competing binary and ternary phases, we predict that the configurational entropy stabilizes a single body-centered cubic (BCC) solid solution at high temperatures, but that it will precipitate a complex intermetallic Laves phase at lower temperatures. Experimentally, we verify formation of the HEA from the melt, confirm that precipitation of the Laves phase at low temperatures, and observe that it reverts to a single BCC phase upon annealing at high temperatures, thus illustrating the entropic stabilization principle of HEAs.

Quaternary alloy system
Consider the four-component alloy system of Cr-Mo-Nb-V. All four elements take bodycentered-cubic (BCC, Pearson type cI2, Strukturbericht A2) structures. Taking elements pairwise, we note that mixing enthalpies for the BCC solid solution are positive for Cr-Mo, Cr-Nb, and Nb-V (see Fig. 1), indicating low-temperature phase separation, while the remaining cases, Cr-V, Mo-Nb, and Mo-V, have negative enthalpies [16,17], indicating stability down to low temperatures. Experimentally, Cr-Nb has a complete miscibility gap at all temperatures below melting, while phase separation has not been seen experimentally in the case of Cr-Mo [18] and Nb-V [19], though it is predicted [20,21]. Cr-Nb and Nb-V are predicted [22] to form Laves phases of types cF24 (Strukturbericht C15) and hP12 (Strukturbericht C14). Both Laves phases are observed experimentally in Cr-Nb [23], while neither has yet been reported in Nb-V. Experimental investigations of the ternaries show that Cr-Mo-Nb phase-separates [24] and Cr-Nb-V forms a C15-Laves phase [25], while Cr-Mo-V and Mo-Nb-V become complete solid solutions [18]. No quaternary phases have been reported experimentally.
Inspection of Fig. 1 offers a clue to why experiments on Cr-Mo and Nb-V fail to observe the predicted miscibility gaps, while separation is seen in Cr-Nb. The mixing enthalpies are relatively low in Cr-Mo and Nb-V, compared with Cr-Nb, so their estimated critical temperatures for phase separation are low (1,213K and 653K, respectively [20,21]). Meanwhile, the melting temperatures of Cr-Mo and Nb-V are relatively high (2,093K and 2,133K, respectively, compared with 1,893K for Cr-Nb). Since atomic diffusion far below the melting temperature can be slow, phase separation will be difficult to observe. The problem is even more severe for complex structures, potentially explaining the failure so far to observe the predicted NbV 2 C14-Laves phase.
Turning to the full quaternary, we will test its stability at high temperatures by calculating its free energy relative to the most-likely competing phases, namely, the coexistence of a Morich BCC medium entropy alloy with a CrNbV-rich Laves phase. A working hypothesis is that the single phase is stabilized by a large configurational entropy of mixing, and that the entropy is maximized for the simple BCC lattice structure in which all sites are equivalent. However, we expect the lower enthalpy of the competing phases to predominate at lower temperatures where the entropy is less relevant, leading to precipitation of a ternary CrNbV Laves phase as temperature drops. At very low temperatures, we anticipate complete phase separation into four coexisting phases.

Free energy calculation
For each phase, we calculate the Helmholtz free energy as a sum of discrete configurational, vibrational, and electronic contributions [22,26], F (x, T, V ) = F c +F v +F e . The composition variable, x, is expressed as a four-component vector in our composition space. Because metallurgical experiments are generally carried out at low pressure (P ≈ 0), the Gibbs free energy, G(x, T ) = min V (F + P V ), of each phase is simply the Helmholtz free energy, F , evaluated at its minimizing volume. If we have a mixture of phases, i, each having a mole fraction γ i , the total free energy is Consider the phase separation of a CrNbV Laves phase from a CrMoNbV BCC solid solution, expressed as where γ ∈ (0, 1) represents the progress of the reaction. The fraction of the Laves phase is 3γ/4, leaving a fraction 1 − 3γ/4 of the BCC phase. The BCC phase becomes increasingly Mo-rich, with the composition of x Mo = 1/(4 − 3γ), while the remaining elements have the composition x α = (1 − γ)/(4 − 3γ) for α = Cr, Nb, or V. Thus, the total free energy for a two-phase mixture of Laves with a partially-transformed BCC structure is and G BCC depends on x Mo . Equation (2) holds for the total Gibbs free energy and also for the separate configurational, vibrational and electronic components.
To calculate the configurational free energy, G c (x, T ), we apply our hybrid Monte Carlo/molecular dynamics method (MCMD [27]) to anneal the short-range chemical order. Over the temperature range of interest (T = 1,000K through melting around 2,133K), the chemical order was found to be moderate and nearly temperature independent. In particular, the configurational entropy, remained within 2% of the ideal value [28]. Thus we take the mixing entropy of the BCC solid solution as ideal, S c = k B α x α log x α . We calculate the configurational enthalpy of formation, relative to pure elements, ∆H c , by quenching our simulated structures to T = 0K through the relaxation of atomic coordinates and lattice parameters. The resulting enthalpies (see Fig. S1a) are nearly independent of the annealing temperature but exhibit a strong variation with respect to composition.
Applying MCMD to the C15-Laves phase of CrNbV reveals that Cr and V readily substitute on the 16d Wyckoff sites but they do not mix to any appreciable degree with Nb, which remains on the 8a sites. Hence, we take the entropy of CrNbV as S c = (2/3)k B log 2. The total Gibbs configurational free energy of formation is the sum of the BCC and Laves terms weighted by their respective fractions, as in Eq. (2).
The vibrational free energy is calculated in the harmonic approximation from the vibrational density of states, g(ν), via the single vibrational mode free energy where g(ν) is obtained by the direct method [29] by means of force constants that we determine using density functional perturbation theory at T = 0K. Figure 2a compares vibrational densities of states for the equiatomic BCC phase and the Laves phase. The Laves phase has 4 some high-frequency phonons that are absent in the BCC phase, while the BCC phase has a small but persistent excess of low-frequency phonons, lowering F v BCC relative to F v Laves ( Fig. 2c and d). These differences reflect the tetrahedral close-packed structure of the Laves phase, compared with the relatively-open BCC structure. To model phase separation, we require the composition dependence of ∆F v . We find that the differences in the BCC vibrational free energy with respect to composition, x Mo , are small, compared with the difference between the BCC and Laves phases, and hence may be treated as composition-independent. The electronic free energy is calculated from the single electron density of states, D(E), predicted at T = 0K by density functional theory (DFT), and the Fermi occupation function, We determine the electronic free energy, Figure 2c shows that the electronic densities of states near E F differ significantly between BCC and Laves phases. This difference has a significant composition dependence, as can be seen within a rigid band model from the irregular shape of D(E) near E F (increasing x Mo increases E F ). Because we utilize the harmonic approximation, thermal expansion is neglected, and the vibrational Helmholtz free energies, F v and F e , can be added directly to the configurational Gibbs free energy, G c . Thus, combining the above results, we have obtained the free energy, G, expressed as a function of transformation, γ, or composition, x Mo , and temperature, T. Figure 3 graphs this function, revealing minima with γ > 0 for T<1,700K, implying phase separation. From T = 1,700K up to melting above T = 2,100K, however, the minima occur at the boundary, γ = 0 (i.e., x Mo = 1/4), indicating that a single BCC phase is stable.

Experimental results
To test this prediction, we formed a sample by arc melting and casting. The as-cast sample was found to be a single BCC phase, with the characteristic of dendrites and interdendrites (Figs. 4 a and d). The single phase remained intact under 21 days annealing at temperatures up to 1,273K, but exhibited precipitation of a C15-Laves phase under annealing of 3 days at 1,473K ( Fig. 4b and d). Subsequent low-temperature annealing failed to restore the single phase, indicating that the phase-separated state is indeed stable at low temperature.
The C15 phase disappeared when the 1,473K-annealed sample was annealed at 1,673K for 12 hours (Fig. 4c and d). Meanwhile, equiaxed grains appeared instead of the original dendrites and interdendrites, indicating the pristine BCC phase was restored at 1,673K. This reversible transformation (BCC + Laves phases at low temperatures, and full BCC phase at high temperatures) is consistent with our theoretical prediction. Further annealing at 1,873K revealed additional impurity-stabilized phases (see Fig. S2b), including a novel CCr 3 Nb 3 carbide of Pearson type cF112. Table 1 reports the compositions of the C15 and BCC phases in the 1,473K-annealed state. Figure S3 shows the distribution of constitutive elements after annealing the 1,473K-annealed sample at 1,673K for 12 hours. No obvious segregation was observed within the BCC phase in the 1,673K-annealed state. Figure S2 presents the microstructures after annealing the alloy at 1,273K and 1,873K. It can be seen that the alloy still keeps the characteristic of dendrites and interdendrites under the 21 days annealing at 1,273K (Fig. S2a). After annealing at 1,873K, multiple phases appeared in the alloy, including BCC, CCr 3 Nb 3 , and an unidentified needle-like phase, (Fig. S2b).
The chemical distribution of the CCr 3 Nb 3 phase in the 1,873K-annealed state was further detected by EDX, as shown in Fig. S4. An Fd3m symmetry was identified by TEM SAED (Fig. 5), suggesting a phase of the structure type of NiTi 2 with a lattice parameter a = 11.7Å. As demonstrated below, this phase is stabilized by carbon impurities and can incorporate silicon substitution in place of Cr.

Discussion
So far we focused on the primary competitor to the HEA, namely, the separation into a Laves phase plus a BCC medium entropy alloy (MEA). In principle, according to the Gibbs phase rule, the system could contain up to four independent phases at low temperature. We employed the Alloy Theoretic Automatic Toolkit cluster expansion [30] to generate a variety of candidate low-temperature binary, ternary, and quaternary configurations based on decorations of the BCC lattice. Based on enthalpies, we predict that the equiatomic CrMoNbV decomposes at low temperatures [22] into four competing phases, Cr 2 V.tI6, CrNbV.cF24, MoNb.oC12, and Mo 4 V 3 .hR7 subject to the stoichiometric relationship 10 CrMoNbV → 3 Cr 2 V + 6 MoNb + 4 CrNbV + 1 Mo 4 V 3 .
Relaxed lattice parameters and Wyckoff coordinates of these structures are given in Table 1.

6
Assuming that all four phases are stoichiometric, we neglect the configurational entropy and set G c = ∆H c . Combining G c with F v and F e , and weighting the contributions according to Eq. (6), we obtain the function displayed in Fig. S5 labeled as G 4 . Notice that G 4 crosses the single-phase equiatomic HEA free energy, G 1 , around T = 800K, indicating a transition from the HEA to the 4-phase mixture at low temperatures. On the same figure, we plot the free energy, G 2 , of the 2-phase mixture, and this value lies below G 1 and G 4 over the range of T = 500-1,700K. Accordingly, the HEA does not transform directly to the 4-phase mixture upon cooling, but instead exhibits a 2-phase region in between. The upper transition at T = 1,700K is continuous, with G 2 joining G 1 tangentially. This trend is reflected in the continuity of composition variables, as illustrated in Fig. 3b of the main text. The lower transition is discontinuous within the present model, with G 4 crossing G 2 transversally at 500K. However, because it corresponds to chemical ordering on the BCC lattice, the actual behavior may be more complicated with multiple intervening phases and possibly continuous transitions. These transitions occur at low temperatures, making it doubtful that they can be observed experimentally.
Because CrNb 2 exhibits a high-temperature C14-Laves phase (Pearson type hP12), and the structure is predicted to be stable at low temperatures in NbV 2 , we checked the free energy of CrNbV as a C14-Laves phase. At this equiatomic composition, the free energy of C15 remains below the free energy of C14, for all temperatures below melting.
As mentioned above, TEM identified a Cr-Nb-rich impurity phase of space group F d3m (group #227), and proposed the structure type of NiTi 2 (Pearson type cF96). This structure proved energetically unfavorable as a Cr-Nb binary. However, carbon impurities at octahedral interstitial sites stabilize the phase in compositions CCr 3 Nb 3 and CCr 2 Nb 4 (both Pearson type cF112). Thus, we predict the existence of a previously unknown stable carbide in the C-Cr-Nb alloy system. Table 2 gives crystallographic details of the CCr 3 Nb 3 .cF112 phase, which we predict to be stable against all known binary and ternary competing phases, with the formation enthalpy of ∆H = −220 meV/atom, relative to pure elements. Substituting a single Si atom at a Cr2 site results in a quaternary that is stable against all known competing phases, with the formation enthalpy of ∆H = −275 /meV/atom. CCr 2 Nb 4 is similar to CCr 3 Nb 3 , but with Nb occupying the 16c site.
In summary, we have predicted the existence of a single-phase refractory HEA. This phase is stabilized at high temperatures by the configurational entropy of mixing despite the immiscibility of constituent elements and a competing complex intermetallic phase, which leads to phase separation at low temperatures. A reversible phase transformation restores the entropically stabilized BCC phase at high temperatures, illustrating one of the foundational principles of HEAs. Our theoretical methods, combining the quantum-mechanical total-energy calculation with statistical mechanics to predict free energies, can be applied to many problems in alloy design. 7

Enthalpy of formation
First-principles calculations are performed using projector augmented wave potentials [31,32] in the generalized gradient approximation [33] as implemented in the Vienna Ab-initio Simulation Package (VASP [34]). Structures were fully relaxed in lattice parameters and internal coordinates at T = 0K using the energy cutoff of 400 eV and k-point densities sufficient to achieve convergence of 1 meV/atom or better. Hybrid Monte Carlo/molecular dynamics simulations (MCMD [27]) generated BCC solid-solution configurations over a series of Mo-rich compositions with an appropriate short-range order. These calculations used default energy cutoffs and a single k-point in systems of size N = 128 atoms. Figure S1a shows the resulting enthalpies of formation, relative to pure elements, which we find fit well to a quadratic polynomial.

Electronic free energy
We utilize the same series of Mo-rich configurations to evaluate electronic densities of states and electronic free energies. The integral in Eq. (5) can be expanded in power series in temperature [22]. Expressing Odd powers drop out, and truncating at T 4 is sufficient even at quite high T . Here we fit coefficients, a and b, to F e (T ) = a(T /1000) 2 + b(T /1000) 4 over the range of T = 0-2,000K. The composition dependence of F e is irregular, as can be inferred from the double minima in D(E) (see Fig. S1). Hence, we model F e by fitting a and b to 6th-order polynomials of x Mo , as illustrated in Fig. S1c.

Vibrational free energy
The vibrational free energy is calculated in the harmonic approximation from the vibrational density of states, g(ν). This is obtained from the force constant matrix, Φ ij = ∂ 2 U/∂R i ∂R j , which VASP determines using density functional perturbation theory at T = 0K. Our cal-culations utilize a plane-wave basis with the energy cutoff of 400 eV, "accurate" precision to avoid wrap-around errors, and an additional support grid for augmentation charges.
The dynamical matrix at the wavevector q is [26,29] We diagonalize D ij (q) to obtain vibrational frequencies, ν(q), and sample the full Brillouin zone to obtain the density of states, g(ν). Representative densities of states, calculated in 27-atom cells for various compositions, are shown in Fig. S6. Notice that g(ν) is relatively insensitive to the composition, but that the BCC DOS differs significantly from the C15-Laves phase. Because we utilize the harmonic approximation, thermal expansion is neglected, and the vibrational Helmholtz free energy, F v , can be added directly to the configurational Gibbs free energy, G c . In the high-temperature range of interest, the free energy varies classically, with quantum effects relevant only for establishing offsets in the energy and entropy. Thus, we can fit F v (T ) with high precision using only a single parameter, an effective Debye temperature, Θ D , chosen to match the quantum-based free energy in the classical regime. We set with D 3 (z) = 1 − 3z/8 + z 2 /20 + · · · the third order Debye function [35]. We set Θ D =393K for the BCC HEAs, and Θ D =420K for the Laves phase, as appropriate values to reproduce the high temperature free energy.

Sample preparation
The CrMoNbV alloy was fabricated by arc-melting the constituent elements (purity > 99.9 weight percent), then drop casting into a water-cooled copper hearth. To achieve a homogeneous distribution of elements in the alloy, the melting and solidification processes were repeated five times. The annealed samples of the CrMoNbV alloy were encapsulated in quartz tubes, filling with high-purity argon after vacuuming several times. The samples were heat treated at 1,273K for 21 days and 1,473K for 3 days, respectively, and then water quenched to validate the stability of the BCC phase. After obtaining the 1,473K-annealed samples, another annealing treatment at 1,673K for 12 hours was performed using a 1,473Kannealed sample to verify the stability of the C15-Laves phase. Further annealing at 1,863K revealed a novel carbide stabilized by impurities.

Microscopy
A LEO Gemini 1525 field emission scanning electron microscopy (SEM) coupled with energydispersive X-ray spectroscopy (EDX) was utilized to characterize the microstructure and composition of this alloy. Transmission-electron microscopy (TEM) was conducted to identify the new phase in the 1,863K-annealed state, using ZEISS LIBRA 200 HT FE MC coupled with EDS. SEM specimens were initially polished with 1200-grit SiC paper and, subsequently, vibration polished using 0.05 µm SiC liquid for the final surface clarification. Focused ion beam (FIB) milling was used for the preparation of TEM specimens, targeting the specific phase. The average compositions of each phase were measured at five different locations to ensure compositional accuracy.

X-ray diffraction
Synchrotron X-ray diffraction experiments were carried out at the Argonne National Laboratory, Advanced Photon Source, using the beamline 11-ID-C. The bulk specimens (∼1 mm thick) after the SEM characterization were measured by a beam energy of 111 keV (0.11798Å) with a beam size of 0.5 × 0.5 mm.