First-Principles Prediction of Incipient Order in Arbitrary High-Entropy Alloys: Exemplified in Ti 0.25CrFeNiAl x

Multi-component solid-solution alloys, including high-entropy alloys, experience segregation or partially-ordering as they are cooled to lower temperatures. For Ti<sub>0.25</sub>CrFeNiAl<sub>x</sub>, experiments suggest a partially-ordered B2 phase, whereas CALculation of PHAse Diagrams (CALPHAD) predicts a region of L2<sub>1</sub>+B2 coexistence. We employ first-principles KKR-CPA electronic-structure to assess stability of phases with arbitrary order having a specific Bravais lattice and composition, as it handles chemical disorder by direct configurational averaging during density-functional theory charge self-consistency. In addition, KKR-CPA linear-response theory is used to predict atomic short-range order (SRO) in the disorder phase, which reveals potentially competing long-range ordered (LRO) phases. With Warren-Cowley pair correlations defined relative to the diffraction lattice, SRO can be analyzed by concentration-waves (site-occupation probabilities) for any partially-ordered cells (including estimated energy gains) that can then be assess directly by KKR-CPA calculations. Our results are in good agreement with experiments and CALPHAD in Al-poor regions (x ≤0.75) and with CALPHAD in Al-rich region ($0.75 ≤x≤1), and they suggest more careful experiments in Al-rich region are needed. Our first-principles KKR-CPA electronic-structure and SRO predictions with additional concentration-wave analysis are shown to be a powerful and fast method to identify and assess competing phases in complex solid-solution alloys.

First-principles prediction of incipient order in arbitrary high-entropy alloys: exemplified in Ti 0.25 CrFeNiAl x

Introduction
Multi-principal-element alloys or complex solid-solution alloys, of which high-entropy alloys (HEAs) are a subset, have established a new paradigm in alloy design [1], but there are many fundamental science questions unanswered. These HEAs have a huge unexplored composition space that in some cases searches have lead to the discovery of novel alloys with unusual properties [1,2,3,4,5,6,7,8,9,10,11,12]. Short-range order (SRO) is one such key property, which, nonetheless, remains less explored in HEAs [13,14,15,17,18,19,20]. The site-pairwise correlation between atoms in HEAs is given by the Warren-Cowley SRO parameters, which are defined with respect to the average diffraction lattice, and indicate the incipient order that competes within the alloy as temperature is decreased. This SRO often indicates the expected partially or fully-ordered state from an order-disorder transistion, where two distant atoms are connected by nonzero long-range order (LRO) correlations as described by LRO parameters [22,23,24,25,26,27,28]. By controlling SRO and LRO, the processing and properties of advanced materials can be manipulated [21]. Thus, predicting SRO in HEAs and assessing its electronic origins (e.g., band-filling, Fermi-surface, atomic-size [i.e., band-width], or charge-transfer) is of great importance.
Generally, SRO assessed in high-temperature, homogeneously disordered HEAs is not biased by an a prior choice of partial LRO structures (unit cells), which may be a problem with more traditional total-energy methods [32,33]. Singh et.al. [19,34,35,36] presented a DFT-based KKR-CPA thermodynamic linear-response theory (a first-principles Landau theory) that predicts all SRO pairs simultaneous (as done for displacement modes, i.e. phonons) and provides the incipient ordering in HEAs [19]. The SRO yields the eigenvectors (chemical modes) and the maximum temperature for the most unstable mode(s) in the HEAs, i.e., the spinodal temperature (T sp ), with respect to the formation of short-wavelength composition waves that characterizes potential ordered structures, i.e., site-probabilities are modulated with a wave-like periodicity [25,26,27,28]. Besides specifying SRO, a great strength of the concentration-wave approach is predicting various competing sub-lattice occupation probabilities for partial LRO phases, which can then be used to assess directly the stabilities of low-temperature competing phases with respect to random alloy [19].
Here, KKR-CPA formation enthalpies and SRO are predicted to assess Ti 0.25 CrFeNiAl x (mole fractions). Our first-principles results for x ≤ 0.75 compare well with the experiments of Liu et.al. [37], both show A2 (bodycentered cubic, bcc) as the most stable phase. However, CALPHAD suggests competing B2 (CsCl) and L2 1 (Heusler) phases in Al richer-regions (0.75 ≤ x ≤ 1). We use the KKR-CPA linear-response SRO and, based on the analytically formulated concentration-wave site occupations, establish potential competing partially-LRO unit cells and then use KKR-CPA to compare directly the partially-LRO states relative to the disordered phase. We find agreement with the CALPHAD for x > 0.75, predicting competing B2 and L2 1 phases, suggesting more careful experiments are needed.

Computational method
DFT-based Approach: KKR-CPA is a Green's function electronic-structure approach to permit charge self-consistency and configurational averaging to be done concomitantly (averaged Green's functions are related directly to observables [38]). Our KKR-CPA code handles core electrons relativistically and semi-core and valence electrons scalar-relativistically (no spin-orbit terms) and the coherent potential approximation (CPA) is used to handle chemical disorder and its configurational averaging [39,40,41]. The LibXC library is used to include generalized gradient approximation (PBE) as the exchange-correlation functional [42]. A variational potential zero v 0 is used to yield kinetic energies and dispersions nearing those of full-potential methods [44]. The Atomic Sphere Approximation (ASA) represents each Voronoi polyhedra (VP) in a unit cell but with periodic boundary conditions enforced and VP integrations are used for charge distributions. For self-consistent densities, a 20 (complex) energy point Gauss-Legendre semi-circular contour integration is used [45], and L max (l, m)=3 in spherical-harmonic basis. Brillouin zone (BZ) is sampled using Monkhorst-Pack method [43] with 12×12×12(6) for A1,A2 (A3) meshes. We used 300 k -points along symmetry lines to visualize dispersion.
Chemical SRO: Formally, the Warren-Cowley SRO parameters (pair correlations), i.e., α µν (k; T ) in Laue units for atom pairs µ,ν, are defined relative to average x-ray scattering lattice, as atomic displacements sum to zero on average (by symmetry, for each spatial direction). For linear-response, the second-order variation of DFT-KKR-CPA free-energy with respect to elemental concentrations {c i µ } at sites i,j is performed [16,17,18,19], and the Warren-Cowley SRO are analytically found for any N-component HEA with concentrations {c µ } at temperature (T) to be given by an where β −1 = k B T with k B is the Boltzmann constant and C µν is a constant matrix element [19]. The arbitrary N th atom is used as "host" due conservation of atoms (N − 1 independent fluctuations, i.e., N µ=1 c µ = 1). The pairwise-interchange energy S (2) µν (k; T ) in linear-response is the chemical stability matrix referenced to the homogeneous HEA [18], which reflects the freeenergy cost for all pair fluctuations with {c i µ (k)} [19,34,35,36]. S (2) µν (k; T ) is a symmetric, (N−1)×(N−1) thermodynamically-averaged functional that is dictated by the electronic dispersion of the HEA. It is not a "pair interaction" in a Hamiltonian. S (2) µν (k; T ) reveals the unstable (Fourier) ordering modes with wavevector k 0 (or clustering if k 0 = (000)) and identifies the origin and associated energy cost for the modes [15,19,17,18,46,47,48,16]. It gives a detailed "fingerprint" for the incipient ordering in an arbitrary HEA.
For completeness, we note that Eq. 1 is exact [18,49], if the configurational average of the functional and its variation are handled exactly. However, the CPA is a mean-field (single-site) approximation to the average, and, although it is often adequate, it is not exact. Nonetheless, while it can be made increasingly accurate via a cluster generalization of the CPA [50], a simple Onsager correction to the single-site CPA is sufficient to remove most of the error [18]. Moreover, this correction eliminates the incorrect topology of mean-field phase diagrams and dramatically improve the transition temperatures [51]. In short, the mean-field S (cpa) µν is corrected by requiring that the optical theorem is obeyed by Eq. 1, i.e., the integral of the SRO intensity is conserved, for details see [51] and [18]. Minimally, a site-diagonal where the SRO is implicitly dependent on Λ µν (T ), and must be corrected at each T > T sp , the spinodal temperature (see below). At large temperatures, Λ µν (∞) = 0 and this coupled set of equations can be solved iteratively, for example, by Newton-Raphson methods [16,17,19]. Spinodal Decomposition & Transition Temperatures: The most unstable SRO mode is where α µν (k o ; T > T sp ) has the largest peak at k o for a specific µ-ν pair [19]. An absolute instability to the k o mode occurs below the spinodal temperature T sp defined at [α −1 (k o ; T sp )] µν = 0. The normal modes are eigenvectors of S (2) (k o , T → T sp ) driving divergence in SRO. These modes are obtained for any arbitrary HEA from S (2) µν (k; T ) using a special oblique coordinate transform in a given Gibbs space, see [19] for details and examples. The free energy cost to establish one of these modes vanishes at T sp . Above T sp , all eigenvalues remains positive, costing energy to establishing the concentration modulations in disordered state. Below T sp , the critical eigenvalues establish the anticipated probability distribution of the LRO state, often, if the dominant fluctuation, even if it is a first-order transition. The KKR-CPA direct evaluation of the formation energy E f of the partially LRO state relative to the disordered state establish order-disorder transition tem -3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 peratures, for example, in ordering systems [52,53] and slightly more complicated in segregating systems [54]. Temperatures are thermodynamically consistent T sp ≤ T od if evaluated correctly.
Concentration-Wave Analysis of SRO: With normal modes of eigenvectors from the SRO Eq. 1, a vector n(r) of probabilities for each element to occupy specific sites in a crystal structure for partially-LRO cell or superlattice, i.e., concentration-waves, can be written as [19] Here, c(r) is an (N − 1)-component vector in site occupation probabilities {c µ } in the Bravais lattice of the N-component homogeneous HEA reference. Whereas, n(r) depends on the type of order and real-space site coordinates, dictated by the LRO parameters η(T ). The sum s runs over the "stars" (inequivalent k that define the order), j s (equivalent k js in the s th -star), and σ (eigenvector branch of the free-energy quadric). The other quantities are LRO parameter η σ s (T ) for the σ th branch and s star; ν σ is (N − 1)component vector of the eigenmode of stability matrix for the σ th branch; and the symmetry coefficient γ σ (k js ) determined by normalization condition and geometry. See Ref. [19] for more details and examples.
The bottom line is that Eq. 3 represents the possible competing types of ordered superlattices (symmetry broken order) that are incipient in the chemical SRO for a fixed Bravais lattice. Each of the anticipated partiallyor fully-ordered cells can be then assessed using KKR-CPA E LRO f relative to E dis f -a direct calculation for state of order rather than estimated using the SRO only. We showcase this in the next section using the SRO and its eigenvectors below to estimate the competing phases and then do direct KKR-CPA calculation for a given partially LRO to confirm. As an aside, the concentration-wave analysis for small-cell type ordering shows that only up to an 8-component HEA can order if at very specific compositions.
To understand the effects of chemical disorder by Al-doping, we calculate Bloch-spectral function (BSF) of A2 Ti 0.25 CrFeNiAli x using KKR-CPA electronic-structure method for x=0,1 cases, see Fig. 2. The BSF is a generalization of band structure when disorder is present. On adding Al to Ti 0.25 CrFeNi, the BSF is smeared out near the Fermi energy (E F ) due to increased disordered scattering. The BSF broadening can directly be related to the inverse of the electron mean free path [36]. Increased k-space smearing at E F indicates a decrease in electron mean-free path for Fig. 2(b) compared to Fig. 2(a) -the shift in dispersion is clearly visible at energies below E F . The Al doping enhances disorder and lowers bonding states [19], stabilizing the A2 phase, as also shown in Fig. 1. For near equiatomic HEA (x = 1 Al mole-fraction), the negative formation energy indicates the favorability for the mixing of alloying elements.
The configurational entropy, dominated by points S pt = k B N µ=1 c µ ln c µ (−k B ln N in equiatomic alloys with c µ = 1/N ) is a key factor for the forma- tion of single-phase HEAs. Point entropy increases with increasing number of alloying components, which suppresses the formation of intermetallic phases [1,4,58], unless enthalpically dominated by favorable chemical interactions amongst pairs of atoms, which grows as 1 2 N (N − 1) ∼ N 2 . So, as N gets larger, enthalpy (N 2 ) can win over entropy (ln N ). Empirically literature suggests a threshold of S pt ∼ 1.5R for an operational definition of high-entropy alloy [59]. Here, the gas constant R is 8.314 Jmol −1 K −1 and k B = R/N A , where N A Avogadro's number.
While S pt in A2 Ti 0.25 CrFeNiAl x is large, a magnetic entropy also has a contribution, see Table 1. For solid solutions, we estimate the chemical entropy by S pt and the magnetic entropy by ∆S mag = c ν ln(1 + µ ν ), here µ ν are the magnetic moment of ν th element in units of gas constant (R) [setting k B =1]. Clearly, increasing Al content increases the chemical entropy, which saturates at equiatomic Al, however, the magnetic entropy slightly decreases with increasing Al as Fe losses its magnetic character, see Table 1.
To reveal the ordering in Ti 0.25 CrFeNiAl x and find out the plausible reasons for disagreement between experiments and CALPHAD, we chose four sets of alloys, i.e., x = 0.25, 0.50, 0.75, 1.00, permitting a one-to-one com- 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 64 65 parison to experimental composition range of Liu et.al. [37]. As SRO is dictated mostly by the electronic-structure of the alloy, the origins of the observed ordering tendencies in Ti 0.25 CrFeNiAl x can be determined, i.e., all the competing effects (e.g., band-filling, Fermi-surface nesting, atomic size, and charge transfer) can be assessed [19,34]. We apply this approach Ti 0.25 CrFeNiAl 0.5 to exemplify SRO predictions and concentration-wave analysis. In Fig. 3, we plot the interchange energies S (2) µν (b), and SRO (a) at 1.15 of the calculated T sp (794 K). As α −1 (H; T )] µν vanishes at T sp , the SRO diverges at k o = H = (111), indicating B2-type (CsCl) ordering instability. At T sp , the Al-Fe pair drives the instability, whereas Al-Ni pair is the most dominant SRO. Clearly, α µν (H) has a dominant SRO peak for Al-Ni pair (followed by Al-Fe and Ti-Ni); however, the instability in S (2) is driven by Al-Fe pairs (followed by Al-Cr). This odd (but correct) result occurs due the probability sum rule and the matrix-inverse relating α µν (k) with S (2) µν (k) [17,18,19,34,36]. The instability at H in A2-Ti 0.25 CrFeNiAl 0.50 occurs when (at least) one of the eigenvalues of the correlation matrix is maximum, i.e., inverse of the corresponding pair-correlation component of the correlation matrix vanishes. The relative polarization of concentration waves is represented by the eigenvector corresponding to the vanishing eigenvalue of aforementioned correlation matrix in the Gibbs space [19,48,56,57]. We extract the eigenvector corresponding to H and P for Ti 0.25 CrFeNiAl 0.50 at T sp to analyze the B2 and L2 1 type ordering. The eigenvectors at T sp helps determine the LRO parameters to solve the concentration-wave Eq. 3 for occupation probabilities [19].
As already discussed, A2 can order into different lower symmetry structures in going from high-temperature (disorder) phase to low-temperature (partially-ordered) phases. At first, on lowering temperature A2 lattice shows regions of B2 superstructure by lowering symmetry along {111}. Further lowering the temperature can further break the A2 symmetry along { 1 2 1 2 1 2 } (secondary ordering) and may order into L2 1 superstructure depending on the material characteristics.
The A2-to-B2 transition is second-order and given by one LRO parameter, while L2 1 has two: η 1 , η 2 . Our KKR-CPA linear-response predicts SRO, which yields correct second-order transitions, but often dictates expected first-order. However, first-order transitions have a discontinuous η(T ) versus T, requiring a thermodynamic simulation to predict directly. L2 1 occurs due to secondary-ordering, a qualitative prediction can be made if we deliberately break the symmetry of Ti 0.25 CrFeniAl 0.50 by ordering formed by wavevectors k 1 = {111} (the B2 ordering) and k 2 = { 1 2 1 2 1 2 } (the Heusler ordering) using the eigenvectors calculated in primary transition, i.e., in A2-B2 ordering, at T = T sp . L2 1 is then represented by a concentration-wave with vector site probabilities n(r) The last term reflects the additional L2 1 ordering over B2 order, enriching minority components along { 1 With a non-stoichiometric Al 0.50 CrFeNiTi 0.25 the above probabilities lead to a partially-ordered L2 1 state. The calculated KKR-CPA energy gain is E A2 − E L2 1 = 10.14 mRy for partially-ordered L2 1 phase with respect to A2 phase, and E B2 − E L2 1 = 6.64 mRy relative to B2 phase, i.e., E A2 − E B2 = 3.48 mRy. For Ti 0.25 CrFeNiAl x , we also found that L2 1 phase is stabilized with respect to A2 and B2 by E A2 − E B2 = 2.98 mRy and E B2 − E L2 1 = 2.49 mRy for x = 0.75; and, for x = 1, by E A2 − E B2 = 8.33 mRy and E B2 − E L2 1 = 6.11 mRy The spinodal decomposition temperature occurs at T sp = (794; 1802; 1190) K for x = (0.50, 0.75; 1.00), respectively. The increase in Al content increases stability with disorder, i.e., disorder phase remains stable over larger temperature range.
x and that from the SRO is easily understood. Upon breaking symmetry with a LRO parameter into B2, the site charges and dispersion change directly responding to the new symmetry and changes in electronic charge density, whereas using the SRO from linear-response, the dispersion is fixed to the random alloy -hence, it is a cheap and fast estimate -albeit a good one, especially for trends. One may consider using this for fast evaluation of HEA before spending time on more accurate calculations, for example. The SRO can be viewed in Gibbs' composition space (a Barycentric coordinate system) and the SRO can be expanded like in finite-element codes, permitting scans of composition space with but a few compositions initially tested. Our calculations of Ti 0.25 CrFeNiAl x predict stable B2 (Al-poor region) and B2+L2 1 (Al-rich region), whereas experiments did not find any signature of L2 1 ordering. In contrast, the CALPHAD phase diagram suggests A2+B2+L2 1 phases, i.e., a A2→B2→L2 1 ordering transformation during slow cooling process [37]. In Fig. 4, we show that there are 8 possible distinct sites in L2 1 superstructure that can be populated without destroying cubic symmetry, i.e., 2π/3 rotation along 111 , and mirror symmetry along 110 . An HEA with more than 8-components, even with stoichiometric compositions, cannot be populated with small-cell cubic order.
Thus, we can infer from Fig. 4 that L2 1 order in quinary systems can exist. So, HEAs often exhibit a series of ordering transitions with one or more partially-ordered phases, and it ultimately reaches to fully-ordered states (if stoichiometric). In general, upon cooling, A2 goes to B2-type order by lowering its symmetry along 111 , e.g., in a binary A2 to B2 with corner (center) sites are perfectly ordered. In a HEA, as temperature is lowered, the B2 phase can only be partially-ordered, then B2 phase breaks symmetry along

Conclusion
For arbitrary HEAs, we presented first-principles KKR-CPA methods for calculating the relative stability and the short-range order (SRO) via thermodynamic linear-response theory (all pair correlations obtained simultaneously) using the electronic-structure of the homogeneous alloy as a reference state. The SRO gives detailed information on the unstable ordering modes inherent in the HEA, as well as their origins.
From a concentration-wave analysis of the SRO, a fast assessment was made for the probability distributions for all possible low-temperature competing partially-ordered states for a HEA in a given Bravais lattice, which directly provides the partially-or fully-ordered unit cells for all competing states. With the relevant competing states identified, we can make an estimate of energy gain for each structure from the SRO parameters (a 'back-ofthe-envelope' estimate that gives correct trends), and we used the partiallyordered unit cells directly in the KKR-CPA to obtain quantitative results for 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 energy difference between all competing states.