**New Advances in High-Entropy Alloys**

Editor

**Yong Zhang**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editor* Yong Zhang University of Science and Technology Beijing China

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Entropy* (ISSN 1099-4300) (available at: https://www.mdpi.com/journal/entropy/special issues/ High-Entropy Alloys).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Article Number*, Page Range.

**ISBN 978-3-03943-619-4 (Pbk) ISBN 978-3-03943-620-0 (PDF)**

c 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**






### **About the Editor**

**Yong Zhang** has been a university professor of Materials Science at USTB since 2004. He has published over 200 papers, including two comprehensive review papers published in the journal of *Progress in Materials Science* (*PMS*). He authored a book entitled "High-Entropy Materials: A Brief Introduction", published by Springer Nature. His contributions include proposing a balance parameter to evaluate the configurational entropy effect over the enthalpy effect in liquid state, which has been verified effective to predict the phase formation of multicomponent materials.

Recent publications:

[1] Zhang, Y.; Yan, X.-H.; Liao, W.-B.; Zhao, K. Effects of Nitrogen Content on the Structure and Mechanical Properties of (Al0.5CrFeNiTi0.25)Nx High-Entropy Films by Reactive Sputtering. Entropy 2018, 20, 624.

[2] Wen-Jie Sheng, Xiao Yang, Jie Zhu, Cong Wang, Yong Zhang, Amorphous phase stability of NbTiAlSiNX high-entropy films. Rare Met. (2018) 37(8):682–689.

[3] Zhang, Y.; Yan, X.H.; Ma, J.; Lu, Z.P.; Zhao, Y.H.; Compositional gradient films constructed by sputtering in a multicomponent Ti–Al–(Cr, Fe, Ni) system. JMR, 2018.

[4] Sheng, W.; Yang, X.; Wang, C.; Zhang, Y. Nano-Crystallization of High-Entropy Amorphous NbTiAlSiWxNyFilms Prepared by Magnetron Sputtering. Entropy 2016, 18, 226.

[5] X.H. Yan, J.S. Li, W.R. Zhang, Y. Zhang, A brief review of high-entropy films, Materials Chemistry and Physics 210 (2018) 12–19.

[6] Zhang, W., Liaw, P.K. & Zhang, Y., Science and technology in high-entropy alloys, Sci. China Mater. (2018) 61: 2.

[7] Rui Xuan Li, P. K. Liaw, Yong Zhang, Synthesis of Al x CoCrFeNi High-entropy Alloys by High-gravity Combustion from Oxides, Materials Science and Engineering A 707, 668-673.

[8] Michael C. Gao, Daniel B. Miracle, David Maurice, Xuehui Yan and Yong Zhang, Jeffrey A. Hawk, High-entropy functional materials, J. Mater. Res., Vol. 0, No. 0, 2018.

[9] Xia S, Lousada CM, Mao H, Maier AC, Korzhavyi PA, Sandstrom R, Wang Y and Zhang Y ¨ (2018) Nonlinear Oxidation Behavior in Pure Ni and Ni-Containing Entropic Alloys. Front. Mater. 5:53. doi: 10.3389/fmats.2018.00053.

[10] Christopher M. Barr, James E. Nathaniel II, Kinga A. Unocic, Junpeng Liu, Yong Zhang,YongqiangWang, Mitra L. Taheri, Exploring radiation induced segregationmechanisms at grain boundaries in equiatomic CoCrFeNiMn high entropy alloy under heavy ion irradiation, Scripta Materialia 156 (2018) 80–84.

[11] Shao Lei, et al., A Low-Cost Lightweight Entropic Alloy with High Strength, Journal of Materials Engineering and Performance, ASM International https://doi.org/10.1007/s11665-018-3720-0.

[12] Xing Q, et al., High-Throughput Screening Solar-Thermal Conversion Films in a Pseudobinary (Cr, Fe, V)–(Ta, W) System, ACS Combinatorial Science, 2018.

[13] Li R., et al., Graded microstructures of Al-Li-Mg-Zn-Cu entropic alloys under supergravity, SCIENCE CHINA Materials, 2019.

[14] Xing Q, et al., Mechanical properties and thermal stability of (NbTiAlSiZi)Nx high-entropy ceramic films at high temperature, JMR, 2018.

[15] Li D., et al., AnnealingeffectfortheAl0.3CoCrFeNihigh-entropyalloyfibers, JACs, 2018.

[16] Y. Zhang\*, J. Liu, S. Chen, X. Xie, P. Liaw, K. Dahmen, J. Qiao, Y. Wang, Serration and noise behaviors in materials. Progress in Materials Science. 2017, 90:358–460.

[17] Y. Zhang, Authored Book, High-entropy Materials: A Brief Introduction. Springer-Nature Publishing. 2019.

[18] X. Yan, J. Ma, Y. Zhang\*, High-throughput screening for biomedical applications in a Ti-Zr-Nb alloy system through masking co-sputtering. Science China: Physics, Mechanics and Astronomy. 2019, 62(9):996111.

### *Editorial* **New Advances in High-Entropy Alloys**

**Yong Zhang 1,2,3,\* and Ruixuan Li <sup>1</sup>**


Received: 6 October 2020; Accepted: 13 October 2020; Published: 15 October 2020

**Keywords:** high-entropy alloys; mechanical behaviors; high-entropy film; low-activation alloys

Exploring new materials is an eternal pursuit in the development of human civilization. In recent years, people have tended to adjust the order/disorder degree to explore new materials. The order/disorder can be measured by entropy, and it can be divided into two parts: the topological ordering and the chemical ordering. The former mainly refers to the ordering in the spatial configuration, e.g. amorphous alloys [1] with short-range ordering but without long-range ordering; the latter mainly refers to the ordering in the chemical occupancy, e.g. high-entropy alloys (HEAs) where components can replace each other. HEAs, in sharp contrast to traditional alloys based on one or two principal elements, have one striking characteristic: the unusually high entropy of mixing. They have not been too much noticed until a review paper entitled: "Microstructure and properties of high entropy alloys" published in 2014 on the journal of Progress in Materials Science [2]. Lots of reports showing they exhibit five recognized performance characteristics, namely strength-plasticity trade-off breaking, irradiation tolerance, corrosion resistance, high impact toughness within a wider temperature range, and high thermal stability [3]. So far, the development of HEAs has gone through three main stages: 1. quinary equal-atomic single-phase solid solution alloys; 2. quaternary or quinary non-equal-atomic multiphase alloys; 3. medium-entropy alloys, high entropy fibers, high entropy films, lightweight HEAs, etc. Nowadays, deeper research on HEAs is urgently needed.

This Special Issue includes forty-two research articles and six review papers dedicated to the frontiers of high entropy materials, from exploring the microstructures by experiment to revealing the structure-performance relationship by simulation. These innovative studies will provide new methods to solve challenging problems and open up new possibilities for emerging research fields.

Configuration entropy plays an important role in the microstructures and properties of HEAs, which is believed to stabilize disordered solid solution phases over intermetallic compounds by lowering the Gibbs free energy. Traditionally, configuration entropy was computed by the empirical formula. In Fernández-Caballero's work [4], a new formalism based on a hybrid combination of the Cluster Expansion 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. Experimentally, Haas et al. [5] used differential scanning calorimetry to determine the thermal entropy and compared it to the configurational entropy. The contributions of entropy and enthalpy to the Gibbs free energy was calculated and examined by them.

In addition to the above basic research on the entropy of the alloy system, research on HEAs mainly focuses on their composition design, preparation and processing technology, mechanical properties, and physical and chemical properties. Other main contents of this issue are listed as follows.

#### **1. Compositional Design**

Because HEAs are characterized by multiple components and high contents of alloying elements, the complexity of their compositional design is greatly increased compared with traditional alloys. Some empirical rules are used early on to realize simple prediction of phase structure. With the development of computer technology and simulation technology, phase diagram and first-principles methods are now widely used to help alloy design. Now, machine learning is also popular by using various algorithms. Here, our Special Issue contains three articles using phase diagram and first principles for compositional design.

A phase diagram is the geometric description of the system under a thermal equilibrium and it is the basis for the study of solidification, phase transformation, crystal growth, and solid-phase transformation. Since the 1970s, the calculation of phase diagrams (CALPHAD) based on the thermodynamic theory has become a new trend. Gorsse et al. [6] and Klaver et al. [7] conducted an in-depth study on the use of CALPHAD methods for composition design. On the other hand, Zheng et al. [8] used first principles to calculate the elastic properties of seventy different compositions in the refractory V–Mo–Nb–Ta–W system. Effects of different elements to the modulus are precisely calculated so as to optimize the composition in this refractory system.

#### **2. Preparation and Processing**

Traditional preparation methods of bulk HEAs mainly include vacuum arc melting, vacuum induction melting, mechanical alloying, and so on. In order to obtain a certain crystal orientation, directional solidification can also be employed. There are also some new preparation methods for bulk HEAs, such as high-gravity combustion synthesis and additive manufacturing. Chen et al. [9] reviewed the application of additive manufacturing methods in the fabrication of HEAs. Compared with the casting counterparts, HEAs prepared by additive manufacturing are found to have a superior yield strength and ductility as a consequence of the fine microstructure formed during the rapid solidification in the fabrication process. As a result, this is an effective method to improve their comprehensive properties.

Different processing methods are closely related to performance. Our Special Issue contains two research papers on the use of annealing to optimize alloys and two review papers on alloys using high pressure and welding.

Annealing is considered to be an effective method to improve the microstructure and properties of alloys, and different annealing temperature and time are closely related to comprehensive properties. Zhuang et al. [10] investigated the effect of annealing temperature on the microstructure, phase constituents, and mechanical properties of Al0.5CoCrFeMoxNi (x = 0, 0.1, 0.2, 0.3, 0.4 and 0.5) HEAs at a fixed annealing time (10 h). They found that the alloys annealed at 80 ◦C showed relatively fine precipitations and microstructures and, therefore, higher hardness and yield stress. Sathiyamoorthi et al. [11] put the high-pressure torsion-processed CoCrNi alloy with a grain size of ~50 nm in different annealing conditions and they explored the optimal processing technology. The sample annealed at 700 ◦C for 15 min exhibited a remarkable combination of tensile strength (~1090 MPa) and strain to failure (~41%).

Pressure, as another fundamental and powerful parameter, has been introduced to the experimental study of HEAs. Many interesting reversible/irreversible phase transitions that were not expected or otherwise invisible before have been observed by applying high pressure. Zhang et al. [12] reviewed recent results in various HEAs obtained using in situ static high-pressure synchrotron radiation X-ray techniques and provided some perspectives for future research.

Welding is an important emerging area with significant potential impact to future application-oriented research and technological developments in HEAs. The selection of feasible welding processes with optimized parameters is essential to enhance the applications of HEAs. Guo et al. [13] reviewed the key recent works on welding of HEAs in detail, focusing on the research of

main HEA systems when applying different welding techniques. They also highlighted the future challenges and main areas to research.

#### **3. Mechanical Behaviors**

Research on the mechanical properties of HEAs is the most extensive. The main research hotspots are three types of alloys: 1, 3D transitional-group-element HEAs; 2, 3D transitional-group-element HEAs with Al or Ti added; 3, refractory high-entropy alloys with excellent high-temperature properties. In addition, other trace alloying elements, such as Mo, Nb, Zr, etc., have also been added to the alloy in order to study the effect of their content on the microstructures and properties of the alloy.

HEAs based on 3D transition metals, such as CrCoNi and CrMnFeCoNi alloys, are revealed to have remarkable mechanical properties. Especially the CoCrFeMnNi alloy is famous as the Cantor alloy for its high plasticity and comparable strength. The tensile creep behavior of the CoCrFeNiMn HEA was systematically investigated by Cao et al. [14] over an intermediate temperature range (500–600 ◦C) and applied stress (140–400 MPa). The alloy exhibited a stress-dependent transition from a low-stress region to a high-stress region, which was characterized by the dynamic recrystallization and the grain-boundary precipitation. This alloy, deformed at high strain rates (900 to 4600 s−1), was investigated by Wang et al. [15]. The yield strength was sensitive to the change of high strain rates, and serration behaviors were also observed on the flow stress curves, which could be predicted by the Zerilli–Armstrong constitutive equation. Bu et al. [16] performed an in situ atomic-scale observation of deformation behaviors in a nano-scaled CoCrCuFeNi alloy. Exceptional strength was realized, and the deformation behaviors, including nano-disturbances and phase transformations, were distinct from those of corresponding bulk HEAs. The same result was confirmed in Mridha's [17] experiments. From the perspective of theory and simulation, Ikeda et al. [18] investigated the impact of compositional fluctuations in the vicinity of stacking faults on CrCoNi and CrMnFeCoNi by employing first-principles calculations. Their research and experimental results complement and perfect each other.

Adding other alloying elements is one of the research hotspots, in order to further improve the strength. The addition of Nb promoted the precipitation of nano-phases and thus increased the strength of the alloy, and Mo addition in the CoCrFeNi alloy effectively helped to increase the corrosion resistance. Detailed work was carried out by Han et al., Wang et al., and Tsau et al. [19–21]. The effect of Zr addition on the CoCrFeNiMn alloy was studied by Zhang et al. [22]. They prepared the alloy by using the ZrH2 powders and a mechanical alloying technique. Multi-phase microstructures formed in the alloys, which can be attributed to the large lattice strain and negative enthalpy of mixing, caused by the addition of Zr. Sun et al. [23] also employed the mechanical alloying technique and prepared the CoCrNiCuZn alloy. Pd addition promoted the local- and long-range lattice distortions in CoCrFeNi alloy and also had effect on the phase stability and phase transformation. This phenomenon was revealed by Zhang et al. [24]. Tan et al. [25] then studied the effect of Mn addition on the microstructures and mechanical properties of this CoCrFeNiPd alloy.

Alloying Al and Ti in 3D transitional Co-Cr-Fe-Ni HEAs shows a strong influence on microstructure and phase composition, as well as the ability to decrease the density. They tend to strengthen the alloy by increasing the lattice distortion, as confirmed by Wu et al. [26]. Löbel et al. [27] studied the influence of Ti in the AlCoCrFeNiTix alloy, and Zhang et al. [28] characterized the partially recrystallized CoCrFeNiTi0.2 alloy. It seems that alloying with Ti leads to an increase in microstructural heterogeneity. Similarly aimed at the Al-Co-Cr-Fe-Ni-Ti alloy after the composition adjustment, Manzoni et al. [29] studied the microstructure and properties of the Al10Co25Cr8Fe15Ni36Ti6 alloy and Haas et al. [30] studied those of the Al10Co25Cr8Fe15Ni36Ti6 alloy. Both kinds of alloy showed beneficial particle-strengthened microstructures. Riva et al. [31] investigated the effect of Sc alloying on the Al2CoCrFeNi, Al0.5CoCrCuFeNi, and AlCoCrCu0.5FeNi HEAs. It caused grain refinement as well as hardness and electrical conductivity increases (up to 20% and 14%, respectively). Furthermore, the solid-state phase transformation kinetics of as-cast and cold rolling deformed Al0.5CoCrFeNi HEAs have been investigated by Wang et al. [32] using the thermal expansion method. Lightweight AlCrMoTi

and AlCrMoTiV HEAs were designed by Kang et al. [33] and they were confirmed as solid solutions with minor ordered B2 phases. They also have superb specific hardness compared to that of commercial alloys.

With the help of their excellent mechanical properties, high-entropy thin film developed by a magnetron sputtering technique has attracted attention, which has exciting potential to make small-structure devices and precision instruments with sizes ranging from nanometers to micrometers. Liao et al. [34] fabricated Al0.3CoCrFeNi film, and Zhang et al. [35] prepared (Al0.5CrFeNiTi0.25)Nx high-entropy films. It showed that the phase structure changes from the amorphous to the face-centered cubie (FCC) structure with increasing N content, which was proved to be related with the atomic size difference in the alloy system. Similarly, (AlCrTiZrV)SixN high entropy film was prepared and characterized by Niu et al. [36]. Then, an Al0.5CoCrCuFeNi coating was successfully prepared on an AZ91D magnesium alloy surface [37], and an AlCoCrFeNi coating was prepared on the stainless steel substrate [38]. The coating obviously had significant high hardness, high wear resistance, and corrosion resistance, and can integrate the properties of the substrate.

The development and preparation of composites have realized the complementary advantages of a variety of different component materials, which can further improve the overall performance of alloys. Zhang et al. [39] prepared a high-entropy ceramic composite by using HfC, Mo2C, TaC, and TiC in pulsed-current processing. Furthermore, Li et al. [40] used a mixture of carbides and oxides in the preparation of an NbMoCrTiAl HEA. With the presence of approximately 7.0 vol. % Al2O3 and 32.2 vol. % TiC-reinforced particles, the compressive fracture strength of the composite reached 1542 MPa, and this was increased by approximately 50% compared with that of the as-cast NbMoCrTiAl HEA. Furthermore, two types of novel HEA/diamond composites were reported, providing a new model for material development [41,42].

It can be found that nano-precipitation phases appear widely in HEAs, which play a vital role in improving their strength and plasticity. In this case, Wang et al. [43] summarized precipitation behavior and precipitation strengthening in HEAs comprehensively, including the morphology evolution of second-phase particles and precipitation strengthening mechanisms. They claimed that the challenge in the future is to design a stable, coherent microstructure in different solid-solution matrices.

Besides room-temperature mechanical properties, HEAs tend to exhibit great high-temperature properties for their high-entropy stabilization effect. The Hf-Nb-Ta-Ti-Zr refractory HEA shows great high-temperature properties and room-temperature properties. Zýka et al. [44] presented investigations of the room-temperature tensile mechanical properties of selected three- and four-element medium entropy alloys derived from the Hf-Nb-Ta-Ti-Zr system, and they found that it is the five-element HEA alloy that exhibits the best combination of strength and elongation. Tseng et al. [45] aimed to reveal the effects of Mo, Nb, Ta, Ti, and Zr on mechanical properties of equiatomic Hf-Mo-Nb-Ta-Ti-Zr alloys. Another VCrFeTaxWx HEA, which has great high-temperature properties and excellent heat-softening resistance, was studied and recommended by Zhang et al. [46].

#### **4. Physical and Chemical Properties**

In addition to mechanical properties, research on the physical and chemical properties of HEAs mainly focuses on soft magnetic properties, radiation resistance, electrical properties, and corrosion resistance.

Because HEAs contain more principal elements, including magnetic elements, such as Fe, Co, Ni, and Mn, and these magnetic elements are mixed with other elements, HEAs exhibit different magnetic properties compared to traditional alloys. Li et al. [47] exhibited the effects of Sn addition on the soft magnetic properties of dual-phase FeCoNi(CuAl)0.8Snx (0 ≤ x ≤ 0.10) HEAs. They showed that saturation magnetization of the alloy increased greatly, while the remanence (Br) decreased after the addition of Sn. Furthermore, thermomagnetic curves indicated that the phases of the alloy will transform from FCC with low Curie temperatures (Tc) to the body-centered cubic (BCC) phase with high Tc at temperatures of 600–700 K. For the FeSiBAlNi alloy, the effects of Co and Gd addition combined

with subsequent annealing on microstructures were investigated by Zhai et al. [48]. FeSiBAlNi-based HEAs also possessed soft magnetism, especially the as-annealed FeSiBAlNiGd alloy, which can be ascribed to the formation of Gd-oxides.

Some HEA systems also have good radiation resistance. Disordered solid-solution phase structure is mainly formed in HEAs. The biggest feature of their structure is the large lattice distortion caused by the difference in atomic size and the high configuration entropy. Therefore, atomic-level stress may be formed. In addition, the sluggish diffusion effect is caused by the interaction of multiple components. As a result, HEAs have shown excellent performance in radiation resistance, which provides new ideas for the development of nuclear materials. CoCrFeCuNi alloys were irradiated with a 100 keV He<sup>+</sup> ion beam by Wang et al. [49]. It was indicated that Cu-rich phases were favorable sites for the nucleation and gathering of He bubbles. At ion doses of 2.5 <sup>×</sup> 10<sup>17</sup> ions/cm2 and 5.0 <sup>×</sup> 10<sup>17</sup> ions/cm2, the HEAs showed obvious hardening, which could be attributed to the formation of large amounts of irradiation defects.

#### **5. Outlook**

The development of HEAs has only been occurring for no more than ten years until now, and people's understanding of this material is still in its infancy. This Special Issue summarizes some of the latest developments in HEAs. Even though not all the questions related to HEAs have been given a reasonable explanation in this issue, it does pay attention to some typical research hotspots, and significant progress can, indeed, be seen. We believe this special issue will guide further research and promote the property breakthrough of HEAs.

**Funding:** Y.Z. would like to thank the financial support from the National Science Foundation of China (Grant No. 51671020), the Fundamental Research Funds for the Central Universities (Grant No. FRF-MP-19-013), and the project of Al-Mg based medium-entropy alloys with Dongguan Eontec company.

**Acknowledgments:** We express our thanks to the authors of the above contributions, and to the journal *Entropy* and MDPI for their support during this work.

**Conflicts of Interest:** The author declare no conflict of interests.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Review* **Deformation Behavior of Bulk Metallic Glasses and High Entropy Alloys under Complex Stress Fields: A Review**

#### **Shunhua Chen 1,2,\*, Jingyuan Wang 1, Lei Xia <sup>3</sup> and Yucheng Wu 2,4,\***


Received: 1 December 2018; Accepted: 9 January 2019; Published: 11 January 2019

**Abstract:** The plastic deformation of bulk metallic glasses (BMGs) depends significantly on applied stress states, and more importantly, in practical applications of BMGs as structural materials, they always deform under complex stress fields. The understanding of deformation behavior of BMGs under complex stress fields is important not only for uncovering the plastic deformation mechanisms of BMGs, but also for developing BMG components with excellent mechanical performance. In this article, we briefly summarize the recent research progress on the deformation behavior of BMGs under complex stress fields, including the formation and propagation of shear bands, tunable macroscopic plasticity, and serrated plastic flows. The effect of complex stress fields on the plastic deformation mechanisms of BMGs is discussed from simple stress gradient to tailored complex stress fields. The deformation behavior of high entropy alloys (HEAs) under complex stress states has also been discussed. Challenges, potential implications and some unresolved issues are proposed.

**Keywords:** bulk metallic glass; complex stress field; shear band; flow serration; deformation mechanism; high entropy alloy

#### **1. Introduction**

As a new class of structural materials, bulk metallic glasses (BMGs) are considered poised for widespread engineering applications [1]. With non-ordered atomic structures, plastic deformation in BMGs is accommodated by initiation and propagation of shear bands, rather than dislocation movements as in crystalline alloys [2–4]. Under applied stress, plastic strain is localized within a thin layer of atoms with a width of about 2–210 nm, forming a shear band [5–10]. Extensive studies have shown that the initiation and propagation of shear bands are significantly dependent on applied stress fields. For example, the burst of a shear band relies on applied shear stress, and a shear band can be initiated only when the applied shear stress exceeds a critical value [5,11–14]. In addition, in practical applications of BMGs as structural materials, they always deform under complex stress states [15]. Studies have shown that BMGs can demonstrate more plastic deformation under complex stress states [16–18]. This implies that although some BMGs are brittle under conventional compression/tension tests with relatively-uniform stress fields, they can still have large capability to deform plastically in practical structural applications. For example, BMG foams can exhibit large nominal plasticity, which is much larger than solid BMG specimens [19–21]. BMG honeycombs and cellular BMG structures also demonstrate similar behavior [22–24]. The understanding of the

deformation behavior of BMGs under complex stress fields is therefore important not only for uncovering the plastic deformation mechanisms of BMGs, but also for developing BMG components with excellent mechanical performance. Although many papers have reviewed the atomic structures, mechanical and physical properties, and functional and structural applications of BMGs [2–4,25–40], a comprehensive review on the plastic deformation behavior of BMGs under complex stress fields has not been reported. In this work, the research progress on the plastic deformation behavior of BMGs under complex stress fields is reviewed, including the initiation and propagation of shear bands, macroscopic plastic deformation behavior, criticality of plastic flows, and transition of deformation modes. The recent contributions on the deformation behavior of high entropy alloys (HEAs) under complex stress fields have also been discussed. The challenges, unresolved issues as well as future research directions are proposed.

#### **2. Initiation of Shear Bands under Complex Stress Fields**

Driven by applied shear stress, the dilatancy of atoms, associated with excess free volume, results in the decrease of viscosity localized in thin layers of shear bands [41,42]. The localization process causes local heating/softening, which was evidenced by experimental observations through fusible coatings [43–45]. The initiation of shear bands in BMGs was studied extensively during the past decades. Greer et al. [4] have summarized the formation mechanisms into three scenarios: the percolation of homogeneous nucleated shear transformation zones (STZs), stemming from the intrinsic structural fluctuations; the nucleation from the sites with extrinsically introduced stress concentrators, such as casting defects; and a two-consecutive-stage formation, including the burst of a viable band from activated STZs at stress-concentrated sites and a following rapid sliding process. Although the shear-banding mechanisms in BMGs are still under debate, all three kinds of scenarios share a common point that the initiation of shear bands is significantly dependent on applied stress fields. Johnson and Samwer [46] have shown that the potential barrier for activating STZs at a given shear stress, *<sup>τ</sup>*, can be expressed as *<sup>W</sup>*<sup>τ</sup> = 4*RG*0*γ*C2[(*τ*<sup>C</sup> − *<sup>τ</sup>*)/*τ*C] 3/2*ζ*Ω, where *τ*<sup>C</sup> is the critical shear stress of the BMG, *γ*<sup>C</sup> is the critical shear strain, *R* is a constant parameter, *G*<sup>0</sup> is the shear modulus of unstressed BMGs at 0 K, Ω is the volume of STZs, and *ζ* is a correction factor. With the change of applied shear stress (*τ*) under complex stress fields, the activation process of STZs will be varied accordingly, resulting in differences in the initiation of shear bands. On the other hand, the annealing of stressed MGs can be used to shape MG films without embrittlement, where the change of free volume in stressed BMGs may also cause different deformation behavior [47]. Hufnagel et al. [48] have indicated that although many previous studies have focused on the deformation behavior of BMGs under uniaxial stress states, the understanding of the deformation behavior as well as structural evolution under more complex stress states is urgently needed to uncover the shear-banding mechanisms, such as the shear localization process. The research progress on the formation mechanisms of shear bands in BMGs has been summarized in many review papers [4,27,30,34,35,48–54], however, the formation of shear bands under complex stress fields has rarely been emphasized. Here, the initiation of shear bands under more complex stress fields is summarized.

With well-known size effect, MGs demonstrate different deformation behavior at submicron scales, for example, larger elastic limit, higher ductility, and more homogeneously plastic deformation with necking effect [55–59]. The introducing of stress concentrators in submicron-sized specimens will facilitate the initiation of shear bands, giving better insight into the shear band formation mechanisms. A simple method to introduce stress concentrators is to create notches in testing specimens. Based on MD simulations of some notched CuZr MGs, Sha et al. [60] have reported that shear bands emanate from the notch root when the plastic zone size ahead of the notches increases beyond a critical value. As shown in Figure 1, the notches can serve as stress concentrators to facilitate the localization of plastic zones, leading to the initiation of shear bands. Following studies have examined the effect of notch geometries and sizes on the formation of shear bands [61–63]. They have shown that the initiation of shear bands can be tailored by introducing notches with appropriate geometries and sizes,

which can serve as stress concentrators for highly localized shear strains. By creating surface roughness in MG nano-pillars can also achieve similar localized shear strains [64]. It can then be concluded that under complex stress fields, stress concentrators serve as origin sites for the initiation of shear bands due to higher localized shear strains. The authors would like to point out that at sub-micron scales, the formation of localized shear bands can also be suppressed till fracture due to the change of deformation modes [65]. This phenomenon has also been observed in notched MG specimens with sub-micron sizes, where the change of notch sizes can tune the deformation modes from localized shear-banding to homogeneous softening and necking [61–63]. A summary on the homogeneous deformation behavior of BMGs under complex stress fields will be given in Section 5 in this paper, and thus not described in detail here.

**Figure 1.** MD simulation results showing the initiation and propagation of shear bands in some notched Cu50Zr50 MG specimens, where the color denotes the localized shear strains. Reprinted from [60] with permission of 2013 AIP Publishing LLC.

Despite the efforts to understand the initiation of shear bands in BMGs under complex stress fields, some problems are still unresolved. For example, how to control the propagation of shear bands through the design of complex stress distributions? With well-known size effect, what's the difference between the mechanisms on the formation of shear bands at submicron scales and macroscopic scales? Further understandings of the initiation and propagation of shear bands under complex stress fields could add more knowledge to the plastic deformation mechanisms of BMGs.

#### **3. Macroscopic Deformation Behavior under Gradient Stress Distribution**

The macroscopic plastic deformation in BMGs is accommodated by the initiation and propagation of shear bands. The rapid propagation of shear bands causes catastrophic failures, while the confinement of the propagation of shear bands can enhance the macroscopic plasticity. Moreover, the impeding of the propagation of shear bands is also beneficial for the initiation and bifurcation of more shear bands, leading to more plastic deformation. The applied complex stress fields in BMGs not only bring to stress concentrators, serving as sites for the initiation of shear bands, but also

influence the propagation of shear bands, even under simple gradient stress distributions. Studies on the deformation behavior of BMGs under gradient stress distributions, for example, the specimens with tailored sample geometries and surface treatments, are discussed here.

#### *3.1. Stress Gradient Resulting from Tailored Sample Geometry (Loading Angle)*

Under compression tests of BMGs, by tailoring a tilted angle between the end surface and loading platen, as shown in Figure 2a inset, significant improvement in the macroscopic plasticity, associated with work-hardening-like behavior, was observed [16]. Based on the scanning electron microscopy (SEM) observations and Finite Element Modeling (FEM) results, Chen et al. [16] have shown that the propagation of shear bands can be stopped in the regions with relatively lower stress concentration orders, leading to the multiplication of more shear bands associated with greatly enhanced macroscopic plasticity. Similar large nominal plasticity was also observed in some other BMG samples with tilted surface angles [17,66,67]. The tilted geometry constrained the initiation of shear bands in the regions with localized strains, and at the same time hindered the propagation of shear bands, leading to the multiplication of more shear bands. Despite the brittleness, some BMG communities still have certain plasticity under compression tests. However, almost all communities of BMGs are brittle under tensile loadings. The investigations on the plastic deformation behavior of the tensile side of bending BMG specimens have indicated that more plastic deformation was found due to the presence of tensile stress gradient [68,69]. The gradient stress distributions may also be possible to be used to improve the plastic deformation behavior of BMGs under tensile loadings. Chen et al. [70] have then introduced gradient stress distribution into some tensile specimens by tailoring tilted angles in Z-shaped BMG specimens, as shown in Figure 2b. It is interesting to find that the tilted BMG specimens have demonstrated more plastic deformation in the stress-concentrated regions, similar to compressive testing results. This implies that although BMGs are brittle in conventional compression/tension tests with relatively uniform stress distributions, they may still demonstrate more plastic deformation in practical engineering applications, where they always deform under complex stress fields. Additionally, with the presence of stress gradients, BMGs can also exhibit attractive performance for engineering applications, such as higher reliability in nominal strains [71] and less dependence on the change of loading rates [72]. The understanding of the deformation behavior of BMGs under stress gradients opens a new window for the practical applications of BMGs.

**Figure 2.** (**a**) Effect of stress gradient on the compressive plastic deformation behavior of BMGs, where the gradient stress distribution was introduced by tailoring a loading angle between the plateau and the specimen surface. A and B are macroscopic nominal stress-strain curves of the specimens without and with stress gradient, respectively (adapted from [16] with permission of 2008 AIP Publishing LLC.); (**b**) Effect of stress gradient on the deformation behavior of BMGs under tensile loadings, where the gradient stress distribution was introduced by designing Z-shaped specimens (adapted from [70] with permission of Elsevier).

#### *3.2. Effect of Surface Residual Stress*

In engineering applications of structural materials, residual stress is usually created to improve the mechanical performance. In BMG communities, the introducing of surface residual stress using techniques, such as shot peening, laser melting and mechanical attrition, has also been investigated to improve the macroscopic plastic deformation behavior. For example, Zhang et al. [73] have treated some Zr-based BMGs using shot peening, and enhanced the bending plasticity significantly (Figure 3). A maximum surface plastic strain of about 0.35% was obtained while the as-cast BMG specimen only has a limited value of about 0.15%. They have shown that the compressive residual stress on the surface of shot-peened specimens can suppress the propagation of shear bands and cracks. Combined with the pre-existed shear bands nucleated underneath the surface, more shear bands will be initiated to achieve larger macroscopic plasticity. By laser surface melting, gradient distribution of both tensile and compressive residual stresses as well as atomic structural changes can also be employed to improve the macroscopic plasticity of BMGs [74,75]. In addition, Wang et al. [76] have obtained an obvious enhancement in the tensile plasticity of BMGs by surface mechanical attrition treatment. Gradient atomic structures were formed during the mechanical attrition process, resulting in a gradient distribution of residual stress with both compressive and tensile stresses. They have shown that the gradient distribution of residual stress can cause the formation of more shear bands, and at the same time, delay the cavitation effect, resulting in a work-hardening mechanism associated with enhanced plasticity [76].

**Figure 3.** Three-point bending behavior of BMGs with surface treatment by shot peening. Reprinted from [73] with permission of Springer Nature.

Despite the use of different methods/techniques to introduce gradient stress distributions to BMGs, it can be concluded that the localized regions with high orders of stress concentration can promote the initiation of more shear bands, and the regions with less stress concentrations impede their propagations. This will further induce the bifurcation and multiplication of shear bands, leading to better macroscopic plasticity. By design of appropriate distributions of more complex stress fields, it is able to tailor the formation and propagation of shear bands, achieving controllable plastic

deformation behavior. In fact, many studies have been devoted to achieving tunable plastic deformation behavior of BMGs by designing complex stress fields, which are given in the following section.

#### **4. Tunable Plastic Deformation Behavior under Tailored Complex Stress Fields**

During the past decade, tunable plastic deformation behavior has been widely achieved by extensive studies based on the design of complex stress fields, for example, to guide the propagation of shear bands and cracks, to tune the criticality of the plastic-flow dynamics, and to obtain large apparent plasticity/elongations. The studies are of significance for giving more insight into the shear-banding mechanisms of BMGs, and obtaining better macroscopic plastic deformation performance for engineering applications.

#### *4.1. Guiding the Propagation of Shear Bands and Cracks*

The initiation of shear bands emanating from notch roots was widely studied to evaluate the fracture behavior of BMGs [77,78]. Tandaiya et al. have shown that by changing the mixity of loading modes (I/II) of bending specimens can tune the plastic zone sizes as well as their distributions [79,80].

Under mixed mode loading conditions, the notch root deforms as that one part sharpens and the other part blunts. The increase of mode II component can enlarge the plastic zone sizes and enhance the localized strain levels ahead of the notch tips, resulting in controllable directions for shear-banding and micro cracks. For some notched BMGs, Yi et al. [81] have demonstrated that by creating pre-existed shear bands can change the plastic deformation behavior as well as the crack propagations around the notch root. To guide and deflect the cracks in those notched BMGs can improve the fracture resistance behavior, which is useful for toughening BMG components with stress concentrations in engineering applications. The controlled shear-banding behavior of notched BMGs was understood by Yang et al. from a perspective of "multiple shear band deformation mechanisms", differing from the conventional materials fracture mechanics [82]. The proposed fracture mechanisms agree well the fracture morphologies of notched BMGs, giving more theoretical understandings on the controlling of shear-banding in notched BMGs. Li et al. [83] have also conducted theoretical analysis on the controlling of shear-banding in BMGs using an instability theory. As shown in Figure 4, the formation of shear bands in notched specimens with mode mixities of 0.5 and 0.75 was successfully predicted by instability analysis, which is highly in line with the FEM results and experimental observations [83]. Such a theory may be further used to predict the shear-banding behavior of BMGs under complex stress fields in practical engineering applications.

**Figure 4.** *Cont*.

**Figure 4.** Shear bands formation in notched specimens with mode mixity (Me) of 0.5 (**a**,**b**) and 0.75 (**c**,**d**), respectively (Me = 1 and 0 stand for mode-I and mode-II loading conditions, respectively). (**a**,**c**) show the predictions by an instability theory, and (**b**,**d**) are the FEM results [83].

#### *4.2. Tunable Criticality in Flow Serrations*

The discontinuous plastic flows in crystalline alloys display serrations or jerky flows in stress-strain curves, which are known as the Portevin-Le Chatelier (PLC) effect [84,85]. Although BMGs have vanished crystalline lattices, the plastic deformation in BMGs is limited to inhomogeneously-localized shear bands, which also results in discontinuous stress-strain behavior similar to the PLC effect in crystalline alloys, i.e., serrated plastic flows. A great consideration of research has shown that the serrated plastic flows of BMGs can evolve to a power-law scaling criticality [40,86–94]. The criticality of the plastic-flow dynamics is usually correlated to larger macroscopic plasticity [86]. The criticality of the flow serrations can be attributed to the multiplication and intersections of shear bands, which suggests that the burst of shear bands may be intrinsically correlated [86,95]. However, for BMGs under tensile loadings, the catastrophic failures lead to limited serrations where the understanding of the plastic-flow dynamics is very challenging. By designing complex stress fields through double-side notches, Chen et al. [96] have reported tunable criticality in the plastic flows of BMGs under mixed mode (I/II) loading conditions. As shown in Figure 5, by the multiplication of shear bands within the regions with complex stress fields, a stable plastic flow stage was observed, resulting in the delay of catastrophic failures [96]. The two-stage plastic flows were also observed under compression tests, where the formation of new shear bands tends to delay catastrophic failures [97]. For the notched specimens with complex stress fields (Figure 5), the flow serrations during the stable plastic flow stages have smaller magnitudes, and evolve to a power-law scaling. The findings suggest that it may be possible to obtain more flow serrations in tensile BMGs by introducing of complex stress fields, where enough serration data could be collected to give more insight into the plastic-flow dynamics of BMGs under tensile loadings.

Thereafter, by tailoring single-side notches with varying radii, complex stress fields with varying stress concentration factors have been introduced to some tensile BMG specimens [76]. Despite different trends in the evolution of amplitudes from compression tests, the flow serrations of the single-side notched BMG specimens also demonstrate power-law criticality within the stable plastic flow stages, similar to the results under compression and nano-indentation tests [76]. The power-law criticality in BMGs might be a universal rule for all kinds of loading conditions. Based on such assumption, the catastrophic failures in BMGs can then be delayed or avoided in practical applications by designing complex stress fields, regardless of loading conditions. Additionally, under complex stress

fields, BMGs can also have high uniformity for accumulating elastic energy during the stress-arising process of the flow serrations [98]. By tailoring of complex stress fields is helpful for uncovering the mechanisms of the plastic flow serrations in BMGs and to achieve better plastic performance.

**Figure 5.** (**a**) Serrated plastic flows of three notched BMG specimens at a loading rate of 0.06 mm/min., where L = 0, 0.2, and 0.4 mm for L00, L20 and L40 specimens, respectively. (**b**) The |d(load)/d(t)| versus time relationships showing the observation of two-stage plastic flows [96].

#### *4.3. Achieving Large Macroscopic Plasticity/Axial Elongation*

The improved macroscopic plastic deformation behavior of BMGs under complex stress fields can be observed in some specimens with casting defects, such as voids [99]. In order to create complex stress fields to block the propagation of shear bands, Zhao et al. [100–103] have introduced notches into BMGs and obtained greatly enhanced nominal strains. An example showing the improvement of nominal strains of notched BMGs was shown in Figure 6a. The introducing of symmetrical notches can result in stable plastic flows confined within the regions between notches. With high orders of stress concentrations, shear bands are easier initiated from the notch bottoms. However, only the stress fields between two symmetrical notches can effectively confine the propagation of shear bands. In the specimens with single notches or two asymmetric notches, significant improvement in nominal plasticity was not observed due to the rapid propagation of shear bands [100–103]. By tailoring the distribution of more notches is useful for achieving larger macroscopic plasticity of BMGs, especially in practical applications of BMGs [104]. Moreover, such a strategy may be particular useful for BMGs due to the unique atomic orders. The notched high-strength steel and ceramic specimens exhibited no obvious improvement or even decreased nominal strains (Figure 6b) [102].

**Figure 6.** Enhanced nominal strain in a notched BMG specimen under compression (**a**), and the comparison of the effect in different materials (**b**). Adapted from [102] with permission of Elsevier.

As compared with the plastic deformation behavior of BMGs under compressive loadings, it is challenging to obtain plasticity under tensile loadings due to catastrophic failures. With vanishing of crystalline lattices and boundaries, the rapid propagation of shear bands cannot be impeded, and leads to brittle failures. Nevertheless, it is fortunately to find that the introducing of complex stress fields can be employed to improve the plastic deformation behavior of the stress-concentrated regions. Qu et al. [105,106] have reported that the complex stress field created by double-side notches can prevent the unstable propagation of shear bands, and result in a stable plastic zone. Li et al. [107] have further examined the effect of notch sizes and shapes on the formation and propagation of shear bands in the stress-concentrated regions, based on FEM simulations, and identified the notch conditions which are beneficial for plastic deformation. By tailoring the plastic deformation in stress-concentrated regions, tunable axial elongations have also been obtained in some curved specimens, as shown in Figure 7 [108]. With both compressive and tensile stresses, the curved segment forms a shear band multiplication stage during tensile loading process [18,108]. Although the large plastic deformation was localized in the stress-concentrated regions, the straightening of the curved segments results in large axial elongations along the loading directions. Despite the brittleness under uniform stress distributions, large axial elongations can still be achieved in BMGs and BMG structures through the geometry design, laying a sound foundation for the practical applications of BMGs as structural materials.

Additionally, tensile ductility in BMGs has also been achieved by introducing complex stress states onto the surface of BMG specimens, such as imprinting method [109,110] and laser surface treatment [111,112]. The enhancement of macroscopic tensile ductility was attributed to the heterogeneities on the specimen surface, such as the mechanical properties and geometries, resulting from tailored complex stress fields [109,110]. The heterogeneities cause the formation of plastic zones with multiple shear bands. The multiplication and intersection of shear bands limit the rapid propagation of shear bands, and therefore enhance the macroscopic ductility. Gao et al. [111] have found that with complex stress fields introduced by laser surface treatment, the shear band propagation was impeded, leading to a relatively homogeneous deformation mode. Dong et al. [112] have further shown that the tensile ductility may be ascribed to the large scale flow, driven by the complex stress

fields introduced by laser engraving. Although the mechanisms on the enhancement of macroscopic ductility of BMGs under complex stress fields are still being debated, in fact, complex stress fields have been widely employed to investigate the deformation mechanisms of BMGs, i.e., under indentation tests. For example, indentation tests have been used to investigate the shear band formation and propagation mechanisms [113,114], and the transition of deformation modes from serrated flows to homogeneous deformation [115,116]. However, since the specimens do not fracture during indentation tests and the studies mainly focused on the corresponding deformation mechanisms, these studies on the indentation tests of BMGs are therefore not discussed in detail here.

**Figure 7.** Tunable large axial elongations in curved BMG specimens under tensile loadings, where the inset image showing the reduced sections of the curved specimens. Adapted from [108] with permission of Elsevier.

Up to date, complex stress fields have been used to characterize the shear band formation and propagation mechanisms, as well as the global deformation behavior, such as the plasticity and fracture. However, theories on how to control the propagation of shear bands and cracks, and subsequently control the macroscopic plastic deformation behavior of BMGs are still urgently needed. In the practical applications of BMGs as structural materials, most of the parts of BMG structures may deform under complex stress fields. The understanding of the macroscopic deformation behavior of BMGs under complex stress fields is also useful for designing BMG structures with better mechanical performance, for example, the development of BMG foams [20], BMG honeycombs [22] and cellular BMGs [117].

#### **5. Transition of Deformation Modes under Complex Stress Fields**

Attributed to the amorphous atomic structures, homogeneous plastic flow in BMGs is usually observed at high temperature [42]. However, at room temperature, the introducing of complex stress fields can also result in the transition of deformation modes from highly localized shear-banding to relatively homogeneous deformation. A pioneered work by Flores and Dauskardte [118] on the strain localization behavior of some notched BMG bars has suggested that the stress states may affect the deformation and failure behavior of BMGs. Want et al. [119] have then examined the plastic deformation behavior of a Zr64.13Cu15.75Ni10.12Al10 (at. %) BMG under multi-axial tensile stress states, introduced by circumferential deep notches. It was surprising to find that the notched BMG specimens demonstrate strain hardening behavior at room temperature, as can be seen in Figure 8. This unusual phenomenon was attributed to the diffusional relaxation driven by multiaxial stress states, where obviously shear-banding behavior was not observed. Further studies have shown that the change of the notch dimensions can cause BMGs to deform plastically through the nucleation and coalescence of

voids/cavies, where shear-banding behavior is suppressed [120]. Since the formation of shear bands occurs from very localized regions into sub-micron scales, the examination of homogeneous plastic deformation at sub-micron scales brings more understandings to the transition of deformation modes from localized shear-banding to homogeneous deformation.

**Figure 8.** (**a**) Strain hardening effect in a notched BMG at room temperature, and (**b**) shows the trace of mircohardness at the notched regions. Adapted from [119] with permission of American Physical Society.

Gu et al. have examined the shear-banding and fracture behavior of some notched Ni-P/Fe-P MGs of about 70 nm in diameter, as shown in Figure 9a [121]. Combined with experimental observations and molecular dynamics (MD) simulations, they showed that plastic deformation in these notched MG specimens initiated from the notch root by forming microscopic voids, and the coalescence of voids resulted in cavitation and final brittle failures. Narayan et al. [122] have studied the deformation mechanisms of some double-side-notched CuZr MGs with varying sharpness (Figure 9b). They have found that the specimens with sharper notches can delay the formation of shear bands, resulting from a higher degree of triaxiality in stress distributions. When the cavitation stress reaches a threshold value, plastic deformation in MGs can transit from shear-banding to microscopic voids coalescence, similar to the MD simulation results reported by Gu et al. [121]. On the other hand, Narayan et al. [122] have also examined the effect of notch depth on the plastic deformation mechanisms of MGs. It was shown that a deeper notch tends to facilitate homogeneously activated STZs and then suppress shear-banding. However, to date the experimental studies on the formation of STZs and shear bands in notched MG specimens are not sufficient enough to make a consensus conclusion on the formation and evolution mechanisms of STZs/shear bands. Since STZs only involve a small cluster of atoms, it is still challenging but necessary to investigate the initiation of STZs under complex stress fields directly through in-situ TEM observations, and how these STZs evolve to shear bands/cracks. MD simulations are therefore more feasible to be employed to investigate the formation of shear bands under complex stress fields.

Based on MD simulations, many previous studies focused on the effect of notch sizes and geometries on the plastic deformation behavior of nanoscaled MG specimens, aiming to shed more light into the transition of deformation modes. Sha et al. have shown that the design of notches with increased depth and sharpness can suppress shear-banding and result in more homogeneous defamation with presence of necking [61]. Similar transition of deformation modes can also be observed in nanopillars with tailored surface roughness [64]. Pan et al. [123] have investigated the physical origin of the homogenous deformation, and found that voids and cavitations were initiated at the notch root when the stress triaxiality exceeds a critical value. On the other hand, Dutta et al. [62] have also observed the transition of deformation modes from shear-banding (Figure 10a,b) to homogeneous deformation with necking (Figure 10c,d). With a relatively sharper notch (Figure 10c), plastic zones first initiated from the notch root due to stress concentrations, and then evolved to incipient shear bands (Figure 10d). However, due to the rapid expansion and coalesce of plastic zones, the formation of shear bands was finally suppressed, showing a necking effect (Figure 10e,f). It is reasonable to see the initiation of some incipient shear bands, where stress concentrators can also serve as sites for the nucleation of STZs, as discussed in Section 2. There may exist a competing process for the formation of localized shear bands and the coalescence of plastic zones, where the rapid coalesce of plastic zones can suppress the initiation of shear bands and result in different plastic deformation mechanisms. This phenomenon is also highly in line with the in-situ TEM observations on the plastic deformation behavior of some nanosized MG specimens [122]. The change of atomic packing in a shear band was observed in a Ni-based MG [124]. More recently, Cui et al. [63] have examined the structural evolution of CuZr MGs during the transition of deformation modes from shear-banding to homogeneously necking in notched MGs. They have shown that the Voronoi volume recovery can be dominant in the localized regions with triaxial stress state, differing from the unnotched specimens. Nevertheless, more effort is still needed for revealing the underlying physical origins of the homogeneous plastic deformation within the localized regions with complex stress fields, such as the movements of atoms, atomic structural evolutions as well as the formation and coalesce of microscopic voids.

**Figure 9.** (**a**) Deformation behavior of single-side-notched Ni-P/Fe-P MGs (Reprinted from [121] with permission of ACS Publications); (**b**) in-situ TEM observation of the deformation behavior of double-side-notched CuZr MG specimens with varying sharpness (adapted from [122] with permission of Elsevier).

**Figure 10.** MD simulation results of the Mises strain plots of unnotched (**a**) and notched (**b**–**f**) specimens, where (**c**–**f**) shows the evolution of plastic zones in a specimen with a sharper notch than (**b**). Adapted from [62] with permission of Elsevier.

#### **6. Deformation Behavior of HEAs under Complex Stress Fields**

HEAs are a new class of alloys, having at least five elements with equal or near equal atomic percentages, where the solvent and solute elements are not easily distinguished [125]. The mixing of multi component in the solution states results in very high entropy associated with unique properties. For example, they can have high strength comparable to BMGs [125]. More importantly, in conventional metallurgical methods, the increase of strength leads to the decrease of ductility, and it seems that the strength and ductility are usually mutually exclusive in conventional alloys [126]. However, some HEAs can exhibit both high strength and ductility [127]. The plastic deformation behavior of HEAs shares similar characteristics as compared with BMGs. For example, differing from conventional alloys, HEAs also have serrated plastic flows, resulting in challenges to accurately predict and control the plastic deformation behavior [40,128]. The serrated plastic flows of HEAs are also significantly dependent on the change of strain rates, temperature and sample dimensions [125,128,129]. For instance, Zou et al. [129] have examined the plastic deformation behavior of some HEAs with Ar+ ion beam-assisted deposition, where flow serrations were obviously observed in the pillar of 70 nm diameter, while the pillars with larger diameters tend to have smooth plastic flows (Figure 11).

**Figure 11.** The compressive stress-strain curves of NbMoTaW HEAs, where IBAD stands for the specimens with Ar+ ion beam-assisted deposition [129].

Recent research has shown that the residual strains may cause the instability of phases and result in the transition of phases [130]. This may further change the mechanical properties of HEAs, for example, transition induced plasticity [127]. Joseph et al. [131] have reported that some laser-fabricated Al0.3CoCrFeNi HEAs exhibited tension/compression asymmetric deformation behavior, without/with mechanical twining. Therefore, the change of applied stress fields could also affect the evolution of microstructures in HEAs, such as the phase transition and mechanical twining, and the resultant flow serrations and mechanical properties. The tailoring of complex stress fields can be used to improve the plasticity and tune the criticality of plastic flows in BMGs, while some research has also shown that the notches can reduce the nominal strains of crystalline high-strength steels [102]. Whether the introducing of complex stress fields is helpful for improving the plasticity and tuning the criticality of flow serrations in HEAs is still a mystery. Up to date, a study focused on the deformation behavior of HEAs under complex stress fields has yet been reported. Regarding that the evolution of microstructures of HEAs is significantly affected by applied stress states, it should be worthy of further investigations on the deformation behavior/mechanisms of HEAs under complex stress fields.

#### **7. Conclusions and Future Directions**

The change of applied stress fields can significantly affect many aspects of the plastic deformation in BMGs, such as the shear band initiation and propagation, evolution of flow serration, and macroscopic plastic performance. Varying complex stress fields have been employed to investigate the mechanisms on the burst of shear bands, guiding the propagation of shear bands/cracks, tunable plastic-flow dynamics, transition of deformation modes and better plastic performance. On one hand, complex stress fields play a significant role for elucidating the mechanisms on the plastic deformation in BMGs. For example, how STZs are initiated and when the formation of STZs is suppressed? On the other hand, complex stress fields can be tailored to achieve controllable plastic deformation behavior for practical engineering applications. A case in point is to design BMG structures/devices with enhanced mechanical properties. However, the understanding on the deformation behavior of BMGs under complex stress fields still faces challenges, and some unresolved issues are summarized and given below:


**Author Contributions:** S.C. prepared the manuscript. S.C. and J.W. collected the data. S.C., L.X. and Y.W. designed the scope of the paper. All authors discussed the conclusions and reviewed the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China grant number 51801049. **Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Review* **Liquid Phase Separation in High-Entropy Alloys—A Review**

#### **Nicholas Derimow and Reza Abbaschian \***

Department of Materials Science and Engineering, University of California, Riverside, CA 92521, USA; nderimow@engr.ucr.edu

**\*** Correspondence: rabba@engr.ucr.edu; Tel.: +1-951-827-6374

Received: 24 October 2018; Accepted: 16 November 2018; Published: 20 November 2018

**Abstract:** It has been 14 years since the discovery of the high-entropy alloys (HEAs), an idea of alloying which has reinvigorated materials scientists to explore unconventional alloy compositions and multicomponent alloy systems. Many authors have referred to these alloys as multi-principal element alloys (MPEAs) or complex concentrated alloys (CCAs) in order to place less restrictions on what constitutes an HEA. Regardless of classification, the research is rooted in the exploration of structure-properties and processing relations in these multicomponent alloys with the aim to surpass the physical properties of conventional materials. More recent studies show that some of these alloys undergo liquid phase separation, a phenomenon largely dictated by low entropy of mixing and positive mixing enthalpy. Studies posit that positive mixing enthalpy of the binary and ternary components contribute substantially to the formation of liquid miscibility gaps. The objective of this review is to bring forth and summarize the findings of the experiments which detail liquid phase separation (LPS) in HEAs, MPEAs, and CCAs and to draw parallels between HEAs and the conventional alloy systems which undergo liquid-liquid separation. Positive mixing enthalpy if not compensated by the entropy of mixing will lead to liquid phase separation. It appears that Co, Ni, and Ti promote miscibility in HEAs/CCAs/MPEAs while Cr, V, and Nb will raise the miscibility gap temperature and increase LPS. Moreover, addition of appropriate amounts of Ni to CoCrCu eliminates immiscibility, such as in cases of dendritically solidifying CoCrCuNi, CoCrCuFeNi, and CoCrCuMnNi.

**Keywords:** high-entropy alloys; liquid phase separation; immiscible alloys; HEAs; multicomponent alloys; miscibility gaps; multi-principal element alloys; MPEAs; complex concentrated alloys; CCAs

#### **1. Introduction**

#### *1.1. Liquid Phase Separation*

Liquid phase separation (LPS), a widely-observed phenomenon in metals, is related directly to the Gibbs free energy of the system, and the most prevailing cases are often two distinct immiscible liquids of varying compositions. Although there is often some degree of solubility between the alloying elements in a metallic system exhibiting LPS, each liquid will have its own equilibrium vapor pressure, such that the vapor pressures of both phases are the same, with a positive deviation from Raoult's law. When positive deviations from Raoult's law are large, phase segregation tends to occur.

The occurrence of liquid phase separation in an alloy can lead to heterogeneous microstructures, which may or may not be desirable depending on the intended application. For example, an alloy exhibiting liquid phase separation would not be suited for use as a structural material due to the heterogeneity of the microstructure; however, it may have potential use as a self-lubricating bearing material, such as the case with Cu-Pb. There have been several comprehensive reviews of

immiscible metal systems of common alloys about the phenomenon [1–5]. Therefore, the scope of this review will focus particularly on the liquid phase separation in the high-entropy alloy (HEA), complex concentrated alloy (CCA), and multi-principal element alloy (MPEA) systems.

#### *1.2. Thermodynamics of Liquid Phase Separation*

Factors such as positive deviations from Raoult's law, positive heat of mixing, and atomic size mismatch in some cases do not overcome the entropy term in the overall Gibbs free energy and cause overall immiscibility in the liquid, as is the case of miscibility between Au-Bi [5]. B. Mott in the late 1950s put together a review of the immiscible liquid metal systems, as well as the corresponding thermodynamic data for each material at the time. The immiscible alloys Mott compiled in the study contained many of the known immiscible binary monotectic alloys of the time [1]. Nearly ten years later, Mott compiled another review detailing the thermodynamics of these metal systems, as well as provided models for predictions of immiscibility in metals [2].

In Table 1, we provide a non-exhaustive table of binary alloys with miscibility gaps in the liquid state. The table expands the tables from Mott's reviews [1,2] by adding immiscible alloys from binary phase diagrams provided by the Centre for Research in Computational Thermochemistry using the FactSage thermochemical software databases [6].


**Table 1.** Binary systems that contain a stable miscibility gap in the liquid state.

The molar Gibbs free energy of a system of stable unmixed liquids is represented additively via the atomic fraction of the free energies of the constituent liquids,

$$G\_{\mathbf{A}+\mathbf{B}+\mathbf{C}\dots}^{\mathbf{L}} = \sum\_{i=\mathbf{A}, \mathbf{B}, \mathbf{C}\dots} x\_i G\_i^{\mathbf{L}} \tag{1}$$

where *xi* <sup>=</sup> A,B,C... are the molar fractions of elements A, B, C, etc. The molar Gibbs free energy of mixing is classically defined as,

$$
\Delta G\_{\rm mix} = \Delta H\_{\rm mix} - T\Delta S\_{\rm mix} \tag{2}
$$

where the entropy of mixing is given as,

$$
\Delta S\_{\text{mix}} = -R \sum\_{i} \mathbf{x}\_{i} \ln \mathbf{x}\_{i} \tag{3}
$$

and the enthalpy of mixing is:

$$
\Delta H\_{\text{mix}} = \sum\_{\substack{i=1, i \neq j}}^n \Delta H\_{x\_i x\_j}^{\text{mix}} \tag{4}
$$

where Δ*H*mix *xi*,*xj* is the interatomic interaction between concentrations of "*i* " and "*j*" elements in the system. Immiscible alloys typically have a positive value of Δ*H*mix, which implies a preference of nearest neighbors of similar atoms as opposed to compound formation with different atoms. Many of the immiscible binary systems can be categorized by their liquid state miscibility gaps and positive enthalpy of mixing, Δ*H*mix, of which extensive thermodynamic treatments are presented in [3,5].

The region of a phase diagram where there is non-mixing of the constituents is defined as a miscibility gap. The liquid miscibility gap in many of the monotectic binary systems assumes a dome-like shape; the shape and location of which may shift with the addition of more alloying elements. For example, one of the most well-studied ternary systems with a stable liquid miscibility gap is the Co-Cu-Fe system [7–17], while with equiatomic additions of Cr and Ni, the CoCrCuFeNi high-entropy alloy solidifies dendritically from a single-phase liquid, as observed by Yeh et al. in 2004 [18].

A generalized equilibrium monotectic phase diagram is presented in Figure 1, where the miscibility gap in the liquid state is present as a dome with label L1 + L2. The size and width of the immiscibility gap varies from system to system; however, the concept is the same. That is, cooling the alloy system from a liquid state in the concentrations that fall within the miscibility gap will lead to the liquid decomposing into two compositionally different liquids, the temperature of which is known as the critical temperature (labeled T*c* in Figure 1).

As the temperature decreases to T1 in Figure 1, the entropy term *T*Δ*S* is smaller than the enthalpy of mixing Δ*H*mix in the free energy of the system (Figure 2); therefore, the free energy of the liquid GL with respect to concentration of the B element in A will also assume a dome shape, presented in Figure 1. If the temperature T1 is held, the equilibrium phases will be L1, L1 + L2, or L2 dependent on composition XB. Cooling the system through T2 until the monotectic temperature, T3, the monotectic reaction will take place, and we will start to see *α* precipitate out of the liquid as L1 is no longer stable until we reach T4, where the remaining equilibrium phases are the solidified *α* and liquid L2. Cooling through the eutectic temperature at T5 to reach T6, we are ultimately left with (*α* + *β*) solid phases.

Due to the lack of experimental data for mixing enthalpies of many binary alloys, a model for generating approximate mixing enthalpies was first developed by Miedema et al. in 1973 [19], which uses the electron density at the Wigner–Seitz cell boundary and the chemical potential of electronic charge of pure metals as input and can be written as Δ*H*mix = ∑*<sup>n</sup> <sup>i</sup>*=1,*i*=*<sup>j</sup>* <sup>Δ</sup>*H*mix *ci*,*cj* . This model was used by Takeuchi et al. in 2005 for the classification of bulk metallic glasses by atomic size difference and heat of mixing [20] and later revisited by Takeuchi in 2010 [21] for mixing enthalpies of

binary alloys, which includes an additional model for sub-regular solutions [22]. The Δ*H*mix of the binary alloys from [21] serve as a starting point for many of the recent calculations of Δ*H*mix for HEAs, MPEAs, and CCAs. Using the calculated binary mixing enthalpies (Δ*H*mix) from Takeuchi et al. [21], Δ*H*mix, much of the values used for determining the mixing enthalpies for HEAs/CCAs/MPEAs were calculated using Equation (4) where Δ*H*mix *xi*,*xj* = 4Ω0*ij xixj* for the *i* th and *j* th elements at A0.50B0.50 concentrations from the tables in [21]. The values for *ci* and *cj* are the normalized atomic concentrations in the multicomponent alloy.

**Figure 1.** Generalized equilibrium monotectic binary phase diagram.

**Figure 2.** Gibbs free energy corresponding to the monotectic phase diagram.

#### *1.3. Metastable Liquid Phase Separation*

Unlike the stable liquid state immiscibility observed in the monotectic binary alloys, there are certain cases where a completely miscible liquid alloy can de-mix in the presence of impurities or when supercooled below the freezing temperature of the alloy, as demonstrated for Co-Cu and Cu-Fe by Nakagawa in 1958 [23]. Since then, there has been an enormous amount of LPS studies on metastable Co-Cu [8,15,16,24–34] and Cu-Fe [15,35–43], as well as the stable LPS that occurs in the combination of all three elements in Co-Cu-Fe [7–17]. Metastable liquid phase separation is defined as the liquid phase separation that occurs when undercooling an alloy such that it enters a miscibility gap that would not have been observed if solidified via conventional methods, presented in the phase diagram from [25] in Figure 3. These studies have shown that when undercooling past freezing, the single-phase alloy liquid will then split into two liquids (L1 + L2), specifically in these cases, into Cu-rich and Cu-lean liquids, which solidify often as spherical globules trapped in the frozen regions of the other liquid (the microstructure of such will be discussed later). This metastable LPS implies that there exists a dome shape similar to the monotectic alloys beneath the liquidus curves in their respective equilibrium phase diagrams.

**Figure 3.** The Co-Cu phase diagram with the dashed line indicating the metastable liquid miscibility gap beneath the liquidus curves.

#### *1.4. High-Entropy Alloys*

The discovery of the high-entropy alloys (HEAs) [18,44–48] has inspired an enormous amount of research into multicomponent alloy design. Since their inception, there have been numerous reviews [49–60] and several books [61–66] that summarize the state-of-the-art for the materials community. These reviews compile and assess the microstructural developments, mechanical properties, crystallography, and thermodynamics of the high-entropy, complex concentrated, and multi-principal element alloy systems.

Many of the ternary alloys that exhibit liquid phase separation contain Cu, as it has a low affinity for mixing with other elements; however, there are a number of other non-Cu-containing ternary alloys with liquid miscibility gaps as well. Table 2 is a non-exhaustive list of the studied ternary alloys with liquid phase miscibility gaps, many of which contain Cu. As was the case with many binary alloys, the Δ*H*mix of these systems are typically positive. One can think of the HEAs/CCAs/MPEAs as the addition of alloying elements to preexisting ternary alloys, some of which may actually contain a stable miscibility gap in the liquid, which is the case with CoCrCu [67] and many of the HEAs that contain Co, Cr, and Cu in equal parts with respect to the other alloying elements in the HEA. It would appear that the increase in the entropy of mixing Δ*S*mix with additional alloying elements stabilizes the solution; however, this strongly depends on the enthalpy of mixing Δ*H*mix, as well as other factors that determine miscibility [58].


**Table 2.** Ternary alloy systems that contain a stable liquid miscibility gap.

One of the first HEAs contained equal parts Co, Cr, Fe, Mn, and Ni, often written alphabetically as CoCrFeMnNi and referred to as the "Cantor alloy" after the alloy's inventor Brian Cantor [44]. Subsequent studies of the Cantor alloy involved substitution of various elements in place of the original five equiatomic elements in the alloy. One of the most popular substitutions is often Cu for Mn [18], as well as the addition of Al to create the widely-studied AlCoCrCuFeNi HEAs [85–87]. Since their inception, the majority of the high entropy alloys that are synthesized consist mostly of the 3D transition metals with other elements substituted intothe well-studied systems in the search for new stable phases, mechanical properties, and reproducible microstructures [58]. There have been several research papers that prescribe methodologies for computationally screening HEAs for those that are likely to create single-phase solid solutions [88–96]. A useful criteria for HEA/CCA/MPEA research is to know whether or not the combination of elements in a proposed HEA system will even mix in the liquid phase, as LPS traditionally is an unwanted phenomena when designing new materials with the goal of enhancing mechanical properties. If a single-phase liquid can be obtained from the proposed combination of elements, then the solidification microstructures will result in more uniformity.

There is much debate as to what exactly constitutes a high-entropy alloy (HEA), complex concentrated alloy (CCA), or multi-principal element alloy (MPEA). In essence, the core idea is very similar behind each definition: multicomponent, "baseless" alloys greater than three elements designed with the goal of surpassing the mechanical properties of traditional alloys. Whether or not which composition will solidify into a single, duplex, or multiple phases is not what this review is concerned with, but rather the observance of liquid phase separation in these multicomponent alloy systems.

Approximately 85% of HEAs in the literature to date are made up of predominantly 3D transition metals, many of which contain Al [58]. As it is impossible to visualize the phase diagram space of multicomponent alloy systems that contain 5+ elements, it can be difficult to know whether these alloy systems will phase separate in the liquid. Many of the well-studied HEA systems, such as CoCrCuFeNi for example [45], contain equiatomic CoCrCu, which has a very large liquid miscibility gap [67,70].

#### **2. Solidification Microstructures**

#### *2.1. Dendritic Microstructure*

Alloy solidification morphology and as-cast microstructure can have many forms depending on the solidification process. The most common microstructures of a solidified alloy can vary from plane front solidification, dendritic morphology, and eutectic microstructures, among others. Many of the solidification microstructures present in HEA literature consist of dendritic growth, which is typically indicative of crystal growth from a liquid with an imposed thermal gradient, as is the case with most arc-melting processes. Dendrites are solid tree-like, branching cellular structures that grow from a liquid phase. The conglomeration of atoms during solidification typically forms a nucleus of spherical shape, which then becomes unstable due to perturbations. The solid shape then begins to express the preferred growth directions of the underlying crystal and consumes atoms from the overall liquid to form a stable solidifying phase [97]. The liquid that is leftover after dendritic solidification is referred to as the interdendritic liquid, which solidifies last, and is referred to as the interdendritic region or interdendrite. The preferred dendritic growth directions for most cubic systems (FCC/BCC) are in the <100> directions, which leads the secondary dendrite arms to grow perpendicular from the primary arm. This is often an easy way to differentiate between the cubic and noncubic crystal structure of the dendritic phase.

The typical dendritic microstructure of an electromagnetically-levitated and solidified alloy is presented in Figure 4. The dendritic morphology usually indicates that the microstructure evolved from a single-phase liquid if the dendritic morphology is uniform; however, if the alloy has a small volume fraction of LPS, the remaining L2 globules may be pushed to the edge of the sample by the growing dendrites. There have been cases where liquid phase separation has occurred in the interdendritic liquid after the growth of primary dendrites [98,99]; however, the general morphology of dendritic microstructures indicates that there was no large-scale liquid-liquid immiscibility between the alloying elements.

**Figure 4.** Backscattered electron image of an as-cast AlMoNi alloy with a dendritic microstructure.

Many HEAs/CCAs/MPEAs solidify with a duplex microstructure, where the dendritic and interdendritic regions have large compositional and crystallographic differences [53,58–60,64]. These have been shown to have interesting mechanical properties; however, they have yet to surpass the mechanical properties of commercial alloys.

#### *2.2. Microstructures Resulting from Liquid Phase Separation*

Alloys that undergo either stable or metastable LPS also have very distinct microstructures that can vary based on the solidification process. Slow cooling rates paired with a static environment can lead to the liquids separating. If the system is a little more dynamic, such as in the case of casting, the process can lead to trapping of the primary liquids in one another, referred to as emulsion (Figure 5). The separated liquids tend to be trapped as spherical globules inside the other liquid, and solidify as such, as is the case of the equiatomic CoCrCu alloy presented in Figure 6. As these liquids can be slightly different in composition than the primary liquid, they are referred to as secondary liquids. Based on morphology alone, one can distinguish the first phase to solidify from the interface between the the two liquids, as the higher melting point liquid will solidify at a higher temperature and will most of the time solidify with protrusions into the other liquid, as is the case with the solidification of CoCrCu presented in Figure 7. The backscattered electron images in Figures 6 and 7 display emulsion of the lighter (Cu-rich) and darker (CoCr-rich) liquids, as well as small protrusions coming from the CoCr-rich secondary liquid in Figure 7.

There have been several efforts to create uniform microstructures of immiscible liquid melts, such that the LPS is evenly distributed. These techniques include free directional and directional solidification [100], rheomixing [101], microgravity experiments [28,102,103], electromagnetic levitation processing [40], and rapid solidification [104]. Much of this work is aimed at the tailoring of the immiscible liquid droplets such that the microstructure has a uniform spread of the immiscible phase [3]. There have also been experiments aimed at the suppression of LPS with additional alloying elements [105].

**Figure 5.** SEM image of liquid phase separation and emulsion of two immiscible liquids L1 and L2 in an undercooled CoCuFe alloy. Image presented with permission from the authors in [16].

**Figure 6.** Backscattered electron image displaying emulsion of CoCr-rich (darker) and Cu-rich (lighter) liquids in an as-cast alloy of CoCrCu.

**Figure 7.** Backscattered electron image displaying emulsion and protrusions of the CoCr-rich (darker) phase into the Cu-rich (lighter) phase in an as-cast alloy of CoCrCu.

Traditionally, observing LPS in metallic systems was done via post-mortem analysis via metallography and microscopy, as metals are not transparent to light and have a very small transparency for X-rays. Recently, through the use of neutron transmission imaging techniques, the direct observation of liquid phase separation in metals was made possible via neutron radiographs taken during heating and cooling of immiscible CoCrCu alloys [72]. These experiments show for the first time an in situ observation of macroscopic LPS in metals and can be applied to any metallic system such that the neutron transmission through the each phase can provide enough contrast between

them. Figure 8 displays neutron radiographs of two stacked CoCrCu arc-melted buttons in a small Al2O3 crucible with an inner diameter of 8 mm. Prior to the in situ testing, the arc-melted CoCrCu buttons underwent stable LPS and solidified with very heterogeneous Cu-rich and Cu-depleted regions, as presented in Figure 9. The resulting as-cast button consists of a non-uniform mix of the two solidified Cu-rich and Cu-depleted phases and can be seen as the lighter (Cu-rich) and darker (Cu-depleted) regions in Figure 8a. During melting (Figure 8b–d), the two liquid phases separate and stack according to density (Cu being the lighter contrast, more dense liquid phase). Note that the solid Cu-rich and CoCr-rich phases were already separated in the arc-melted buttons, but at a finer scale. The stable liquid phases then agglomerated and separated at the macro-scale. A full sequence of images in the form of a movie can be found in [72].

**Figure 8.** Melting and liquid phase separation of stacked CoCrCu samples. (**a**) During initial heating, the two as-cast buttons are intact; (**b**) the Cu-rich phase melts first between 1075 and 1100 °C and (**c**) pools at the bottom of the crucible; (**d**) the Cu-lean phase fully melts upon heating to 1500 °C and stacks based on density due to the influence of gravity. Images displayed with permission from the authors in [72].

**Figure 9.** Neutron radiograph of two stacked arc-melted CoCrCu buttons in an alumina crucible with brighter regions corresponding to Cu-rich phases and darker regions corresponding to CoCr-rich phases. Image displayed with permission from the authors in [72].

#### **3. High-Entropy Alloys Exhibiting Liquid Phase Separation**

#### *3.1. HEAs Containing Cu*

The first occurrence of LPS in HEAs was observed by Hsu et al. in 2007 with a study of the alloying behavior of AlCoCrCuNi-based HEAs with additions of Fe, Ag, and Au [106]. The addition of Ag to the AlCoCrCuNi HEA to create AgAlCoCrCuNi was found to phase separate in the liquid, which resulted in the solidification microstructure consisting of Cu-rich globules embedded in Cu-depleted phases, contrary to the typical dendritic solidification microstructures observed for AlCoCrCuNi. Hsu suggested that in order to achieve effective mixing in the liquid, the Δ*H*mix for atom pairs should not exceed 10 kJ/mol and that "mutual interaction between elements, based on their mixing enthalpies, should be taken into account when designing high-entropy alloys" [106].

These alloys were then revisited by Munitz et al. in 2013 where AgAlCoCrCuNi and AgAlCoCrCuFeNi were synthesized to study the melt separation behavior [107]. It was observed that the Cu-rich immiscible liquid tended to flow to the bottom of the buttons during arc-melting, as well as residual Cu-rich liquid being trapped in the interdendritic region of the Cu-depleted dendritic phase. Undercooling experiments were also carried out for a similar alloy of Al1.8CoCrCu3.5FeNi; however, no metastable liquid miscibility gap was found at the undercoolings obtained in the study (∼150 K) [107].

A similar alloy composition of AlCoCrCuFeNiSi0.5 doped with Y2O3 was synthesized via laser cladding with the intention to form a core-shell structure in HEAs, inspired by the LPS observed in HEAs and binary monotectics. The undoped AlCoCrCuFeNiSi0.5 did not undergo LPS, while the addition of 1 wt.% nanosized Y2O3 caused the liquids to separate into egg-like globules of Cu-rich liquids inside the Cu-depleted liquid [108].

Recent studies into Co-free Al2.2CrCuFeNi2 revealed what the authors referred to as anomalous "sunflower-like" solidification microstructures [109], where it was suggested that LPS occurs in the no longer stable-depleted interdendritic liquid, occurring due to changes in the composition. Munitz et al. suggested that the liquid phase separation is due to constitutional changes and not temperature changes, where the authors referred to this phenomena as "constitutional LPS" (CLPS) [99]. For the Al2.2CrCuFeNi2 alloy, constitutional LPS occurred in the interdendritic liquid, where the interdendritic liquid decomposed into a CrFe-rich L1 and a Cu-rich L2. As the dendritic skeleton was already formed, the heavier Cu-rich liquid accumulated in the interdendritic region and the cast bottom, while the CrFe-rich spheres underwent solidification emulsified in the Cu-rich liquid.

A large study of several HEAs by Munitz et al. in 2017 was undertaken to explore the effects of Al, Co, Cr, Ni, Ti, and V on the miscibility gap temperature of several HEA systems. It was shown that Al, Co, Ni, and Ti lowered the miscibility gap temperature, while Cr, V, and Nb raised the miscibility gap temperature and increased LPS in these systems, the alloys of which are found in Table 3. Many of the HEAs studied by Munitz et al. contained equiatomic CoCrCu, which was experimentally determined to have a large liquid miscibility gap [67]. It is peculiar that systems such as CoCrCuFeNi will solidify dendritically [18], while similar alloys of CoCrCu [67,71,72] and CoCuFe [7–12,14–17,73] have been shown to have large liquid miscibility gaps. A recent study by Derimow et al. investigated the solidification microstructures of equiatomic CoCrCu with added Fe, Mn, Ni, V, FeMn, FeNi, FeV, MnNi, MnV, and NiV to the composition. It was found that only three of the alloys solidified dendritically (CoCrCuNi, CoCrCuFeNi, and CoCrCuMnNi), while the remaining combinations underwent stable LPS [71]. Derimow et al. also suggested that the positive mixing enthalpy of each of the systems was responsible for the LPS and presented a tree diagram for approximating the likelihood of which elements will cluster together in the melt, presented in Figure 10.

The tree diagram in Figure 10 indicates that AB atoms are more likely to cluster than ABC, AC, or BC, thereby rejecting the C element. This can be seen in case of the ternary CoCrCu, where the LPS consists of CoCr-rich and Cu-rich liquids [67,71] and solidifies with the microstructure shown in Figure 6.

**Figure 10.** Tree diagram representing the probability of clustering based on the mixing enthalpies of binary combinations of elements ABC.

Table 3 compiles the multicomponent alloys found in the literature that have been shown to undergo LPS. Each of these systems have a solidification morphology similar to the microstructure presented in Figure 6. From this table, the common component in all but one of these systems is equiatomic Cu. This is in part due to the positive mixing enthalpy Cu has with many of the alloying elements in the system. It should be noted that although Ni appears to promote miscibility in Cu-containing HEAs, this may be ineffective if the repulsion between Cu and the majority of the alloying elements is too great.


**Table 3.** Multicomponent alloy systems with reported liquid miscibility gaps.

#### *3.2. CoCrCuFeNi*

One of the seminal HEAs synthesized by Yeh et al. was the CoCrCuFeNi HEA [18]. Along with the Cantor alloy (CoCrFeMnNi) [44], these alloys served as many starting points for the addition and subtraction of alloying elements. The first study to note the liquid phase separation in HEAs in a similar composition of AgAlCoCrCuNi attributed the presence of LPS to the positive mixing enthalpies between Ag and the rest of the alloying elements [106].

An induction melting study by Wu et al. involving CoCrCuFeNi demonstrated the first occurrences of LPS in this alloy [116]. In Wu's study, several combinations of the CoCrCuFeNi alloy with varied Fe and Ni were studied to investigate the effects on the microstructure and crystallography of the system. When Fe and Ni were varied to create CoCrCuFe0.5Ni and CoCrCuFeNi0.5, spherical Cu-rich separations were observed by post-mortem analysis of the solidified samples [116]. The authors rationalized that the two alloy compositions with positive mixing enthalpies were too

great to be overcome by the entropy of the system, which would thereby lower the overall Gibbs energy of the system in the molten state.

Previous studies of supercooling and rapid solidification via electromagnetic levitation melting by Elder et al. characterized large undercoolings of 150 K for Cu-Fe and 75 K for Co-Cu that produce the LPS microstructures for these metastable alloys [120]. Elder et al. listed several techniques to achieve high undercoolings such as melt emulsification, melting in molten slag or fused silica (glass fluxing), free fall in a drop tube, and electromagnetic levitation techniques to achieve undercooled temperatures [120]. Using the molten fused silica technique, rapid solidification studies of CoCrCuFeNi by Liu et al. were carried out to investigate rapid solidification effects on the microstructure and phase stability of CoCrCuFe*x*Ni HEAs (where x = 1, 1.5, and 2.0) [115]. It was found that LPS occurred in all three compositions below the critical undercooling temperature, <sup>Δ</sup>*T*crit, where <sup>Δ</sup>*Tx*=1.0 crit = 160 K, <sup>Δ</sup>*Tx*=1.5 crit = 190 K, and <sup>Δ</sup>*Tx*=2.0 crit = 293 K. When <sup>Δ</sup>*<sup>T</sup>* <sup>&</sup>gt; <sup>Δ</sup>*T*crit, the microstructure is consistent with that of LPS, such that there were Cu-rich spheres present throughout the material. Due to the dendritic solidification behavior of CoCrCuFeNi through regular solidification routes [18], this alloy is classified as having a metastable liquid phase miscibility gap that is present when undercooled past Δ*T*crit [115]. The metastable liquid miscibility gap was also confirmed by Wang et al. for CoCrCuFeNi, where the authors also achieved an exceptionally high degree of undercooling of 381 K (0.23T*m*) for the alloy [117]. The Δ*T*crit for equiatomic CoCrCuFeNi was also studied by Guo et al., and it was found that metastable LPS occurs when Δ*T* > Δ*T*crit = 100 K [118]. The authors also show that the yield strength and elongation of equiatomic CoCrCuFeNi significantly decrease when the alloy undergoes liquid phase separation due to the non-uniformity of the resultant microstructure [118].

Recently, Wang et al. showed that with the addition of 3 at.% Sn to CoCrCuFeNi, the alloy undergoes the same characteristic liquid phase separation when undercooled past Δ*T*crit = 100 K [121]. The study showed that the LPS produced an increase in hardness of the Cu-depleted phases due to the separation of Cu and Sn in the liquid [121].

Previous studies on similar alloys have shown that Mo improves the strength in the AlCoCrFeNi and AlCoCrCuFeNi HEAs [122–124]. However, it has been shown that with varied Cu concentrations in CoCrCu*x*FeMoNi where x ≥ 0.5, LPS occurs in a similar fashion to the other Cu-containing HEAs [114]. The addition of Mo to the CoCrCuFeNi alloy was investigated by Wu et al. to elucidate the solidification process in these alloys, as there had been a lack of thorough studies on the solidification microstructures of these alloys [114]. Due to the Cu-rich sphere emulsion in Cu-depleted phases, likely due to the positive mixing enthalpy between Cu and the remaining alloying elements, Wu et al. suggested that the Δ*H*mix criteria for the prediction of single phase formation be amended to include the possibility of liquid phase separation in the liquid when Δ*H*mix > 0. In order to further study the Mo-containing HEAs, Peng et al. synthesized a Co-depleted CrCu*x*FeMo*y*Ni HEA to further elucidate the effects of the large positive Δ*H*mix between Cu and Mo on the solidification process and microstructure. It was found that Cu-rich and Cu-depleted LPS occurs in the CrCu*x*FeMo*y*Ni when x and y = 0.5 and 1, attributed to Δ*H*mix > 0 for these alloy combinations [119].

#### **4. Closing**

The field of high-entropy alloys, complex concentrated alloys, and multi-principal element alloys continues to grow, with new studies producing valuable insights for the materials community with the overarching goal of creating new alloys that exceed the properties of conventional materials. This relatively new class of material is not much different from the conventional alloys, being that they are still subject to the same thermodynamic rules that are imposed on them. The main caveats are that with the increase of alloying elements, orthogonal element phase diagram visualization becomes impossible; therefore, creative ideas are warranted to help understand the nature of the solidification of these alloys. Positive mixing enthalpy, if not compensated by the entropy of mixing, will cause liquid phase separation. It appears that Co, Ni, and Ti promote miscibility in multicomponent alloys, while Cr, V, and Nb will raise the miscibility gap temperature and increase LPS. Moreover, for equiatomic

CoCrCu, which has a large liquid miscibility gap, the addition of appropriate amounts of Ni eliminates immiscibility. The indication of such an example is the CoCrCuFeNi alloy, which will solidify dendritically, while similar alloys of CoCrCu and CoCuFe show strong immiscibility. Moreover, when Fe, Mn, Ni, V, FeMn, FeNi, FeV, MnNi, MnV, and NiV are added to to equiatomic CoCrCu, only three of the alloys solidify dendritically (CoCrCuNi, CoCrCuFeNi, and CoCrCuMnNi), while the remaining combinations undergo stable LPS. In the case of CoCrCuNiV, it appears that the addition of Ni in equiatomic amounts was not enough to overcome the positive mixing enthalpy interaction between Cu and V, as the CoCrCuV alloy also exhibits stable LPS. From the table of listed multicomponent alloys that undergo LPS, Cu is found in all but one of these combinations, which indicates that Cu containing HEAs may contain a metastable liquid miscibility gap such as the case with CoCrCuFeNi.

**Author Contributions:** Conceptualization, N.D. and R.A.; investigation, N.D. and R.A.; resources, R.A.; writing—original draft preparation, N.D. and R.A.; writing—review and editing, N.D. and R.A.; supervision, R.A.; funding acquisition, R.A.

**Funding:** This research was funded by the 'WINSTON CHUNG Professor in Sustainability' endowment at UC Riverside.

**Acknowledgments:** The authors would like to acknowledge Dr. Abraham Munitz for his lifetime of work on liquid phase separation.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

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

**Antonio Fernández-Caballero 1,2, Mark Fedorov 3, Jan S. Wróbel 3, Paul M. Mummery <sup>1</sup> and Duc Nguyen-Manh 2,\***


Received: 7 November 2018; Accepted: 11 January 2019; Published: 15 January 2019

**Abstract:** 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).

**Keywords:** multicomponent; ab initio; configuration entropy; matrix formulation; cluster expansion; cluster variation method; monte carlo; thermodynamic integration

#### **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 *Sconfig*/*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–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 *Sconfig*/*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 Al1.3CoCrCuFeNi 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–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*<sup>28</sup> 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 60Co 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 Cr25Fe25Mn25Ni25 and also the composition Cr18Fe27Mn27*Ni*<sup>18</sup> 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 *T* = 0 K of the alloy can be found from enthalpies of mixing, Δ*HMixing DFT* [*σ*], of derivative alloy structures [16,25] denoted by varying length arrays of spin-like variables*σ<sup>i</sup>* 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*[*σ*]*<sup>p</sup>* = 0. For each member structure, the reference energy in the calculation of the enthalpy of mixing is the pure element total energy, *Etotal*[*p*], with the same crystallographic symmetry. The alloy enthalpy of mixing is defined as follows:

$$\Delta H\_{DFT}^{\text{Mixing}}[\vec{\sigma} \equiv \{\sigma\_1, \sigma\_2, \dots \}] = E\_{\text{total}}[\vec{\sigma} \equiv \{\sigma\_1, \sigma\_2, \dots \}] - \sum\_{p=1,\dots,K} \mathbf{x}\_p[\vec{\sigma} \equiv \{\sigma\_1, \sigma\_2, \dots \}] E\_{\text{total}}[p]. \tag{1}$$

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

$$
\Delta H\_{CE}^{Mixing}[\vec{\sigma} \equiv \{\sigma\_1, \sigma\_2, \dots \}] = \sum\_{\omega, n\_\star(s)} m\_{\omega, n}^{(s)} l\_{\omega, n}^{(s)} \langle \Gamma\_{\omega', \mu'}^{(s')}[\vec{\sigma} \equiv \{\sigma\_1, \sigma\_2, \dots \}] \rangle\_{\omega, n, (s)}.\tag{2}$$

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 *ji* taking values 0, ..., *K* − 1. For each reference cluster, there is an associated effective cluster interaction *J* (*s*) *<sup>ω</sup>*,*n*, and multiplicity per lattice site, *<sup>m</sup>*(*s*) *<sup>ω</sup>*,*n*. The term Γ(*<sup>s</sup>* ) *<sup>ω</sup>*,*<sup>n</sup>*[*σ*] *<sup>ω</sup>*,*n*,(*s*) used in the definition of enthalpy of mixing in Equation (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*(*s*) *<sup>ω</sup>*,*<sup>n</sup> J* (*s*) *ω*,*n* Γ(*<sup>s</sup>* ) *<sup>ω</sup>*,*<sup>n</sup>*[*σ*] *<sup>ω</sup>*,*n*,(*s*) 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:

$$\begin{split} \langle \Omega\_{\omega',n'}^{(s')}[\vec{\sigma}] \rangle\_{\omega,\eta,(s)} &= \langle \gamma\_{j\_1,\mathbb{K}}[\sigma\_1] \gamma\_{j\_2,\mathbb{K}}[\sigma\_2] \cdots \gamma\_{j\_{\omega,\mathbb{K}}}[\sigma\_{\omega}] \rangle \\ &= \frac{1}{\Omega[\omega,n]} \sum\_{u=1}^{\Omega[\omega,n]} \gamma\_{j\_1,\mathbb{K}}[\sigma\_{\overleftarrow{(\omega,n)}\_{u\_1}}] \gamma\_{j\_2,\mathbb{K}}[\sigma\_{\overleftarrow{(\omega,n)}\_{u\_2}}] \cdots \gamma\_{j\_{\omega,\mathbb{K}}}[\sigma\_{\overleftarrow{(\omega,n)}\_{u\_{\omega}}}], \end{split} \tag{3}$$

where the point function is written as:

$$\gamma\_{j,K}[\sigma\_i] = \begin{cases} 1, & \text{if } j = 0 \\ -\cos\left(2\pi \lceil \frac{j}{2} \rceil \frac{\sigma\_i}{K} \right), & \text{if } j > 0 \text{ and odd,} \\ -\sin\left(2\pi \lceil \frac{j}{2} \rceil \frac{\sigma\_i}{K} \right), & \text{if } j > 0 \text{ and even,} \end{cases} \tag{4}$$

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, { ←−−→ (*ω*, *<sup>n</sup>*)*u*}*u*=1,2,··· ,Ω[*ω*,*n*] where each cluster labeled by *<sup>u</sup>* contains *<sup>ω</sup>* sites, each site denoted by the vector, −−−→ (*ω*, *<sup>n</sup>*)*ui* . The *<sup>i</sup>*-union of the *<sup>ω</sup>* vectors, −−−→ (*ω*, *<sup>n</sup>*)*ui* , composes the cluster of vectors, −−−→ (*ω*, *<sup>n</sup>*)*u*<sup>1</sup> , −−−→ (*ω*, *<sup>n</sup>*)*u*<sup>2</sup> , ··· , −−−→ (*ω*, *<sup>n</sup>*)*u<sup>ω</sup>* <sup>≡</sup> ←−−→ (*ω*, *n*)*u*. The *u*1, *u*2, ··· , *u<sup>ω</sup>* 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 { ←−−→ (*ω*, *<sup>n</sup>*)*u*}. The clusters ←−−→ (*ω*, *<sup>n</sup>*)*<sup>u</sup>* are symmetrically equivalent to the reference cluster (*ω*, *<sup>n</sup>*) sites by a symmetry operation and are enumerated by an index taking values 1, 2, ··· , Ω[*ω*, *n*]. Ω[*ω*, *n*]

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 Γ(*<sup>s</sup>* ) *<sup>ω</sup>*,*<sup>n</sup>*[*σ*] *<sup>ω</sup>*=2,*n*,(*s*) and *<sup>K</sup>* 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 *<sup>K</sup>*−<sup>1</sup> ) <sup>=</sup> (*<sup>ω</sup>* <sup>+</sup> *<sup>K</sup>* <sup>−</sup> <sup>1</sup>)! (*<sup>K</sup>* <sup>−</sup> <sup>1</sup>)!(*ω*)! 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 Γ(*<sup>s</sup>* ) *<sup>ω</sup>*,*<sup>n</sup>*[*σ*] *ω*,*n*,(*s*) . For convenience in notation, from now on, we will use Γ(*s*) *<sup>ω</sup>*,*n*[*σ*] = Γ(*<sup>s</sup>* ) *<sup>ω</sup>*,*<sup>n</sup>*[*σ*] *ω*,*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 sites {*τ*1, *τ*2, ··· , *τω*} ⊂ (*ω*, *n*) into {*τg*[1],*τg*[2], ··· ,*τg*[*ω*]} ⊂ (*ω*, *n*) as the permutation connecting the decorations (*s*)=(1, 2, ··· , *ω*) ⊂ {(*s*)(*<sup>ω</sup>*,*n*)} and (*s* )=(*g*[1], *g*[2], ··· *g*[*ω*]) ⊂ {(*s*)(*<sup>ω</sup>*,*n*)}. This is given by the following equation:

$$\mathbf{g}\{\vec{\pi}\_1, \vec{\pi}\_2, \dots, \vec{\pi}\_{\omega}\} = \{\vec{\pi}\_{\mathbf{g}[1]}, \vec{\pi}\_{\mathbf{g}[2]}, \dots, \vec{\pi}\_{\mathbf{g}[\omega]}\}.\tag{5}$$

With the symmetry relations between decorations (*s*) and (*s* ) in any given cluster, we are able to retrieve all *K<sup>ω</sup>* 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], *τ<sup>K</sup>*

$$(\tau\_{\mathbb{K}})\_{\vec{\jmath}i} \equiv \gamma\_{\mathbb{J},\mathbb{K}}[\sigma\_{i}]\_{\mathsf{\prime}} \tag{6}$$

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 *τ<sup>K</sup>* matrix can be obtained from the complex conjugate matrix of *τ<sup>K</sup>* by taking the real and imaginary part and riffling their rows. From Equation (4), the following expression results for a system with K components:

$$(\tau\_{\mathcal{K}}^{-1})\_{ij} = \begin{cases} \frac{1}{\mathcal{K}\_{\mathcal{I}}} & \text{if } j = 0 \\ -\frac{2}{\mathcal{K}} \cos \left( 2\pi \lceil \frac{j}{2} \rceil \frac{\upsilon\_{i}}{\mathcal{K}} \right), & \text{if } j > 0 \text{ and } j - 1 < \mathcal{K} \text{ and } j \text{ odd}, \\ -\frac{2}{\mathcal{K}} \sin \left( 2\pi \lceil \frac{j}{2} \rceil \frac{\upsilon\_{i}}{\mathcal{K}} \right), & \text{if } j > 0 \text{ and } j \text{ even}, \\ -\frac{1}{\mathcal{K}} \cos \left( 2\pi \lceil \frac{j}{2} \rceil \frac{\upsilon\_{i}}{\mathcal{K}} \right), & \text{if } j - 1 = \mathcal{K} \text{ and } j \text{ odd}. \end{cases} \tag{7}$$

To the best of our knowledge, Equation (7) represents a new formulation for the inverse of *τ<sup>K</sup>* matrix to ensure that the basis set defined by Equation (4) is rigorously orthonormal. The size of *τ<sup>K</sup>* matrix is *K*x*K*, where *K* is the number of components. It is convenient to perform the matrix multiplications with all *K<sup>ω</sup>* 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 *τ*−<sup>1</sup> *<sup>K</sup>* :

$$
\pi\_4 = \begin{pmatrix} 1 & 1 & 1 & 1 \\ -1 & 0 & 1 & 0 \\ 0 & -1 & 0 & 1 \\ -1 & 1 & -1 & 1 \end{pmatrix} \qquad \pi\_4^{-1} = \frac{1}{4} \begin{pmatrix} 1 & -2 & 0 & -1 \\ 1 & 0 & -2 & 1 \\ 1 & 2 & 0 & -1 \\ 1 & 0 & 2 & 1 \end{pmatrix} . \tag{8}
$$

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:

$$\langle \Gamma\_{\omega,n}^{(ij\cdots)}[\vec{\sigma}] \rangle = \sum\_{\forall [(pq\cdots)]}^{K^{\omega}} \gamma\_{i,\mathcal{K}}[\vec{\sigma}\_p] \gamma\_{j,\mathcal{K}}[\vec{\sigma}\_q] \cdots \jmath\_{\omega,n}^{(pq\cdots)}[\vec{\sigma}] \equiv \overbrace{(\tau\_K \otimes \cdots \otimes \tau\_K)}^{\omega}\_{ij\cdots,pq\cdots} \jmath\_{\omega,n}^{(pq\cdots)}[\vec{\sigma}].\tag{9}$$

With the aid of the matrix formulation from Equation (9) and the generalized form of the inverse matrix *τ*−<sup>1</sup> *<sup>K</sup>* in Equation (7), one can express the *ω*-cluster probabilities into a matrix form:

$$\langle y\_{\omega,n}^{(pq\cdots)} \vert \vec{\sigma} \rangle = \overbrace{\left(\tau\_K^{-1} \otimes \cdots \otimes \tau\_K^{-1}\right)}^{\omega}\_{pq\cdots,j\cdots,j\cdots,\ell} \langle \Gamma\_{\omega,n}^{(ij\cdots)} \vert \vec{\sigma} \rangle \rangle\_{\prime} \tag{10}$$

where we have used the notation for direct product of matrices (*τ*−<sup>1</sup> *<sup>K</sup>* ⊗···⊗ *<sup>τ</sup>*−<sup>1</sup> *<sup>K</sup>* )*pq*··· ,*ij*··· <sup>=</sup> (*τK*)*p*,*i*(*τK*)*q*,*<sup>j</sup>* ··· . 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 *<sup>ω</sup>*, *<sup>K</sup>ω*x*Kω*; in particular, for *<sup>ω</sup>* <sup>=</sup> 4, the matrices have 256 <sup>×</sup> 256 entries. The cluster probabilities are normalized as expressed in Equation (11)

$$\sum\_{\substack{\vec{\sigma}\in[(pq\cdots)]}}^{k^{\omega}} y\_{\omega,\mathfrak{u}}^{(pq\cdots)}[\vec{\sigma}] = 1.\tag{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 (*ω*, *n*) but not the *i*-th sub-cluster (*ωi*, *ni*):

$$\sum\_{(p\eta\cdots\cdots)\in\{\lceil (p\eta\cdots)\_{\omega\_i\eta\_i}\rceil\}} y\_{\omega\_i\mu}^{(p\eta\cdots)}[\vec{\sigma}] = y\_{\omega\_i\mu\_i}^{(p\eta\cdots)}[\vec{\sigma}].\tag{12}$$

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, *α* (*pq*) 2,*<sup>n</sup>* [*σ*]( *p* = *q*), 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 *xp*[*σ*]. For the four component Cr-Fe-Mn-Ni system, there are six chemically distinct SRO parameters: *α*(01) 2,*<sup>n</sup>* [*σ*] for Cr-Fe; *<sup>α</sup>*(02) 2,*<sup>n</sup>* [*σ*] for Cr-Mn; *<sup>α</sup>*(03) 2,*<sup>n</sup>* [*σ*] for Cr-Ni; *α*(12) 2,*<sup>n</sup>* [*σ*] for Fe-Mn; *<sup>α</sup>*(13) 2,*<sup>n</sup>* [*σ*] for Fe-Ni; and *<sup>α</sup>*(23) 2,*<sup>n</sup>* [*σ*] 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*:

$$y\_{2,n}^{(pq)}[\vec{\sigma}] = \mathbf{x}\_p[\vec{\sigma}]\mathbf{x}\_q[\vec{\sigma}](1 - a\_{2,n}^{(pq)}[\vec{\sigma}]).\tag{13}$$

#### *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 Δ*HMixing CE* [*σ*] is described by a set of symmetry unique probability distributions *y* (*pq*···) *<sup>ω</sup>*,*<sup>n</sup>* [*σ*] characterized by decorations of a chosen maximal (*ω*, *<sup>n</sup>*) 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*<sup>1</sup> , *ηω*2,*n*<sup>2</sup> ··· , *ηωs*[*ω*,*n*],*ns*[*ω*,*n*] of the probability distributions of the sub-clusters (*ωi*, *ni*) ⊆ (*ω*, *n*); *i* = 1, ··· , [*ω*, *n*] [19] with decorations (*pq* ···)*ωi*,*ni* corresponding to the components (*pq* ···)*ωi*,*ni* of all decorations (*pq* ···) ∈ {(*pq* ···)(*<sup>ω</sup>*,*n*)} occupied by sites of the (*ωi*, *ni*) sub-cluster

$$y\_{\omega\mu}^{(pq\cdots)}[\vec{\sigma}] = (y\_{\omega\_1,n\_1}^{(pq\cdots)}[\vec{\sigma}])^{\eta\_{\omega\_1\mu\_1}}(y\_{\omega\_2,n\_2}^{(pq\cdots)}[\vec{\sigma}])^{\eta\_{\omega\_2\mu\_2}}\cdots \left(y\_{\omega\_{[\omega\_{[\omega]}n]}n\_{s[\omega]}}^{(pq\cdots)}[\vec{\sigma}]\right)^{\eta\_{\omega\_s[\omega]}n\_s[\omega\_s]}[\vec{\sigma}])^{\eta\_{\omega\_s[\omega]}n\_s[\omega\_s\mu]}.\tag{14}$$

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, <sup>∑</sup>[*ω*,*n*] *<sup>i</sup>*=<sup>1</sup> *ωiηωi*,*ni* = −1. 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:

$$\chi\_{\omega\_{i},n\_{i}}^{(pq\cdots)}[\vec{\sigma}] \underset{T \to \infty}{\to} \prod\_{j=1}^{\omega\_{i}} y\_{1,1}^{((pq\cdots)\_{j})\_{\omega\_{j},n\_{i}}}[\vec{\sigma}] = \prod\_{j=1}^{\omega\_{i}} \mathbf{x}\_{((pq\cdots)\_{j})\_{\omega\_{i},n\_{i}}}[\vec{\sigma}]\_{\prime} \tag{15}$$

where *y* ((*pq*···)*j*)*ωi*,*ni* 1,1 [*σ*] = *<sup>x</sup>*(*pq*···*j*)*ωi*,*ni* [*σ*] is the concentration of the (*pq* ···)*<sup>j</sup>* 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

$$\mathcal{S}\_{\omega,n}[\vec{\sigma}] \equiv \sum\_{i=1}^{s[\omega,n]} \eta\_{\omega\_{i},n\_{i}} \widetilde{\mathcal{S}}\_{\omega\_{i},n\_{i}}[\vec{\sigma}] = \sum\_{i=1}^{s[\omega,n]} \eta\_{\omega\_{i},n\_{i}} \sum\_{\forall \left( (p\eta \cdots)\_{\omega\_{i}n\_{i}} \right)}^{\mathbf{s}[\omega,n]} k y\_{\omega\_{i},n\_{i}}^{(pq\cdots)\_{\omega\_{j}n\_{i}}}[\vec{\sigma}] \ln(y\_{\omega\_{i},n\_{i}}^{(pq\cdots)\_{\omega\_{j}n\_{i}}}[\vec{\sigma}]),\tag{16}$$

where <sup>∼</sup> *Sωi*,*ni* , is the entropy contribution to cluster (*ω*, *n*) from the sub-cluster (*ωi*, *ni*):

$$\widetilde{S}\_{\omega\_{i},n\_{i}}[\vec{\sigma}] = \sum\_{\forall (p\eta\cdots)\_{\omega\_{i}n\_{i}}}^{K^{\omega\_{i}}} k\_{B} y\_{\omega\_{i},n\_{i}}^{(p\eta\cdots)\_{\omega\_{i}n\_{i}}}[\vec{\sigma}] \ln(y\_{\omega\_{i},n\_{i}}^{(p\eta\cdots)\_{\omega\_{i}n\_{i}}}[\vec{\sigma}]).\tag{17}$$

The set of integers *ηω*1,*n*<sup>1</sup> , *ηω*2,*n*<sup>2</sup> ··· , *ηωs*[*ω*,*n*],*ns*[*ω*,*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 *<sup>N</sup>ωi*,*ni* and the sub-cluster multiplicity *<sup>N</sup><sup>β</sup> ωi*,*ni* . The site multiplicity can be determined by calculating the number of symmetry operators, N*ωi*,*ni* , that stabilize the *ωi*, *ni* cluster i.e., *g* ∈ *G* such that the application of *g* into the set of cluster positions {*τ*1, *<sup>τ</sup>*2, ··· , *<sup>τ</sup>ω*} results in a permutation of the set. Then, *<sup>N</sup>ωi*,*ni* <sup>=</sup> <sup>|</sup>*G*<sup>|</sup> |N*ωi*,*ni* | , where |*G*| is the order of the point group associated with the space group *G* and |N*ωi*,*ni* | is the order of the group N*ωi*,*ni* . The sub-cluster multiplicity *N<sup>β</sup> <sup>ω</sup>i*,*ni* is just the frequency of the cluster *β* that is contained in the cluster (*ωi*, *ni*)

$$\eta\_{\omega\_{\bar{\nu}}\imath\_{\bar{\nu}}} = -\aleph\_{\omega\_{\bar{\nu}}\imath\_{\bar{\nu}}} - \sum\_{(\omega\_{\bar{\nu}}\imath\_{\bar{\nu}}) \subsetneq \pounds\_{\bar{\nu}}(\omega,n)} \aleph\_{\omega\_{\bar{\nu}}\imath\_{\bar{\nu}}}^{\notin} \eta\_{\beta}.\tag{18}$$

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=11,*n*1=<sup>1</sup> = −1 and for a 2-body cluster (*ω* = 2, *n*) in the *n*th coordination shell, we

get *ηω*2,*n*<sup>2</sup> = −*N*2,*<sup>n</sup>* and *ηω*1,*n*<sup>1</sup> = 2*N*2,*<sup>n</sup>* − 1, where *N*2,*<sup>n</sup>* 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, *Cconf* derived from the fluctuations of enthalpy of mixing at temperature values in a fine grid of temperature values:

$$S\_{conf}[T] = \int\_0^T \frac{C\_{conf}}{T'} dT' = \int\_0^T \frac{\left<\Delta H\_{CE}^{Mixing}\right>^2 - \left<\Delta H\_{CE}^{Mixing}\right>^2}{T'^3} dT' \tag{19}$$

where Δ*HMixing CE* 2 and Δ*HMixing 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.

#### *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–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]3d74s1, [Ar]3d54s1, [Ar]3d64s1 and [Ar]3d94s1, 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−<sup>6</sup> eV/cell, and force components were relaxed to 10−<sup>3</sup> 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, Δ*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 Cr18Fe27Mn27Ni28 and equiatomic Cr25Fe25Mn25Ni25. 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 *kBln*(4) for equiatomic and −*kB*{0.18 ∗ *ln*(0.18) + 2 ∗ 0.27 ∗ *ln*(0.27) + 0.28 ∗ *ln*(0.28)} for Cr18Fe27Mn27Ni28 [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 *FeNi*<sup>3</sup> in *L*12 structure and anti-ferromagnetic *MnNi* in *L*10 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*<sup>4</sup> 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).

#### *3.2. Full Set of Cluster Decorations*

In general, the set of temperature dependent decorated cluster correlation functions, Γ(*s*) *<sup>ω</sup>*,*n*[*σ*], constitute a set of *K<sup>ω</sup>* quantities for a given maximal cluster (*ω*, *n*) in the CE. From each of the *K<sup>ω</sup>* 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*, *ni*) ∈ (*ω*, *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*, *ni*) is included in the maximal cluster, (*ωi*, *ni*) ∈ (*ω*, *n*), then the decorations from the sub-cluster (*s*)(*<sup>ω</sup>i*,*ni*) can be transferred into the decorations of the maximal cluster, (*s*)(*<sup>ω</sup>*,*n*), by using the intrinsic space group symmetry of the parent lattice. In general, any given sub-cluster (*ωi*, *ni*) 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*) *<sup>ω</sup>*,*n*[*σ*], corresponding to the maximal cluster (*ω*, *n*), is generated by studying which of the decorations (*s*)*ωi*,*ni* form the sub-clusters (*ωi*, *ni*) in (*ω*, *n*). The decorations can have empty cluster, i.e., at least one integer entry in 0 ⊆ (*s*)*ωi*,*ni* .

For the case of the maximal cluster (*ω* = 2, *n* = 1) consisting of sites {(1, 1, 1),(1, 3/2, 3/2)}, see Table 1, we have a full set of 4<sup>2</sup> = 16 2-tuples of decorations, each associated with a temperature dependent decorated cluster correlation function. It can be noted that the cluster {(1, 1, 1),(1, 3/2, 3/2)} contains two sub-clusters: the point cluster {(1, 1, 1)} that has decorations {(0),(1),(2),(3)}; and the cluster itself which has decorations {(1, 1),(2, 1),(3, 1),(2, 2),(2, 3),(3, 3)}. For each of these decorations, there are corresponding decorated cluster correlation functions. In this case, the point cluster is contained twice: once in (1, 1, 1) ⊂ {(1, 1, 1),(1, 3/2, 3/2)} and also in (1, 3/2, 3/2) ⊂ {(1, 1, 1),(1, 3/2, 3/2)}. Similarly the cluster itself {(1, 1, 1),(1, 3/2, 3/2)} is contained once. Thus a point cluster decorated by (*s*)=(3) is transfered to *ω* = 2-tuples notation in the maximal cluster as (3, 0) or (0, 3). This can be captured by the permutation operators {(1, *X*),(*X*, 1)} (see Table 2) where *X* stands for an empty cluster where a 0 should be placed. Similarly, the decorations denoted by (*s*)=(1, 2) in the 2-tuples notation are equivalent to the symmetry equivalent decorations (*s*)=(1, 2) and (*s* )=(2, 1) ≡ (1, 2). The symmetry effect can again be captured by permutation operators {(1, 2),(2, 1)}, where *X* no longer appears, since the *ω<sup>i</sup>* = 2-tuples have the same order (*ω* = *ωi*) as the *ω* = 2-tuples from the maximal cluster (*ω* = 2, *n* = 1). The analysis outlined above for the 2-body cluster {(1, 1, 1),(1, 3/2, 3/2)} with *K* = 4 components can be extended to the remaining (*ω* = 2, *n* > 1) maximal clusters. The permutation operators that carry out the appropriate decorations (*s*)*ωi*,*ni* of sub-clusters (*ωi*, *ni*) into (*ω*, *n*) 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 4 + 6 = 10 symmetrically unique correlation functions <sup>Γ</sup>(*ij*) *<sup>ω</sup>*,*n*[*σ*].

Referring to the maximal cluster (*ω* = 3, *n* = 1), from Table 2, there are three contained sub-clusters (*ω*<sup>1</sup> = 1, *n*<sup>1</sup> = 1), (*ω*<sup>2</sup> = 2, *n*<sup>2</sup> = 1) and (*ω*<sup>3</sup> = 3, *n*<sup>3</sup> = 1) within the maximal cluster (*ω* = 3, *n* = 1). The point sub-cluster (*ω*<sup>1</sup> = 1, *n*<sup>1</sup> = 1) has, according to Table 1, four decorations associated with it; the 2-body cluster (*ω*<sup>2</sup> = 2, *n*<sup>2</sup> = 1) has six decorations associated with itself; and finally the sub-cluster (*ω*<sup>3</sup> = 3, *n*<sup>3</sup> = 1) has 10 decorations associated with itself. The total number of symmetry unique decorated cluster correlation functions, which is smaller than *K<sup>ω</sup>* = 43 = 64, but can nevertheless fully describe cluster (*ω* = 3, *n* = 1) with *K* = 4 components, is therefore given by 4 + 6 + 10 = 20 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 (*ω*<sup>1</sup> = 1, *n*<sup>1</sup> = 1) to (*ω* = 3, *n* = 1), is 3; from (*ω*<sup>2</sup> = 2, *n*<sup>2</sup> = 1), six permutation operators, and, from the cluster itself (*ω*<sup>3</sup> = 3, *n*<sup>3</sup> = 1), there are three permutation operators. The total number of permutation operators is therefore 3 + 6 + 6 = 15. The remaining decorations, which are not listed in Table 1, add up to 64 − 20 = 44 and are obtained by using appropriately the 15 permutation operators corresponding to the maximal cluster (*ω* = 3, *n* = 1). Similarly, the maximal cluster (*ω* = 3, *n* = 2) has 34 symmetrically unique decorations and 12 permutation operators; and (*ω* = 4, *n* = 1) has 35 symmetry unique decorations with 64 permutation operators (all possible permutations for *K* = 4).


**Table 1.** *ω*, highest coordination shell *n*, decoration (*s*) 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 (*s*) is the same as the sequence of points in the relevant cluster. The canonical order for decoration indexes, (*s*)*i*, 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.

#### *3.3. Four-Body Probability Functions from Monte Carlo Simulations*

After performing semi-canonical Monte Carlo simulations for alloy compositions Cr18Fe27Mn27Ni28 and equiatomic Cr25Fe25Mn25Ni25, 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 *ω*, *n* and temperature, we generate the corresponding *τ*−<sup>1</sup> matrices with dimension *Kω*x*K<sup>ω</sup>* and a vector formed by all the decorated correlation cluster functions, <sup>Γ</sup>(*ij*···) *<sup>ω</sup>*,*<sup>n</sup>* [*σ*] with dimension *<sup>K</sup>ω*. The Monte Carlo simulations

output only the symmetrically unique correlation functions (see Table 1), which are less than the *K<sup>ω</sup>* required for matrix operation. In order to generate the full set of *K<sup>ω</sup>* from the symmetrically unique correlation functions, we devise a set of permutation operators *gi* (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, *y* (*s*) *<sup>ω</sup>*,*n*, for any possible decoration, (*s*), 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 ((1, 1, 1),(3/2, 3/2, 1),(3/2, 1, 1/2),(1, 3/2, 1/2)) 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 Cr25Fe25Mn25Ni25 as a function of temperature; (**b**) the same as in (**a**) but for the composition Cr18Fe27Mn27Ni28.

The most important 4-body bonding configurations in the temperature range from 0 to 3000 K for the two alloy systems (the equiatomic Cr25Fe25Mn25Ni25 and Cr18Fe27Mn27Ni28) 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 *FeCr*2*MnNi*<sup>4</sup> 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.

For the system with average composition Cr18Fe27Mn27Ni28, 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.

**Table 2.** List of permutation operators for generating the full set of decorations represented by *K<sup>ω</sup>* dimensional integer arrays with entries taking values from 0 to *K* − 1. 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 (*s*) column in Table 1 by permuting the entries in (*s*) 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, *<sup>g</sup>* <sup>∈</sup> *<sup>O</sup>*<sup>5</sup> *<sup>h</sup>*, is implicitly assumed in order to convolute the symmetry unique into the full set of decorations. The last four columns represent (*ωi*, *ni*) *<sup>i</sup>*th sub-cluster of the maximal cluster (*ω*, *<sup>n</sup>*); *<sup>N</sup>ωi*,*ni* site multiplicity of (*ωi*, *ni*); *<sup>N</sup>ωi*,*ni <sup>ω</sup>*,*<sup>n</sup>* sub-cluster multiplicity of the cluster (*ωi*, *ni*) ∈ (*ω*, *n*); and *ηωi*,*ni* , the sub-cluster (*ωi*, *ni*) contribution to configuration entropy expression corresponding to the maximal cluster (*ω*, *n*)


The prediction of the high probabilities for Mn-Mn-Ni-Ni and Cr-Fe-Fe-Fe clusters in Cr18Fe27Mn27Ni28 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 L10 structure in Strukturbericht notation with Mn and Ni atoms, whereas the composition Cr-Fe-Fe-Fe is related the L12 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 *MnNi*-*L*10 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 *CrFe*3-*L*12 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 *FeCrNi* system persists until temperatures of 600 K, and that for the atomic compositions between CrFe2 and CrFe3 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 Cr25Fe25Mn25Ni25 and Cr18Fe27Mn27Ni28 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**) L12 phase corresponding to CrFe3 with Cr and Fe; (**c**) L10 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 (*<sup>ω</sup>* <sup>=</sup> 4, *<sup>n</sup>* <sup>=</sup> 1) are used in Equation (17) to obtain the quantities <sup>∼</sup> *Sωi*,*ni* [*σ*]. There is one quantity <sup>∼</sup> *Sωi*,*ni* [*σ*] and a corresponding weight *ηωi*,*ni* for each of the (*ωi*, *ni*) sub-clusters : (*ω* = 1, *n* = 1), (*ω* = 2, *n* = 1); (*ω* = 3, *n* = 1); and (*ω* = 4, *n* = 1) included in the maximal cluster (*ω* = 4, *n* = 1). The total configuration entropy, *<sup>S</sup>ω*,*n*[*σ*], is calculated by adding the *ηωi*,*ni* weighted quantities <sup>∼</sup> *Sωi*,*ni* [*σ*]. 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 *kB*, 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 (−*kB*∑<sup>4</sup> *<sup>p</sup> xp*[*σ*]*ln*(*xp*[*σ*])) 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 Cr25Fe25Mn25Ni25.

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 L10 towards the disorder configuration *A*1 in Figure 2. Similarly, in the alloy with composition Cr18Fe27Mn27Ni28, 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 L12-CrFe3 and L10-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.

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

**Author Contributions:** A.F.-C. and D.N.-M. developed theoretical formalism and applied the formulation to the CrFeMnNi system. A.F.-C., D.N.-M. and J.S.W. wrote the manuscript and P.M.M. checked grammar. J.S.W. and M.F. provided the results of DFT-based semi canonical Monte Carlo simulations for the CrFeMnNi system. All the authors made comments on the final version before its submission.

*Entropy* **2019**, *21*, 68

**Funding:** This research was funded by the Euratom research and training programme 2014–2018 and 2019–2020 under Grant Agreement No. 633053 and by the Research Council UK (RCUK) Energy Programme (Grant Number EP/P012450/1). The views and opinions expressed herein do not necessarily reflect those of the European Commission. A.F.-C. was funded by the EPSRC grant (EP/L01680X/1) through the Materials for Demanding Environments Center for Doctoral Training. M.F. and J.S.W. were funded by the Foundation for Polish Science grant HOMING (no. Homing/2016-1/12). The HOMING programme is co-financed by the European Union under the European Regional Development Fund.

**Acknowledgments:** The simulations were partially carried out by M.F. with the support of the Interdisciplinary Centre for Mathematical and Computational Modelling (ICM), University of Warsaw, under Grant No. GA69-30. D.N.-M. and J.S.W. acknowledges the support from high-performance computing facility MARCONI (Bologna, Italy) provided by EUROfusion. D.N.-M. also acknowledges the support from the Institute of Materials Science (IMS) at Los Alamos (NM, USA) for the IMS Rapid Response 2018 visit to the Los Alamos National Laboratory (LANL).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
