Configurational Entropy in Multicomponent Alloys: Matrix Formulation from Ab Initio Based Hamiltonian and Application to the FCC Cr-Fe-Mn-Ni System

Configuration entropy is believed to stabilize disordered solid solution phases in multicomponent systems at elevated temperatures over intermetallic compounds by lowering the Gibbs free energy. Traditionally, the increment of configuration entropy with temperature was computed by time-consuming thermodynamic integration methods. In this work, a new formalism based on a hybrid combination of the Cluster Expansion (CE) Hamiltonian and Monte Carlo simulations is developed to predict the configuration entropy as a function of temperature from multi-body cluster probability in a multi-component system with arbitrary average composition. The multi-body probabilities are worked out by explicit inversion and direct product of a matrix formulation within orthonomal sets of point functions in the clusters obtained from symmetry independent correlation functions. The matrix quantities are determined from semi canonical Monte Carlo simulations with Effective Cluster Interactions (ECIs) derived from Density Functional Theory (DFT) calculations. The formalism is applied to analyze the 4-body cluster probabilities for the quaternary system Cr-Fe-Mn-Ni as a function of temperature and alloy concentration. It is shown that, for two specific compositions (Cr25Fe25Mn25Ni25 and Cr18Fe27Mn27Ni28), the high value of probabilities for Cr-Fe-Fe-Fe and Mn-Mn-Ni-Ni are strongly correlated with the presence of the ordered phases L12-CrFe3 and L10-MnNi, respectively. These results are in an excellent agreement with predictions of these ground state structures by ab initio calculations. The general formalism is used to investigate the configuration entropy as a function of temperature and for 285 different alloy compositions. It is found that our matrix formulation of cluster probabilities provides an efficient tool to compute configuration entropy in multi-component alloys in a comparison with the result obtained by the thermodynamic integration method. At high temperatures, it is shown that many-body cluster correlations still play an important role in understanding the configuration entropy before reaching the solid solution limit of high-entroy alloys (HEAs).


Introduction
Multicomponent systems, called High Entropy Alloys (HEAs), are crystalline solids that form predominantly in a single phase. These systems were brought to wider attention through the work of Cantor et al. [1] and Yeh et al. [2]. At or near equiatomic ratio of compositions, the configuration entropy is maximized and given by the formula S con f ig /R = lnK where R is the ideal gas constant, and K is the total number of different components. It has been proposed to define K = 4 as the minimum number of components at or near equiatomic composition for these systems to be called HEAs [3][4][5].
There are various sources of entropy in a system in addition to configuration entropy including vibration, electronic and magnetic contributions. The co-existence of multiple phases in the equilibrium state of a chemical system with a given overall composition at a specified temperature is found by minimizing the Gibbs free energy of the whole system. When the system undergoes a phase transition of either order-disorder or spinodal decomposition, the relative differences of configurational entropy are dominant in comparison to those contributions from vibration and magnetic terms in FeCoCrNi [6] where the relative difference of configuration entropy dominates or is of similar magnitude to vibration entropy. These multicomponent alloys have transition metals Co, Cr, Fe, and Ni that exhibit magnetic behavior and have magnetic phase transitions characterized by Curie temperatures sensitive to concentrations [7].
In general, the formulation of configuration entropy given by the expression S con f ig /R = lnK only holds for disordered alloys with equiatomic composition at temperatures near their melting point. On the other hand, at low temperatures, the lowest free energy structures can be ordered intermetallic phases or partially ordered structures. In the intermediate temperature range, there exists short range order, which is related to the nature of the chemical environment of each atomic species, containing ordering or segregation preferences specially at low temperatures, because random disordering is favored at the high temperatures. The importance of short-range order (SRO) was manifested in diffuse X-ray scattering measurements of alloys containing ordered superstructure domains. The SRO effect is also seen to significantly affect electrical resistivity properties [8,9], and to influence the alloy strength by hindering dislocation motion [10] in concentrated alloys. Recently, it has been demonstrated how short-range order can be predicted from ab initio based Hamiltonian in combination with Monte Carlo simulations [11] in multi-component Mo-Nb-Ta-V-W systems. The study predicted that a strong SRO parameter may lead to the formation of Mo-Ta ordering in the B2 structure after quenching down from temperatures as high as 3000 K. HEAs such as Al 1.3 CoCrCuFeNi have been reported to be susceptible to complex phase transitions including segregation, precipitation, chemical ordering and spinodal decomposition into a complex microstructure containing regions of Body-Centered Cubic (BCC), Face-Centered Cubic (FCC) and B2 phases [12]. Physical models have been successfully applied to predict the formation of these phases on the basis of the average number of valence electrons [13].
The prediction of equilibrium thermodynamic properties (free energies, and phase diagrams) is one of the goals of computational materials science. Lattice statistical models involving an Ising-like Hamiltonian developed from ab initio enthalpies of mixing have become an important tool in the computation of the thermodynamic properties of alloy systems. In particular, the cluster expansion method [14][15][16] expands the Hamiltonian into contributions from an optimized set of clusters, each term weighted by Effective Cluster Interactions. The thermodynamic properties at temperature are obtained from the ECIs by computing the free energies from semi-canonical Monte Carlo simulations in combination with the thermodynamic integration technique [17].
The Cluster Variation Method (CVM) formalism [18,19] expresses the free energy in terms of enthalpies of mixing and configuration entropies as a function of the temperature dependence for a specific set of clusters by minimizing the free energy from variational principles. The correlation function parameters in terms of which configuration entropy and enthalpy of mixing are expressed grow exponentially with component number K and size of the maximal cluster ω. Due to this, the CVM has been commonly applied to clusters consisting of four sites in a regular tetrahedron or the so-called tetrahedron-octahedron in the FCC lattice [20,21], or the four sites in a irregular tetrahedron in the BCC lattice [22]. Within the tetrahedron approximation, the CVM calculations for the FCC involve the empty lattice, isolated points, first nearest neighbor pair, triangle and tetrahedron [20]. It is generally accepted that the integrated Monte Carlo method is more accurate, but the calculations are time consuming due to the need for many passes to obtain free energies at any given temperature from the disordered high temperature configuration. A hybrid approach taking the Monte Carlo calculated correlation functions at temperature and computing the configuration entropy value from analytic CVM expression has been used in a binary FCC lattice model [23]. It is shown that accurate free energies can be obtained for ordered and disordered phases at arbitrary chemical concentration and temperature without thermodynamic integration provided that use is made of high-order CVM entropy expressions. However, the clusters optimized from the Cluster Expansion (CE) method reported in the literature do not often match those used for the CVM. More importantly, calculations of the configuration entropy for multi-component systems represent a very serious challenge for computational materials science. In this work, we close this gap by developing a new methodology based on matrix formulation to calculate analytically the cluster probabilities for arbitrary K-component alloys from the correlation functions obtained by the hybrid Monte Carlo and CE Hamiltonian. We apply our method to the FCC four-component CrFeMnNi system for investigating configuration entropy as a function of temperature and alloy composition. Two specific compositions are chosen to illustrate the importance of many-body cluster probability functions in multi-component systems. The first one is equatomic composition Cr 25 Fe 25 Mn 25 Ni 25 , which is of a great interest for the HEA community. The second one is Cr 18 Fe 27 Mn 27 Ni 28 which has been used recently to design radiation tolerant materials for advanced nuclear reactor systems. This alloy has been studied in preference to Cantor's one [1] due to the removal of Co, which can cause activation by transmutation of 60 Co isotope in nuclear reactors [24]. This paper is organized as follows. In Section 2, a new matrix formulation for a K-component system based on the orthonormal sets of cluster expansion initially introduced within the Alloy Theoretic Automated Toolkit (ATAT) program [15,16] is presented. This formulation provides a rigorous link of cluster correlation functions obtained from Monte Carlo simulations with multi-body probabilities. Notation is developed to treat as general as possible arbitrary cluster sizes, ω, number of components, K, and temperature. In Section 3, the hybrid method is illustrated to calculate 4-body cluster probabilities as a function of temperature for the two specific compositions of equiatomic Cr 25 Fe 25 Mn 25 Ni 25 and also the composition Cr 18 Fe 27 Mn 27 Ni 18 from [24]. Here, the connection with the tetrahedron approximation in CVM is discussed due to the presence of a first nearest-neighbor 4-body cluster interaction within the investigated CE Hamiltonian. In Section 4, configuration entropy is discussed at two temperatures 1000 and 3000 K for different cluster decorations obtained from the CE method and compared with the integrated Monte Carlo results. The changes in configuration entropy are attributed to the presence of ordered phases that are more stable at low temperatures and the complementary tendency towards disordered random solution of the alloy at the given average composition. The main conclusions of this work are given in Section 5.

Matrix Formulation of Cluster Expansion
Let us consider an alloy system with arbitrary number, K, of components and crystalline lattice symmetry, in which the disordered phase is described by space group, G. The CE Hamiltonian at T = 0 K of the alloy can be found from enthalpies of mixing, ∆H Mixing DFT [ σ], of derivative alloy structures [16,25] denoted by varying length arrays of spin-like variables σ i taking values from 0 to K − 1. In general, a derivative structure of the alloy has non-zero atomic concentrations for each of p = 1, ..., K different chemical elements forming it i.e., x[ σ] p = 0. For each member structure, the reference energy in the calculation of the enthalpy of mixing is the pure element total energy, E total [p], with the same crystallographic symmetry. The alloy enthalpy of mixing is defined as follows: The enthalpies of mixing are calculated from the total energies E total [ σ] and E total [p] of the alloy and the pure element systems, respectively, by density functional theory for the underlying lattice. In the CE, the ∆H mixing DFT [ σ] of an arbitrary structure σ is expanded into a sum of reference clusters, ∆H

Mixing CE
[ σ] and can be written in the following formula [15][16][17]: (s) . ( In Equation (2), each reference cluster is characterized by three labels ω, n and (s). The label ω denotes the total number sites; n is an auxiliary label that refers to highest order coordination shell contained in the reference cluster; and finally the label (s) = (j 1 , j 2 , · · · , j ω ) denotes decoration of the cluster by point functions with dimension equal to ω and j i taking values 0, ..., K − 1. For each reference cluster, there is an associated effective cluster interaction J  (2) is the thermally averaged cluster correlation function over all clusters ω , n , (s ) which are equivalent to a space group symmetry element of the disordered phase within the reference decorated cluster defined by ω, n, (s). Overall, the triple product of multiplicity, ECI and correlation function, m gives the energetic contribution per lattice site of the reference cluster ω, n and (s) to the enthalpy of mixing of the particular structure configuration given by σ.
For multicomponent systems, the choice of a set of basis functions is important for the matrix formulation of the CE method [14]. The simplest set consists of the successive K − 1 powers of the pseudo-spin configuration variable {σ} as was originally suggested by Taggart [26]. In this work, the correlation function is defined as the product of point functions initially proposed in the ATAT program [15,16]. Configurational average of the correlation functions is then given by the following formula: where the point function is written as: In Equation (3), the correlation function of the system is averaged over the arbitrary crystal structure of the alloy system with configuration σ over all the set of Ω[ω, n] clusters, is defined as the number of times the reference cluster (ω, n) is contained in a structural configuration which can be obtained from Monte Carlo simulation. The sites in the cluster { ←−−→ (ω, n) u } can allocate all of the possible decoration values ∈ (s) which we express as (s) = (j 1 , j 2 , · · · , j ω ). These integer indexes corresponding to the decorated cluster are parameters for the point functions in Equation (4).
For an arbitrary ω-sites cluster and K components, the number of ω-tuples formed with integer entries running from 0 · · · K − 1 can be calculated. For ω = 2, the formula calculates the number of symmetrically unique decorations (s) for a two-point cluster Γ ω=2,n, (s) and K component alloy system reduces to (K + 1)K/2 [11], where K is the number of components. For higher order clusters, the total number of decorations depend on the cluster coordinates and the space group symmetry G of the high temperature disordered phase i.e., FCC or BCC and can not be simply expressed in terms of K (in general there could be less than ( ω+K−1 number of symmetrically unique correlation functions). ATAT [15,16] numerically works out all the number of symmetrically unique decorated clusters ω , n , (s ) equivalent to ω, n, (s) and uses one correlation function per set of equivalent decorated clusters Γ . For convenience in notation, from now on, we ω,n, (s) . In general, two decorations (s) and (s ) are symmetrically equivalent if there is at least one element in the space group symmetry g ∈ G that transforms the }. This is given by the following equation: With the symmetry relations between decorations (s) and (s ) in any given cluster, we are able to retrieve all K ω possible ω-tuples from the set of unique symmetries corresponding to an ω cluster. The point functions γ j,K [σ i ] used to define the general correlation functions are related to the multi-body cluster probabilities (see below in Equation (9)) by direct products of a linear transformation [27], τ K where the new matrix τ is constructed from these point functions,γ j,K [σ i ]. It can be trivially shown that the inverse of Van der Monde matrices with complex entries equal to roots of unity is its complex conjugate. The inverse of the τ K matrix can be obtained from the complex conjugate matrix of τ K by taking the real and imaginary part and riffling their rows. From Equation (4), the following expression results for a system with K components: To the best of our knowledge, Equation (7) represents a new formulation for the inverse of τ K matrix to ensure that the basis set defined by Equation (4) is rigorously orthonormal. The size of τ K matrix is KxK, where K is the number of components. It is convenient to perform the matrix multiplications with all K ω decorations formed from ω-tuples with integer elements running from 0 to K − 1. In particular, for the case of four component K = 4, these matrices τ become by applying Equations (4) and (7) for the inverse τ −1 K : We duplicate the symmetrically unique decorations (s) whenever two decorations are connected by symmetry of the disordered structure. For relating two equivalent decorations, we find it convenient to use a permutation representation of the space group operator as permutations of ω-site tuples. We use the property of invariance of the cluster expansion to obtain probability distributions from correlation functions [28]. As a consequence of the compact formalism and by using Equation (4) and (7), the expression of correlation functions can be rewritten into a matrix form: With the aid of the matrix formulation from Equation (9) and the generalized form of the inverse matrix τ −1 K in Equation (7), one can express the ω-cluster probabilities into a matrix form: where we have used the notation for direct product of matrices In addition, we also implied summation over repeated indexes on the right-hand side of Equation (10). Note that the size of these matrices increases exponentially with cluster size ω, K ω xK ω ; in particular, for ω = 4, the matrices have 256 × 256 entries. The cluster probabilities are normalized as expressed in Equation (11) K As a consequence of the normalization of the cluster probabilities, it is possible to separate probabilities of the decorated sub-clusters from the probabilities of the maximal cluster by partial summations over all possible decorations for the sites that belong to the maximal cluster (ω, n) but not the i-th sub-cluster (ω i , n i ): From Equation (10) and for the case with ω = 2, it follows that the generalized expression for SRO of species p and q at the nth shell, α , can be interpreted as the tendency to order or segregate species p and q with respect the disordered random probability given by the product of their elemental, p, bulk concentration for Mn-Ni. The SRO allows a quantitative description of the interactions between atoms as a function of temperature to predict order-disorder transition temperatures [11,29]. In particular, the matrix formalism from Equations (7) and (10) allows one to generalize the SRO treatment for an arbitrary number of components, K:

Configuration Entropy in the Matrix Formulation
In general, a thermodynamical system in state σ and with enthalpy of mixing given by the CE Hamiltonian ∆H

Mixing CE
[ σ] is described by a set of symmetry unique probability distributions characterized by decorations of a chosen maximal (ω, n) cluster. In practice, the chosen maximal cluster (ω, n) contains few points. As a mean field approximation [27], the CVM can be rationalized as a factorization of the probability distributions of the (ω, n) maximal cluster into integer powers η ω 1 ,n 1 , η ω 2 ,n 2 · · · , η ω s[ω,n] ,n s[ω,n] of the probability distributions of the sub-clusters (ω i , n i ) ⊆ (ω, n); i = 1, · · · , [ω, n] [19] with decorations (pq · · · ) ω i ,n i corresponding to the components (pq · · · ) ω i ,n i of all decorations (pq · · · ) ∈ {(pq · · · ) (ω,n) } occupied by sites of the (ω i , n i ) sub-cluster The following expression, a consequence of the CVM factorization scheme, can be derived from [19] and is the reason why the disordered configuration at high temperature is reproduced from This occurs because, at the high temperature limit, the CVM probabilities tend to products of composition of the species involved in the decorations of the cluster: where y is the concentration of the (pq · · · ) j equal to one of the integers 0, 1, · · · , K − 1 in the disordered state of the alloy. A natural consequence of the factorization scheme chosen for the CVM multi-body probabilities is that the configuration entropy assumes the formulation where ∼ S ω i ,n i , is the entropy contribution to cluster (ω, n) from the sub-cluster (ω i , n i ): The set of integers η ω 1 ,n 1 , η ω 2 ,n 2 · · · , η ω s[ω,n] ,n s[ω,n] are the mean field integer coefficients associated with the partition function of the alloy system. In the theory of regular mixtures, the coefficients can be found from the recursive heuristic expression after Kikuchi, [30], Barker [31] and whose formulation was explicitly derived by using group theoretic methods by Gratias et al. [32]. It requires the determination of two quantities: the site multiplicity N ω i ,n i and the sub-cluster multiplicity N β ω i ,n i . The site multiplicity can be determined by calculating the number of symmetry operators, N ω i ,n i , that stabilize the ω i , n i cluster i.e., g ∈ G such that the application of g into the set of cluster positions , where |G| is the order of the point group associated with the space group G and |N ω i ,n i | is the order of the group N ω i ,n i . The sub-cluster multiplicity N β ω i ,n i is just the frequency of the cluster β that is contained in the cluster (ω i , n i ) In particular, the formulation applied to a point cluster retrieves the Bragg-Williams approximation for the maximal cluster (ω = 1, n = 1), i.e., a site cluster, giving the entropy weight of η ω 1 =1 1 ,n 1 =1 = −1 and for a 2-body cluster (ω = 2, n) in the nth coordination shell, we get η ω 2 ,n 2 = −N 2,n and η ω 1 ,n 1 = 2N 2,n − 1, where N 2,n is site multiplicity of the cluster (2, n) calculated from the sites in this cluster and the space group G.
In this work, the above matrix formulation is applied in the hybrid CE-Monte Carlo method which performs the free energy minimization from the CE Hamiltonian in a combination with Monte Carlo simulations. Within the process, the Monte Carlo method produces the correlation functions for the equilibrium configurations found at each of the temperatures investigated. The hybrid approach uses these correlation functions in the analytic expressions for configuration entropy.
There is an alternative to the hybrid approach for the entropy calculation where the thermodynamic integration method can be used. Here, entropy is calculated from the configuration contribution to the specific heat at constant volume, C con f derived from the fluctuations of enthalpy of mixing at temperature values in a fine grid of temperature values: where ∆H

and ∆H
Mixing CE 2 are the square of the mean and mean square enthalpies of mixing, respectively, calculated by averaging over all the MC steps at the accumulation stage for a given temperature. The accuracy of evaluation of configuration entropy depends on the size of temperature integration step and the number of MC steps performed at the accumulation stage [33]. The integration of specific heat is performed from 0 K to the temperature T. In order to calculate the configuration entropy at a given temperature, the value for specific heats at lower temperatures is thus required. For example, assuming that the chosen temperature step is equal to 5 K, to evaluate configuration entropy numerically at 3000 K would require computing the specific heat at 600 smaller temperatures. In contrast, using Equation (17), the configuration entropy can be computed analytically from the correlation functions at any given temperature and alloy composition. Our experiences show that the computational time using the hybrid method can be of two orders faster than those by the thermodynamic integration.
In the semi canonical Monte Carlo calculations performed, the temperature range and temperature steps are important. The accuracy of the thermodynamic integration method to calculate configuration entropy generally requires a smaller temperature integration step, ∆T = 5 K and also depends on the number of Monte Carlo passes [33]. In particular, we use a cell containing 2048 atoms distributed into a 8 × 8 × 8 primitive unit cell, and average compositions for the ensemble given by Cr 18 [42]. By quenching down systematically with the temperature step of ∆ T = 5 K, various equilibrium configurations were obtained at lower temperatures. For thermodynamic integration calculation of configuration entropy, we integrated numerically the specific heat at constant volume using the theoretical formula (19) starting from the lowest temperature value 0 K to 3000 K.

Cluster Probability Functions in FCC Cr-Fe-Mn-Ni Alloys
We apply the matrix formulation of cluster expansion outlined in Section 2.1 to the FCC CrFeMnNi system for investigating the temperature and composition dependent cluster probability distribution functions and the configuration entropy.

Cluster Expansion Hamiltonian for FCC CrFeMnNi
The DFT enthalpies of mixing for the FCC CrFeMnNi system were used to map iteratively into the cluster expansion Hamiltonian given by Equation (2) using the ATAT package [16]. The mapping has been performed systematically from the six binary and four ternary constituent subsystems of the considered quaternary. The database of structures for the cluster expansion consisted of 835 structures categorized by the difference of local environments in binaries (structures with two chemical elements: 58 CrFe, 55 CrMn, 77 CrNi, 58 FeMn, 54 FeNi and 52 MnNi), ternaries (structures with three chemical elements: 89 CrFeMn, 85 CrFeNi, 46 FeMnNi, and 66 CrMnNi); and 191 quaternaries CrFeMnNi. More information about the type of binary and ternary structures used in our DFT database has been detailed in our previous work [17]. Structures are typically ordered ones with their composition ranging for each constituent element from 5% to 95%. It is important to stress here that, different to other studies of HEAs, our DFT database for constructing the CE Hamiltonian didn't include randomly distributed structures such as the so-called special quasi-random structures (SQSs). The latter, however, can be generated within the present approach from Monte-Carlo simulations after obtaining the reliable ECIs. The set of clusters which have minimized the cross-validation score of 12.95 meV/atom, consists of six 2-body, two 3-body and one 4-body ECIs. In the present work, the clusters with the same sizes and relative positions have been included consistently in the considered subsystems. The prediction of the corresponding ground-state intermetallic phases is in a good agreement with the experimental binary phase diagrams available (for example, the ferromagnetic FeNi 3 in L1 2 structure and anti-ferromagnetic MnNi in L1 0 structure) as well as with the previous theoretical study for the ternary and ferrimagnetic CrFe 2 Ni phase in NiCu 2 Zn structure [17]. A new intermetallic phase FeCr 2 MnNi 4 is predicted for the quaternary system and full results of magnetic properties in CrFeMnNi systems will be discussed in a separate work.
The CE Hamiltonian for the quaternary system CrFeMnNi consists, in total, of 83 different decorated clusters distributed among 10 different non-decorated clusters: four decorated clusters for the point cluster, (ω = 1, n = 1); six decorated clusters for each pair cluster (ω = 2, n = 1, · · · , 6); 10 decorated clusters for the non-decorated cluster (ω = 3, n = 1); 18 decorated clusters for the non-decorated cluster (ω = 3, n = 2); and 15 decorated clusters for the non-decorated cluster (ω = 4, n = 1). The decoration labels, (s) required to specify the clusters, are listed in Table 1. Each non-decorated cluster is defined by the coordinates (specified with respect to the standard Cartesian coordinate system in units of lattice spacing) of the lattice sites that it includes. The decorations define the chemical species allocated to each site in strictly the same order i.e., for (ω = 2, n = 1) cluster the decoration (2, 3) means that species 2 is allocated for site with coordinates (1, 1, 1) and species 3 is allocated for site with coordinates (1/2, 3/2, 3/2).

Full Set of Cluster Decorations
In general, the set of temperature dependent decorated cluster correlation functions, Γ (s) ω,n [ σ] , constitute a set of K ω quantities for a given maximal cluster (ω, n) in the CE. From each of the K ω temperature dependent decorated cluster correlation functions, the matrix formalism described in Section 2.1 generates the set of temperature dependent multi-body probability functions describing the temperature dependent behavior associated with the maximal cluster (ω, n). The symmetry unique decorations for all of the clusters in the CE, ∀(ω, n), are reported in the ATAT clusters output file and here they are listed in Table 1. If one cluster is included into another, it becomes a sub-cluster (ω i , n i ) ∈ (ω, n); these sub-clusters have been classified according to their inclusion into maximal clusters, and are listed in Table 2. If the sub-cluster (ω i , n i ) is included in the maximal cluster, (ω i , n i ) ∈ (ω, n), then the decorations from the sub-cluster (s) (ω i ,n i ) can be transferred into the decorations of the maximal cluster, (s) (ω,n) , by using the intrinsic space group symmetry of the parent lattice. In general, any given sub-cluster (ω i , n i ) can be found several times within the cluster (ω, n). For convenience, symmetry of the cluster decorations is implemented by means of permutation operators in Table 2 including X for the empty site. The full set of decorated cluster correlation functions, Γ (s) ω,n [ σ] , corresponding to the maximal cluster (ω, n), is generated by studying which of the decorations (s) ω i ,n i form the sub-clusters (ω i , n i ) in (ω, n). The decorations can have empty cluster, i.e., at least one integer entry in 0 ⊆ (s) ω i ,n i .
As an important case of study, we chose the maximal cluster given by (ω = 4, n = 1) to obtain the 4-body configuration probabilities. Figure 1a,b show the plots of 35 symmetrically unique probabilities for the 4-body maximal cluster defined by the set of lattice sites ((1, 1, 1), (3/2, 3/2, 1), (3/2, 1, 1/2), (1, 3/2, 1/2)) as a function of temperature.  Figure 1a,b. It should be noted that, for both compositions, the cluster configuration referred to as Cr-Cr-Cr-Cr appears as the least probable of all the cluster configurations. This finding is consistent with the fact that Cr has the BCC ground-state and therefore the probability of finding a Cr cluster in the FCC lattice is negligible, in particular in the low temperature region.
For the equiatomic composition, the probability of Cr-Fe-Mn-Ni is particularly very high at temperatures between 0-900 K. The presence of 4-body Cr-Fe-Mn-Ni clusters at low temperatures demonstrates the relationship with our DFT/CE prediction of the new ordered phase FeCr 2 MnNi 4 in the quaternary system. Furthermore, in the temperature range below 900 K, the cluster configuration Cr-Fe-Fe-Ni shows that it is the second most probable configuration. This configuration is directly correlated with the ordered Fe 2 CrNi structure predicted in our earlier study of phase stability in the ternary Fe-Cr-Ni system [17]. In the temperature region between 900 K and 1200 K, an increase of probability of Mn-Mn-Ni-Ni cluster is significantly important. Beyond 1200 K, all of the cluster configuration probabilities tend to the solid solution or random configuration.
The total configuration entropy, S ω,n [ σ], is calculated by adding the η ω i ,n i weighted quantities By applying the same arguments developed for the cluster (ω = 4, n = 1), the total configuration entropy for each of the 10 different maximal clusters can be calculated. The sub-clusters and the permutation operators contained in each of the 10 maximal clusters appearing in the cluster expansion are detailed in Table 2.
Composition dependent entropies at fixed temperatures 1000 K and 3000 K are shown in Figure 3a,b. At each of these temperature values, any of the configuration entropy expressions from Table 2 provides approximately the same value. At 3000 K, the entropy is maximized in the center of the tetrahedron, namely at the equiatomic composition to 1.37 k B , and it is decreased upon lowering the temperature at all of the composition points. The variation of configuration entropy is too complex to be captured in the fixed temperature plots, therefore we use specifically the equiatomic composition and calculate the configuration entropy as a function of temperature in Figure 4. Regarding the high temperature limit value of entropy in Figure 4, we use as maximal cluster that from the 4-body, 3-body and 2-body in the nearest neighbor cluster probabilities. Above 1200 K, the cluster approximation for the configuration entropy approaches the high temperature limit of disordered solute solution (−k B ∑ 4 p x p [ σ]ln(x p [ σ])) much more rapidly than the thermodynamic integration. The reason the high temperature limit is preserved can be clearly seen from the factorization of the CE expression that resulted from analytical derivation in Equation (15). It is important to stress that, within the high-temperature limit, our theoretical prediction recovers the ideal configuration entropy of mixing for solution phase as it has been originally proposed by the phenomenological guidelines in designing high-entropy alloys using thermodynamic and topological parameters of the constituent elements [44].  It is worth mentioning that the entropy increase ( Figure 4) as a function of temperature for the equiatomic alloy composition can thus be understood by the corresponding decrease in probabilities (Figure 1) of the ordered state of the structures L1 0 towards the disorder configuration A1 in Figure 2. Similarly, in the alloy with composition Cr 18 Fe 27 Mn 27 Ni 28 , the increase of configuration entropy is correlated with the decrease in high probability of the configurations of Cr-Fe-Fe-Fe and Mn-Mn-Ni-Ni clusters originally related to the ordered FCC-like L1 2 -CrFe 3 and L1 0 -MnNi structure, respectively.
For the cluster (ω = 4, n = 1), the calculated configuration entropy in Figure 4 appears to be negative at low temperature with respect to the positive values correctly predicted by the thermodynamic integration method. This demonstrates clearly that this maximum cluster within the tetrahedron approximation conventionally adopted within the CVM is only valid to describe the configuration entropy for the four-component Cr-Fe-Mn-Ni system in the high-temperature limit. It is shown from Figure 4 that the first nearest-neighbor triangle cluster (ω = 3, n = 1), which plays an important role within the tetrahedron approximation, also gives physically incorrect and negative entropy contribution at low temperature. As it has been discussed in the introduction, in difference from the CVM tetrahedron approximation, besides the first nearest-neighbor cluster contributions, the present CE results also include other pair clusters up to the sixth nearest neighbors and the second nearest-neighbor triangle cluster contributions. For example, in Figure 4, the third nearest-neighbor pair cluster (ω = 2, n = 3) gives the significantly positive contribution to configuration entropy in the entire range of temperature. These additional cluster contributions, in turn, ensure the correct behavior of the configuration entropy obtained from the thermodynamic integration. Therefore, from a statistical physics point of view, the hybrid technique combining Monte Carlo simulations with the CE method can be considered to be more advanced than the mean-field approach advocated within the CVM.

Conclusions
In this work, we develop a matrix formalism to study multi-body ordering probabilities beyond pair approximation previously used for investigating the SRO and configuration entropy in multi-component alloys by using a hybrid combination of CE and Monte Carlo methods. The cluster probabilities are worked out by explicit inversion within the orthonormal sets of the point functions adopted in the ATAT package and a direct product of a matrix formulation obtained from symmetrically independent correlation functions. The correlation functions are determined from semi-canonical Monte Carlo simulations and ECIs derived from DFT calculations. We apply our method to the quaternary FCC Cr-Fe-Mn-Ni system by considering 285 different alloy compositions covering the compositional space of the quaternary alloy as a function of temperature. To further assess our formulated expressions for configuration entropy, focus is put on two alloy compositions, equiatomic Cr 25 Fe 25 Mn 25 Ni 25 and Cr 18 Fe 27 Mn 27 Ni 28 , to obtain cluster probabilities and understand the variation in configuration entropy with temperature due to the lowering of the ordering probability corresponding to the ordered configuration. The cluster probability plots against temperature show that, for the composition Cr 18 Fe 27 Mn 27 Ni 28 , there is a high probability of formation of the L1 0 MnNi phase at low temperatures below 1300 K and for the L1 2 CrFe 3 phase in the temperature range 500-1200 K. Similarly, for the equiatomic composition, the L1 0 MnNi phase appears stable in the temperature range 900-1200 K, while the lower temperature region is preferred for the configuration Cr-Fe-Mn-Ni. The configuration Cr-Cr-Cr-Cr was found to be the least probable configuration at all temperature ranges for both compositions of the Cr-Fe-Mn-Ni system. Furthermore, the configuration entropy as a function of temperature was derived from these probabilities: the high-temperature limit is in accordance with random solid solution approximation, but, at low temperatures, the entropy is seen to be reduced due to ordering or segregation tendencies which are in turn determined by multi-body probability functions including chemical short-range order.
We believe that the present study will help to promote further understanding of derivative phases as a function of temperature and multi-component alloy composition from disordered solid solutions in the phenomenological description of HEAs. By applying the formalism to the Cr-Fe-Mn-Ni system, it will also serve as a benchmarking example in designing radiation tolerant materials for advanced nuclear reactor systems by using the ab initio based CE method. The composition Cr 18 Fe 27 Mn 27 Ni 28 multicomponent system has been studied for advanced nuclear applications due to promising irradiation resistance with regards to void swelling.