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 (Cr 25Fe 25Mn 25Ni 25 and Cr 18Fe 27Mn 27Ni 28), 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 L1 2 -CrFe 3 and L1 0-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).

1. 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 where R is the ideal gas constant, and K is the total number of different components. It has been proposed to define 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 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 AlCoCrCuFeNi 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 , which is of a great interest for the HEA community. The second one is 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 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 CrFeMnNi and also the composition CrFeMn 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.

2. Methods

2.1. 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 K of the alloy can be found from enthalpies of mixing, , of derivative alloy structures [16,25] denoted by varying length arrays of spin-like variables taking values from 0 to . In general, a derivative structure of the alloy has non-zero atomic concentrations for each of different chemical elements forming it i.e., . For each member structure, the reference energy in the calculation of the enthalpy of mixing is the pure element total energy, , with the same crystallographic symmetry. The alloy enthalpy of mixing is defined as follows:

The enthalpies of mixing are calculated from the total energies and of the alloy and the pure element systems, respectively, by density functional theory for the underlying lattice. In the CE, the of an arbitrary structure is expanded into a sum of reference clusters, and can be written in the following formula [15,16,17]:

In Equation (2), each reference cluster is characterized by three labels and . 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 = () denotes decoration of the cluster by point functions with dimension equal to and taking values . For each reference cluster, there is an associated effective cluster interaction , and multiplicity per lattice site, . The term used in the definition of enthalpy of mixing in Equation (2) is the thermally averaged cluster correlation function over all clusters which are equivalent to a space group symmetry element of the disordered phase within the reference decorated cluster defined by . Overall, the triple product of multiplicity, ECI and correlation function, gives the energetic contribution per lattice site of the reference cluster and 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 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 clusters, where each cluster labeled by u contains sites, each site denoted by the vector, . The i-union of the vectors, , composes the cluster of vectors, . The are numeric labels referring to the values taken by variable u in Equation (3) that serve to enumerate the vectors making up the u cluster of vectors . The clusters are symmetrically equivalent to the reference cluster sites by a symmetry operation and are enumerated by an index taking values . is defined as the number of times the reference cluster is contained in a structural configuration which can be obtained from Monte Carlo simulation. The sites in the cluster can allocate all of the possible decoration values which we express as . 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 can be calculated. For , the formula calculates the number of symmetrically unique decorations (s) for a two-point cluster and K component alloy system reduces to [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 number of symmetrically unique correlation functions). ATAT [15,16] numerically works out all the number of symmetrically unique decorated clusters equivalent to and uses one correlation function per set of equivalent decorated clusters . For convenience in notation, from now on, we will use . In general, two decorations and are symmetrically equivalent if there is at least one element in the space group symmetry that transforms the sites into as the permutation connecting the decorations and . This is given by the following equation:

With the symmetry relations between decorations and in any given cluster, we are able to retrieve all possible -tuples from the set of unique symmetries corresponding to an cluster. The point functions 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], where the new matrix is constructed from these point functions,. 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 matrix can be obtained from the complex conjugate matrix of 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 matrix to ensure that the basis set defined by Equation (4) is rigorously orthonormal. The size of matrix is KxK, where K is the number of components. It is convenient to perform the matrix multiplications with all decorations formed from -tuples with integer elements running from 0 to . In particular, for the case of four component , these matrices become by applying Equations (4) and (7) for the inverse :

We duplicate the symmetrically unique decorations () 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 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 , x; in particular, for , the matrices have 256 × 256 entries. The cluster probabilities are normalized as expressed in Equation (11)

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 () but not the -th sub-cluster ():

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 the four component Cr-Fe-Mn-Ni system, there are six chemically distinct SRO parameters: for Cr-Fe; for Cr-Mn; for Cr-Ni; for Fe-Mn; for Fe-Ni; and 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:

2.2. Configuration Entropy in the Matrix Formulation

In general, a thermodynamical system in state and with enthalpy of mixing given by the CE Hamiltonian is described by a set of symmetry unique probability distributions characterized by decorations of a chosen maximal () cluster. In practice, the chosen maximal cluster () contains few points. As a mean field approximation [27], the CVM can be rationalized as a factorization of the probability distributions of the () maximal cluster into integer powers of the probability distributions of the sub-clusters [19] with decorations corresponding to the components of all decorations occupied by sites of the () 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 multi-body probabilities, . 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 is the concentration of the equal to one of the integers 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 , is the entropy contribution to cluster from the sub-cluster :

The set of integers 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 and the sub-cluster multiplicity . The site multiplicity can be determined by calculating the number of symmetry operators, , that stabilize the cluster i.e., such that the application of g into the set of cluster positions results in a permutation of the set. Then, where is the order of the point group associated with the space group G and is the order of the group . The sub-cluster multiplicity is just the frequency of the cluster that is contained in the cluster ()

In particular, the formulation applied to a point cluster retrieves the Bragg–Williams approximation for the maximal cluster , i.e., a site cluster, giving the entropy weight of and for a 2-body cluster in the nth coordination shell, we get and , where is site multiplicity of the cluster 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, derived from the fluctuations of enthalpy of mixing at temperature values in a fine grid of temperature values: where and 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.

2.3. Computational Details

The DFT enthalpies of mixing given by Equation (1) were calculated from fully relaxed spin-polarized DFT total energy calculations performed using the projector augmented wave (PAW) method [34] implemented in Vienna Ab initio Simulation Package (VASP) [35,36,37,38,39]. Exchange and correlation were treated in the generalized gradient approximation (GGA) and the Perdew-Burke-Ernzerhof (PBE) functional [40]. The core configurations of Fe, Cr, Mn and Ni in PAW potentials are [Ar]3d4s, [Ar]3d4s, [Ar]3d4s and [Ar]3d4s, respectively. The total energies were calculated using the Monkhorst–Pack [41] mesh of 12 × 12 × 12 k-point for a four-atom FCC cubic cell. The plane wave cut-off energy used in the calculations was 400 eV. The total energy convergence criterion was set to 10 eV/cell, and force components were relaxed to 10 eV/nm.

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, 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 CrFeMnNi and equiatomic CrFeMnNi. The Monte Carlo simulations were performed from random configurations at high temperature (3000 K) where the configuration entropy at the high-temperature limit is given by for equiatomic and for CrFeMnNi [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.

3. 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.

3.1. 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 in structure and anti-ferromagnetic in structure) as well as with the previous theoretical study for the ternary and ferrimagnetic phase in structure [17]. A new intermetallic phase 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, (); six decorated clusters for each pair cluster (); 10 decorated clusters for the non-decorated cluster (); 18 decorated clusters for the non-decorated cluster (); and 15 decorated clusters for the non-decorated cluster (). The decoration labels, 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 () cluster the decoration means that species 2 is allocated for site with coordinates and species 3 is allocated for site with coordinates .

Table 1

, highest coordination shell n, decoration and coordinates of points in the relevant clusters on the Face-Centered Cubic (FCC) lattice. The coordinates are referred to the simple cubic Bravais lattice. Index is the same as the sequence of points in the relevant cluster. The canonical order for decoration indexes, , is 0, 1, 2 and 3 is Cr, Fe, Mn and Ni. All values of the Effective Cluster Interactions (ECIs) obtained from the present CE study are shown in the last column.

ω	n	(s)	Coordinates	ECI (meV/atom)	ω	n	(s)	Coordinates	ECI (meV/atom)
1	1	(0)	(1,1,1)	+0.00	3	1	(1,1,1)	(1,1,1)	+0.20
		(1)		+0.11			(2,1,1)	(3/2,1,1/2)	+0.90
		(2)		−0.04			(3,1,1)	(1,3/2,1/2)	+1.60
		(3)		−0.01			(2,2,1)		−3.40
2	1	(1,1)	(1,1,1)	+9.40			(3,2,1)		−0.50
		(2,1)	(1,3/2,3/2)	−0.10			(3,3,1)		+1.20
		(3,1)		+3.40			(2,2,2)		+0.20
		(2,2)		+0.40			(3,2,2)		+2.00
		(3,2)		+1.30			(3,3,2)		−0.50
		(3,3)		+6.00			(3,3,3)		+0.00
2	2	(1,1)	(1,1,1)	−9.20	3	2	(1,1,1)	(1,1,1)	−0.60
		(2,1)	(1,1,0)	+0.40			(2,1,1)	(1,3/2,1/2)	+1.00
		(3,1)		−4.40			(3,1,1)	(1,1,0)	+0.90
		(2,2)		−11.60			(1,2,1)		−1.80
		(3,2)		−3.50			(2,2,1)		+2.30
		(3,3)		−8.80			(3,2,1)		−0.70
2	3	(1,1)	(1,1,1)	0.90			(1,3,1)		−2.10
		(2,1)	(2,3/2,3/2)	2.60			(2,3,1)		−0.30
		(3,1)		3.50			(3,3,1)		−0.60
		(2,2)		1.60			(2,1,2)		−6.30
		(3,2)		−0.30			(3,1,2)		−1.20
		(3,3)		0.10			(2,2,2)		-0.60
2	4	(1,1)	(1,1,1)	−0.40			(3,2,2)		−0.10
		(2,1)	(2,1,2)	2.40			(2,3,2)		+1.60
		(3,1)		1.20			(3,3,2)		+0.60
		(2,2)		0.50			(3,1,3)		−2.00
		(3,2)		0.60			(3,2,3)		+0.30
		(3,3)		−0.80			(3,3,3)		−1.40
2	5	(1,1)	(1,1,1)	−1.00	4	1	(1,1,1,1)	(1,1,1)	−3.30
		(2,1)	(1,3/2,−1/2)	−3.50			(2,1,1,1)	(3/2,3/2,1)	+2.00
		(3,1)		−2.00			(3,1,1,1)	(3/2,1,1/2)	+0.70
		(2,2)		0.20			(2,2,1,1)	(1,3/2,1/2)	−2.90
		(3,2)		0.90			(3,2,1,1)		+0.60
		(3,3)		0.40			(3,3,1,1)		+0.70
2	6	(1,1)	(1,1,1)	0.80			(2,2,2,1)		−0.60
		(2,1)	(2,2,0)	2.10			(3,2,2,1)		−1.10
		(3,1)		1.00			(3,3,2,1)		+1.30
		(2,2)		−2.70			(3,3,3,1)		+2.60
		(3,2)		−1.30			(2,2,2,2)		−0.50
		(3,3)		0.50			(3,2,2,2)		+4.90
							(3,3,2,2)		+1.00
							(3,3,3,2)		−1.60
							(3,3,3,3)		−1.70

3.2. Full Set of Cluster Decorations

In general, the set of temperature dependent decorated cluster correlation functions, , constitute a set of quantities for a given maximal cluster in the CE. From each of the 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 . The symmetry unique decorations for all of the clusters in the CE, , 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 ; these sub-clusters have been classified according to their inclusion into maximal clusters, and are listed in Table 2. If the sub-cluster is included in the maximal cluster, , then the decorations from the sub-cluster can be transferred into the decorations of the maximal cluster, , by using the intrinsic space group symmetry of the parent lattice. In general, any given sub-cluster can be found several times within the cluster . 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, , corresponding to the maximal cluster , is generated by studying which of the decorations form the sub-clusters in . The decorations can have empty cluster, i.e., at least one integer entry in .

Table 2

List of permutation operators for generating the full set of decorations represented by dimensional integer arrays with entries taking values from 0 to . The symmetry operators (see Equation (5)) represented here in permutation form in column 2 act on the set of 83 symmetry unique decorations indicated under the column in Table 1 by permuting the entries in or by introducing the empty cluster X for sub-clusters belong to a given cluster; examples with discussion are provided in Section 3.2. The space group of the disordered FCC structure, , is implicitly assumed in order to convolute the symmetry unique into the full set of decorations. The last four columns represent ith sub-cluster of the maximal cluster ; site multiplicity of (); sub-cluster multiplicity of the cluster ; and , the sub-cluster () contribution to configuration entropy expression corresponding to the maximal cluster ()

Maximal Cluster	Permutation Operators {g_1,g_2,⋯}	Sub-Cluster {ω_i,n_i}	N_ω_i,n_i	Nω,nω_i,n_i	η_ω_i,n_i
ω=1,n=1	{(1)}	ω_1=1,n_1=1	1	1	−1
ω=2,n=1	{(1,X),(X,1)}	ω_1=1,n_1=1	1	2	11
	{(1,2),(2,1)}	ω_2=2,n_2=1	6	1	−6
ω=2,n=2	{(1,X),(X,1)}	ω_1=1,n_1=1	1	2	5
	{(1,2),(2,1)}	ω_2=2,n_2=2	3	1	−3
ω=2,n=3	{(1,X),(X,1)}	ω_1=1,n_1=1	1	2	23
	{(1,2),(2,1)}	ω_2=2,n_2=3	12	1	−12
ω=2,n=4	{(1,X),(X,1)}	ω_1=1,n_1=1	1	2	11
	{(1,2),(2,1)}	ω_2=2,n_2=4	6	1	−6
ω=2,n=5	{(1,X),(X,1)}	ω_1=1,n_1=1	1	2	23
	{(1,2),(2,1)}	ω_2=2,n_2=5	12	1	−12
ω=2,n=6	{(1,X),(X,1)}	ω_1=1,n_1=1	1	2	7
	{(1,2),(2,1)}	ω_2=2,n_2=6	4	1	−4
ω=3,n=1	{(1,X,X),(X,X,1),(X,1,X)}	ω_1=1,n_1=1	1	3	−13
	{(1,2,X),(2,1,X),(2,X,1),(1,X,2),(X,1,2),(X,2,1)}	ω_2=2,n_2=1	6	3	18
	{(1,3,2),(3,2,1),(2,1,3),(3,1,2),(2,3,1),(1,2,3)}	ω_3=3,n_3=1	8	1	−8
ω=3,n=2	{(1,X,X),(X,X,1),(X,1,X)}	ω_1=1,n_1=1	1	3	−19
	{(2,1,X),(1,2,X)}	ω_2=2,n_2=2	1	1	9
	{(1,X,2),(2,X,1),(X,1,2),(X,2,1)}	ω_3=2,n_3=1	6	2	18
	{(3,1,2),(1,3,2),(1,2,3)}	ω_4=3,n_4=2	12	1	−12
ω=4,n=1	{(1,X,X,X),(X,X,X,1),(X,X,1,X),(X,1,X,X)}	ω_1=1,n_1=1	1	4	−5
	{(X,X,1,2),(X,X,2,1),(1,2,X,X),(2,1,X,X)	ω_2=2,n_2=1	6	6	6
	(X,1,X,2),(X,2,X,1),(2,X,1,X),(1,X,2,X)	ω_2=2,n_2=1
	(1,X,X,2),(2,X,X,1),(X,1,2,X),(X,2,1,X)}	ω_2=2,n_2=1
	{(1,3,2,X),(3,1,X,2),(2,X,1,3),(X,2,3,1)	ω_3=3,n_3=1	8	4	0
	(3,2,1,X),(1,X,3,2),(X,1,2,3),(2,3,X,1)	ω_3=3,n_3=1
	(2,1,3,X),(X,3,1,2),(1,2,X,3),(3,X,2,1)	ω_3=3,n_3=1
	(2,X,3,1),(X,2,1,3),(1,3,X,2),(3,1,2,X)	ω_3=3,n_3=1
	(1,X,2,3),(3,2,X,1),(2,3,1,X),(X,1,3,2)	ω_3=3,n_3=1
	(3,X,1,2),(1,2,3,X),(X,3,2,1),(2,1,X,3)}	ω_3=3,n_3=1
	{(1,4,3,2),(4,1,2,3),(3,2,1,4),(2,3,4,1)	ω_4=4,n_4=1	2	1	−2
	(4,3,1,2),(1,2,4,3),(2,1,3,4),(3,4,2,1)	ω_4=4,n_4=1
	(3,1,4,2),(2,4,1,3),(1,3,2,4),(4,2,3,1)	ω_4=4,n_4=1
	(3,2,4,1),(2,3,1,4),(1,4,2,3),(4,1,3,2)	ω_4=4,n_4=1
	(1,2,3,4),(4,3,2,1),(3,4,1,2),(2,1,4,3)	ω_4=4,n_4=1
	(4,2,1,3),(1,3,4,2),(2,4,3,1),(3,1,2,4)}	ω_4=4,n_4=1

For the case of the maximal cluster consisting of sites , see Table 1, we have a full set of 2-tuples of decorations, each associated with a temperature dependent decorated cluster correlation function. It can be noted that the cluster contains two sub-clusters: the point cluster that has decorations ; and the cluster itself which has decorations . For each of these decorations, there are corresponding decorated cluster correlation functions. In this case, the point cluster is contained twice: once in and also in . Similarly the cluster itself is contained once. Thus a point cluster decorated by is transfered to -tuples notation in the maximal cluster as or . This can be captured by the permutation operators (see Table 2) where X stands for an empty cluster where a 0 should be placed. Similarly, the decorations denoted by in the 2-tuples notation are equivalent to the symmetry equivalent decorations and . The symmetry effect can again be captured by permutation operators , where X no longer appears, since the -tuples have the same order () as the -tuples from the maximal cluster . The analysis outlined above for the 2-body cluster with components can be extended to the remaining maximal clusters. The permutation operators that carry out the appropriate decorations of sub-clusters into cluster -tuples notation are detailed in Table 2. It should be noted that, for any 2-body cluster, the point and corresponding pair sub-clusters add up, forming symmetrically unique correlation functions .

Referring to the maximal cluster , from Table 2, there are three contained sub-clusters , and within the maximal cluster . The point sub-cluster has, according to Table 1, four decorations associated with it; the 2-body cluster has six decorations associated with itself; and finally the sub-cluster has 10 decorations associated with itself. The total number of symmetry unique decorated cluster correlation functions, which is smaller than , but can nevertheless fully describe cluster with components, is therefore given by decorated cluster correlation functions. From these 20 symmetry unique decorated cluster correlations, it is possible to generate a total of 64 decorated cluster correlation functions by using the permutation operators in Table 2. The number of permutation operators to transfer decorations from sub-cluster to , is 3; from , six permutation operators, and, from the cluster itself there are three permutation operators. The total number of permutation operators is therefore . The remaining decorations, which are not listed in Table 1, add up to and are obtained by using appropriately the 15 permutation operators corresponding to the maximal cluster . Similarly, the maximal cluster has 34 symmetrically unique decorations and 12 permutation operators; and has 35 symmetry unique decorations with 64 permutation operators (all possible permutations for ).

3.3. Four-Body Probability Functions from Monte Carlo Simulations

After performing semi-canonical Monte Carlo simulations for alloy compositions CrFeMnNi and equiatomic CrFeMnNi, we use the correlation functions calculated at each of the temperatures to study the temperature dependent probabilities of decorated clusters corresponding to the equilibrium configuration of the Monte Carlo super-cell following the formalism described in the Methods section. For each non-decorated cluster and temperature, we generate the corresponding matrices with dimension x and a vector formed by all the decorated correlation cluster functions, with dimension . The Monte Carlo simulations output only the symmetrically unique correlation functions (see Table 1), which are less than the required for matrix operation. In order to generate the full set of from the symmetrically unique correlation functions, we devise a set of permutation operators (see Table 2) that indicate how an arbitrary decoration -tuple of integers is obtained from the symmetrically unique decorations belonging to the non-decorated cluster or one of its sub-clusters (see the third column of Table 2 for identifying the cluster-sub-cluster relation). Table 2 contains all the operators necessary to generate cluster probabilities, , for any possible decoration, , of 10 non-decorated clusters employed in the cluster expansion Hamiltonian.

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 as a function of temperature.

Figure 1

4-body probabilities obtained from the hybrid Cluster Expansion (CE)-Monte Carlo calculations. (a) all the 4-body probabilities for the equiatomic composition CrFeMnNi as a function of temperature; (b) the same as in (a) but for the composition CrFeMnNi.

The most important 4-body bonding configurations in the temperature range from 0 to 3000 K for the two alloy systems (the equiatomic CrFeMnNi and CrFeMnNi) are highlighted in 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 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 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.

For the system with average composition CrFeMnNi, the cluster probability of decorations given by Mn-Mn-Ni-Ni appears to dominate until 1300 K, where the solid solution or random configuration begins. Furthermore, the second most probable cluster configuration in the temperature range 500–1200 K appears to be Cr-Fe-Fe-Fe. Again, these findings are remarkable and the origin of these clusters would need a more detailed discussion.

The prediction of the high probabilities for Mn-Mn-Ni-Ni and Cr-Fe-Fe-Fe clusters in CrFeMnNi alloy composition shown in Figure 1b can be explained by the decoration of the first nearest neighbor four-atom cluster interaction obtained from the CE Hamiltonian. It is important to stress that, in this case, the CE method reproduces the well-known result from the tetrahedron approximation in the CVM [18,19]. The configuration Mn-Mn-Ni-Ni is understood to be related to L1 structure in Strukturbericht notation with Mn and Ni atoms, whereas the composition Cr-Fe-Fe-Fe is related the L1 structure in Strukturbericht notation with Cr and Fe atoms. Both of these structures are depicted in Figure 2 and they are in a full agreement with our first-principles investigations. Indeed, the - structure is predicted by both the DFT and CE Hamiltonian to be one of the ground-state structures not only for the binary Mn-Ni system but also for the quaternary Fe-Cr-Mn-Ni one. The - structure has also been found as the lowest-energy binary structure in the previous study of the phase stability of Fe-Cr binary in the FCC lattice (see Figure 2a from [17]). In particular, the latter structure has also been found to be stable due to the strong anti-ferromagnetic interaction between Cr and Fe from the magnetic cluster expansion (MCE) performed for the ternary CrFeNi system [43]. In that work, Monte Carlo simulations using MCE indicated that ordered magnetic structures in the Ni-rich corner (ferromagnetic) of the system persists until temperatures of 600 K, and that for the atomic compositions between CrFe and CrFe the anti-ferromagnetic order was retained beyond 500 K.

Figure 2

(a) most probable phase at high temperature (disordered structure); (b,c): two most probable ordered phases at low temperature in the equiatomic CrFeMnNi and CrFeMnNi HEAs compositions. Cr, Fe, Mn and Ni are illustrated in blue, red, yellow and green respectively. (a) A1 phase, sites are occupied by Cr, Mn, Fe, and Ni in probabilities determined by their average concentration in the system; (b) L1 phase corresponding to CrFe with Cr and Fe; (c) L1 phase corresponding to MnNi with Mn and Ni.

4. Configuration Entropy in a Cr-Fe-Mn-Ni System

For each temperature value in the range 0–3000 K, the configuration probabilities corresponding to cluster () are used in Equation (17) to obtain the quantities . There is one quantity and a corresponding weight for each of the () sub-clusters: (), (); (); and () included in the maximal cluster (). The total configuration entropy, , is calculated by adding the weighted quantities . By applying the same arguments developed for the cluster (), 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 , 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 () 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].

Figure 3

Composition dependent entropies obtained from Monte Carlo simulations in CE. (a) Composition dependent entropy at fixed temperature 1000 K; (b) Composition dependent entropy at fixed temperature 3000 K.

Figure 4

Temperature dependence of configuration entropy evaluated at various levels of cluster approxinations and compared with the thermodynamic integration result at the equiatomic composition CrFeMnNi.

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 towards the disorder configuration in Figure 2. Similarly, in the alloy with composition CrFeMnNi, 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-CrFe and L1-MnNi structure, respectively.

For the cluster (), 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 (), 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 () 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.

5. 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 CrFeMnNi and CrFeMnNi, 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 CrFeMnNi, there is a high probability of formation of the L1 MnNi phase at low temperatures below 1300 K and for the L1 CrFe phase in the temperature range 500–1200 K. Similarly, for the equiatomic composition, the L1 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 CrFeMnNi multicomponent system has been studied for advanced nuclear applications due to promising irradiation resistance with regards to void swelling.